MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pumi.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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"
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
31using namespace std;
32
33namespace mfem
34{
35
36static 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
59static 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
76static 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)
93static 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
104static 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
116static 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
130static 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
163static 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
200static 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
223static 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
245PumiMesh::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
253void 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
276void 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
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;
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
342void 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);
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;
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.
443ParPumiMesh::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);
552 while ((ent=apf_mesh->iterate(itr)))
553 {
554 apf::ModelEntity* mdEnt = apf_mesh->toModel(ent);
555 if (apf_mesh->getModelType(mdEnt) == BcDim)
556 {
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 }
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 }
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 }
793 for (int i = 0; i < sedges.Size(); i++)
794 {
795 group_sedge.AddConnection(sedges[i].one, i);
796 }
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 }
806 for (int i = 0; i < svert_group.Size(); i++)
807 {
808 group_svert.AddConnection(svert_group[i], i);
809 }
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
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
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.
956void ParPumiMesh::UpdateMesh(const ParMesh* AdaptedpMesh)
957{
958 // Destroy the ParMesh data fields.
959 delete pncmesh;
960 pncmesh = NULL;
961
963
964 for (int i = 0; i < shared_edges.Size(); i++)
965 {
967 }
968 shared_quads.DeleteAll();
969 shared_trias.DeleteAll();
970 shared_edges.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++;
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 // Duplicate the faces and faces_info.
1032 faces.SetSize(AdaptedpMesh->faces.Size());
1033 for (int i = 0; i < faces.Size(); i++)
1034 {
1035 Element *face = AdaptedpMesh->faces[i]; // in 1D the faces are NULL
1036 faces[i] = (face) ? face->Duplicate(this) : NULL;
1037 }
1038 AdaptedpMesh->faces_info.Copy(faces_info);
1039
1040 // Do NOT copy the element-to-element Table, el_to_el
1041 el_to_el = NULL;
1042
1043 // Do NOT copy the face-to-edge Table, face_edge
1044 face_edge = NULL;
1045
1046 // Copy the edge-to-vertex Table, edge_vertex
1047 edge_vertex = (AdaptedpMesh->edge_vertex) ?
1048 new Table(*(AdaptedpMesh->edge_vertex)) : NULL;
1049
1050 // Copy the attributes and bdr_attributes
1051 AdaptedpMesh->attributes.Copy(attributes);
1052 AdaptedpMesh->bdr_attributes.Copy(bdr_attributes);
1053
1054 // PUMI meshes cannot use NURBS meshes.
1055 MFEM_VERIFY(AdaptedpMesh->NURBSext == NULL,
1056 "invalid adapted mesh: it is a NURBS mesh");
1057 NURBSext = NULL;
1058
1059 // PUMI meshes cannot use NCMesh/ParNCMesh.
1060 MFEM_VERIFY(AdaptedpMesh->ncmesh == NULL && AdaptedpMesh->pncmesh == NULL,
1061 "invalid adapted mesh: it is a non-conforming mesh");
1062 ncmesh = NULL;
1063 pncmesh = NULL;
1064
1065 // Parallel Implications
1066 AdaptedpMesh->group_svert.Copy(group_svert);
1067 AdaptedpMesh->group_sedge.Copy(group_sedge);
1068 group_stria = AdaptedpMesh->group_stria;
1069 group_squad = AdaptedpMesh->group_squad;
1070 AdaptedpMesh->gtopo.Copy(gtopo);
1071
1072 MyComm = AdaptedpMesh->MyComm;
1073 NRanks = AdaptedpMesh->NRanks;
1074 MyRank = AdaptedpMesh->MyRank;
1075
1076 // Duplicate the shared_edges
1077 shared_edges.SetSize(AdaptedpMesh->shared_edges.Size());
1078 for (int i = 0; i < shared_edges.Size(); i++)
1079 {
1080 shared_edges[i] = AdaptedpMesh->shared_edges[i]->Duplicate(this);
1081 }
1082
1083 // Duplicate the shared_trias and shared_quads
1084 shared_trias = AdaptedpMesh->shared_trias;
1085 shared_quads = AdaptedpMesh->shared_quads;
1086
1087 // Copy the shared-to-local index Arrays
1088 AdaptedpMesh->svert_lvert.Copy(svert_lvert);
1089 AdaptedpMesh->sedge_ledge.Copy(sedge_ledge);
1090 AdaptedpMesh->sface_lface.Copy(sface_lface);
1091
1092 // Do not copy face-neighbor data (can be generated if needed)
1093 have_face_nbr_data = false;
1094
1095 // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
1096 // and the FiniteElementSpace (as a ParFiniteElementSpace)
1097 if (AdaptedpMesh->Nodes)
1098 {
1099 FiniteElementSpace *fes = AdaptedpMesh->Nodes->FESpace();
1100 const FiniteElementCollection *fec = fes->FEColl();
1101 FiniteElementCollection *fec_copy =
1103 ParFiniteElementSpace *pfes_copy =
1104 new ParFiniteElementSpace(this, fec_copy, fes->GetVDim(),
1105 fes->GetOrdering());
1106 Nodes = new ParGridFunction(pfes_copy);
1107 Nodes->MakeOwner(fec_copy);
1108 *Nodes = *(AdaptedpMesh->Nodes);
1109 own_nodes = 1;
1110 }
1111}
1112
1113int ParPumiMesh::RotationPUMItoMFEM(apf::Mesh2* apf_mesh,
1114 apf::MeshEntity* ent,
1115 int elemId)
1116{
1117 int type = apf_mesh->getType(ent);
1118 MFEM_CONTRACT_VAR(type);
1119 MFEM_ASSERT(apf::isSimplex(type),
1120 "only implemented for simplex entity types");
1121 // get downward vertices of PUMI element
1122 apf::Downward vs;
1123 int nv = apf_mesh->getDownward(ent,0,vs);
1124 int pumi_vid[12];
1125 for (int i = 0; i < nv; i++)
1126 {
1127 pumi_vid[i] = apf::getNumber(v_num_loc, vs[i], 0, 0);
1128 }
1129
1130 // get downward vertices of MFEM element
1131 mfem::Array<int> mfem_vid;
1132 this->GetElementVertices(elemId, mfem_vid);
1133
1134 // get rotated indices of PUMI element
1135 int pumi_vid_rot[12];
1136 for (int i = 0; i < nv; i++)
1137 {
1138 pumi_vid_rot[i] = mfem_vid.Find(pumi_vid[i]);
1139 }
1140 apf::Downward vs_rot;
1141 for (int i = 0; i < nv; i++)
1142 {
1143 vs_rot[i] = vs[pumi_vid_rot[i]];
1144 }
1145 return findSimplexRotation(apf_mesh, ent, vs_rot);
1146}
1147
1148// Convert parent coordinate form a PUMI tet to an MFEM tet
1150 apf::MeshEntity* tet,
1151 int elemId,
1152 apf::NewArray<apf::Vector3>& pumi_xi,
1153 bool checkOrientation)
1154{
1155 int type = apf_mesh->getType(tet);
1156 MFEM_ASSERT(apf::isSimplex(type),
1157 "only implemented for simplex entity types");
1158 int num_nodes = pumi_xi.size();
1159 IntegrationRule mfem_xi(num_nodes);
1160 int rotation = checkOrientation ? RotationPUMItoMFEM(apf_mesh, tet, elemId):0;
1161 for (int i = 0; i < num_nodes; i++)
1162 {
1163 // for non zero "rotation", rotate the xi
1164 if (rotation)
1165 {
1166 unrotateSimplexXi(pumi_xi[i], apf::Mesh::typeDimension[type], rotation);
1167 }
1168 IntegrationPoint& ip = mfem_xi.IntPoint(i);
1169 double tmp_xi[3];
1170 pumi_xi[i].toArray(tmp_xi);
1171 ip.Set(tmp_xi,3);
1172 }
1173 return mfem_xi;
1174}
1175
1176// Convert parent coordinate from MFEM tet to PUMI tet
1177void ParPumiMesh::ParentXisMFEMtoPUMI(apf::Mesh2* apf_mesh,
1178 int elemId,
1179 apf::MeshEntity* tet,
1180 const IntegrationRule& mfem_xi,
1181 apf::NewArray<apf::Vector3>& pumi_xi,
1182 bool checkOrientation)
1183{
1184 int type = apf_mesh->getType(tet);
1185 MFEM_ASSERT(apf::isSimplex(type),
1186 "only implemented for simplex entity types");
1187 int num_nodes = mfem_xi.Size();
1188 if (!pumi_xi.allocated())
1189 {
1190 pumi_xi.allocate(num_nodes);
1191 }
1192 else
1193 {
1194 pumi_xi.resize(num_nodes);
1195 }
1196
1197 int rotation = checkOrientation ? RotationPUMItoMFEM(apf_mesh, tet, elemId):0;
1198 for (int i = 0; i < num_nodes; i++)
1199 {
1200 IntegrationPoint ip = mfem_xi.IntPoint(i);
1201 pumi_xi[i] = apf::Vector3(ip.x, ip.y, ip.z);
1202
1203 // for non zero "rotation", un-rotate the xi
1204 if (rotation)
1205 {
1206 rotateSimplexXi(pumi_xi[i], apf::Mesh::typeDimension[type], rotation);
1207 }
1208 }
1209}
1210
1211
1212// Transfer a mixed vector-scalar field (i.e. velocity,pressure) and the
1213// magnitude of the vector field to use for mesh adaptation.
1214void ParPumiMesh::FieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
1215 ParGridFunction* grid_vel,
1216 ParGridFunction* grid_pr,
1217 apf::Field* vel_field,
1218 apf::Field* pr_field,
1219 apf::Field* vel_mag_field)
1220{
1221 apf::FieldShape* field_shape = getShape(vel_field);
1222 int dim = apf_mesh->getDimension();
1223
1224 apf::MeshEntity* ent;
1225 apf::MeshIterator* itr = apf_mesh->begin(dim);
1226 int iel = 0;
1227 while ((ent = apf_mesh->iterate(itr)))
1228 {
1229 apf::NewArray<apf::Vector3> pumi_nodes;
1230 apf::getElementNodeXis(field_shape, apf_mesh, ent, pumi_nodes);
1232 apf_mesh, ent, iel, pumi_nodes, true);
1233 // Get the solution
1236 grid_vel->GetVectorValues(*eltr, mfem_nodes, vel);
1237 Vector pr;
1238 grid_pr->GetValues(iel, mfem_nodes, pr, 1);
1239
1240 int non = 0;
1241 for (int d = 0; d <= dim; d++)
1242 {
1243 if (!field_shape->hasNodesIn(d)) { continue; }
1244 apf::Downward a;
1245 int na = apf_mesh->getDownward(ent,d,a);
1246 for (int i = 0; i < na; i++)
1247 {
1248 int type = apf_mesh->getType(a[i]);
1249 int nan = field_shape->countNodesOn(type);
1250 for (int n = 0; n < nan; n++)
1251 {
1252 apf::Vector3 v(vel.GetColumn(non));
1253 apf::setVector(vel_field, a[i], n, v);
1254 apf::setScalar(pr_field, a[i], n, pr[non]);
1255 apf::setScalar(vel_mag_field, a[i], n, v.getLength());
1256 non++;
1257 }
1258 }
1259 }
1260 iel++;
1261 }
1262 apf_mesh->end(itr);
1263}
1264
1265// Transfer a scalar field its magnitude to use for mesh adaptation.
1266void ParPumiMesh::FieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
1267 ParGridFunction* grid_pr,
1268 apf::Field* pr_field,
1269 apf::Field* pr_mag_field)
1270{
1271 apf::FieldShape* field_shape = getShape(pr_field);
1272 int dim = apf_mesh->getDimension();
1273
1274 apf::MeshEntity* ent;
1275 apf::MeshIterator* itr = apf_mesh->begin(dim);
1276 int iel = 0;
1277 while ((ent = apf_mesh->iterate(itr)))
1278 {
1279 apf::NewArray<apf::Vector3> pumi_nodes;
1280 apf::getElementNodeXis(field_shape, apf_mesh, ent, pumi_nodes);
1282 apf_mesh, ent, iel, pumi_nodes, true);
1283 // Get the solution
1284 Vector vals;
1285 grid_pr->GetValues(iel, mfem_nodes, vals, 1);
1286
1287 int non = 0;
1288 for (int d = 0; d <= dim; d++)
1289 {
1290 if (!field_shape->hasNodesIn(d)) { continue; }
1291 apf::Downward a;
1292 int na = apf_mesh->getDownward(ent,d,a);
1293 for (int i = 0; i < na; i++)
1294 {
1295 int type = apf_mesh->getType(a[i]);
1296 int nan = field_shape->countNodesOn(type);
1297 for (int n = 0; n < nan; n++)
1298 {
1299 double pr = vals[non];
1300 double pr_mag = pr >= 0 ? pr : -pr;
1301 apf::setScalar(pr_field, a[i], n, pr);
1302 apf::setScalar(pr_mag_field, a[i], n, pr_mag);
1303 non++;
1304 }
1305 }
1306 }
1307 iel++;
1308 }
1309 apf_mesh->end(itr);
1310}
1311
1312// Transfer a vector field and the magnitude of the vector field to use for mesh
1313// adaptation
1314void ParPumiMesh::VectorFieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
1315 ParGridFunction* grid_vel,
1316 apf::Field* vel_field,
1317 apf::Field* vel_mag_field)
1318{
1319 apf::FieldShape* field_shape = getShape(vel_field);
1320 int dim = apf_mesh->getDimension();
1321
1322 apf::MeshEntity* ent;
1323 apf::MeshIterator* itr = apf_mesh->begin(dim);
1324 int iel = 0;
1325 while ((ent = apf_mesh->iterate(itr)))
1326 {
1327 apf::NewArray<apf::Vector3> pumi_nodes;
1328 apf::getElementNodeXis(field_shape, apf_mesh, ent, pumi_nodes);
1330 apf_mesh, ent, iel, pumi_nodes, true);
1331 // Get the solution
1334 grid_vel->GetVectorValues(*eltr, mfem_nodes, vel);
1335
1336 int non = 0;
1337 for (int d = 0; d <= dim; d++)
1338 {
1339 if (!field_shape->hasNodesIn(d)) { continue; }
1340 apf::Downward a;
1341 int na = apf_mesh->getDownward(ent,d,a);
1342 for (int i = 0; i < na; i++)
1343 {
1344 int type = apf_mesh->getType(a[i]);
1345 int nan = field_shape->countNodesOn(type);
1346 for (int n = 0; n < nan; n++)
1347 {
1348 apf::Vector3 v(vel.GetColumn(non));
1349 apf::setScalar(vel_mag_field, a[i], n, v.getLength());
1350 apf::setVector(vel_field, a[i], n, v);
1351 non++;
1352 }
1353 }
1354 }
1355 iel++;
1356 }
1357 apf_mesh->end(itr);
1358}
1359
1360void ParPumiMesh::NedelecFieldMFEMtoPUMI(apf::Mesh2* apf_mesh,
1361 ParGridFunction* gf,
1362 apf::Field* nedelec_field)
1363{
1364 apf::FieldShape* nedelecFieldShape = nedelec_field->getShape();
1365 int dim = apf_mesh->getDimension();
1366
1367 // loop over all elements
1368 size_t elemNo = 0;
1369 apf::MeshEntity* ent;
1370 apf::MeshIterator* it = apf_mesh->begin(dim);
1371 while ( (ent = apf_mesh->iterate(it)) )
1372 {
1373 // get all the pumi nodes and rotate them
1374 apf::NewArray<apf::Vector3> pumi_nodes;
1375 apf::getElementNodeXis(nedelecFieldShape, apf_mesh, ent, pumi_nodes);
1377 apf_mesh, ent, elemNo, pumi_nodes, true);
1378 // evaluate the vector field on the mfem nodes
1379 ElementTransformation* eltr = this->GetElementTransformation(elemNo);
1380 DenseMatrix mfem_field_vals;
1381 gf->GetVectorValues(*eltr, mfem_nodes, mfem_field_vals);
1382
1383 // compute and store dofs on ND field
1384 int non = 0;
1385 for (int d = 0; d <= dim; d++)
1386 {
1387 if (!nedelecFieldShape->hasNodesIn(d)) { continue; }
1388 apf::Downward a;
1389 int na = apf_mesh->getDownward(ent,d,a);
1390 for (int i = 0; i < na; i++)
1391 {
1392 int type = apf_mesh->getType(a[i]);
1393 int nan = nedelecFieldShape->countNodesOn(type);
1394 apf::MeshElement* me = apf::createMeshElement(apf_mesh, a[i]);
1395 for (int n = 0; n < nan; n++)
1396 {
1397 apf::Vector3 xi, tangent;
1398 nedelecFieldShape->getNodeXi(type, n, xi);
1399 nedelecFieldShape->getNodeTangent(type, n, tangent);
1400 apf::Vector3 pumi_field_vector(mfem_field_vals.GetColumn(non));
1401 apf::Matrix3x3 J;
1402 apf::getJacobian(me, xi, J);
1403 double dof = (J * pumi_field_vector) * tangent;
1404 apf::setScalar(nedelec_field, a[i], n, dof);
1405 non++;
1406 }
1407 apf::destroyMeshElement(me);
1408 }
1409 }
1410 elemNo++;
1411 }
1412 apf_mesh->end(it); // end loop over all elements
1413}
1414
1415void ParPumiMesh::FieldPUMItoMFEM(apf::Mesh2* apf_mesh,
1416 apf::Field* field,
1417 ParGridFunction* grid)
1418{
1419 int nc = apf::countComponents(field);
1420 ParFiniteElementSpace* fes = grid->ParFESpace();
1421 ParMesh* pmesh = fes->GetParMesh();
1422
1423 int dim = apf_mesh->getDimension();
1424
1425 apf::MeshIterator* it = apf_mesh->begin(dim);
1426 for (int i = 0; i < pmesh->GetNE(); i++)
1427 {
1428 const FiniteElement* mfem_elem = fes->GetFE(i);
1429 const IntegrationRule &mfem_xi = mfem_elem->GetNodes();
1430 int non = mfem_xi.Size();
1431 apf::MeshEntity* ent = apf_mesh->iterate(it);
1432 apf::NewArray<apf::Vector3> pumi_xi(non);
1433 ParentXisMFEMtoPUMI(apf_mesh,
1434 i,
1435 ent,
1436 mfem_xi,
1437 pumi_xi,
1438 true);
1439 Array<int> vdofs;
1440 fes->GetElementVDofs(i, vdofs);
1441 apf::MeshElement* me = apf::createMeshElement(apf_mesh, ent);
1442 apf::Element* el = apf::createElement(field, me);
1443 for (int j = 0; j < non; j++)
1444 {
1445 apf::DynamicVector values(nc);
1446 apf::getComponents(el, pumi_xi[j], &values[0]);
1447 // Fill the nodes list
1448 for (int c = 0; c < nc; c++)
1449 {
1450 int dof_loc = j + c * non;
1451 (grid->GetData())[vdofs[dof_loc]] = values[c];
1452 }
1453 }
1454 apf::destroyElement(el);
1455 apf::destroyMeshElement(me);
1456 }
1457 apf_mesh->end(it);
1458}
1459
1460}
1461
1462#endif // MFEM_USE_MPI
1463#endif // MFEM_USE_SCOREC
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition array.hpp:261
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void DeleteAll()
Delete the whole array.
Definition array.hpp:864
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
Definition array.hpp:828
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition array.hpp:874
T * begin()
STL-like begin. Returns pointer to the first element of the array.
Definition array.hpp:302
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void GetColumn(int c, Vector &col) const
Abstract data type element.
Definition element.hpp:29
virtual Element * Duplicate(Mesh *m) const =0
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name.
Definition fe_coll.cpp:126
virtual const char * Name() const
Definition fe_coll.hpp:79
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:287
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:725
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
Abstract class for all finite elements.
Definition fe_base.hpp:239
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
Class for PUMI grid functions.
Definition pumi.hpp:152
GridFunctionPumi(Mesh *m, apf::Mesh2 *PumiM, apf::Numbering *v_num_loc, const int mesh_order)
Construct a GridFunction from a PUMI mesh.
Definition pumi.cpp:887
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition gridfunc.cpp:521
FiniteElementCollection * fec
Used when the grid function is read from a file. It can also be set explicitly, see MakeOwner().
Definition gridfunc.hpp:40
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
Definition gridfunc.hpp:122
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
Definition gridfunc.hpp:34
int VectorDim() const
Definition gridfunc.cpp:323
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
Definition gridfunc.cpp:395
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition gridfunc.cpp:714
void SetComm(MPI_Comm comm)
Set the MPI communicator to 'comm'.
void Create(ListOfIntegerSets &groups, int mpitag)
Set up the group topology given the list of sets of shared entities.
void Copy(GroupTopology &copy) const
Copy the internal data to the external 'copy'.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
A set of integers.
Definition sets.hpp:24
void Recreate(const int n, const int *p)
Create an integer set from C-array 'p' of 'n' integers. Overwrites any existing set data.
Definition sets.cpp:68
Class for integration point with weight.
Definition intrules.hpp:35
void Set(const real_t *p, const int dim)
Definition intrules.hpp:46
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
List of integer sets.
Definition sets.hpp:63
int Size()
Return the number of integer sets in the list.
Definition sets.hpp:70
int Insert(IntegerSet &s)
Check to see if set 's' is in the list. If not append it to the end of the list. Returns the index of...
Definition sets.cpp:91
Mesh data type.
Definition mesh.hpp:56
Array< Vertex > vertices
Definition mesh.hpp:97
int meshgen
Definition mesh.hpp:80
Element * NewElement(int geom)
Definition mesh.cpp:4401
int GetNEdges() const
Return the number of edges.
Definition mesh.hpp:1232
Array< FaceInfo > faces_info
Definition mesh.hpp:224
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
int NumOfBdrElements
Definition mesh.hpp:71
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1438
Array< Element * > faces
Definition mesh.hpp:99
int Dim
Definition mesh.hpp:68
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition mesh.hpp:1235
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition mesh.cpp:3135
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
Table * el_to_face
Definition mesh.hpp:228
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
long sequence
Definition mesh.hpp:87
Table * bel_to_edge
Definition mesh.hpp:232
Table * el_to_edge
Definition mesh.hpp:227
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
Table * edge_vertex
Definition mesh.hpp:239
int NumOfVertices
Definition mesh.hpp:71
Array< int > be_to_face
Definition mesh.hpp:230
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition mesh.hpp:1223
GridFunction * Nodes
Definition mesh.hpp:252
Table * el_to_el
Definition mesh.hpp:229
Table * face_edge
Definition mesh.hpp:238
int NumOfElements
Definition mesh.hpp:71
int NumOfFaces
Definition mesh.hpp:72
int spaceDim
Definition mesh.hpp:69
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition mesh.cpp:3241
int own_nodes
Definition mesh.hpp:253
void FreeElement(Element *E)
Definition mesh.cpp:13056
Array< Element * > boundary
Definition mesh.hpp:98
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition mesh.hpp:291
int NumOfEdges
Definition mesh.hpp:72
Operation last_operation
Definition mesh.hpp:302
Array< Element * > elements
Definition mesh.hpp:92
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
A pair of objects.
Abstract parallel finite element space.
Definition pfespace.hpp:29
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
const FiniteElement * GetFE(int i) const override
Definition pfespace.cpp:534
Class for parallel grid function.
Definition pgridfunc.hpp:33
ParFiniteElementSpace * ParFESpace() const
Class for parallel meshes.
Definition pmesh.hpp:34
Array< Element * > shared_edges
Definition pmesh.hpp:69
Array< int > sface_lface
Definition pmesh.hpp:86
Table group_sedge
Definition pmesh.hpp:78
friend class ParPumiMesh
Definition pmesh.hpp:38
Table group_svert
Shared objects in each group.
Definition pmesh.hpp:77
bool have_face_nbr_data
Definition pmesh.hpp:429
MPI_Comm MyComm
Definition pmesh.hpp:45
Array< Vert4 > shared_quads
Definition pmesh.hpp:74
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition pmesh.cpp:1524
Table group_squad
Definition pmesh.hpp:80
Array< int > svert_lvert
Shared to local index mapping.
Definition pmesh.hpp:83
Array< int > sedge_ledge
Definition pmesh.hpp:84
GroupTopology gtopo
Definition pmesh.hpp:426
void Destroy()
Definition pmesh.cpp:6725
void DeleteFaceNbrData()
Definition pmesh.cpp:1986
void FinalizeParTopo()
Definition pmesh.cpp:898
Array< Vert3 > shared_trias
Definition pmesh.hpp:73
ParNCMesh * pncmesh
Definition pmesh.hpp:439
Table group_stria
Definition pmesh.hpp:79
void FieldPUMItoMFEM(apf::Mesh2 *apf_mesh, apf::Field *field, ParGridFunction *grid)
Transfer a field from PUMI to MFEM after mesh adapt [Scalar and Vector].
Definition pumi.cpp:1415
void VectorFieldMFEMtoPUMI(apf::Mesh2 *apf_mesh, ParGridFunction *grid_vel, apf::Field *vel_field, apf::Field *vel_mag_field)
Transfer field from MFEM mesh to PUMI mesh [Vector].
Definition pumi.cpp:1314
void FieldMFEMtoPUMI(apf::Mesh2 *apf_mesh, ParGridFunction *grid_vel, ParGridFunction *grid_pr, apf::Field *vel_field, apf::Field *pr_field, apf::Field *vel_mag_field)
Transfer field from MFEM mesh to PUMI mesh [Mixed].
Definition pumi.cpp:1214
void ParentXisMFEMtoPUMI(apf::Mesh2 *apf_mesh, int elemId, apf::MeshEntity *tet, const IntegrationRule &mfem_xi, apf::NewArray< apf::Vector3 > &pumi_xi, bool checkOrientation=true)
Convert the parent coordinate from MFEM to PUMI.
Definition pumi.cpp:1177
int RotationPUMItoMFEM(apf::Mesh2 *apf_mesh, apf::MeshEntity *tet, int elemId)
Returns the PUMI-to-MFEM permutation (aka rotation, aka orientation)
Definition pumi.cpp:1113
void UpdateMesh(const ParMesh *AdaptedpMesh)
Update the mesh after adaptation.
Definition pumi.cpp:956
IntegrationRule ParentXisPUMItoMFEM(apf::Mesh2 *apf_mesh, apf::MeshEntity *tet, int elemId, apf::NewArray< apf::Vector3 > &pumi_xi, bool checkOrientation=true)
Convert the parent coordinate from PUMI to MFEM.
Definition pumi.cpp:1149
void NedelecFieldMFEMtoPUMI(apf::Mesh2 *apf_mesh, ParGridFunction *gf, apf::Field *nedelec_field)
Transfer Nedelec field from MFEM mesh to PUMI mesh [Vector].
Definition pumi.cpp:1360
void ReadSCORECMesh(apf::Mesh2 *apf_mesh, apf::Numbering *v_num_loc, const int curved)
Definition pumi.cpp:342
PumiMesh(apf::Mesh2 *apf_mesh, int generate_edges=0, int refine=1, bool fix_orientation=true)
Generate an MFEM mesh from a PUMI mesh.
Definition pumi.cpp:245
void CountBoundaryEntity(apf::Mesh2 *apf_mesh, const int BcDim, int &NumBC)
Definition pumi.cpp:253
void Load(apf::Mesh2 *apf_mesh, int generate_edges=0, int refine=1, bool fix_orientation=true)
Load a PUMI mesh (following the steps in the MFEM Load function).
Definition pumi.cpp:276
Data type line segment element.
Definition segment.hpp:23
void ShiftUpI()
Definition table.cpp:115
void Clear()
Definition table.cpp:381
void AddConnection(int r, int c)
Definition table.hpp:80
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition table.cpp:81
void Copy(Table &copy) const
Definition table.cpp:390
void MakeJ()
Definition table.cpp:91
void AddAColumnInRow(int r)
Definition table.hpp:77
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
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
void vel(const Vector &x, real_t t, Vector &u)