MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pumi.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "pumi.hpp"
13 
14 #ifdef MFEM_USE_PUMI
15 #ifdef MFEM_USE_MPI
16 
17 #include "mesh_headers.hpp"
18 #include "../fem/fem.hpp"
19 #include "../general/sort_pairs.hpp"
20 #include "../general/text.hpp"
21 #include "../general/sets.hpp"
22 
23 #include <iostream>
24 #include <sstream>
25 #include <fstream>
26 #include <limits>
27 #include <cmath>
28 #include <cstring>
29 #include <ctime>
30 
31 using namespace std;
32 
33 namespace mfem
34 {
35 
36 PumiMesh::PumiMesh(apf::Mesh2* apf_mesh, int generate_edges, int refine,
37  bool fix_orientation)
38 {
39  Load(apf_mesh, generate_edges, refine, fix_orientation);
40 }
41 
42 Element *PumiMesh::ReadElement(apf::MeshEntity* Ent, const int geom,
43  apf::Downward Verts,
44  const int Attr, apf::Numbering* vert_num)
45 {
46  Element *el;
47  int nv, *v;
48 
49  // Create element in MFEM
50  el = NewElement(geom);
51  nv = el->GetNVertices();
52  v = el->GetVertices();
53 
54  // Fill the connectivity
55  for (int i = 0; i < nv; ++i)
56  {
57  v[i] = apf::getNumber(vert_num, Verts[i], 0, 0);
58  }
59 
60  // Assign attribute
61  el->SetAttribute(Attr);
62 
63  return el;
64 }
65 
66 void PumiMesh::CountBoundaryEntity(apf::Mesh2* apf_mesh, const int BcDim,
67  int &NumBc)
68 {
69  apf::MeshEntity* ent;
70  apf::MeshIterator* itr = apf_mesh->begin(BcDim);
71 
72  while ((ent=apf_mesh->iterate(itr)))
73  {
74  apf::ModelEntity* mdEnt = apf_mesh->toModel(ent);
75  if (apf_mesh->getModelType(mdEnt) == BcDim)
76  {
77  NumBc++;
78  }
79  }
80  apf_mesh->end(itr);
81 
82  // Check if any boundary is detected
83  if (NumBc==0)
84  {
85  MFEM_ABORT("no boundary detected!");
86  }
87 }
88 
89 void PumiMesh::Load(apf::Mesh2* apf_mesh, int generate_edges, int refine,
90  bool fix_orientation)
91 {
92  int curved = 0, read_gf = 1;
93 
94  // Add a check on apf_mesh just in case
95  Clear();
96 
97  // First number vertices
98  apf::Field* apf_field_crd = apf_mesh->getCoordinateField();
99  apf::FieldShape* crd_shape = apf::getShape(apf_field_crd);
100  apf::Numbering* v_num_loc = apf::createNumbering(apf_mesh, "VertexNumbering",
101  crd_shape, 1);
102  // Check if it is a curved mesh
103  curved = (crd_shape->getOrder() > 1) ? 1 : 0;
104 
105  // Read mesh
106  ReadSCORECMesh(apf_mesh, v_num_loc, curved);
107 #ifdef MFEM_DEBUG
108  mfem::out << "After ReadSCORECMesh" << endl;
109 #endif
110  // at this point the following should be defined:
111  // 1) Dim
112  // 2) NumOfElements, elements
113  // 3) NumOfBdrElements, boundary
114  // 4) NumOfVertices, with allocated space in vertices
115  // 5) curved
116  // 5a) if curved == 0, vertices must be defined
117  // 5b) if curved != 0 and read_gf != 0,
118  // 'input' must point to a GridFunction
119  // 5c) if curved != 0 and read_gf == 0,
120  // vertices and Nodes must be defined
121 
122  // FinalizeTopology() will:
123  // - assume that generate_edges is true
124  // - assume that refine is false
125  // - does not check the orientation of regular and boundary elements
126  FinalizeTopology();
127 
128  if (curved && read_gf)
129  {
130  // Check it to be only Quadratic if higher order
131  Nodes = new GridFunctionPumi(this, apf_mesh, v_num_loc,
132  crd_shape->getOrder());
133  edge_vertex = NULL;
134  own_nodes = 1;
135  spaceDim = Nodes->VectorDim();
136 
137  // Set the 'vertices' from the 'Nodes'
138  for (int i = 0; i < spaceDim; i++)
139  {
140  Vector vert_val;
141  Nodes->GetNodalValues(vert_val, i+1);
142  for (int j = 0; j < NumOfVertices; j++)
143  {
144  vertices[j](i) = vert_val(j);
145  }
146  }
147  }
148 
149  // Delete numbering
150  apf::destroyNumbering(v_num_loc);
151 
152  Finalize(refine, fix_orientation);
153 }
154 
155 void PumiMesh::ReadSCORECMesh(apf::Mesh2* apf_mesh, apf::Numbering* v_num_loc,
156  const int curved)
157 {
158  // Here fill the element table from SCOREC MESH
159  // The vector of element pointers is generated with attr and connectivity
160 
161  apf::MeshIterator* itr = apf_mesh->begin(0);
162  apf::MeshEntity* ent;
163  NumOfVertices = 0;
164  while ((ent = apf_mesh->iterate(itr)))
165  {
166  // IDs start from 0
167  apf::number(v_num_loc, ent, 0, 0, NumOfVertices);
168  NumOfVertices++;
169  }
170  apf_mesh->end(itr);
171 
172  Dim = apf_mesh->getDimension();
173  NumOfElements = countOwned(apf_mesh,Dim);
174  elements.SetSize(NumOfElements);
175 
176  // Read elements from SCOREC Mesh
177  itr = apf_mesh->begin(Dim);
178  unsigned int j=0;
179  while ((ent = apf_mesh->iterate(itr)))
180  {
181  // Get vertices
182  apf::Downward verts;
183  apf_mesh->getDownward(ent,0,verts); // num_vert
184  // Get attribute Tag vs Geometry
185  int attr = 1;
186 
187  int geom_type = apf_mesh->getType(ent);
188  elements[j] = ReadElement(ent, geom_type, verts, attr, v_num_loc);
189  j++;
190  }
191  // End iterator
192  apf_mesh->end(itr);
193 
194  // Read Boundaries from SCOREC Mesh
195  // First we need to count them
196  int BCdim = Dim - 1;
197  NumOfBdrElements = 0;
198  CountBoundaryEntity(apf_mesh, BCdim, NumOfBdrElements);
199  boundary.SetSize(NumOfBdrElements);
200  j=0;
201 
202  // Read boundary from SCOREC mesh
203  itr = apf_mesh->begin(BCdim);
204  while ((ent = apf_mesh->iterate(itr)))
205  {
206  // Check if this mesh entity is on the model boundary
207  apf::ModelEntity* mdEnt = apf_mesh->toModel(ent);
208  if (apf_mesh->getModelType(mdEnt) == BCdim)
209  {
210  apf::Downward verts;
211  apf_mesh->getDownward(ent, 0, verts);
212  int attr = 1;
213  int geom_type = apf_mesh->getType(ent);
214  boundary[j] = ReadElement( ent, geom_type, verts, attr, v_num_loc);
215  j++;
216  }
217  }
218  apf_mesh->end(itr);
219 
220  // Fill vertices
221  vertices.SetSize(NumOfVertices);
222 
223  if (!curved)
224  {
225  apf::MeshIterator* itr = apf_mesh->begin(0);
226  spaceDim = Dim;
227 
228  while ((ent = apf_mesh->iterate(itr)))
229  {
230  unsigned int id = apf::getNumber(v_num_loc, ent, 0, 0);
231  apf::Vector3 Crds;
232  apf_mesh->getPoint(ent,0,Crds);
233 
234  for (int ii=0; ii<spaceDim; ii++)
235  {
236  vertices[id](ii) = Crds[ii];
237  }
238  }
239  apf_mesh->end(itr);
240  }
241 }
242 
243 // ParPumiMesh implementation
244 Element *ParPumiMesh::ReadElement(apf::MeshEntity* Ent, const int geom,
245  apf::Downward Verts,
246  const int Attr, apf::Numbering* vert_num)
247 {
248  Element *el;
249  int nv, *v;
250 
251  // Create element in MFEM
252  el = NewElement(geom);
253  nv = el->GetNVertices();
254  v = el->GetVertices();
255 
256  // Fill the connectivity
257  for (int i = 0; i < nv; ++i)
258  {
259  v[i] = apf::getNumber(vert_num, Verts[i], 0, 0);
260  }
261 
262  // Assign attribute
263  el->SetAttribute(Attr);
264 
265  return el;
266 }
267 
268 // This function loads a parallel PUMI mesh and returns the parallel MFEM mesh
269 // corresponding to it.
270 ParPumiMesh::ParPumiMesh(MPI_Comm comm, apf::Mesh2* apf_mesh)
271 {
272  // Set the communicator for gtopo
273  gtopo.SetComm(comm);
274 
275  MyComm = comm;
276  MPI_Comm_size(MyComm, &NRanks);
277  MPI_Comm_rank(MyComm, &MyRank);
278 
279  Dim = apf_mesh->getDimension();
280  spaceDim = Dim;// mesh.spaceDim;
281 
282  apf::MeshIterator* itr;
283  apf::MeshEntity* ent;
284 
285  // Global numbering of vertices. This is necessary to build a local numbering
286  // that has the same ordering in each process.
287  apf::Numbering* vLocNum =
288  apf::numberOwnedDimension(apf_mesh, "AuxVertexNumbering", 0);
289  apf::GlobalNumbering* VertexNumbering = apf::makeGlobal(vLocNum, true);
290  apf::synchronize(VertexNumbering);
291 
292  // Take this process global vertex IDs and sort
293  Array<Pair<long,int> > thisVertIds(apf_mesh->count(0));
294  itr = apf_mesh->begin(0);
295  for (int i = 0; (ent = apf_mesh->iterate(itr)); i++)
296  {
297  long id = apf::getNumber(VertexNumbering, ent, 0, 0);
298  thisVertIds[i] = Pair<long,int>(id, i);
299  }
300  apf_mesh->end(itr);
301  apf::destroyGlobalNumbering(VertexNumbering);
302  thisVertIds.Sort();
303  // Set thisVertIds[i].one = j where j is such that thisVertIds[j].two = i.
304  // Thus, the mapping i -> thisVertIds[i].one is the inverse of the mapping
305  // j -> thisVertIds[j].two.
306  for (int j = 0; j < thisVertIds.Size(); j++)
307  {
308  const int i = thisVertIds[j].two;
309  thisVertIds[i].one = j;
310  }
311 
312  // Create local numbering that respects the global ordering
313  apf::Field* apf_field_crd = apf_mesh->getCoordinateField();
314  apf::FieldShape* crd_shape = apf::getShape(apf_field_crd);
315  apf::Numbering* v_num_loc = apf::createNumbering(apf_mesh,
316  "LocalVertexNumbering",
317  crd_shape, 1);
318 
319  // Construct the numbering v_loc_num and set the coordinates of the vertices.
320  NumOfVertices = thisVertIds.Size();
321  vertices.SetSize(NumOfVertices);
322  itr = apf_mesh->begin(0);
323  for (int i = 0; (ent = apf_mesh->iterate(itr)); i++)
324  {
325  const int id = thisVertIds[i].one;
326  // Assign as local number
327  apf::number(v_num_loc, ent, 0, 0, id);
328 
329  apf::Vector3 Crds;
330  apf_mesh->getPoint(ent,0,Crds);
331 
332  for (int ii=0; ii<spaceDim; ii++)
333  {
334  vertices[id](ii) = Crds[ii]; // Assuming the IDs are ordered and from 0
335  }
336  }
337  apf_mesh->end(itr);
338  thisVertIds.DeleteAll();
339 
340  // Fill the elements
341  NumOfElements = countOwned(apf_mesh,Dim);
342  elements.SetSize(NumOfElements);
343 
344  // Read elements from SCOREC Mesh
345  itr = apf_mesh->begin(Dim);
346  for (int j = 0; (ent = apf_mesh->iterate(itr)); j++)
347  {
348  // Get vertices
349  apf::Downward verts;
350  apf_mesh->getDownward(ent,0,verts);
351 
352  // Get attribute Tag vs Geometry
353  int attr = 1;
354  int geom_type = apf_mesh->getType(ent);
355  elements[j] = ReadElement(ent, geom_type, verts, attr, v_num_loc);
356  }
357  // End iterator
358  apf_mesh->end(itr);
359 
360  // Count number of boundaries by classification
361  int BcDim = Dim - 1;
362  itr = apf_mesh->begin(BcDim);
363  NumOfBdrElements = 0;
364  while ((ent=apf_mesh->iterate(itr)))
365  {
366  apf::ModelEntity* mdEnt = apf_mesh->toModel(ent);
367  if (apf_mesh->getModelType(mdEnt) == BcDim)
368  {
369  NumOfBdrElements++;
370  }
371  }
372  apf_mesh->end(itr);
373 
374  boundary.SetSize(NumOfBdrElements);
375  // Read boundary from SCOREC mesh
376  itr = apf_mesh->begin(BcDim);
377  for (int bdr_ctr = 0; (ent = apf_mesh->iterate(itr)); )
378  {
379  // Check if this mesh entity is on the model boundary
380  apf::ModelEntity* mdEnt = apf_mesh->toModel(ent);
381  if (apf_mesh->getModelType(mdEnt) == BcDim)
382  {
383  apf::Downward verts;
384  apf_mesh->getDownward(ent, 0, verts);
385  int attr = 1 ;
386  int geom_type = apf_mesh->getType(ent);
387  boundary[bdr_ctr++] = ReadElement(ent, geom_type, verts, attr,
388  v_num_loc);
389  }
390  }
391  apf_mesh->end(itr);
392 
393  // The next two methods are called by FinalizeTopology() called below:
394  // Mesh::SetMeshGen();
395  // Mesh::SetAttributes();
396 
397  // This is called by the default Mesh constructor
398  // Mesh::InitTables();
399 
400  this->FinalizeTopology();
401 
402  ListOfIntegerSets groups;
403  IntegerSet group;
404 
405  // The first group is the local one
406  group.Recreate(1, &MyRank);
407  groups.Insert(group);
408 
409  MFEM_ASSERT(Dim >= 3 || GetNFaces() == 0,
410  "[proc " << MyRank << "]: invalid state");
411 
412  // Determine shared faces
414  // Initially sfaces[i].one holds the global face id.
415  // Then it is replaced by the group id of the shared face.
416  if (Dim > 2)
417  {
418  // Number the faces globally and enumerate the local shared faces
419  // following the global enumeration. This way we ensure that the ordering
420  // of the shared faces within each group (of processors) is the same in
421  // each processor in the group.
422  apf::Numbering* AuxFaceNum =
423  apf::numberOwnedDimension(apf_mesh, "AuxFaceNumbering", 2);
424  apf::GlobalNumbering* GlobalFaceNum = apf::makeGlobal(AuxFaceNum, true);
425  apf::synchronize(GlobalFaceNum);
426 
427  itr = apf_mesh->begin(2);
428  while ((ent = apf_mesh->iterate(itr)))
429  {
430  if (apf_mesh->isShared(ent))
431  {
432  long id = apf::getNumber(GlobalFaceNum, ent, 0, 0);
433  sfaces.Append(Pair<long,apf::MeshEntity*>(id, ent));
434  }
435  }
436  apf_mesh->end(itr);
437  sfaces.Sort();
438  apf::destroyGlobalNumbering(GlobalFaceNum);
439 
440  // Replace the global face id in sfaces[i].one with group id.
441  for (int i = 0; i < sfaces.Size(); i++)
442  {
443  ent = sfaces[i].two;
444 
445  const int thisNumAdjs = 2;
446  int eleRanks[thisNumAdjs];
447 
448  // Get the IDs
449  apf::Parts res;
450  apf_mesh->getResidence(ent, res);
451  int kk = 0;
452  for (std::set<int>::iterator itr = res.begin();
453  itr != res.end(); ++itr)
454  {
455  eleRanks[kk++] = *itr;
456  }
457 
458  group.Recreate(2, eleRanks);
459  sfaces[i].one = groups.Insert(group) - 1;
460  }
461  }
462 
463  // Determine shared edges
465  // Initially sedges[i].one holds the global edge id.
466  // Then it is replaced by the group id of the shared edge.
467  if (Dim > 1)
468  {
469  // Number the edges globally and enumerate the local shared edges
470  // following the global enumeration. This way we ensure that the ordering
471  // of the shared edges within each group (of processors) is the same in
472  // each processor in the group.
473  apf::Numbering* AuxEdgeNum =
474  apf::numberOwnedDimension(apf_mesh, "EdgeNumbering", 1);
475  apf::GlobalNumbering* GlobalEdgeNum = apf::makeGlobal(AuxEdgeNum, true);
476  apf::synchronize(GlobalEdgeNum);
477 
478  itr = apf_mesh->begin(1);
479  while ((ent = apf_mesh->iterate(itr)))
480  {
481  if (apf_mesh->isShared(ent))
482  {
483  long id = apf::getNumber(GlobalEdgeNum, ent, 0, 0);
484  sedges.Append(Pair<long,apf::MeshEntity*>(id, ent));
485  }
486  }
487  apf_mesh->end(itr);
488  sedges.Sort();
489  apf::destroyGlobalNumbering(GlobalEdgeNum);
490 
491  // Replace the global edge id in sedges[i].one with group id.
492  Array<int> eleRanks;
493  for (int i = 0; i < sedges.Size(); i++)
494  {
495  ent = sedges[i].two;
496 
497  // Number of adjacent element
498  apf::Parts res;
499  apf_mesh->getResidence(ent, res);
500  eleRanks.SetSize(res.size());
501 
502  // Get the IDs
503  int kk = 0;
504  for (std::set<int>::iterator itr = res.begin();
505  itr != res.end(); itr++)
506  {
507  eleRanks[kk++] = *itr;
508  }
509 
510  // Generate the group
511  group.Recreate(eleRanks.Size(), eleRanks);
512  sedges[i].one = groups.Insert(group) - 1;
513  }
514  }
515 
516  // Determine shared vertices
518  // The entries sverts[i].one hold the local vertex ids.
519  Array<int> svert_group;
520  {
521  itr = apf_mesh->begin(0);
522  while ((ent = apf_mesh->iterate(itr)))
523  {
524  if (apf_mesh->isShared(ent))
525  {
526  int vtId = apf::getNumber(v_num_loc, ent, 0, 0);
527  sverts.Append(Pair<int,apf::MeshEntity*>(vtId, ent));
528  }
529  }
530  apf_mesh->end(itr);
531  sverts.Sort();
532 
533  // Determine svert_group
534  svert_group.SetSize(sverts.Size());
535  Array<int> eleRanks;
536  for (int i = 0; i < sverts.Size(); i++)
537  {
538  ent = sverts[i].two;
539 
540  // Number of adjacent element
541  apf::Parts res;
542  apf_mesh->getResidence(ent, res);
543  eleRanks.SetSize(res.size());
544 
545  // Get the IDs
546  int kk = 0;
547  for (std::set<int>::iterator itr = res.begin();
548  itr != res.end(); itr++)
549  {
550  eleRanks[kk++] = *itr;
551  }
552 
553  group.Recreate(eleRanks.Size(), eleRanks);
554  svert_group[i] = groups.Insert(group) - 1;
555  }
556  }
557 
558  // Build group_stria and group_squad.
559  // Also allocate shared_trias, shared_quads, and sface_lface.
560  group_stria.MakeI(groups.Size()-1);
561  group_squad.MakeI(groups.Size()-1);
562  for (int i = 0; i < sfaces.Size(); i++)
563  {
564  apf::Mesh::Type ftype = apf_mesh->getType(sfaces[i].two);
565  if (ftype == apf::Mesh::TRIANGLE)
566  {
567  group_stria.AddAColumnInRow(sfaces[i].one);
568  }
569  else if (ftype == apf::Mesh::QUAD)
570  {
571  group_squad.AddAColumnInRow(sfaces[i].one);
572  }
573  }
574  group_stria.MakeJ();
575  group_squad.MakeJ();
576  {
577  int nst = 0;
578  for (int i = 0; i < sfaces.Size(); i++)
579  {
580  apf::Mesh::Type ftype = apf_mesh->getType(sfaces[i].two);
581  if (ftype == apf::Mesh::TRIANGLE)
582  {
583  group_stria.AddConnection(sfaces[i].one, nst++);
584  }
585  else if (ftype == apf::Mesh::QUAD)
586  {
587  group_squad.AddConnection(sfaces[i].one, i-nst);
588  }
589  }
590  shared_trias.SetSize(nst);
591  shared_quads.SetSize(sfaces.Size()-nst);
592  sface_lface.SetSize(sfaces.Size());
593  }
594  group_stria.ShiftUpI();
595  group_squad.ShiftUpI();
596 
597  // Build group_sedge
598  group_sedge.MakeI(groups.Size()-1);
599  for (int i = 0; i < sedges.Size(); i++)
600  {
601  group_sedge.AddAColumnInRow(sedges[i].one);
602  }
603  group_sedge.MakeJ();
604  for (int i = 0; i < sedges.Size(); i++)
605  {
606  group_sedge.AddConnection(sedges[i].one, i);
607  }
608  group_sedge.ShiftUpI();
609 
610  // Build group_svert
611  group_svert.MakeI(groups.Size()-1);
612  for (int i = 0; i < svert_group.Size(); i++)
613  {
614  group_svert.AddAColumnInRow(svert_group[i]);
615  }
616  group_svert.MakeJ();
617  for (int i = 0; i < svert_group.Size(); i++)
618  {
619  group_svert.AddConnection(svert_group[i], i);
620  }
621  group_svert.ShiftUpI();
622 
623  // Build shared_trias and shared_quads. They are allocated above.
624  {
625  int nst = 0;
626  for (int i = 0; i < sfaces.Size(); i++)
627  {
628  ent = sfaces[i].two;
629 
630  apf::Downward verts;
631  apf_mesh->getDownward(ent,0,verts);
632 
633  int *v, nv = 0;
634  apf::Mesh::Type ftype = apf_mesh->getType(ent);
635  if (ftype == apf::Mesh::TRIANGLE)
636  {
637  v = shared_trias[nst++].v;
638  nv = 3;
639  }
640  else if (ftype == apf::Mesh::QUAD)
641  {
642  v = shared_quads[i-nst].v;
643  nv = 4;
644  }
645  for (int j = 0; j < nv; ++j)
646  {
647  v[j] = apf::getNumber(v_num_loc, verts[j], 0, 0);
648  }
649  }
650  }
651 
652  // Build shared_edges and allocate sedge_ledge
653  shared_edges.SetSize(sedges.Size());
654  sedge_ledge. SetSize(sedges.Size());
655  for (int i = 0; i < sedges.Size(); i++)
656  {
657  ent = sedges[i].two;
658 
659  apf::Downward verts;
660  apf_mesh->getDownward(ent, 0, verts);
661  int id1, id2;
662  id1 = apf::getNumber(v_num_loc, verts[0], 0, 0);
663  id2 = apf::getNumber(v_num_loc, verts[1], 0, 0);
664  if (id1 > id2) { swap(id1,id2); }
665 
666  shared_edges[i] = new Segment(id1, id2, 1);
667  }
668 
669  // Build svert_lvert
670  svert_lvert.SetSize(sverts.Size());
671  for (int i = 0; i < sverts.Size(); i++)
672  {
673  svert_lvert[i] = sverts[i].one;
674  }
675 
676  // Build the group communication topology
677  gtopo.Create(groups, 822);
678 
679  // Determine sedge_ledge and sface_lface
680  FinalizeParTopo();
681 
682  // Set nodes for higher order mesh
683  int curved = (crd_shape->getOrder() > 1) ? 1 : 0;
684  if (curved) // curved mesh
685  {
686  GridFunctionPumi auxNodes(this, apf_mesh, v_num_loc,
687  crd_shape->getOrder());
688  Nodes = new ParGridFunction(this, &auxNodes);
689  Nodes->Vector::Swap(auxNodes);
690  this->edge_vertex = NULL;
691  own_nodes = 1;
692  }
693 }
694 
695 
696 // GridFunctionPumi Implementation needed for high order meshes
697 GridFunctionPumi::GridFunctionPumi(Mesh* m, apf::Mesh2* PumiM,
698  apf::Numbering* v_num_loc,
699  const int mesh_order)
700 {
701  int spDim = m->SpaceDimension();
702  // Note: default BasisType for 'fec' is GaussLobatto.
703  fec = new H1_FECollection(mesh_order, m->Dimension());
704  int ordering = Ordering::byVDIM; // x1y1z1/x2y2z2/...
705  fes = new FiniteElementSpace(m, fec, spDim, ordering);
706  int data_size = fes->GetVSize();
707 
708  // Read PUMI mesh data
709  this->SetSize(data_size);
710  double* PumiData = this->GetData();
711 
712  apf::MeshEntity* ent;
713  apf::MeshIterator* itr;
714 
715  // Assume all element type are the same i.e. tetrahedral
716  const FiniteElement* H1_elem = fes->GetFE(0);
717  const IntegrationRule &All_nodes = H1_elem->GetNodes();
718  int nnodes = All_nodes.Size();
719 
720  // Loop over elements
721  apf::Field* crd_field = PumiM->getCoordinateField();
722 
723  int nc = apf::countComponents(crd_field);
724  int iel = 0;
725  itr = PumiM->begin(m->Dimension());
726  while ((ent = PumiM->iterate(itr)))
727  {
728  Array<int> vdofs;
729  fes->GetElementVDofs(iel, vdofs);
730 
731  // Create PUMI element to interpolate
732  apf::MeshElement* mE = apf::createMeshElement(PumiM, ent);
733  apf::Element* elem = apf::createElement(crd_field, mE);
734 
735  // Vertices are already interpolated
736  for (int ip = 0; ip < nnodes; ip++)
737  {
738  // Take parametric coordinates of the node
739  apf::Vector3 param;
740  param[0] = All_nodes.IntPoint(ip).x;
741  param[1] = All_nodes.IntPoint(ip).y;
742  param[2] = All_nodes.IntPoint(ip).z;
743 
744  // Compute the interpolating coordinates
745  apf::DynamicVector phCrd(nc);
746  apf::getComponents(elem, param, &phCrd[0]);
747 
748  // Fill the nodes list
749  for (int kk = 0; kk < spDim; ++kk)
750  {
751  int dof_ctr = ip + kk * nnodes;
752  PumiData[vdofs[dof_ctr]] = phCrd[kk];
753  }
754  }
755  iel++;
756  apf::destroyElement(elem);
757  apf::destroyMeshElement(mE);
758  }
759  PumiM->end(itr);
760 
761  sequence = 0;
762 }
763 
764 // Copy the adapted mesh to the original mesh and increase the sequence to be
765 // able to Call Update() methods of FESpace, Linear and Bilinear forms.
766 void ParPumiMesh::UpdateMesh(const ParMesh* AdaptedpMesh)
767 {
768  // Destroy the ParMesh data fields.
769  delete pncmesh;
770  pncmesh = NULL;
771 
772  DeleteFaceNbrData();
773 
774  for (int i = 0; i < shared_edges.Size(); i++)
775  {
776  FreeElement(shared_edges[i]);
777  }
778  shared_quads.DeleteAll();
779  shared_trias.DeleteAll();
780  shared_edges.DeleteAll();
781  group_svert.Clear();
782  group_sedge.Clear();
783  group_stria.Clear();
784  group_squad.Clear();
785  svert_lvert.DeleteAll();
786  sedge_ledge.DeleteAll();
787  sface_lface.DeleteAll();
788 
789  // Destroy the Mesh data fields.
790  Destroy();
791 
792  // Assuming Dim, spaceDim, geom type is unchanged
793  MFEM_ASSERT(Dim == AdaptedpMesh->Dim, "");
794  MFEM_ASSERT(spaceDim == AdaptedpMesh->spaceDim, "");
795  MFEM_ASSERT(meshgen == AdaptedpMesh->meshgen, "");
796 
797  NumOfVertices = AdaptedpMesh->GetNV();
798  NumOfElements = AdaptedpMesh->GetNE();
799  NumOfBdrElements = AdaptedpMesh->GetNBE();
800  NumOfEdges = AdaptedpMesh->GetNEdges();
801  NumOfFaces = AdaptedpMesh->GetNFaces();
802 
803  meshgen = AdaptedpMesh->meshgen;
804 
805  // Sequence is increased by one to trigger update in FEspace etc.
806  sequence++;
807  last_operation = Mesh::NONE;
808 
809  // Duplicate the elements
810  elements.SetSize(NumOfElements);
811  for (int i = 0; i < NumOfElements; i++)
812  {
813  elements[i] = AdaptedpMesh->GetElement(i)->Duplicate(this);
814  }
815 
816  // Copy the vertices
817  AdaptedpMesh->vertices.Copy(vertices);
818 
819  // Duplicate the boundary
820  boundary.SetSize(NumOfBdrElements);
821  for (int i = 0; i < NumOfBdrElements; i++)
822  {
823  boundary[i] = AdaptedpMesh->GetBdrElement(i)->Duplicate(this);
824  }
825 
826  // Copy the element-to-face Table, el_to_face
827  el_to_face = (AdaptedpMesh->el_to_face) ?
828  new Table(*(AdaptedpMesh->el_to_face)) : NULL;
829 
830  // Copy the boundary-to-face Array, be_to_face.
831  AdaptedpMesh->be_to_face.Copy(be_to_face);
832 
833  // Copy the element-to-edge Table, el_to_edge
834  el_to_edge = (AdaptedpMesh->el_to_edge) ?
835  new Table(*(AdaptedpMesh->el_to_edge)) : NULL;
836 
837  // Copy the boudary-to-edge Table, bel_to_edge (3D)
838  bel_to_edge = (AdaptedpMesh->bel_to_edge) ?
839  new Table(*(AdaptedpMesh->bel_to_edge)) : NULL;
840 
841  // Copy the boudary-to-edge Array, be_to_edge (2D)
842  AdaptedpMesh->be_to_edge.Copy(be_to_edge);
843 
844  // Duplicate the faces and faces_info.
845  faces.SetSize(AdaptedpMesh->faces.Size());
846  for (int i = 0; i < faces.Size(); i++)
847  {
848  Element *face = AdaptedpMesh->faces[i]; // in 1D the faces are NULL
849  faces[i] = (face) ? face->Duplicate(this) : NULL;
850  }
851  AdaptedpMesh->faces_info.Copy(faces_info);
852 
853  // Do NOT copy the element-to-element Table, el_to_el
854  el_to_el = NULL;
855 
856  // Do NOT copy the face-to-edge Table, face_edge
857  face_edge = NULL;
858 
859  // Copy the edge-to-vertex Table, edge_vertex
860  edge_vertex = (AdaptedpMesh->edge_vertex) ?
861  new Table(*(AdaptedpMesh->edge_vertex)) : NULL;
862 
863  // Copy the attributes and bdr_attributes
864  AdaptedpMesh->attributes.Copy(attributes);
865  AdaptedpMesh->bdr_attributes.Copy(bdr_attributes);
866 
867  // PUMI meshes cannot use NURBS meshes.
868  MFEM_VERIFY(AdaptedpMesh->NURBSext == NULL,
869  "invalid adapted mesh: it is a NURBS mesh");
870  NURBSext = NULL;
871 
872  // PUMI meshes cannot use NCMesh/ParNCMesh.
873  MFEM_VERIFY(AdaptedpMesh->ncmesh == NULL && AdaptedpMesh->pncmesh == NULL,
874  "invalid adapted mesh: it is a non-conforming mesh");
875  ncmesh = NULL;
876  pncmesh = NULL;
877 
878  // Parallel Implications
879  AdaptedpMesh->group_svert.Copy(group_svert);
880  AdaptedpMesh->group_sedge.Copy(group_sedge);
881  group_stria = AdaptedpMesh->group_stria;
882  group_squad = AdaptedpMesh->group_squad;
883  AdaptedpMesh->gtopo.Copy(gtopo);
884 
885  MyComm = AdaptedpMesh->MyComm;
886  NRanks = AdaptedpMesh->NRanks;
887  MyRank = AdaptedpMesh->MyRank;
888 
889  // Duplicate the shared_edges
890  shared_edges.SetSize(AdaptedpMesh->shared_edges.Size());
891  for (int i = 0; i < shared_edges.Size(); i++)
892  {
893  shared_edges[i] = AdaptedpMesh->shared_edges[i]->Duplicate(this);
894  }
895 
896  // Duplicate the shared_trias and shared_quads
897  shared_trias = AdaptedpMesh->shared_trias;
898  shared_quads = AdaptedpMesh->shared_quads;
899 
900  // Copy the shared-to-local index Arrays
901  AdaptedpMesh->svert_lvert.Copy(svert_lvert);
902  AdaptedpMesh->sedge_ledge.Copy(sedge_ledge);
903  AdaptedpMesh->sface_lface.Copy(sface_lface);
904 
905  // Do not copy face-neighbor data (can be generated if needed)
906  have_face_nbr_data = false;
907 
908  // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
909  // and the FiniteElementSpace (as a ParFiniteElementSpace)
910  if (AdaptedpMesh->Nodes)
911  {
912  FiniteElementSpace *fes = AdaptedpMesh->Nodes->FESpace();
913  const FiniteElementCollection *fec = fes->FEColl();
914  FiniteElementCollection *fec_copy =
915  FiniteElementCollection::New(fec->Name());
916  ParFiniteElementSpace *pfes_copy =
917  new ParFiniteElementSpace(this, fec_copy, fes->GetVDim(),
918  fes->GetOrdering());
919  Nodes = new ParGridFunction(pfes_copy);
920  Nodes->MakeOwner(fec_copy);
921  *Nodes = *(AdaptedpMesh->Nodes);
922  own_nodes = 1;
923  }
924 }
925 
926 // Transfer a mixed vector-scalar field (i.e. velocity,pressure) and the
927 // magnitude of the vector field to use for mesh adaptation.
928 void ParPumiMesh::FieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
929  ParGridFunction* grid_vel,
930  ParGridFunction* grid_pr,
931  apf::Field* VelField,
932  apf::Field* PrField,
933  apf::Field* VelMagField)
934 {
935  apf::FieldShape* VelFieldShape = getShape(VelField);
936  int num_nodes = 4 * VelFieldShape->countNodesOn(0) + // Vertex
937  6 * VelFieldShape->countNodesOn(1) + // Edge
938  4 * VelFieldShape->countNodesOn(2) + // Triangle
939  VelFieldShape->countNodesOn(4); // Tetrahedron
940 
941  // Define integration points
942  IntegrationRule pumi_nodes(num_nodes);
943  int ip_cnt = 0;
944  apf::Vector3 xi_crd(0.,0.,0.);
945 
946  // Create a template of dof holders coordinates in parametric coordinates.
947  // The ordering is taken care of when the field is transferred to PUMI.
948 
949  // Dofs on Vertices
950  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
951  double pt_crd[3] = {0., 0., 0.};
952  ip.Set(pt_crd, 3);
953  for (int kk = 0; kk < 3; kk++)
954  {
955  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
956  double pt_crd[3] = {0.,0.,0.};
957  pt_crd[kk] = 1.0;
958  ip.Set(pt_crd, 3);
959  }
960  // Dofs on Edges
961  if (VelFieldShape->hasNodesIn(apf::Mesh::EDGE))
962  {
963  const int nn = VelFieldShape->countNodesOn(apf::Mesh::EDGE);
964  for (int ii = 0; ii < 6; ii++)
965  {
966  for (int jj = 0; jj < nn; jj++)
967  {
968  VelFieldShape->getNodeXi(apf::Mesh::EDGE, jj, xi_crd);
969  xi_crd[0] = 0.5 * (xi_crd[0] + 1.);// from (-1,1) to (0,1)
970  double pt_crd[3] = {0., 0., 0.};
971  switch (ii)
972  {
973  case 0:
974  pt_crd[0] = xi_crd[0];
975  break;
976  case 1:
977  pt_crd[0] = 1. - xi_crd[0];
978  pt_crd[1] = xi_crd[0];
979  break;
980  case 2:
981  pt_crd[1] = xi_crd[0];
982  break;
983  case 3:
984  pt_crd[2] = xi_crd[0];
985  break;
986  case 4:
987  pt_crd[0] = 1. - xi_crd[0];
988  pt_crd[2] = xi_crd[0];
989  break;
990  case 5:
991  pt_crd[1] = 1. - xi_crd[0];
992  pt_crd[2] = xi_crd[0];
993  break;
994  }
995  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
996  ip.Set(pt_crd, 3);
997  }
998  }
999  }
1000  // Dofs on Faces
1001  if (VelFieldShape->hasNodesIn(apf::Mesh::TRIANGLE))
1002  {
1003  const int nn = VelFieldShape->countNodesOn(apf::Mesh::TRIANGLE);
1004  for (int ii = 0; ii < 4; ii++)
1005  {
1006  for (int jj = 0; jj < nn; jj++)
1007  {
1008  VelFieldShape->getNodeXi(apf::Mesh::TRIANGLE, jj, xi_crd);
1009  double pt_crd[3] = {0., 0., 0.};
1010  switch (ii)
1011  {
1012  case 0:
1013  pt_crd[0] = xi_crd[0];
1014  pt_crd[1] = xi_crd[1];
1015  break;
1016  case 1:
1017  pt_crd[0] = xi_crd[0];
1018  pt_crd[2] = xi_crd[2];
1019  break;
1020  case 2:
1021  pt_crd[0] = xi_crd[0];
1022  pt_crd[1] = xi_crd[1];
1023  pt_crd[2] = xi_crd[2];
1024  break;
1025  case 3:
1026  pt_crd[1] = xi_crd[0];
1027  pt_crd[2] = xi_crd[1];
1028  break;
1029  }
1030  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1031  ip.Set(pt_crd, 3);
1032  }
1033  }
1034  }
1035  MFEM_ASSERT(ip_cnt == num_nodes, "");
1036 
1037  // Other dofs
1038  apf::MeshEntity* ent;
1039  apf::MeshIterator* itr = apf_mesh->begin(3);
1040  int iel = 0;
1041  while ((ent = apf_mesh->iterate(itr)))
1042  {
1043  // Get the solution
1044  Vector u_vel, v_vel, w_vel;
1045  grid_vel->GetValues(iel, pumi_nodes, u_vel, 1);
1046  grid_vel->GetValues(iel, pumi_nodes, v_vel, 2);
1047  grid_vel->GetValues(iel, pumi_nodes, w_vel, 3);
1048 
1049  Vector pr;
1050  grid_pr->GetValues(iel, pumi_nodes, pr, 1);
1051 
1052  // Transfer
1053  apf::Downward vtxs;
1054  int num_vts = apf_mesh->getDownward(ent, 0, vtxs);
1055  for (int kk = 0; kk < num_vts; kk++)
1056  {
1057  double mag = u_vel[kk] * u_vel[kk] + v_vel[kk] * v_vel[kk] +
1058  w_vel[kk] * w_vel[kk];
1059  mag = sqrt(mag);
1060  apf::setScalar(VelMagField, vtxs[kk], 0, mag);
1061  // Set vel
1062  double vels[3] = {u_vel[kk], v_vel[kk], w_vel[kk]};
1063  apf::setComponents(VelField, vtxs[kk], 0, vels);
1064 
1065  // Set Pr
1066  apf::setScalar(PrField, vtxs[kk], 0, pr[kk]);
1067  }
1068 
1069  int dofId = num_vts;
1070 
1071  apf::EntityShape* es = VelFieldShape->getEntityShape(apf::Mesh::TET);
1072  // Edge Dofs
1073  if (VelFieldShape->hasNodesIn(apf::Mesh::EDGE))
1074  {
1075  int ndOnEdge = VelFieldShape->countNodesOn(apf::Mesh::EDGE);
1076  Array<int> order(ndOnEdge);
1077 
1078  apf::Downward edges;
1079  int num_edge = apf_mesh->getDownward(ent, apf::Mesh::EDGE, edges);
1080  for (int ii = 0 ; ii < num_edge; ++ii)
1081  {
1082  es->alignSharedNodes(apf_mesh, ent, edges[ii], order);
1083  for (int jj = 0; jj < ndOnEdge; jj++)
1084  {
1085  int cnt = dofId + order[jj];
1086  double mag = u_vel[cnt] * u_vel[cnt] +
1087  v_vel[cnt] * v_vel[cnt] +
1088  w_vel[cnt] * w_vel[cnt];
1089  mag = sqrt(mag);
1090  apf::setScalar(VelMagField, edges[ii], jj, mag);
1091 
1092  // Set vel
1093  double vels[3] = {u_vel[cnt], v_vel[cnt], w_vel[cnt]};
1094  apf::setComponents(VelField, edges[ii], jj, vels);
1095 
1096  // Set Pr
1097  apf::setScalar(PrField, edges[ii], jj, pr[cnt]);
1098 
1099  }
1100  // Counter
1101  dofId += ndOnEdge;
1102  }
1103  }
1104  // Face Dofs
1105  if (VelFieldShape->hasNodesIn(apf::Mesh::TRIANGLE))
1106  {
1107  int ndOnFace = VelFieldShape->countNodesOn(apf::Mesh::TRIANGLE);
1108  Array<int> order(ndOnFace);
1109 
1110  apf::Downward faces;
1111  int num_face = apf_mesh->getDownward(ent, apf::Mesh::TRIANGLE, faces);
1112  for (int ii = 0; ii < num_face; ii++)
1113  {
1114  if ( ndOnFace > 1)
1115  {
1116  es->alignSharedNodes(apf_mesh, ent, faces[ii], order);
1117  }
1118  else
1119  {
1120  order[0] = 0;
1121  }
1122  for (int jj = 0; jj < ndOnFace; jj++)
1123  {
1124  int cnt = dofId + order[jj];
1125  double mag = u_vel[cnt] * u_vel[cnt] +
1126  v_vel[cnt] * v_vel[cnt] +
1127  w_vel[cnt] * w_vel[cnt];
1128  mag = sqrt(mag);
1129  apf::setScalar(VelMagField, faces[ii], jj, mag);
1130 
1131  // Set vel
1132  double vels[3] = {u_vel[cnt], v_vel[cnt], w_vel[cnt]};
1133  apf::setComponents(VelField, faces[ii], jj, vels);
1134 
1135  // Set Pr
1136  apf::setScalar(PrField, faces[ii], jj, pr[cnt]);
1137  }
1138  // Counter
1139  dofId += ndOnFace;
1140  }
1141  }
1142 
1143  iel++;
1144  }
1145  apf_mesh->end(itr);
1146 }
1147 
1148 // Transfer a scalar field its magnitude to use for mesh adaptation.
1149 void ParPumiMesh::FieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
1150  ParGridFunction* grid_pr,
1151  apf::Field* PrField,
1152  apf::Field* PrMagField)
1153 {
1154  apf::FieldShape* PrFieldShape = getShape(PrField);
1155  int num_nodes = 4 * PrFieldShape->countNodesOn(0) + // Vertex
1156  6 * PrFieldShape->countNodesOn(1) + // Edge
1157  4 * PrFieldShape->countNodesOn(2) + // Triangle
1158  PrFieldShape->countNodesOn(4); // Tetrahedron
1159 
1160  // Define integration points
1161  IntegrationRule pumi_nodes(num_nodes);
1162  int ip_cnt = 0;
1163  apf::Vector3 xi_crd(0.,0.,0.);
1164 
1165  // Create a template of dof holders coordinates in parametric coordinates.
1166  // The ordering is taken care of when the field is transferred to PUMI.
1167 
1168  // Dofs on Vertices
1169  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1170  double pt_crd[3] = {0., 0., 0.};
1171  ip.Set(pt_crd, 3);
1172  for (int kk = 0; kk < 3; kk++)
1173  {
1174  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1175  double pt_crd[3] = {0.,0.,0.};
1176  pt_crd[kk] = 1.0;
1177  ip.Set(pt_crd, 3);
1178  }
1179  // Dofs on Edges
1180  if (PrFieldShape->hasNodesIn(apf::Mesh::EDGE))
1181  {
1182  const int nn = PrFieldShape->countNodesOn(apf::Mesh::EDGE);
1183  for (int ii = 0; ii < 6; ii++)
1184  {
1185  for (int jj = 0; jj < nn; jj++)
1186  {
1187  PrFieldShape->getNodeXi(apf::Mesh::EDGE, jj, xi_crd);
1188  xi_crd[0] = 0.5 * (xi_crd[0] + 1.); // from (-1,1) to (0,1)
1189  double pt_crd[3] = {0., 0., 0.};
1190  switch (ii)
1191  {
1192  case 0:
1193  pt_crd[0] = xi_crd[0];
1194  break;
1195  case 1:
1196  pt_crd[0] = 1. - xi_crd[0];
1197  pt_crd[1] = xi_crd[0];
1198  break;
1199  case 2:
1200  pt_crd[1] = xi_crd[0];
1201  break;
1202  case 3:
1203  pt_crd[2] = xi_crd[0];
1204  break;
1205  case 4:
1206  pt_crd[0] = 1. - xi_crd[0];
1207  pt_crd[2] = xi_crd[0];
1208  break;
1209  case 5:
1210  pt_crd[1] = 1. - xi_crd[0];
1211  pt_crd[2] = xi_crd[0];
1212  break;
1213  }
1214  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1215  ip.Set(pt_crd, 3);
1216  }
1217  }
1218  }
1219  // Dofs on Faces
1220  if (PrFieldShape->hasNodesIn(apf::Mesh::TRIANGLE))
1221  {
1222  const int nn = PrFieldShape->countNodesOn(apf::Mesh::TRIANGLE);
1223  for (int ii = 0; ii < 4; ii++)
1224  {
1225  for (int jj = 0; jj < nn; jj++)
1226  {
1227  PrFieldShape->getNodeXi(apf::Mesh::TRIANGLE, jj, xi_crd);
1228  double pt_crd[3] = {0., 0., 0.};
1229  switch (ii)
1230  {
1231  case 0:
1232  pt_crd[0] = xi_crd[0];
1233  pt_crd[1] = xi_crd[1];
1234  break;
1235  case 1:
1236  pt_crd[0] = xi_crd[0];
1237  pt_crd[2] = xi_crd[2];
1238  break;
1239  case 2:
1240  pt_crd[0] = xi_crd[0];
1241  pt_crd[1] = xi_crd[1];
1242  pt_crd[2] = xi_crd[2];
1243  break;
1244  case 3:
1245  pt_crd[1] = xi_crd[0];
1246  pt_crd[2] = xi_crd[1];
1247  break;
1248  }
1249  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1250  ip.Set(pt_crd, 3);
1251  }
1252  }
1253  }
1254  MFEM_ASSERT(ip_cnt == num_nodes, "");
1255 
1256  // Other dofs
1257  apf::MeshEntity* ent;
1258  apf::MeshIterator* itr = apf_mesh->begin(3);
1259  int iel = 0;
1260  while ((ent = apf_mesh->iterate(itr)))
1261  {
1262  // Get the solution
1263  Vector pr;
1264  grid_pr->GetValues(iel, pumi_nodes, pr, 1);
1265 
1266  // Transfer
1267  apf::Downward vtxs;
1268  int num_vts = apf_mesh->getDownward(ent, 0, vtxs);
1269  for (int kk = 0; kk < num_vts; kk++)
1270  {
1271  double mag;
1272  (pr[kk] >= 0. ? mag = pr[kk] : mag = -pr[kk]);
1273  apf::setScalar(PrMagField, vtxs[kk], 0, mag);
1274 
1275  // Set Pr
1276  apf::setScalar(PrField, vtxs[kk], 0, pr[kk]);
1277  }
1278 
1279  int dofId = num_vts;
1280 
1281  apf::EntityShape* es = PrFieldShape->getEntityShape(apf::Mesh::TET);
1282  // Edge Dofs
1283  if (PrFieldShape->hasNodesIn(apf::Mesh::EDGE))
1284  {
1285  int ndOnEdge = PrFieldShape->countNodesOn(apf::Mesh::EDGE);
1286  Array<int> order(ndOnEdge);
1287 
1288  apf::Downward edges;
1289  int num_edge = apf_mesh->getDownward(ent, apf::Mesh::EDGE, edges);
1290  for (int ii = 0 ; ii < num_edge; ++ii)
1291  {
1292  es->alignSharedNodes(apf_mesh, ent, edges[ii], order);
1293  for (int jj = 0; jj < ndOnEdge; jj++)
1294  {
1295  int cnt = dofId + order[jj];
1296  double mag;
1297  (pr[cnt] >= 0. ? mag = pr[cnt] : mag = -pr[cnt]);
1298  apf::setScalar(PrMagField, edges[ii], jj, mag);
1299 
1300  // Set Pr
1301  apf::setScalar(PrField, edges[ii], jj, pr[cnt]);
1302 
1303  }
1304  // Counter
1305  dofId += ndOnEdge;
1306  }
1307  }
1308 
1309  // Face Dofs
1310  if (PrFieldShape->hasNodesIn(apf::Mesh::TRIANGLE))
1311  {
1312  int ndOnFace = PrFieldShape->countNodesOn(apf::Mesh::TRIANGLE);
1313  Array<int> order(ndOnFace);
1314 
1315  apf::Downward faces;
1316  int num_face = apf_mesh->getDownward(ent, apf::Mesh::TRIANGLE, faces);
1317  for (int ii = 0; ii < num_face; ii++)
1318  {
1319  if ( ndOnFace > 1)
1320  {
1321  es->alignSharedNodes(apf_mesh, ent, faces[ii], order);
1322  }
1323  else
1324  {
1325  order[0] = 0;
1326  }
1327  for (int jj = 0; jj < ndOnFace; jj++)
1328  {
1329  int cnt = dofId + order[jj];
1330  double mag;
1331  (pr[cnt] >= 0. ? mag = pr[cnt] : mag = -pr[cnt]);
1332  apf::setScalar(PrMagField, faces[ii], jj, mag);
1333 
1334  // Set Pr
1335  apf::setScalar(PrField, faces[ii], jj, pr[cnt]);
1336  }
1337  // Counter
1338  dofId += ndOnFace;
1339  }
1340  }
1341 
1342  iel++;
1343  }
1344  apf_mesh->end(itr);
1345 }
1346 
1347 // Transfer a vector field and the magnitude of the vector field to use for mesh
1348 // adaptation
1349 void ParPumiMesh::VectorFieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
1350  ParGridFunction* grid_vel,
1351  apf::Field* VelField,
1352  apf::Field* VelMagField)
1353 {
1354  apf::FieldShape* VelFieldShape = getShape(VelField);
1355  int num_nodes = 4 * VelFieldShape->countNodesOn(0) + // Vertex
1356  6 * VelFieldShape->countNodesOn(1) + // Edge
1357  4 * VelFieldShape->countNodesOn(2) + // Triangle
1358  VelFieldShape->countNodesOn(4);// Tetrahedron
1359 
1360  // Define integration points
1361  IntegrationRule pumi_nodes(num_nodes);
1362  int ip_cnt = 0;
1363  apf::Vector3 xi_crd(0.,0.,0.);
1364 
1365  // Create a template of dof holders coordinates in parametric coordinates.
1366  // The ordering is taken care of when the field is transferred to PUMI.
1367 
1368  // Dofs on Vertices
1369  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1370  double pt_crd[3] = {0., 0., 0.};
1371  ip.Set(pt_crd, 3);
1372  for (int kk = 0; kk < 3; kk++)
1373  {
1374  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1375  double pt_crd[3] = {0.,0.,0.};
1376  pt_crd[kk] = 1.0;
1377  ip.Set(pt_crd, 3);
1378  }
1379  // Dofs on Edges
1380  if (VelFieldShape->hasNodesIn(apf::Mesh::EDGE))
1381  {
1382  const int nn = VelFieldShape->countNodesOn(apf::Mesh::EDGE);
1383  for (int ii = 0; ii < 6; ii++)
1384  {
1385  for (int jj = 0; jj < nn; jj++)
1386  {
1387  VelFieldShape->getNodeXi(apf::Mesh::EDGE, jj, xi_crd);
1388  xi_crd[0] = 0.5 * (xi_crd[0] + 1.); // from (-1,1) to (0,1)
1389  double pt_crd[3] = {0., 0., 0.};
1390  switch (ii)
1391  {
1392  case 0:
1393  pt_crd[0] = xi_crd[0];
1394  break;
1395  case 1:
1396  pt_crd[0] = 1. - xi_crd[0];
1397  pt_crd[1] = xi_crd[0];
1398  break;
1399  case 2:
1400  pt_crd[1] = xi_crd[0];
1401  break;
1402  case 3:
1403  pt_crd[2] = xi_crd[0];
1404  break;
1405  case 4:
1406  pt_crd[0] = 1. - xi_crd[0];
1407  pt_crd[2] = xi_crd[0];
1408  break;
1409  case 5:
1410  pt_crd[1] = 1. - xi_crd[0];
1411  pt_crd[2] = xi_crd[0];
1412  break;
1413  }
1414  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1415  ip.Set(pt_crd, 3);
1416  }
1417  }
1418  }
1419  // Dofs on Faces
1420  if (VelFieldShape->hasNodesIn(apf::Mesh::TRIANGLE))
1421  {
1422  const int nn = VelFieldShape->countNodesOn(apf::Mesh::TRIANGLE);
1423  for (int ii = 0; ii < 4; ii++)
1424  {
1425  for (int jj = 0; jj < nn; jj++)
1426  {
1427  VelFieldShape->getNodeXi(apf::Mesh::TRIANGLE, jj, xi_crd);
1428  double pt_crd[3] = {0., 0., 0.};
1429  switch (ii)
1430  {
1431  case 0:
1432  pt_crd[0] = xi_crd[0];
1433  pt_crd[1] = xi_crd[1];
1434  break;
1435  case 1:
1436  pt_crd[0] = xi_crd[0];
1437  pt_crd[2] = xi_crd[2];
1438  break;
1439  case 2:
1440  pt_crd[0] = xi_crd[0];
1441  pt_crd[1] = xi_crd[1];
1442  pt_crd[2] = xi_crd[2];
1443  break;
1444  case 3:
1445  pt_crd[1] = xi_crd[0];
1446  pt_crd[2] = xi_crd[1];
1447  break;
1448  }
1449  IntegrationPoint& ip = pumi_nodes.IntPoint(ip_cnt++);
1450  ip.Set(pt_crd, 3);
1451  }
1452  }
1453  }
1454  MFEM_ASSERT(ip_cnt == num_nodes, "");
1455 
1456  // Other dofs
1457  apf::MeshEntity* ent;
1458  apf::MeshIterator* itr = apf_mesh->begin(3);
1459  int iel = 0;
1460  while ((ent = apf_mesh->iterate(itr)))
1461  {
1462  // Get the solution
1463  Vector u_vel, v_vel, w_vel;
1464  grid_vel->GetValues(iel, pumi_nodes, u_vel, 1);
1465  grid_vel->GetValues(iel, pumi_nodes, v_vel, 2);
1466  grid_vel->GetValues(iel, pumi_nodes, w_vel, 3);
1467 
1468  // Transfer
1469  apf::Downward vtxs;
1470  int num_vts = apf_mesh->getDownward(ent, 0, vtxs);
1471  for (int kk = 0; kk < num_vts; kk++)
1472  {
1473  double mag = u_vel[kk] * u_vel[kk] + v_vel[kk] * v_vel[kk] +
1474  w_vel[kk] * w_vel[kk];
1475  mag = sqrt(mag);
1476  apf::setScalar(VelMagField, vtxs[kk], 0, mag);
1477  // Set vel
1478  double vels[3] = {u_vel[kk], v_vel[kk], w_vel[kk]};
1479  apf::setComponents(VelField, vtxs[kk], 0, vels);
1480  }
1481 
1482  int dofId = num_vts;
1483 
1484  apf::EntityShape* es = VelFieldShape->getEntityShape(apf::Mesh::TET);
1485  // Edge Dofs
1486  if (VelFieldShape->hasNodesIn(apf::Mesh::EDGE))
1487  {
1488  int ndOnEdge = VelFieldShape->countNodesOn(apf::Mesh::EDGE);
1489  Array<int> order(ndOnEdge);
1490 
1491  apf::Downward edges;
1492  int num_edge = apf_mesh->getDownward(ent, apf::Mesh::EDGE, edges);
1493  for (int ii = 0 ; ii < num_edge; ++ii)
1494  {
1495  es->alignSharedNodes(apf_mesh, ent, edges[ii], order);
1496  for (int jj = 0; jj < ndOnEdge; jj++)
1497  {
1498  int cnt = dofId + order[jj];
1499  double mag = u_vel[cnt] * u_vel[cnt] +
1500  v_vel[cnt] * v_vel[cnt] +
1501  w_vel[cnt] * w_vel[cnt];
1502  mag = sqrt(mag);
1503  apf::setScalar(VelMagField, edges[ii], jj, mag);
1504 
1505  // Set vel
1506  double vels[3] = {u_vel[cnt], v_vel[cnt], w_vel[cnt]};
1507  apf::setComponents(VelField, edges[ii], jj, vels);
1508  }
1509  // Counter
1510  dofId += ndOnEdge;
1511  }
1512  }
1513 
1514  // Face Dofs
1515  if (VelFieldShape->hasNodesIn(apf::Mesh::TRIANGLE))
1516  {
1517  int ndOnFace = VelFieldShape->countNodesOn(apf::Mesh::TRIANGLE);
1518  Array<int> order(ndOnFace);
1519 
1520  apf::Downward faces;
1521  int num_face = apf_mesh->getDownward(ent, apf::Mesh::TRIANGLE, faces);
1522  for (int ii = 0; ii < num_face; ii++)
1523  {
1524  if ( ndOnFace > 1)
1525  {
1526  es->alignSharedNodes(apf_mesh, ent, faces[ii], order);
1527  }
1528  else
1529  {
1530  order[0] = 0;
1531  }
1532  for (int jj = 0; jj < ndOnFace; jj++)
1533  {
1534  int cnt = dofId + order[jj];
1535  double mag = u_vel[cnt] * u_vel[cnt] +
1536  v_vel[cnt] * v_vel[cnt] +
1537  w_vel[cnt] * w_vel[cnt];
1538  mag = sqrt(mag);
1539  apf::setScalar(VelMagField, faces[ii], jj, mag);
1540 
1541  // Set vel
1542  double vels[3] = {u_vel[cnt], v_vel[cnt], w_vel[cnt]};
1543  apf::setComponents(VelField, faces[ii], jj, vels);
1544  }
1545  // Counter
1546  dofId += ndOnFace;
1547  }
1548  }
1549 
1550  iel++;
1551  }
1552  apf_mesh->end(itr);
1553 }
1554 
1555 void ParPumiMesh::FieldPUMItoMFEM(apf::Mesh2* apf_mesh,
1556  apf::Field* ScalarField,
1557  ParGridFunction* Pr)
1558 {
1559  // Pr->Update();
1560  // Find local numbering
1561  v_num_loc = apf_mesh->findNumbering("LocalVertexNumbering");
1562 
1563  // Loop over field to copy
1564  getShape(ScalarField);
1565  apf::MeshEntity* ent;
1566  apf::MeshIterator* itr = apf_mesh->begin(0);
1567  while ((ent = apf_mesh->iterate(itr)))
1568  {
1569  unsigned int id = apf::getNumber(v_num_loc, ent, 0, 0);
1570  double fieldVal = apf::getScalar(ScalarField, ent, 0);
1571 
1572  (Pr->GetData())[id] = fieldVal;
1573  }
1574  apf_mesh->end(itr);
1575 
1576  // Check for higher order
1577  getShape(ScalarField);
1578  if ( Pr->FESpace()->GetOrder(1) > 1 )
1579  {
1580  // Assume all element type are the same i.e. tetrahedral
1581  const FiniteElement* H1_elem = Pr->FESpace()->GetFE(1);
1582  const IntegrationRule &All_nodes = H1_elem->GetNodes();
1583  int nnodes = All_nodes.Size();
1584 
1585  // Loop over elements
1586  int nc = apf::countComponents(ScalarField);
1587  int iel = 0;
1588  itr = apf_mesh->begin(3);
1589  while ((ent = apf_mesh->iterate(itr)))
1590  {
1591  Array<int> vdofs;
1592  Pr->FESpace()->GetElementVDofs(iel, vdofs);
1593 
1594  // Create PUMI element to interpolate
1595  apf::MeshElement* mE = apf::createMeshElement(apf_mesh, ent);
1596  apf::Element* elem = apf::createElement(ScalarField, mE);
1597 
1598  // Vertices are already interpolated
1599  for (int ip = 0; ip < nnodes; ip++) //num_vert
1600  {
1601  // Take parametric coordinates of the node
1602  apf::Vector3 param;
1603  param[0] = All_nodes.IntPoint(ip).x;
1604  param[1] = All_nodes.IntPoint(ip).y;
1605  param[2] = All_nodes.IntPoint(ip).z;
1606 
1607  // Compute the interpolating coordinates
1608  apf::DynamicVector phCrd(nc);
1609  apf::getComponents(elem, param, &phCrd[0]);
1610 
1611  // Fill the nodes list
1612  for (int kk = 0; kk < nc; ++kk)
1613  {
1614  int dof_ctr = ip + kk * nnodes;
1615  (Pr->GetData())[vdofs[dof_ctr]] = phCrd[kk];
1616  }
1617  }
1618  iel++;
1619  apf::destroyElement(elem);
1620  apf::destroyMeshElement(mE);
1621  }
1622  apf_mesh->end(itr);
1623  }
1624 }
1625 
1626 }
1627 
1628 #endif // MFEM_USE_MPI
1629 #endif // MFEM_USE_SCOREC
Abstract class for Finite Elements.
Definition: fe.hpp:232
void Set(const double *p, const int dim)
Definition: intrules.hpp:37
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:393
int Size() const
Logical size of the array.
Definition: array.hpp:124
void Recreate(const int n, const int *p)
Definition: sets.cpp:59
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
T * end()
Definition: array.hpp:266
int NRanks
Definition: pmesh.hpp:39
virtual Element * Duplicate(Mesh *m) const =0
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:696
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:172
const Geometry::Type geom
Definition: ex1.cpp:40
Array< int > sface_lface
Definition: pmesh.hpp:79
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:812
Class for PUMI grid functions.
Definition: pumi.hpp:118
int VectorDim() const
Definition: gridfunc.cpp:302
GridFunction * Nodes
Definition: mesh.hpp:164
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Array< Vert3 > shared_trias
Definition: pmesh.hpp:66
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:166
Array< Element * > faces
Definition: mesh.hpp:85
Array< int > sedge_ledge
Definition: pmesh.hpp:77
Table group_stria
Definition: pmesh.hpp:72
Table * el_to_face
Definition: mesh.hpp:145
ParNCMesh * pncmesh
Definition: pmesh.hpp:247
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
MPI_Comm MyComm
Definition: pmesh.hpp:38
Array< Element * > shared_edges
Definition: pmesh.hpp:62
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:707
A pair of objects.
Definition: sort_pairs.hpp:23
const IntegrationRule & GetNodes() const
Definition: fe.hpp:367
int Insert(IntegerSet &s)
Definition: sets.cpp:82
T * begin()
Definition: array.hpp:265
const Element * GetElement(int i) const
Definition: mesh.hpp:775
void Sort()
Sorts the array. This requires operator&lt; to be defined for T.
Definition: array.hpp:229
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:442
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:442
int Dimension() const
Definition: mesh.hpp:744
Table * el_to_edge
Definition: mesh.hpp:144
Array< Vert4 > shared_quads
Definition: pmesh.hpp:67
void Copy(GroupTopology &copy) const
Copy.
int SpaceDimension() const
Definition: mesh.hpp:745
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
Array< Vertex > vertices
Definition: mesh.hpp:83
int GetOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:95
int meshgen
Definition: mesh.hpp:70
Table * bel_to_edge
Definition: mesh.hpp:148
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:690
Array< int > be_to_edge
Definition: mesh.hpp:147
virtual const char * Name() const
Definition: fe_coll.hpp:53
Class for integration point with weight.
Definition: intrules.hpp:25
Table group_sedge
Definition: pmesh.hpp:71
A set of integers.
Definition: sets.hpp:23
Table group_svert
Shared objects in each group.
Definition: pmesh.hpp:70
Table * edge_vertex
Definition: mesh.hpp:151
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
Array< FaceInfo > faces_info
Definition: mesh.hpp:141
virtual int GetNVertices() const =0
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:192
int Dim
Definition: mesh.hpp:59
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:699
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:702
Vector data type.
Definition: vector.hpp:48
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:395
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Array< int > svert_lvert
Shared to local index mapping.
Definition: pmesh.hpp:76
int MyRank
Definition: pmesh.hpp:39
int spaceDim
Definition: mesh.hpp:60
Class for parallel grid function.
Definition: pgridfunc.hpp:32
List of integer sets.
Definition: sets.hpp:54
Array< int > be_to_face
Definition: mesh.hpp:149
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
GroupTopology gtopo
Definition: pmesh.hpp:234
Class for parallel meshes.
Definition: pmesh.hpp:32
Abstract data type element.
Definition: element.hpp:28
Data type line segment element.
Definition: segment.hpp:22
Table group_squad
Definition: pmesh.hpp:73
void Copy(Table &copy) const
Definition: table.cpp:390
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:187
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:779