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