MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pmesh.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 "../config/config.hpp"
13
14#ifdef MFEM_USE_MPI
15
16#include "mesh_headers.hpp"
17#include "../fem/fem.hpp"
18#include "../general/sets.hpp"
20#include "../general/text.hpp"
22
23#include <iostream>
24#include <fstream>
25
26using namespace std;
27
28namespace mfem
29{
30
31ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
32 : Mesh(pmesh, false),
33 group_svert(pmesh.group_svert),
34 group_sedge(pmesh.group_sedge),
35 group_stria(pmesh.group_stria),
36 group_squad(pmesh.group_squad),
37 glob_elem_offset(-1),
38 glob_offset_sequence(-1),
39 gtopo(pmesh.gtopo)
40{
41 MyComm = pmesh.MyComm;
42 NRanks = pmesh.NRanks;
43 MyRank = pmesh.MyRank;
44
45 // Duplicate the shared_edges
46 shared_edges.SetSize(pmesh.shared_edges.Size());
47 for (int i = 0; i < shared_edges.Size(); i++)
48 {
49 shared_edges[i] = pmesh.shared_edges[i]->Duplicate(this);
50 }
51
54
55 // Copy the shared-to-local index Arrays
59
60 // Do not copy face-neighbor data (can be generated if needed)
61 have_face_nbr_data = false;
62
63 // If pmesh has a ParNURBSExtension, it was copied by the Mesh copy ctor, so
64 // there is no need to do anything here.
65
66 // Copy ParNCMesh, if present
67 if (pmesh.pncmesh)
68 {
69 pncmesh = new ParNCMesh(*pmesh.pncmesh);
71 }
72 else
73 {
74 pncmesh = NULL;
75 }
77
78 // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
79 // and the FiniteElementSpace (as a ParFiniteElementSpace)
80 if (pmesh.Nodes && copy_nodes)
81 {
82 FiniteElementSpace *fes = pmesh.Nodes->FESpace();
83 const FiniteElementCollection *fec = fes->FEColl();
84 FiniteElementCollection *fec_copy =
86 ParFiniteElementSpace *pfes_copy =
87 new ParFiniteElementSpace(*fes, *this, fec_copy);
88 Nodes = new ParGridFunction(pfes_copy);
89 Nodes->MakeOwner(fec_copy);
90 *Nodes = *pmesh.Nodes;
91 own_nodes = 1;
92 }
93}
94
96{
97 Swap(mesh);
98}
99
101{
102 Swap(mesh);
103 return *this;
104}
105
106ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
107 int part_method)
108 : glob_elem_offset(-1)
109 , glob_offset_sequence(-1)
110 , gtopo(comm)
111{
112 int *partitioning = NULL;
113 Array<bool> activeBdrElem;
114
115 MyComm = comm;
116 MPI_Comm_size(MyComm, &NRanks);
117 MPI_Comm_rank(MyComm, &MyRank);
118
119 if (mesh.Nonconforming())
120 {
121 if (partitioning_)
122 {
123 partitioning = partitioning_;
124 }
125 ncmesh = pncmesh = new ParNCMesh(comm, *mesh.ncmesh, partitioning);
126 if (!partitioning)
127 {
128 partitioning = new int[mesh.GetNE()];
129 for (int i = 0; i < mesh.GetNE(); i++)
130 {
131 partitioning[i] = pncmesh->InitialPartition(i);
132 }
133 }
134
135 pncmesh->Prune();
136
138 pncmesh->OnMeshUpdated(this);
139
141
142 // SetMeshGen(); // called by Mesh::InitFromNCMesh(...) above
143 meshgen = mesh.meshgen; // copy the global 'meshgen'
144
147
148 // Copy attribute and bdr_attribute names
151
153 }
154 else // mesh.Conforming()
155 {
156 Dim = mesh.Dim;
157 spaceDim = mesh.spaceDim;
158
159 ncmesh = pncmesh = NULL;
160
161 if (partitioning_)
162 {
163 partitioning = partitioning_;
164 }
165 else
166 {
167 partitioning = mesh.GeneratePartitioning(NRanks, part_method);
168 }
169
170 // re-enumerate the partitions to better map to actual processor
171 // interconnect topology !?
172
173 Array<int> vert_global_local;
174 NumOfVertices = BuildLocalVertices(mesh, partitioning, vert_global_local);
175 NumOfElements = BuildLocalElements(mesh, partitioning, vert_global_local);
176
177 Table *edge_element = NULL;
178 NumOfBdrElements = BuildLocalBoundary(mesh, partitioning,
179 vert_global_local,
180 activeBdrElem, edge_element);
181
182 SetMeshGen();
183 meshgen = mesh.meshgen; // copy the global 'meshgen'
184
187
188 // Copy attribute and bdr_attribute names
191
193
194 if (Dim > 1)
195 {
196 el_to_edge = new Table;
198 }
199
200 STable3D *faces_tbl = NULL;
201 if (Dim == 3)
202 {
203 faces_tbl = GetElementToFaceTable(1);
204 }
205
207
208 // Make sure the be_to_face array is initialized.
209 // In 2D, it will be set in the above call to Mesh::GetElementToEdgeTable.
210 // In 3D, it will be set in GetElementToFaceTable.
211 // In 1D, we need to set it manually.
212 if (Dim == 1)
213 {
215 for (int i = 0; i < NumOfBdrElements; ++i)
216 {
217 be_to_face[i] = boundary[i]->GetVertices()[0];
218 }
219 }
220
221 ListOfIntegerSets groups;
222 {
223 // the first group is the local one
224 IntegerSet group;
225 group.Recreate(1, &MyRank);
226 groups.Insert(group);
227 }
228
229 MFEM_ASSERT(mesh.GetNFaces() == 0 || Dim >= 3, "");
230
231 Array<int> face_group(mesh.GetNFaces());
232 Table *vert_element = mesh.GetVertexToElementTable(); // we must delete this
233
234 FindSharedFaces(mesh, partitioning, face_group, groups);
235 int nsedges = FindSharedEdges(mesh, partitioning, edge_element, groups);
236 int nsvert = FindSharedVertices(partitioning, vert_element, groups);
237
238 // build the group communication topology
239 gtopo.Create(groups, 822);
240
241 // fill out group_sface, group_sedge, group_svert
242 int ngroups = groups.Size()-1, nstris, nsquads;
243 BuildFaceGroup(ngroups, mesh, face_group, nstris, nsquads);
244 BuildEdgeGroup(ngroups, *edge_element);
245 BuildVertexGroup(ngroups, *vert_element);
246
247 // build shared_faces and sface_lface mapping
248 BuildSharedFaceElems(nstris, nsquads, mesh, partitioning, faces_tbl,
249 face_group, vert_global_local);
250 delete faces_tbl;
251
252 // build shared_edges and sedge_ledge mapping
253 BuildSharedEdgeElems(nsedges, mesh, vert_global_local, edge_element);
254 delete edge_element;
255
256 // build svert_lvert mapping
257 BuildSharedVertMapping(nsvert, vert_element, vert_global_local);
258 delete vert_element;
259 }
260
261 if (mesh.NURBSext)
262 {
263 MFEM_ASSERT(mesh.GetNodes() &&
264 mesh.GetNodes()->FESpace()->GetNURBSext() == mesh.NURBSext,
265 "invalid NURBS mesh");
266 NURBSext = new ParNURBSExtension(comm, mesh.NURBSext, partitioning,
267 activeBdrElem);
268 }
269
270 if (mesh.GetNodes()) // curved mesh
271 {
272 if (!NURBSext)
273 {
274 Nodes = new ParGridFunction(this, mesh.GetNodes());
275 }
276 else
277 {
278 const FiniteElementSpace *glob_fes = mesh.GetNodes()->FESpace();
282 new ParFiniteElementSpace(this, nfec, glob_fes->GetVDim(),
283 glob_fes->GetOrdering());
284 Nodes = new ParGridFunction(pfes);
285 Nodes->MakeOwner(nfec); // Nodes will own nfec and pfes
286 }
287 own_nodes = 1;
288
289 Array<int> gvdofs, lvdofs;
290 Vector lnodes;
291 int element_counter = 0;
292 for (int i = 0; i < mesh.GetNE(); i++)
293 {
294 if (partitioning[i] == MyRank)
295 {
296 Nodes->FESpace()->GetElementVDofs(element_counter, lvdofs);
297 mesh.GetNodes()->FESpace()->GetElementVDofs(i, gvdofs);
298 mesh.GetNodes()->GetSubVector(gvdofs, lnodes);
299 Nodes->SetSubVector(lvdofs, lnodes);
300 element_counter++;
301 }
302 }
303
304 // set meaningful values to 'vertices' even though we have Nodes,
305 // for compatibility (e.g., Mesh::GetVertex())
307 }
308
309 if (partitioning != partitioning_)
310 {
311 delete [] partitioning;
312 }
313
314 have_face_nbr_data = false;
315}
316
317
319 const int* partitioning,
320 Array<int> &vert_global_local)
321{
322 vert_global_local.SetSize(mesh.GetNV());
323 vert_global_local = -1;
324
325 int vert_counter = 0;
326 for (int i = 0; i < mesh.GetNE(); i++)
327 {
328 if (partitioning[i] == MyRank)
329 {
330 Array<int> vert;
331 mesh.GetElementVertices(i, vert);
332 for (int j = 0; j < vert.Size(); j++)
333 {
334 if (vert_global_local[vert[j]] < 0)
335 {
336 vert_global_local[vert[j]] = vert_counter++;
337 }
338 }
339 }
340 }
341
342 // re-enumerate the local vertices to preserve the global ordering
343 vert_counter = 0;
344 for (int i = 0; i < vert_global_local.Size(); i++)
345 {
346 if (vert_global_local[i] >= 0)
347 {
348 vert_global_local[i] = vert_counter++;
349 }
350 }
351
352 vertices.SetSize(vert_counter);
353
354 for (int i = 0; i < vert_global_local.Size(); i++)
355 {
356 if (vert_global_local[i] >= 0)
357 {
358 vertices[vert_global_local[i]].SetCoords(mesh.SpaceDimension(),
359 mesh.GetVertex(i));
360 }
361 }
362
363 return vert_counter;
364}
365
366int ParMesh::BuildLocalElements(const Mesh& mesh, const int* partitioning,
367 const Array<int>& vert_global_local)
368{
369 const int nelems = std::count_if(partitioning,
370 partitioning + mesh.GetNE(), [this](int i) { return i == MyRank;});
371
372 elements.SetSize(nelems);
373
374 // Determine elements, enumerating the local elements to preserve the global
375 // order. This is used, e.g. by the ParGridFunction ctor that takes a global
376 // GridFunction as input parameter.
377 int element_counter = 0;
378 for (int i = 0; i < mesh.GetNE(); i++)
379 {
380 if (partitioning[i] == MyRank)
381 {
382 elements[element_counter] = mesh.GetElement(i)->Duplicate(this);
383 int *v = elements[element_counter]->GetVertices();
384 int nv = elements[element_counter]->GetNVertices();
385 for (int j = 0; j < nv; j++)
386 {
387 v[j] = vert_global_local[v[j]];
388 }
389 ++element_counter;
390 }
391 }
392
393 return element_counter;
394}
395
396int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
397 const Array<int>& vert_global_local,
398 Array<bool>& activeBdrElem,
399 Table*& edge_element)
400{
401 int nbdry = 0;
402 if (mesh.NURBSext)
403 {
404 activeBdrElem.SetSize(mesh.GetNBE());
405 activeBdrElem = false;
406 }
407 // build boundary elements
408 if (Dim == 3)
409 {
410 for (int i = 0; i < mesh.GetNBE(); i++)
411 {
412 int face, o, el1, el2;
413 mesh.GetBdrElementFace(i, &face, &o);
414 mesh.GetFaceElements(face, &el1, &el2);
415 if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
416 {
417 nbdry++;
418 if (mesh.NURBSext)
419 {
420 activeBdrElem[i] = true;
421 }
422 }
423 }
424
425 int bdrelem_counter = 0;
426 boundary.SetSize(nbdry);
427 for (int i = 0; i < mesh.GetNBE(); i++)
428 {
429 int face, o, el1, el2;
430 mesh.GetBdrElementFace(i, &face, &o);
431 mesh.GetFaceElements(face, &el1, &el2);
432 if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] == MyRank)
433 {
434 boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
435 int *v = boundary[bdrelem_counter]->GetVertices();
436 int nv = boundary[bdrelem_counter]->GetNVertices();
437 for (int j = 0; j < nv; j++)
438 {
439 v[j] = vert_global_local[v[j]];
440 }
441 bdrelem_counter++;
442 }
443 }
444 }
445 else if (Dim == 2)
446 {
447 edge_element = new Table;
448 Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
449
450 for (int i = 0; i < mesh.GetNBE(); i++)
451 {
452 int edge = mesh.GetBdrElementFaceIndex(i);
453 int el1 = edge_element->GetRow(edge)[0];
454 if (partitioning[el1] == MyRank)
455 {
456 nbdry++;
457 if (mesh.NURBSext)
458 {
459 activeBdrElem[i] = true;
460 }
461 }
462 }
463
464 int bdrelem_counter = 0;
465 boundary.SetSize(nbdry);
466 for (int i = 0; i < mesh.GetNBE(); i++)
467 {
468 int edge = mesh.GetBdrElementFaceIndex(i);
469 int el1 = edge_element->GetRow(edge)[0];
470 if (partitioning[el1] == MyRank)
471 {
472 boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
473 int *v = boundary[bdrelem_counter]->GetVertices();
474 int nv = boundary[bdrelem_counter]->GetNVertices();
475 for (int j = 0; j < nv; j++)
476 {
477 v[j] = vert_global_local[v[j]];
478 }
479 bdrelem_counter++;
480 }
481 }
482 }
483 else if (Dim == 1)
484 {
485 for (int i = 0; i < mesh.GetNBE(); i++)
486 {
487 int vert = mesh.boundary[i]->GetVertices()[0];
488 int el1, el2;
489 mesh.GetFaceElements(vert, &el1, &el2);
490 if (partitioning[el1] == MyRank)
491 {
492 nbdry++;
493 }
494 }
495
496 int bdrelem_counter = 0;
497 boundary.SetSize(nbdry);
498 for (int i = 0; i < mesh.GetNBE(); i++)
499 {
500 int vert = mesh.boundary[i]->GetVertices()[0];
501 int el1, el2;
502 mesh.GetFaceElements(vert, &el1, &el2);
503 if (partitioning[el1] == MyRank)
504 {
505 boundary[bdrelem_counter] = mesh.GetBdrElement(i)->Duplicate(this);
506 int *v = boundary[bdrelem_counter]->GetVertices();
507 v[0] = vert_global_local[v[0]];
508 bdrelem_counter++;
509 }
510 }
511 }
512
513 return nbdry;
514}
515
516void ParMesh::FindSharedFaces(const Mesh &mesh, const int *partitioning,
517 Array<int> &face_group,
518 ListOfIntegerSets &groups)
519{
520 IntegerSet group;
521
522 // determine shared faces
523 face_group.SetSize(mesh.GetNFaces());
524 for (int i = 0; i < face_group.Size(); i++)
525 {
526 int el[2];
527 face_group[i] = -1;
528 mesh.GetFaceElements(i, &el[0], &el[1]);
529 if (el[1] >= 0)
530 {
531 el[0] = partitioning[el[0]];
532 el[1] = partitioning[el[1]];
533 if ((el[0] == MyRank && el[1] != MyRank) ||
534 (el[0] != MyRank && el[1] == MyRank))
535 {
536 group.Recreate(2, el);
537 face_group[i] = groups.Insert(group) - 1;
538 }
539 }
540 }
541}
542
543int ParMesh::FindSharedEdges(const Mesh &mesh, const int *partitioning,
544 Table*& edge_element,
545 ListOfIntegerSets& groups)
546{
547 IntegerSet group;
548
549 // determine shared edges
550 int sedge_counter = 0;
551 if (!edge_element)
552 {
553 edge_element = new Table;
554 if (Dim == 1)
555 {
556 edge_element->SetDims(0,0);
557 }
558 else
559 {
560 Transpose(mesh.ElementToEdgeTable(), *edge_element, mesh.GetNEdges());
561 }
562 }
563
564 for (int i = 0; i < edge_element->Size(); i++)
565 {
566 int me = 0, others = 0;
567 for (int j = edge_element->GetI()[i]; j < edge_element->GetI()[i+1]; j++)
568 {
569 int k = edge_element->GetJ()[j];
570 int rank = partitioning[k];
571 edge_element->GetJ()[j] = rank;
572 if (rank == MyRank)
573 {
574 me = 1;
575 }
576 else
577 {
578 others = 1;
579 }
580 }
581
582 if (me && others)
583 {
584 sedge_counter++;
585 group.Recreate(edge_element->RowSize(i), edge_element->GetRow(i));
586 edge_element->GetRow(i)[0] = groups.Insert(group) - 1;
587 }
588 else
589 {
590 edge_element->GetRow(i)[0] = -1;
591 }
592 }
593
594 return sedge_counter;
595}
596
597int ParMesh::FindSharedVertices(const int *partitioning, Table *vert_element,
598 ListOfIntegerSets &groups)
599{
600 IntegerSet group;
601
602 // determine shared vertices
603 int svert_counter = 0;
604 for (int i = 0; i < vert_element->Size(); i++)
605 {
606 int me = 0, others = 0;
607 for (int j = vert_element->GetI()[i]; j < vert_element->GetI()[i+1]; j++)
608 {
609 vert_element->GetJ()[j] = partitioning[vert_element->GetJ()[j]];
610 if (vert_element->GetJ()[j] == MyRank)
611 {
612 me = 1;
613 }
614 else
615 {
616 others = 1;
617 }
618 }
619
620 if (me && others)
621 {
622 svert_counter++;
623 group.Recreate(vert_element->RowSize(i), vert_element->GetRow(i));
624 vert_element->GetI()[i] = groups.Insert(group) - 1;
625 }
626 else
627 {
628 vert_element->GetI()[i] = -1;
629 }
630 }
631 return svert_counter;
632}
633
634void ParMesh::BuildFaceGroup(int ngroups, const Mesh &mesh,
635 const Array<int> &face_group,
636 int &nstria, int &nsquad)
637{
638 // build group_stria and group_squad
639 group_stria.MakeI(ngroups);
640 group_squad.MakeI(ngroups);
641
642 for (int i = 0; i < face_group.Size(); i++)
643 {
644 if (face_group[i] >= 0)
645 {
646 if (mesh.GetFace(i)->GetType() == Element::TRIANGLE)
647 {
648 group_stria.AddAColumnInRow(face_group[i]);
649 }
650 else
651 {
652 group_squad.AddAColumnInRow(face_group[i]);
653 }
654 }
655 }
656
659
660 nstria = nsquad = 0;
661 for (int i = 0; i < face_group.Size(); i++)
662 {
663 if (face_group[i] >= 0)
664 {
665 if (mesh.GetFace(i)->GetType() == Element::TRIANGLE)
666 {
667 group_stria.AddConnection(face_group[i], nstria++);
668 }
669 else
670 {
671 group_squad.AddConnection(face_group[i], nsquad++);
672 }
673 }
674 }
675
678}
679
680void ParMesh::BuildEdgeGroup(int ngroups, const Table &edge_element)
681{
682 group_sedge.MakeI(ngroups);
683
684 for (int i = 0; i < edge_element.Size(); i++)
685 {
686 if (edge_element.GetRow(i)[0] >= 0)
687 {
688 group_sedge.AddAColumnInRow(edge_element.GetRow(i)[0]);
689 }
690 }
691
693
694 int sedge_counter = 0;
695 for (int i = 0; i < edge_element.Size(); i++)
696 {
697 if (edge_element.GetRow(i)[0] >= 0)
698 {
699 group_sedge.AddConnection(edge_element.GetRow(i)[0], sedge_counter++);
700 }
701 }
702
704}
705
706void ParMesh::BuildVertexGroup(int ngroups, const Table &vert_element)
707{
708 group_svert.MakeI(ngroups);
709
710 for (int i = 0; i < vert_element.Size(); i++)
711 {
712 if (vert_element.GetI()[i] >= 0)
713 {
714 group_svert.AddAColumnInRow(vert_element.GetI()[i]);
715 }
716 }
717
719
720 int svert_counter = 0;
721 for (int i = 0; i < vert_element.Size(); i++)
722 {
723 if (vert_element.GetI()[i] >= 0)
724 {
725 group_svert.AddConnection(vert_element.GetI()[i], svert_counter++);
726 }
727 }
728
730}
731
732void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
733 const Mesh& mesh, int *partitioning,
734 const STable3D *faces_tbl,
735 const Array<int> &face_group,
736 const Array<int> &vert_global_local)
737{
738 shared_trias.SetSize(ntri_faces);
739 shared_quads.SetSize(nquad_faces);
740 sface_lface. SetSize(ntri_faces + nquad_faces);
741
742 if (Dim < 3) { return; }
743
744 int stria_counter = 0;
745 int squad_counter = 0;
746 for (int i = 0; i < face_group.Size(); i++)
747 {
748 if (face_group[i] < 0) { continue; }
749
750 const Element *face = mesh.GetFace(i);
751 const int *fv = face->GetVertices();
752 switch (face->GetType())
753 {
755 {
756 shared_trias[stria_counter].Set(fv);
757 int *v = shared_trias[stria_counter].v;
758 for (int j = 0; j < 3; j++)
759 {
760 v[j] = vert_global_local[v[j]];
761 }
762 const int lface = (*faces_tbl)(v[0], v[1], v[2]);
763 sface_lface[stria_counter] = lface;
764 if (meshgen == 1) // Tet-only mesh
765 {
766 Tetrahedron *tet = dynamic_cast<Tetrahedron *>
767 (elements[faces_info[lface].Elem1No]);
768 // mark the shared face for refinement by reorienting
769 // it according to the refinement flag in the tetrahedron
770 // to which this shared face belongs to.
771 if (tet->GetRefinementFlag())
772 {
773 tet->GetMarkedFace(faces_info[lface].Elem1Inf/64, v);
774 // flip the shared face in the processor that owns the
775 // second element (in 'mesh')
776 int gl_el1, gl_el2;
777 mesh.GetFaceElements(i, &gl_el1, &gl_el2);
778 if (MyRank == partitioning[gl_el2])
779 {
780 std::swap(v[0], v[1]);
781 }
782 }
783 }
784 stria_counter++;
785 break;
786 }
787
789 {
790 shared_quads[squad_counter].Set(fv);
791 int *v = shared_quads[squad_counter].v;
792 for (int j = 0; j < 4; j++)
793 {
794 v[j] = vert_global_local[v[j]];
795 }
796 sface_lface[shared_trias.Size() + squad_counter] =
797 (*faces_tbl)(v[0], v[1], v[2], v[3]);
798 squad_counter++;
799 break;
800 }
801
802 default:
803 MFEM_ABORT("unknown face type: " << face->GetType());
804 break;
805 }
806 }
807}
808
809void ParMesh::BuildSharedEdgeElems(int nedges, Mesh& mesh,
810 const Array<int>& vert_global_local,
811 const Table* edge_element)
812{
813 // The passed in mesh is still the global mesh. "this" mesh is the
814 // local partitioned mesh.
815
816 shared_edges.SetSize(nedges);
817 sedge_ledge. SetSize(nedges);
818
819 {
820 DSTable v_to_v(NumOfVertices);
822
823 int sedge_counter = 0;
824 for (int i = 0; i < edge_element->Size(); i++)
825 {
826 if (edge_element->GetRow(i)[0] >= 0)
827 {
828 Array<int> vert;
829 mesh.GetEdgeVertices(i, vert);
830
831 shared_edges[sedge_counter] =
832 new Segment(vert_global_local[vert[0]],
833 vert_global_local[vert[1]], 1);
834
835 sedge_ledge[sedge_counter] = v_to_v(vert_global_local[vert[0]],
836 vert_global_local[vert[1]]);
837
838 MFEM_VERIFY(sedge_ledge[sedge_counter] >= 0, "Error in v_to_v.");
839
840 sedge_counter++;
841 }
842 }
843 }
844}
845
847 const mfem::Table *vert_element,
848 const Array<int> &vert_global_local)
849{
850 // build svert_lvert
851 svert_lvert.SetSize(nvert);
852
853 int svert_counter = 0;
854 for (int i = 0; i < vert_element->Size(); i++)
855 {
856 if (vert_element->GetI()[i] >= 0)
857 {
858 svert_lvert[svert_counter++] = vert_global_local[i];
859 }
860 }
861}
862
863
864// protected method, used by Nonconforming(De)Refinement and Rebalance
866 : MyComm(pncmesh.MyComm)
867 , NRanks(pncmesh.NRanks)
868 , MyRank(pncmesh.MyRank)
869 , glob_elem_offset(-1)
870 , glob_offset_sequence(-1)
871 , gtopo(MyComm)
872 , pncmesh(NULL)
873{
876 have_face_nbr_data = false;
877}
878
880{
881 if (glob_offset_sequence != sequence) // mesh has changed
882 {
883 long long local_elems = NumOfElements;
884 MPI_Scan(&local_elems, &glob_elem_offset, 1, MPI_LONG_LONG, MPI_SUM,
885 MyComm);
886 glob_elem_offset -= local_elems;
887
888 glob_offset_sequence = sequence; // don't recalculate until refinement etc.
889 }
890}
891
893{
894 int loc_meshgen = meshgen;
895 MPI_Allreduce(&loc_meshgen, &meshgen, 1, MPI_INT, MPI_BOR, MyComm);
896}
897
899{
900 // Determine sedge_ledge
902 if (shared_edges.Size())
903 {
904 DSTable v_to_v(NumOfVertices);
906 for (int se = 0; se < shared_edges.Size(); se++)
907 {
908 const int *v = shared_edges[se]->GetVertices();
909 const int l_edge = v_to_v(v[0], v[1]);
910 MFEM_ASSERT(l_edge >= 0, "invalid shared edge");
911 sedge_ledge[se] = l_edge;
912 }
913 }
914
915 // Determine sface_lface
916 const int nst = shared_trias.Size();
917 sface_lface.SetSize(nst + shared_quads.Size());
918 if (sface_lface.Size())
919 {
920 auto faces_tbl = std::unique_ptr<STable3D>(GetFacesTable());
921 for (int st = 0; st < nst; st++)
922 {
923 const int *v = shared_trias[st].v;
924 sface_lface[st] = (*faces_tbl)(v[0], v[1], v[2]);
925 }
926 for (int sq = 0; sq < shared_quads.Size(); sq++)
927 {
928 const int *v = shared_quads[sq].v;
929 sface_lface[nst+sq] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
930 }
931 }
932}
933
934ParMesh::ParMesh(MPI_Comm comm, istream &input, bool refine, int generate_edges,
935 bool fix_orientation)
936 : glob_elem_offset(-1)
937 , glob_offset_sequence(-1)
938 , gtopo(comm)
939{
940 MyComm = comm;
941 MPI_Comm_size(MyComm, &NRanks);
942 MPI_Comm_rank(MyComm, &MyRank);
943
944 have_face_nbr_data = false;
945 pncmesh = NULL;
946
947 Load(input, generate_edges, refine, fix_orientation);
948}
949
950void ParMesh::Load(istream &input, int generate_edges, int refine,
951 bool fix_orientation)
952{
954
955 // Tell Loader() to read up to 'mfem_serial_mesh_end' instead of
956 // 'mfem_mesh_end', as we have additional parallel mesh data to load in from
957 // the stream.
958 Loader(input, generate_edges, "mfem_serial_mesh_end");
959
960 ReduceMeshGen(); // determine the global 'meshgen'
961
962 if (Conforming())
963 {
964 LoadSharedEntities(input);
965 }
966 else
967 {
968 // the ParNCMesh instance was already constructed in 'Loader'
969 pncmesh = dynamic_cast<ParNCMesh*>(ncmesh);
970 MFEM_ASSERT(pncmesh, "internal error");
971
972 // in the NC case we don't need to load extra data from the file,
973 // as the shared entities can be constructed from the ghost layer
975 }
976
977 Finalize(refine, fix_orientation);
978
980
981 // note: attributes and bdr_attributes are local lists
982
983 // TODO: NURBS meshes?
984}
985
986void ParMesh::LoadSharedEntities(istream &input)
987{
988 string ident;
989 skip_comment_lines(input, '#');
990
991 // read the group topology
992 input >> ident;
993 MFEM_VERIFY(ident == "communication_groups",
994 "input stream is not a parallel MFEM mesh");
995 gtopo.Load(input);
996
997 skip_comment_lines(input, '#');
998
999 // read and set the sizes of svert_lvert, group_svert
1000 {
1001 int num_sverts;
1002 input >> ident >> num_sverts;
1003 MFEM_VERIFY(ident == "total_shared_vertices", "invalid mesh file");
1004 svert_lvert.SetSize(num_sverts);
1005 group_svert.SetDims(GetNGroups()-1, num_sverts);
1006 }
1007 // read and set the sizes of sedge_ledge, group_sedge
1008 if (Dim >= 2)
1009 {
1010 skip_comment_lines(input, '#');
1011 int num_sedges;
1012 input >> ident >> num_sedges;
1013 MFEM_VERIFY(ident == "total_shared_edges", "invalid mesh file");
1014 sedge_ledge.SetSize(num_sedges);
1015 shared_edges.SetSize(num_sedges);
1016 group_sedge.SetDims(GetNGroups()-1, num_sedges);
1017 }
1018 else
1019 {
1020 group_sedge.SetSize(GetNGroups()-1, 0); // create empty group_sedge
1021 }
1022 // read and set the sizes of sface_lface, group_{stria,squad}
1023 if (Dim >= 3)
1024 {
1025 skip_comment_lines(input, '#');
1026 int num_sface;
1027 input >> ident >> num_sface;
1028 MFEM_VERIFY(ident == "total_shared_faces", "invalid mesh file");
1029 sface_lface.SetSize(num_sface);
1032 }
1033 else
1034 {
1035 group_stria.SetSize(GetNGroups()-1, 0); // create empty group_stria
1036 group_squad.SetSize(GetNGroups()-1, 0); // create empty group_squad
1037 }
1038
1039 // read, group by group, the contents of group_svert, svert_lvert,
1040 // group_sedge, shared_edges, group_{stria,squad}, shared_{trias,quads}
1041 int svert_counter = 0, sedge_counter = 0;
1042 for (int gr = 1; gr < GetNGroups(); gr++)
1043 {
1044 skip_comment_lines(input, '#');
1045#if 0
1046 // implementation prior to prism-dev merge
1047 int g;
1048 input >> ident >> g; // group
1049 if (g != gr)
1050 {
1051 mfem::err << "ParMesh::ParMesh : expecting group " << gr
1052 << ", read group " << g << endl;
1053 mfem_error();
1054 }
1055#endif
1056 {
1057 int nv;
1058 input >> ident >> nv; // shared_vertices (in this group)
1059 MFEM_VERIFY(ident == "shared_vertices", "invalid mesh file");
1060 nv += svert_counter;
1061 MFEM_VERIFY(nv <= group_svert.Size_of_connections(),
1062 "incorrect number of total_shared_vertices");
1063 group_svert.GetI()[gr] = nv;
1064 for ( ; svert_counter < nv; svert_counter++)
1065 {
1066 group_svert.GetJ()[svert_counter] = svert_counter;
1067 input >> svert_lvert[svert_counter];
1068 }
1069 }
1070 if (Dim >= 2)
1071 {
1072 int ne, v[2];
1073 input >> ident >> ne; // shared_edges (in this group)
1074 MFEM_VERIFY(ident == "shared_edges", "invalid mesh file");
1075 ne += sedge_counter;
1076 MFEM_VERIFY(ne <= group_sedge.Size_of_connections(),
1077 "incorrect number of total_shared_edges");
1078 group_sedge.GetI()[gr] = ne;
1079 for ( ; sedge_counter < ne; sedge_counter++)
1080 {
1081 group_sedge.GetJ()[sedge_counter] = sedge_counter;
1082 input >> v[0] >> v[1];
1083 shared_edges[sedge_counter] = new Segment(v[0], v[1], 1);
1084 }
1085 }
1086 if (Dim >= 3)
1087 {
1088 int nf, tstart = shared_trias.Size(), qstart = shared_quads.Size();
1089 input >> ident >> nf; // shared_faces (in this group)
1090 for (int i = 0; i < nf; i++)
1091 {
1092 int geom, *v;
1093 input >> geom;
1094 switch (geom)
1095 {
1096 case Geometry::TRIANGLE:
1097 shared_trias.SetSize(shared_trias.Size()+1);
1098 v = shared_trias.Last().v;
1099 for (int ii = 0; ii < 3; ii++) { input >> v[ii]; }
1100 break;
1101 case Geometry::SQUARE:
1102 shared_quads.SetSize(shared_quads.Size()+1);
1103 v = shared_quads.Last().v;
1104 for (int ii = 0; ii < 4; ii++) { input >> v[ii]; }
1105 break;
1106 default:
1107 MFEM_ABORT("invalid shared face geometry: " << geom);
1108 }
1109 }
1110 group_stria.AddColumnsInRow(gr-1, shared_trias.Size()-tstart);
1111 group_squad.AddColumnsInRow(gr-1, shared_quads.Size()-qstart);
1112 }
1113 }
1114 if (Dim >= 3)
1115 {
1116 MFEM_VERIFY(shared_trias.Size() + shared_quads.Size()
1117 == sface_lface.Size(),
1118 "incorrect number of total_shared_faces");
1119 // Define the J arrays of group_stria and group_squad -- they just contain
1120 // consecutive numbers starting from 0 up to shared_trias.Size()-1 and
1121 // shared_quads.Size()-1, respectively.
1123 for (int i = 0; i < shared_trias.Size(); i++)
1124 {
1125 group_stria.GetJ()[i] = i;
1126 }
1128 for (int i = 0; i < shared_quads.Size(); i++)
1129 {
1130 group_squad.GetJ()[i] = i;
1131 }
1132 }
1133}
1134
1135ParMesh::ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type)
1136{
1137 MakeRefined_(*orig_mesh, ref_factor, ref_type);
1138}
1139
1140void ParMesh::MakeRefined_(ParMesh &orig_mesh, int ref_factor, int ref_type)
1141{
1142 MyComm = orig_mesh.GetComm();
1143 NRanks = orig_mesh.GetNRanks();
1144 MyRank = orig_mesh.GetMyRank();
1145 face_nbr_el_to_face = nullptr;
1146 glob_elem_offset = -1;
1148 gtopo = orig_mesh.gtopo;
1149 have_face_nbr_data = false;
1150 pncmesh = NULL;
1151
1152 Array<int> ref_factors(orig_mesh.GetNE());
1153 ref_factors = ref_factor;
1154 Mesh::MakeRefined_(orig_mesh, ref_factors, ref_type);
1155
1156 // Need to initialize:
1157 // - shared_edges, shared_{trias,quads}
1158 // - group_svert, group_sedge, group_{stria,squad}
1159 // - svert_lvert, sedge_ledge, sface_lface
1160
1161 meshgen = orig_mesh.meshgen; // copy the global 'meshgen'
1162
1163 H1_FECollection rfec(ref_factor, Dim, ref_type);
1164 ParFiniteElementSpace rfes(&orig_mesh, &rfec);
1165
1166 // count the number of entries in each row of group_s{vert,edge,face}
1167 group_svert.MakeI(GetNGroups()-1); // exclude the local group 0
1171 for (int gr = 1; gr < GetNGroups(); gr++)
1172 {
1173 // orig vertex -> vertex
1174 group_svert.AddColumnsInRow(gr-1, orig_mesh.GroupNVertices(gr));
1175 // orig edge -> (ref_factor-1) vertices and (ref_factor) edges
1176 const int orig_ne = orig_mesh.GroupNEdges(gr);
1177 group_svert.AddColumnsInRow(gr-1, (ref_factor-1)*orig_ne);
1178 group_sedge.AddColumnsInRow(gr-1, ref_factor*orig_ne);
1179 // orig face -> (?) vertices, (?) edges, and (?) faces
1180 const int orig_nt = orig_mesh.GroupNTriangles(gr);
1181 if (orig_nt > 0)
1182 {
1184 const int nvert = Geometry::NumVerts[geom];
1185 RefinedGeometry &RG =
1186 *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1187
1188 // count internal vertices
1189 group_svert.AddColumnsInRow(gr-1, orig_nt*rfec.DofForGeometry(geom));
1190 // count internal edges
1191 group_sedge.AddColumnsInRow(gr-1, orig_nt*(RG.RefEdges.Size()/2-
1192 RG.NumBdrEdges));
1193 // count refined faces
1194 group_stria.AddColumnsInRow(gr-1, orig_nt*(RG.RefGeoms.Size()/nvert));
1195 }
1196 const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1197 if (orig_nq > 0)
1198 {
1199 const Geometry::Type geom = Geometry::SQUARE;
1200 const int nvert = Geometry::NumVerts[geom];
1201 RefinedGeometry &RG =
1202 *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1203
1204 // count internal vertices
1205 group_svert.AddColumnsInRow(gr-1, orig_nq*rfec.DofForGeometry(geom));
1206 // count internal edges
1207 group_sedge.AddColumnsInRow(gr-1, orig_nq*(RG.RefEdges.Size()/2-
1208 RG.NumBdrEdges));
1209 // count refined faces
1210 group_squad.AddColumnsInRow(gr-1, orig_nq*(RG.RefGeoms.Size()/nvert));
1211 }
1212 }
1213
1216
1220
1226
1227 Array<int> rdofs;
1228 for (int gr = 1; gr < GetNGroups(); gr++)
1229 {
1230 // add shared vertices from original shared vertices
1231 const int orig_n_verts = orig_mesh.GroupNVertices(gr);
1232 for (int j = 0; j < orig_n_verts; j++)
1233 {
1234 rfes.GetVertexDofs(orig_mesh.GroupVertex(gr, j), rdofs);
1235 group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[0])-1);
1236 }
1237
1238 // add refined shared edges; add shared vertices from refined shared edges
1239 const int orig_n_edges = orig_mesh.GroupNEdges(gr);
1240 if (orig_n_edges > 0)
1241 {
1242 const Geometry::Type geom = Geometry::SEGMENT;
1243 const int nvert = Geometry::NumVerts[geom];
1244 RefinedGeometry &RG = *GlobGeometryRefiner.Refine(geom, ref_factor);
1245 const int *c2h_map = rfec.GetDofMap(geom, ref_factor); // FIXME hp
1246
1247 for (int e = 0; e < orig_n_edges; e++)
1248 {
1249 rfes.GetSharedEdgeDofs(gr, e, rdofs);
1250 MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
1251 // add the internal edge 'rdofs' as shared vertices
1252 for (int j = 2; j < rdofs.Size(); j++)
1253 {
1254 group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[j])-1);
1255 }
1256 for (int j = 0; j < RG.RefGeoms.Size(); j += nvert)
1257 {
1258 Element *elem = NewElement(geom);
1259 int *v = elem->GetVertices();
1260 for (int k = 0; k < nvert; k++)
1261 {
1262 int cid = RG.RefGeoms[j+k]; // local Cartesian index
1263 v[k] = rdofs[c2h_map[cid]];
1264 }
1265 group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
1266 }
1267 }
1268 }
1269 // add refined shared faces; add shared edges and shared vertices from
1270 // refined shared faces
1271 const int orig_nt = orig_mesh.GroupNTriangles(gr);
1272 if (orig_nt > 0)
1273 {
1275 const int nvert = Geometry::NumVerts[geom];
1276 RefinedGeometry &RG =
1277 *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1278 const int num_int_verts = rfec.DofForGeometry(geom);
1279 const int *c2h_map = rfec.GetDofMap(geom, ref_factor); // FIXME hp
1280
1281 for (int f = 0; f < orig_nt; f++)
1282 {
1283 rfes.GetSharedTriangleDofs(gr, f, rdofs);
1284 MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
1285 // add the internal face 'rdofs' as shared vertices
1286 for (int j = rdofs.Size()-num_int_verts; j < rdofs.Size(); j++)
1287 {
1288 group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[j])-1);
1289 }
1290 // add the internal (for the shared face) edges as shared edges
1291 for (int j = 2*RG.NumBdrEdges; j < RG.RefEdges.Size(); j += 2)
1292 {
1294 int *v = elem->GetVertices();
1295 for (int k = 0; k < 2; k++)
1296 {
1297 v[k] = rdofs[c2h_map[RG.RefEdges[j+k]]];
1298 }
1299 group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
1300 }
1301 // add refined shared faces
1302 for (int j = 0; j < RG.RefGeoms.Size(); j += nvert)
1303 {
1304 shared_trias.SetSize(shared_trias.Size()+1);
1305 int *v = shared_trias.Last().v;
1306 for (int k = 0; k < nvert; k++)
1307 {
1308 int cid = RG.RefGeoms[j+k]; // local Cartesian index
1309 v[k] = rdofs[c2h_map[cid]];
1310 }
1311 group_stria.AddConnection(gr-1, shared_trias.Size()-1);
1312 }
1313 }
1314 }
1315 const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1316 if (orig_nq > 0)
1317 {
1318 const Geometry::Type geom = Geometry::SQUARE;
1319 const int nvert = Geometry::NumVerts[geom];
1320 RefinedGeometry &RG =
1321 *GlobGeometryRefiner.Refine(geom, ref_factor, ref_factor);
1322 const int num_int_verts = rfec.DofForGeometry(geom);
1323 const int *c2h_map = rfec.GetDofMap(geom, ref_factor); // FIXME hp
1324
1325 for (int f = 0; f < orig_nq; f++)
1326 {
1327 rfes.GetSharedQuadrilateralDofs(gr, f, rdofs);
1328 MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
1329 // add the internal face 'rdofs' as shared vertices
1330 for (int j = rdofs.Size()-num_int_verts; j < rdofs.Size(); j++)
1331 {
1332 group_svert.AddConnection(gr-1, svert_lvert.Append(rdofs[j])-1);
1333 }
1334 // add the internal (for the shared face) edges as shared edges
1335 for (int j = 2*RG.NumBdrEdges; j < RG.RefEdges.Size(); j += 2)
1336 {
1338 int *v = elem->GetVertices();
1339 for (int k = 0; k < 2; k++)
1340 {
1341 v[k] = rdofs[c2h_map[RG.RefEdges[j+k]]];
1342 }
1343 group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
1344 }
1345 // add refined shared faces
1346 for (int j = 0; j < RG.RefGeoms.Size(); j += nvert)
1347 {
1348 shared_quads.SetSize(shared_quads.Size()+1);
1349 int *v = shared_quads.Last().v;
1350 for (int k = 0; k < nvert; k++)
1351 {
1352 int cid = RG.RefGeoms[j+k]; // local Cartesian index
1353 v[k] = rdofs[c2h_map[cid]];
1354 }
1355 group_squad.AddConnection(gr-1, shared_quads.Size()-1);
1356 }
1357 }
1358 }
1359 }
1364
1366
1367 if (Nodes != NULL)
1368 {
1369 // This call will turn the Nodes into a ParGridFunction
1370 SetCurvature(1, GetNodalFESpace()->IsDGSpace(), spaceDim,
1371 GetNodalFESpace()->GetOrdering());
1372 }
1373}
1374
1375ParMesh ParMesh::MakeRefined(ParMesh &orig_mesh, int ref_factor, int ref_type)
1376{
1377 ParMesh mesh;
1378 mesh.MakeRefined_(orig_mesh, ref_factor, ref_type);
1379 return mesh;
1380}
1381
1383{
1384 ParMesh mesh;
1385
1386 mesh.MyComm = orig_mesh.GetComm();
1387 mesh.NRanks = orig_mesh.GetNRanks();
1388 mesh.MyRank = orig_mesh.GetMyRank();
1389 mesh.glob_elem_offset = -1;
1390 mesh.glob_offset_sequence = -1;
1391 mesh.gtopo = orig_mesh.gtopo;
1392 mesh.have_face_nbr_data = false;
1393 mesh.pncmesh = NULL;
1394 mesh.meshgen = orig_mesh.meshgen;
1395
1396 H1_FECollection fec(1, orig_mesh.Dimension());
1397 ParFiniteElementSpace fes(&orig_mesh, &fec);
1398
1399 Array<int> vglobal(orig_mesh.GetNV());
1400 for (int iv=0; iv<orig_mesh.GetNV(); ++iv)
1401 {
1402 vglobal[iv] = fes.GetGlobalTDofNumber(iv);
1403 }
1404 mesh.MakeSimplicial_(orig_mesh, vglobal);
1405
1406 // count the number of entries in each row of group_s{vert,edge,face}
1407 mesh.group_svert.MakeI(mesh.GetNGroups()-1); // exclude the local group 0
1408 mesh.group_sedge.MakeI(mesh.GetNGroups()-1);
1409 mesh.group_stria.MakeI(mesh.GetNGroups()-1);
1410 mesh.group_squad.MakeI(mesh.GetNGroups()-1);
1411 for (int gr = 1; gr < mesh.GetNGroups(); gr++)
1412 {
1413 mesh.group_svert.AddColumnsInRow(gr-1, orig_mesh.GroupNVertices(gr));
1414 mesh.group_sedge.AddColumnsInRow(gr-1, orig_mesh.GroupNEdges(gr));
1415 // Every quad gives an extra edge
1416 const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1417 mesh.group_sedge.AddColumnsInRow(gr-1, orig_nq);
1418 // Every quad is subdivided into two triangles
1419 mesh.group_stria.AddColumnsInRow(gr-1, 2*orig_nq);
1420 // Existing triangles remain unchanged
1421 const int orig_nt = orig_mesh.GroupNTriangles(gr);
1422 mesh.group_stria.AddColumnsInRow(gr-1, orig_nt);
1423 }
1424 mesh.group_svert.MakeJ();
1426
1427 mesh.group_sedge.MakeJ();
1428 mesh.shared_edges.Reserve(mesh.group_sedge.Size_of_connections());
1430
1431 mesh.group_stria.MakeJ();
1432 mesh.shared_trias.Reserve(mesh.group_stria.Size_of_connections());
1433 mesh.sface_lface.SetSize(mesh.shared_trias.Size());
1434
1435 mesh.group_squad.MakeJ();
1436
1437 constexpr int ntris = 2, nv_tri = 3, nv_quad = 4;
1438
1439 Array<int> dofs;
1440 for (int gr = 1; gr < mesh.GetNGroups(); gr++)
1441 {
1442 // add shared vertices from original shared vertices
1443 const int orig_n_verts = orig_mesh.GroupNVertices(gr);
1444 for (int j = 0; j < orig_n_verts; j++)
1445 {
1446 fes.GetVertexDofs(orig_mesh.GroupVertex(gr, j), dofs);
1447 mesh.group_svert.AddConnection(gr-1, mesh.svert_lvert.Append(dofs[0])-1);
1448 }
1449
1450 // add original shared edges
1451 const int orig_n_edges = orig_mesh.GroupNEdges(gr);
1452 for (int e = 0; e < orig_n_edges; e++)
1453 {
1454 int iedge, o;
1455 orig_mesh.GroupEdge(gr, e, iedge, o);
1456 Element *elem = mesh.NewElement(Geometry::SEGMENT);
1457 Array<int> edge_verts;
1458 orig_mesh.GetEdgeVertices(iedge, edge_verts);
1459 elem->SetVertices(edge_verts);
1460 mesh.group_sedge.AddConnection(gr-1, mesh.shared_edges.Append(elem)-1);
1461 }
1462 // add original shared triangles
1463 const int orig_nt = orig_mesh.GroupNTriangles(gr);
1464 for (int e = 0; e < orig_nt; e++)
1465 {
1466 int itri, o;
1467 orig_mesh.GroupTriangle(gr, e, itri, o);
1468 const int *v = orig_mesh.GetFace(itri)->GetVertices();
1469 mesh.shared_trias.SetSize(mesh.shared_trias.Size()+1);
1470 int *v2 = mesh.shared_trias.Last().v;
1471 for (int iv=0; iv<nv_tri; ++iv) { v2[iv] = v[iv]; }
1472 mesh.group_stria.AddConnection(gr-1, mesh.shared_trias.Size()-1);
1473 }
1474 // add triangles from split quads and add resulting diagonal edge
1475 const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
1476 if (orig_nq > 0)
1477 {
1478 static const int trimap[12] =
1479 {
1480 0, 0, 0, 1,
1481 1, 2, 1, 2,
1482 2, 3, 3, 3
1483 };
1484 static const int diagmap[4] = { 0, 2, 1, 3 };
1485 for (int f = 0; f < orig_nq; ++f)
1486 {
1487 int iquad, o;
1488 orig_mesh.GroupQuadrilateral(gr, f, iquad, o);
1489 const int *v = orig_mesh.GetFace(iquad)->GetVertices();
1490 // Split quad according the smallest (global) vertex
1491 int vg[nv_quad];
1492 for (int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
1493 int iv_min = std::min_element(vg, vg+nv_quad) - vg;
1494 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
1495 // Add diagonal
1496 Element *diag = mesh.NewElement(Geometry::SEGMENT);
1497 int *v_diag = diag->GetVertices();
1498 v_diag[0] = v[diagmap[0 + isplit*2]];
1499 v_diag[1] = v[diagmap[1 + isplit*2]];
1500 mesh.group_sedge.AddConnection(gr-1, mesh.shared_edges.Append(diag)-1);
1501 // Add two new triangles
1502 for (int itri=0; itri<ntris; ++itri)
1503 {
1504 mesh.shared_trias.SetSize(mesh.shared_trias.Size()+1);
1505 int *v2 = mesh.shared_trias.Last().v;
1506 for (int iv=0; iv<nv_tri; ++iv)
1507 {
1508 v2[iv] = v[trimap[itri + isplit*2 + iv*ntris*2]];
1509 }
1510 mesh.group_stria.AddConnection(gr-1, mesh.shared_trias.Size()-1);
1511 }
1512 }
1513 }
1514 }
1515 mesh.group_svert.ShiftUpI();
1516 mesh.group_sedge.ShiftUpI();
1517 mesh.group_stria.ShiftUpI();
1518
1519 mesh.FinalizeParTopo();
1520
1521 return mesh;
1522}
1523
1524void ParMesh::Finalize(bool refine, bool fix_orientation)
1525{
1526 const int meshgen_save = meshgen; // Mesh::Finalize() may call SetMeshGen()
1527 // 'mesh_geoms' is local, so there's no need to save and restore it.
1528
1529 Mesh::Finalize(refine, fix_orientation);
1530
1531 meshgen = meshgen_save;
1532 // Note: if Mesh::Finalize() calls MarkTetMeshForRefinement() then the
1533 // shared_trias have been rotated as necessary.
1534
1535 // Setup secondary parallel mesh data: sedge_ledge, sface_lface
1537}
1538
1539int ParMesh::GetLocalElementNum(long long global_element_num) const
1540{
1542 long long local = global_element_num - glob_elem_offset;
1543 if (local < 0 || local >= NumOfElements) { return -1; }
1544 return local;
1545}
1546
1547long long ParMesh::GetGlobalElementNum(int local_element_num) const
1548{
1550 return glob_elem_offset + local_element_num;
1551}
1552
1554{
1555 // Determine the largest attribute number across all processors
1556 int max_attr = attr.Size() ? attr.Max() : 1 /*allow empty ranks*/;
1557 int glb_max_attr = -1;
1558 MPI_Allreduce(&max_attr, &glb_max_attr, 1, MPI_INT, MPI_MAX, MyComm);
1559
1560 // Create marker arrays to indicate which attributes are present
1561 // assuming attribute numbers are in the range [1,glb_max_attr].
1562 bool *attr_marker = new bool[glb_max_attr];
1563 bool *glb_attr_marker = new bool[glb_max_attr];
1564 for (int i = 0; i < glb_max_attr; i++)
1565 {
1566 attr_marker[i] = false;
1567 }
1568 for (int i = 0; i < attr.Size(); i++)
1569 {
1570 attr_marker[attr[i] - 1] = true;
1571 }
1572 MPI_Allreduce(attr_marker, glb_attr_marker, glb_max_attr,
1573 MPI_C_BOOL, MPI_LOR, MyComm);
1574 delete [] attr_marker;
1575
1576 // Translate from the marker array to a unique, sorted list of attributes
1577 attr.SetSize(0);
1578 attr.Reserve(glb_max_attr);
1579 for (int i = 0; i < glb_max_attr; i++)
1580 {
1581 if (glb_attr_marker[i])
1582 {
1583 attr.Append(i + 1);
1584 }
1585 }
1586 delete [] glb_attr_marker;
1587}
1588
1590{
1591 // Determine the attributes occurring in local interior and boundary elements
1593
1595 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1596 {
1597 MFEM_WARNING("Non-positive boundary element attributes found!");
1598 }
1599
1601 if (attributes.Size() > 0 && attributes[0] <= 0)
1602 {
1603 MFEM_WARNING("Non-positive element attributes found!");
1604 }
1605}
1606
1608{
1609 // maximum number of boundary elements over all ranks
1610 int maxNumOfBdrElements;
1611 MPI_Allreduce(&NumOfBdrElements, &maxNumOfBdrElements, 1,
1612 MPI_INT, MPI_MAX, MyComm);
1613 return (maxNumOfBdrElements > 0);
1614}
1615
1616void ParMesh::GroupEdge(int group, int i, int &edge, int &o) const
1617{
1618 int sedge = group_sedge.GetRow(group-1)[i];
1619 edge = sedge_ledge[sedge];
1620 int *v = shared_edges[sedge]->GetVertices();
1621 o = (v[0] < v[1]) ? (+1) : (-1);
1622}
1623
1624void ParMesh::GroupTriangle(int group, int i, int &face, int &o) const
1625{
1626 int stria = group_stria.GetRow(group-1)[i];
1627 face = sface_lface[stria];
1628 // face gives the base orientation
1629 MFEM_ASSERT(faces[face]->GetType() == Element::TRIANGLE,
1630 "Expecting a triangular face.");
1631
1632 o = GetTriOrientation(faces[face]->GetVertices(), shared_trias[stria].v);
1633}
1634
1635void ParMesh::GroupQuadrilateral(int group, int i, int &face, int &o) const
1636{
1637 int squad = group_squad.GetRow(group-1)[i];
1638 face = sface_lface[shared_trias.Size()+squad];
1639 // face gives the base orientation
1640 MFEM_ASSERT(faces[face]->GetType() == Element::QUADRILATERAL,
1641 "Expecting a quadrilateral face.");
1642
1643 o = GetQuadOrientation(faces[face]->GetVertices(), shared_quads[squad].v);
1644}
1645
1647 GroupCommunicator& sedge_comm) const
1648{
1649 Table &gr_sedge = sedge_comm.GroupLDofTable();
1650 gr_sedge.SetDims(GetNGroups(), shared_edges.Size());
1651 gr_sedge.GetI()[0] = 0;
1652 for (int gr = 1; gr <= GetNGroups(); gr++)
1653 {
1654 gr_sedge.GetI()[gr] = group_sedge.GetI()[gr-1];
1655 }
1656 for (int k = 0; k < shared_edges.Size(); k++)
1657 {
1658 if (ordering == 1)
1659 {
1660 gr_sedge.GetJ()[k] =k;
1661 }
1662 else
1663 {
1664 gr_sedge.GetJ()[k] = group_sedge.GetJ()[k];
1665 }
1666 }
1667 sedge_comm.Finalize();
1668}
1669
1671 GroupCommunicator& svert_comm) const
1672{
1673 Table &gr_svert = svert_comm.GroupLDofTable();
1674 gr_svert.SetDims(GetNGroups(), svert_lvert.Size());
1675 gr_svert.GetI()[0] = 0;
1676 for (int gr = 1; gr <= GetNGroups(); gr++)
1677 {
1678 gr_svert.GetI()[gr] = group_svert.GetI()[gr-1];
1679 }
1680 for (int k = 0; k < svert_lvert.Size(); k++)
1681 {
1682 if (ordering == 1)
1683 {
1684 gr_svert.GetJ()[k] = k;
1685 }
1686 else
1687 {
1688 gr_svert.GetJ()[k] = group_svert.GetJ()[k];
1689 }
1690 }
1691 svert_comm.Finalize();
1692}
1693
1695 GroupCommunicator& squad_comm) const
1696{
1697 Table &gr_squad = squad_comm.GroupLDofTable();
1698 gr_squad.SetDims(GetNGroups(), shared_quads.Size());
1699 gr_squad.GetI()[0] = 0;
1700 for (int gr = 1; gr <= GetNGroups(); gr++)
1701 {
1702 gr_squad.GetI()[gr] = group_squad.GetI()[gr-1];
1703 }
1704 for (int k = 0; k < shared_quads.Size(); k++)
1705 {
1706 if (ordering == 1)
1707 {
1708 gr_squad.GetJ()[k] = k;
1709 }
1710 else
1711 {
1712 gr_squad.GetJ()[k] = group_squad.GetJ()[k];
1713 }
1714 }
1715 squad_comm.Finalize();
1716}
1717
1719 GroupCommunicator& stria_comm) const
1720{
1721 Table &gr_stria = stria_comm.GroupLDofTable();
1722 gr_stria.SetDims(GetNGroups(), shared_trias.Size());
1723 gr_stria.GetI()[0] = 0;
1724 for (int gr = 1; gr <= GetNGroups(); gr++)
1725 {
1726 gr_stria.GetI()[gr] = group_stria.GetI()[gr-1];
1727 }
1728 for (int k = 0; k < shared_trias.Size(); k++)
1729 {
1730 if (ordering == 1)
1731 {
1732 gr_stria.GetJ()[k] = k;
1733 }
1734 else
1735 {
1736 gr_stria.GetJ()[k] = group_stria.GetJ()[k];
1737 }
1738 }
1739 stria_comm.Finalize();
1740}
1741
1743{
1744 Array<int> order;
1745 GetEdgeOrdering(v_to_v, order); // local edge ordering
1746
1747 // create a GroupCommunicator on the shared edges
1748 GroupCommunicator sedge_comm(gtopo);
1749 GetSharedEdgeCommunicator(0, sedge_comm);
1750
1751 Array<int> sedge_ord(shared_edges.Size());
1752 Array<Pair<int,int> > sedge_ord_map(shared_edges.Size());
1753 for (int k = 0; k < shared_edges.Size(); k++)
1754 {
1755 // sedge_ledge may be undefined -- use shared_edges and v_to_v instead
1756 const int sedge = group_sedge.GetJ()[k];
1757 const int *v = shared_edges[sedge]->GetVertices();
1758 sedge_ord[k] = order[v_to_v(v[0], v[1])];
1759 }
1760
1761 sedge_comm.Bcast<int>(sedge_ord, 1);
1762
1763 for (int k = 0, gr = 1; gr < GetNGroups(); gr++)
1764 {
1765 const int n = group_sedge.RowSize(gr-1);
1766 if (n == 0) { continue; }
1767 sedge_ord_map.SetSize(n);
1768 for (int j = 0; j < n; j++)
1769 {
1770 sedge_ord_map[j].one = sedge_ord[k+j];
1771 sedge_ord_map[j].two = j;
1772 }
1773 SortPairs<int, int>(sedge_ord_map, n);
1774 for (int j = 0; j < n; j++)
1775 {
1776 const int sedge_from = group_sedge.GetJ()[k+j];
1777 const int *v = shared_edges[sedge_from]->GetVertices();
1778 sedge_ord[k+j] = order[v_to_v(v[0], v[1])];
1779 }
1780 std::sort(&sedge_ord[k], &sedge_ord[k] + n);
1781 for (int j = 0; j < n; j++)
1782 {
1783 const int sedge_to = group_sedge.GetJ()[k+sedge_ord_map[j].two];
1784 const int *v = shared_edges[sedge_to]->GetVertices();
1785 order[v_to_v(v[0], v[1])] = sedge_ord[k+j];
1786 }
1787 k += n;
1788 }
1789
1790#ifdef MFEM_DEBUG
1791 {
1792 Array<Pair<int, real_t> > ilen_len(order.Size());
1793
1794 for (int i = 0; i < NumOfVertices; i++)
1795 {
1796 for (DSTable::RowIterator it(v_to_v, i); !it; ++it)
1797 {
1798 int j = it.Index();
1799 ilen_len[j].one = order[j];
1800 ilen_len[j].two = GetLength(i, it.Column());
1801 }
1802 }
1803
1804 SortPairs<int, real_t>(ilen_len, order.Size());
1805
1806 real_t d_max = 0.;
1807 for (int i = 1; i < order.Size(); i++)
1808 {
1809 d_max = std::max(d_max, ilen_len[i-1].two-ilen_len[i].two);
1810 }
1811
1812#if 0
1813 // Debug message from every MPI rank.
1814 mfem::out << "proc. " << MyRank << '/' << NRanks << ": d_max = " << d_max
1815 << endl;
1816#else
1817 // Debug message just from rank 0.
1818 real_t glob_d_max;
1819 MPI_Reduce(&d_max, &glob_d_max, 1, MPITypeMap<real_t>::mpi_type, MPI_MAX, 0,
1820 MyComm);
1821 if (MyRank == 0)
1822 {
1823 mfem::out << "glob_d_max = " << glob_d_max << endl;
1824 }
1825#endif
1826 }
1827#endif
1828
1829 // use 'order' to mark the tets, the boundary triangles, and the shared
1830 // triangle faces
1831 for (int i = 0; i < NumOfElements; i++)
1832 {
1833 if (elements[i]->GetType() == Element::TETRAHEDRON)
1834 {
1835 elements[i]->MarkEdge(v_to_v, order);
1836 }
1837 }
1838
1839 for (int i = 0; i < NumOfBdrElements; i++)
1840 {
1841 if (boundary[i]->GetType() == Element::TRIANGLE)
1842 {
1843 boundary[i]->MarkEdge(v_to_v, order);
1844 }
1845 }
1846
1847 for (int i = 0; i < shared_trias.Size(); i++)
1848 {
1849 Triangle::MarkEdge(shared_trias[i].v, v_to_v, order);
1850 }
1851}
1852
1853// For a line segment with vertices v[0] and v[1], return a number with
1854// the following meaning:
1855// 0 - the edge was not refined
1856// 1 - the edge e was refined once by splitting v[0],v[1]
1858 int *middle)
1859{
1860 int m, *v = edge->GetVertices();
1861
1862 if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
1863 {
1864 return 1;
1865 }
1866 else
1867 {
1868 return 0;
1869 }
1870}
1871
1872void ParMesh::GetFaceSplittings(const int *fv, const HashTable<Hashed2> &v_to_v,
1873 Array<unsigned> &codes)
1874{
1875 typedef Triple<int,int,int> face_t;
1876 Array<face_t> face_stack;
1877
1878 unsigned code = 0;
1879 face_stack.Append(face_t(fv[0], fv[1], fv[2]));
1880 for (unsigned bit = 0; face_stack.Size() > 0; bit++)
1881 {
1882 if (bit == 8*sizeof(unsigned))
1883 {
1884 codes.Append(code);
1885 code = bit = 0;
1886 }
1887
1888 const face_t &f = face_stack.Last();
1889 int mid = v_to_v.FindId(f.one, f.two);
1890 if (mid == -1)
1891 {
1892 // leave a 0 at bit 'bit'
1893 face_stack.DeleteLast();
1894 }
1895 else
1896 {
1897 code += (1 << bit); // set bit 'bit' to 1
1898 mid += NumOfVertices;
1899 face_stack.Append(face_t(f.three, f.one, mid));
1900 face_t &r = face_stack[face_stack.Size()-2];
1901 r = face_t(r.two, r.three, mid);
1902 }
1903 }
1904 codes.Append(code);
1905}
1906
1908 const Array<unsigned> &codes, int &pos)
1909{
1910 typedef Triple<int,int,int> face_t;
1911 Array<face_t> face_stack;
1912
1913 bool need_refinement = 0;
1914 face_stack.Append(face_t(v[0], v[1], v[2]));
1915 for (unsigned bit = 0, code = codes[pos++]; face_stack.Size() > 0; bit++)
1916 {
1917 if (bit == 8*sizeof(unsigned))
1918 {
1919 code = codes[pos++];
1920 bit = 0;
1921 }
1922
1923 if ((code & (1 << bit)) == 0) { face_stack.DeleteLast(); continue; }
1924
1925 const face_t &f = face_stack.Last();
1926 int mid = v_to_v.FindId(f.one, f.two);
1927 if (mid == -1)
1928 {
1929 mid = v_to_v.GetId(f.one, f.two);
1930 int ind[2] = { f.one, f.two };
1931 vertices.Append(Vertex());
1932 AverageVertices(ind, 2, vertices.Size()-1);
1933 need_refinement = 1;
1934 }
1935 mid += NumOfVertices;
1936 face_stack.Append(face_t(f.three, f.one, mid));
1937 face_t &r = face_stack[face_stack.Size()-2];
1938 r = face_t(r.two, r.three, mid);
1939 }
1940 return need_refinement;
1941}
1942
1944 Array<HYPRE_BigInt> *offsets[]) const
1945{
1946 if (HYPRE_AssumedPartitionCheck())
1947 {
1948 Array<HYPRE_BigInt> temp(N);
1949 MPI_Scan(loc_sizes, temp.GetData(), N, HYPRE_MPI_BIG_INT, MPI_SUM, MyComm);
1950 for (int i = 0; i < N; i++)
1951 {
1952 offsets[i]->SetSize(3);
1953 (*offsets[i])[0] = temp[i] - loc_sizes[i];
1954 (*offsets[i])[1] = temp[i];
1955 }
1956 MPI_Bcast(temp.GetData(), N, HYPRE_MPI_BIG_INT, NRanks-1, MyComm);
1957 for (int i = 0; i < N; i++)
1958 {
1959 (*offsets[i])[2] = temp[i];
1960 // check for overflow
1961 MFEM_VERIFY((*offsets[i])[0] >= 0 && (*offsets[i])[1] >= 0,
1962 "overflow in offsets");
1963 }
1964 }
1965 else
1966 {
1968 MPI_Allgather(loc_sizes, N, HYPRE_MPI_BIG_INT, temp.GetData(), N,
1969 HYPRE_MPI_BIG_INT, MyComm);
1970 for (int i = 0; i < N; i++)
1971 {
1972 Array<HYPRE_BigInt> &offs = *offsets[i];
1973 offs.SetSize(NRanks+1);
1974 offs[0] = 0;
1975 for (int j = 0; j < NRanks; j++)
1976 {
1977 offs[j+1] = offs[j] + temp[i+N*j];
1978 }
1979 // Check for overflow
1980 MFEM_VERIFY(offs[MyRank] >= 0 && offs[MyRank+1] >= 0,
1981 "overflow in offsets");
1982 }
1983 }
1984}
1985
1987{
1988 if (!have_face_nbr_data)
1989 {
1990 return;
1991 }
1992
1993 have_face_nbr_data = false;
1997 for (int i = 0; i < face_nbr_elements.Size(); i++)
1998 {
2000 }
2001 face_nbr_elements.DeleteAll();
2002 face_nbr_vertices.DeleteAll();
2005}
2006
2007void ParMesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
2008{
2010 space_dim = (space_dim == -1) ? spaceDim : space_dim;
2012 if (discont)
2013 {
2014 nfec = new L2_FECollection(order, Dim, BasisType::GaussLobatto);
2015 }
2016 else
2017 {
2018 nfec = new H1_FECollection(order, Dim);
2019 }
2020 ParFiniteElementSpace* nfes = new ParFiniteElementSpace(this, nfec, space_dim,
2021 ordering);
2022 auto pnodes = new ParGridFunction(nfes);
2023 GetNodes(*pnodes);
2024 NewNodes(*pnodes, true);
2025 Nodes->MakeOwner(nfec);
2026}
2027
2029{
2031 ParFiniteElementSpace *npfes = dynamic_cast<ParFiniteElementSpace*>(nfes);
2032 if (npfes)
2033 {
2034 SetNodalFESpace(npfes);
2035 }
2036 else
2037 {
2039 }
2040}
2041
2043{
2045 ParGridFunction *nodes = new ParGridFunction(npfes);
2046 SetNodalGridFunction(nodes, true);
2047}
2048
2050{
2051 if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
2052 {
2054 ParFiniteElementSpace *pfes =
2055 new ParFiniteElementSpace(*Nodes->FESpace(), *this);
2056 ParGridFunction *new_nodes = new ParGridFunction(pfes);
2057 *new_nodes = *Nodes;
2058 if (Nodes->OwnFEC())
2059 {
2060 new_nodes->MakeOwner(Nodes->OwnFEC());
2061 Nodes->MakeOwner(NULL); // takes away ownership of 'fec' and 'fes'
2062 delete Nodes->FESpace();
2063 }
2064 delete Nodes;
2065 Nodes = new_nodes;
2066 }
2067}
2068
2070{
2072 {
2073 return;
2074 }
2075
2076 if (Nonconforming())
2077 {
2078 // With ParNCMesh we can set up face neighbors mostly without communication.
2079 pncmesh->GetFaceNeighbors(*this);
2080 have_face_nbr_data = true;
2081
2083 return;
2084 }
2085
2086 Table *gr_sface;
2087 int *s2l_face;
2088 bool del_tables = false;
2089 if (Dim == 1)
2090 {
2091 gr_sface = &group_svert;
2092 s2l_face = svert_lvert;
2093 }
2094 else if (Dim == 2)
2095 {
2096 gr_sface = &group_sedge;
2097 s2l_face = sedge_ledge;
2098 }
2099 else
2100 {
2101 s2l_face = sface_lface;
2102 if (shared_trias.Size() == sface_lface.Size())
2103 {
2104 // All shared faces are Triangular
2105 gr_sface = &group_stria;
2106 }
2107 else if (shared_quads.Size() == sface_lface.Size())
2108 {
2109 // All shared faced are Quadrilateral
2110 gr_sface = &group_squad;
2111 }
2112 else
2113 {
2114 // Shared faces contain a mixture of triangles and quads
2115 gr_sface = new Table;
2116 del_tables = true;
2117
2118 // Merge the Tables group_stria and group_squad
2119 gr_sface->MakeI(group_stria.Size());
2120 for (int gr=0; gr<group_stria.Size(); gr++)
2121 {
2122 gr_sface->AddColumnsInRow(gr,
2123 group_stria.RowSize(gr) +
2124 group_squad.RowSize(gr));
2125 }
2126 gr_sface->MakeJ();
2127 const int nst = shared_trias.Size();
2128 for (int gr=0; gr<group_stria.Size(); gr++)
2129 {
2130 gr_sface->AddConnections(gr,
2131 group_stria.GetRow(gr),
2132 group_stria.RowSize(gr));
2133 for (int c=0; c<group_squad.RowSize(gr); c++)
2134 {
2135 gr_sface->AddConnection(gr,
2136 nst + group_squad.GetRow(gr)[c]);
2137 }
2138 }
2139 gr_sface->ShiftUpI();
2140 }
2141 }
2142
2143 ExchangeFaceNbrData(gr_sface, s2l_face);
2144
2145 if (Dim == 3)
2146 {
2148 }
2149
2150 if (del_tables) { delete gr_sface; }
2151
2152 if ( have_face_nbr_data ) { return; }
2153
2154 have_face_nbr_data = true;
2155
2157}
2158
2159void ParMesh::ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
2160{
2161 int num_face_nbrs = 0;
2162 for (int g = 1; g < GetNGroups(); g++)
2163 {
2164 if (gr_sface->RowSize(g-1) > 0)
2165 {
2166 num_face_nbrs++;
2167 }
2168 }
2169
2170 face_nbr_group.SetSize(num_face_nbrs);
2171
2172 if (num_face_nbrs == 0)
2173 {
2174 have_face_nbr_data = true;
2175 return;
2176 }
2177
2178 {
2179 // sort face-neighbors by processor rank
2180 Array<Pair<int, int> > rank_group(num_face_nbrs);
2181
2182 for (int g = 1, counter = 0; g < GetNGroups(); g++)
2183 {
2184 if (gr_sface->RowSize(g-1) > 0)
2185 {
2186 MFEM_ASSERT(gtopo.GetGroupSize(g) == 2, "group size is not 2!");
2187
2188 const int *nbs = gtopo.GetGroup(g);
2189 int lproc = (nbs[0]) ? nbs[0] : nbs[1];
2190 rank_group[counter].one = gtopo.GetNeighborRank(lproc);
2191 rank_group[counter].two = g;
2192 counter++;
2193 }
2194 }
2195
2196 SortPairs<int, int>(rank_group, rank_group.Size());
2197
2198 for (int fn = 0; fn < num_face_nbrs; fn++)
2199 {
2200 face_nbr_group[fn] = rank_group[fn].two;
2201 }
2202 }
2203
2204 MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
2205 MPI_Request *send_requests = requests;
2206 MPI_Request *recv_requests = requests + num_face_nbrs;
2207 MPI_Status *statuses = new MPI_Status[num_face_nbrs];
2208
2209 int *nbr_data = new int[6*num_face_nbrs];
2210 int *nbr_send_data = nbr_data;
2211 int *nbr_recv_data = nbr_data + 3*num_face_nbrs;
2212
2213 Array<int> el_marker(GetNE());
2214 Array<int> vertex_marker(GetNV());
2215 el_marker = -1;
2216 vertex_marker = -1;
2217
2218 Array<int> fcs, cor;
2219
2220 Table send_face_nbr_elemdata, send_face_nbr_facedata;
2221
2222 send_face_nbr_elements.MakeI(num_face_nbrs);
2223 send_face_nbr_vertices.MakeI(num_face_nbrs);
2224 send_face_nbr_elemdata.MakeI(num_face_nbrs);
2225 send_face_nbr_facedata.MakeI(num_face_nbrs);
2226 for (int fn = 0; fn < num_face_nbrs; fn++)
2227 {
2228 int nbr_group = face_nbr_group[fn];
2229 int num_sfaces = gr_sface->RowSize(nbr_group-1);
2230 int *sface = gr_sface->GetRow(nbr_group-1);
2231 for (int i = 0; i < num_sfaces; i++)
2232 {
2233 int lface = s2l_face[sface[i]];
2234 int el = faces_info[lface].Elem1No;
2235 if (el_marker[el] != fn)
2236 {
2237 el_marker[el] = fn;
2239
2240 const int nv = elements[el]->GetNVertices();
2241 const int *v = elements[el]->GetVertices();
2242 for (int j = 0; j < nv; j++)
2243 if (vertex_marker[v[j]] != fn)
2244 {
2245 vertex_marker[v[j]] = fn;
2247 }
2248
2249 const int nf = elements[el]->GetNFaces();
2250
2251 send_face_nbr_elemdata.AddColumnsInRow(fn, nv + nf + 2);
2252 }
2253 }
2254 send_face_nbr_facedata.AddColumnsInRow(fn, 2*num_sfaces);
2255
2256 nbr_send_data[3*fn ] = send_face_nbr_elements.GetI()[fn];
2257 nbr_send_data[3*fn+1] = send_face_nbr_vertices.GetI()[fn];
2258 nbr_send_data[3*fn+2] = send_face_nbr_elemdata.GetI()[fn];
2259
2260 int nbr_rank = GetFaceNbrRank(fn);
2261 int tag = 0;
2262
2263 MPI_Isend(&nbr_send_data[3*fn], 3, MPI_INT, nbr_rank, tag, MyComm,
2264 &send_requests[fn]);
2265 MPI_Irecv(&nbr_recv_data[3*fn], 3, MPI_INT, nbr_rank, tag, MyComm,
2266 &recv_requests[fn]);
2267 }
2270 send_face_nbr_elemdata.MakeJ();
2271 send_face_nbr_facedata.MakeJ();
2272 el_marker = -1;
2273 vertex_marker = -1;
2274 const int nst = shared_trias.Size();
2275 for (int fn = 0; fn < num_face_nbrs; fn++)
2276 {
2277 int nbr_group = face_nbr_group[fn];
2278 int num_sfaces = gr_sface->RowSize(nbr_group-1);
2279 int *sface = gr_sface->GetRow(nbr_group-1);
2280 for (int i = 0; i < num_sfaces; i++)
2281 {
2282 const int sf = sface[i];
2283 int lface = s2l_face[sf];
2284 int el = faces_info[lface].Elem1No;
2285 if (el_marker[el] != fn)
2286 {
2287 el_marker[el] = fn;
2289
2290 const int nv = elements[el]->GetNVertices();
2291 const int *v = elements[el]->GetVertices();
2292 for (int j = 0; j < nv; j++)
2293 if (vertex_marker[v[j]] != fn)
2294 {
2295 vertex_marker[v[j]] = fn;
2297 }
2298
2299 send_face_nbr_elemdata.AddConnection(fn, GetAttribute(el));
2300 send_face_nbr_elemdata.AddConnection(
2301 fn, GetElementBaseGeometry(el));
2302 send_face_nbr_elemdata.AddConnections(fn, v, nv);
2303
2304 if (Dim == 3)
2305 {
2306 const int nf = elements[el]->GetNFaces();
2307 GetElementFaces(el, fcs, cor);
2308 send_face_nbr_elemdata.AddConnections(fn, cor, nf);
2309 }
2310 }
2311 send_face_nbr_facedata.AddConnection(fn, el);
2312 int info = faces_info[lface].Elem1Inf;
2313 // change the orientation in info to be relative to the shared face
2314 // in 1D and 2D keep the orientation equal to 0
2315 if (Dim == 3)
2316 {
2317 const int *lf_v = faces[lface]->GetVertices();
2318 if (sf < nst) // triangle shared face
2319 {
2320 info += GetTriOrientation(shared_trias[sf].v, lf_v);
2321 }
2322 else // quad shared face
2323 {
2324 info += GetQuadOrientation(shared_quads[sf-nst].v, lf_v);
2325 }
2326 }
2327 send_face_nbr_facedata.AddConnection(fn, info);
2328 }
2329 }
2332 send_face_nbr_elemdata.ShiftUpI();
2333 send_face_nbr_facedata.ShiftUpI();
2334
2335 // convert the vertex indices in send_face_nbr_elemdata
2336 // convert the element indices in send_face_nbr_facedata
2337 for (int fn = 0; fn < num_face_nbrs; fn++)
2338 {
2339 int num_elems = send_face_nbr_elements.RowSize(fn);
2340 int *elems = send_face_nbr_elements.GetRow(fn);
2341 int num_verts = send_face_nbr_vertices.RowSize(fn);
2342 int *verts = send_face_nbr_vertices.GetRow(fn);
2343 int *elemdata = send_face_nbr_elemdata.GetRow(fn);
2344 int num_sfaces = send_face_nbr_facedata.RowSize(fn)/2;
2345 int *facedata = send_face_nbr_facedata.GetRow(fn);
2346
2347 for (int i = 0; i < num_verts; i++)
2348 {
2349 vertex_marker[verts[i]] = i;
2350 }
2351
2352 for (int el = 0; el < num_elems; el++)
2353 {
2354 const int nv = elements[elems[el]]->GetNVertices();
2355 const int nf = (Dim == 3) ? elements[elems[el]]->GetNFaces() : 0;
2356 elemdata += 2; // skip the attribute and the geometry type
2357 for (int j = 0; j < nv; j++)
2358 {
2359 elemdata[j] = vertex_marker[elemdata[j]];
2360 }
2361 elemdata += nv + nf;
2362
2363 el_marker[elems[el]] = el;
2364 }
2365
2366 for (int i = 0; i < num_sfaces; i++)
2367 {
2368 facedata[2*i] = el_marker[facedata[2*i]];
2369 }
2370 }
2371
2372 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
2373
2374 Array<int> recv_face_nbr_facedata;
2375 Table recv_face_nbr_elemdata;
2376
2377 // fill-in face_nbr_elements_offset, face_nbr_vertices_offset
2378 face_nbr_elements_offset.SetSize(num_face_nbrs + 1);
2379 face_nbr_vertices_offset.SetSize(num_face_nbrs + 1);
2380 recv_face_nbr_elemdata.MakeI(num_face_nbrs);
2383 for (int fn = 0; fn < num_face_nbrs; fn++)
2384 {
2386 face_nbr_elements_offset[fn] + nbr_recv_data[3*fn];
2388 face_nbr_vertices_offset[fn] + nbr_recv_data[3*fn+1];
2389 recv_face_nbr_elemdata.AddColumnsInRow(fn, nbr_recv_data[3*fn+2]);
2390 }
2391 recv_face_nbr_elemdata.MakeJ();
2392
2393 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2394
2395 // send and receive the element data
2396 for (int fn = 0; fn < num_face_nbrs; fn++)
2397 {
2398 int nbr_rank = GetFaceNbrRank(fn);
2399 int tag = 0;
2400
2401 MPI_Isend(send_face_nbr_elemdata.GetRow(fn),
2402 send_face_nbr_elemdata.RowSize(fn),
2403 MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
2404
2405 MPI_Irecv(recv_face_nbr_elemdata.GetRow(fn),
2406 recv_face_nbr_elemdata.RowSize(fn),
2407 MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
2408 }
2409
2410 // convert the element data into face_nbr_elements
2411 face_nbr_elements.SetSize(face_nbr_elements_offset[num_face_nbrs]);
2412 face_nbr_el_ori.reset(new Table(face_nbr_elements_offset[num_face_nbrs], 6));
2413 while (true)
2414 {
2415 int fn;
2416 MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
2417
2418 if (fn == MPI_UNDEFINED)
2419 {
2420 break;
2421 }
2422
2423 int vert_off = face_nbr_vertices_offset[fn];
2424 int elem_off = face_nbr_elements_offset[fn];
2425 int num_elems = face_nbr_elements_offset[fn+1] - elem_off;
2426 int *recv_elemdata = recv_face_nbr_elemdata.GetRow(fn);
2427
2428 for (int i = 0; i < num_elems; i++)
2429 {
2430 Element *el = NewElement(recv_elemdata[1]);
2431 el->SetAttribute(recv_elemdata[0]);
2432 recv_elemdata += 2;
2433 int nv = el->GetNVertices();
2434 for (int j = 0; j < nv; j++)
2435 {
2436 recv_elemdata[j] += vert_off;
2437 }
2438 el->SetVertices(recv_elemdata);
2439 recv_elemdata += nv;
2440 if (Dim == 3)
2441 {
2442 int nf = el->GetNFaces();
2443 int * fn_ori = face_nbr_el_ori->GetRow(elem_off);
2444 for (int j = 0; j < nf; j++)
2445 {
2446 fn_ori[j] = recv_elemdata[j];
2447 }
2448 recv_elemdata += nf;
2449 }
2450 face_nbr_elements[elem_off++] = el;
2451 }
2452 }
2453 face_nbr_el_ori->Finalize();
2454
2455 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2456
2457 // send and receive the face data
2458 recv_face_nbr_facedata.SetSize(
2459 send_face_nbr_facedata.Size_of_connections());
2460 for (int fn = 0; fn < num_face_nbrs; fn++)
2461 {
2462 int nbr_rank = GetFaceNbrRank(fn);
2463 int tag = 0;
2464
2465 MPI_Isend(send_face_nbr_facedata.GetRow(fn),
2466 send_face_nbr_facedata.RowSize(fn),
2467 MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
2468
2469 // the size of the send and receive face data is the same
2470 MPI_Irecv(&recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]],
2471 send_face_nbr_facedata.RowSize(fn),
2472 MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
2473 }
2474
2475 // transfer the received face data into faces_info
2476 while (true)
2477 {
2478 int fn;
2479 MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
2480
2481 if (fn == MPI_UNDEFINED)
2482 {
2483 break;
2484 }
2485
2486 int elem_off = face_nbr_elements_offset[fn];
2487 int nbr_group = face_nbr_group[fn];
2488 int num_sfaces = gr_sface->RowSize(nbr_group-1);
2489 int *sface = gr_sface->GetRow(nbr_group-1);
2490 int *facedata =
2491 &recv_face_nbr_facedata[send_face_nbr_facedata.GetI()[fn]];
2492
2493 for (int i = 0; i < num_sfaces; i++)
2494 {
2495 const int sf = sface[i];
2496 int lface = s2l_face[sf];
2497 FaceInfo &face_info = faces_info[lface];
2498 face_info.Elem2No = -1 - (facedata[2*i] + elem_off);
2499 int info = facedata[2*i+1];
2500 // change the orientation in info to be relative to the local face
2501 if (Dim < 3)
2502 {
2503 info++; // orientation 0 --> orientation 1
2504 }
2505 else
2506 {
2507 int nbr_ori = info%64, nbr_v[4];
2508 const int *lf_v = faces[lface]->GetVertices();
2509
2510 if (sf < nst) // triangle shared face
2511 {
2512 // apply the nbr_ori to sf_v to get nbr_v
2513 const int *perm = tri_t::Orient[nbr_ori];
2514 const int *sf_v = shared_trias[sf].v;
2515 for (int j = 0; j < 3; j++)
2516 {
2517 nbr_v[perm[j]] = sf_v[j];
2518 }
2519 // get the orientation of nbr_v w.r.t. the local face
2520 nbr_ori = GetTriOrientation(lf_v, nbr_v);
2521 }
2522 else // quad shared face
2523 {
2524 // apply the nbr_ori to sf_v to get nbr_v
2525 const int *perm = quad_t::Orient[nbr_ori];
2526 const int *sf_v = shared_quads[sf-nst].v;
2527 for (int j = 0; j < 4; j++)
2528 {
2529 nbr_v[perm[j]] = sf_v[j];
2530 }
2531 // get the orientation of nbr_v w.r.t. the local face
2532 nbr_ori = GetQuadOrientation(lf_v, nbr_v);
2533 }
2534
2535 info = 64*(info/64) + nbr_ori;
2536 }
2537 face_info.Elem2Inf = info;
2538 }
2539 }
2540
2541 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2542
2543 // allocate the face_nbr_vertices
2544 face_nbr_vertices.SetSize(face_nbr_vertices_offset[num_face_nbrs]);
2545
2546 delete [] nbr_data;
2547
2548 delete [] statuses;
2549 delete [] requests;
2550}
2551
2553{
2554 if (!have_face_nbr_data)
2555 {
2556 ExchangeFaceNbrData(); // calls this method at the end
2557 }
2558 else if (Nodes == NULL)
2559 {
2560 if (Nonconforming())
2561 {
2562 // with ParNCMesh we already have the vertices
2563 return;
2564 }
2565
2566 int num_face_nbrs = GetNFaceNeighbors();
2567
2568 if (!num_face_nbrs) { return; }
2569
2570 MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
2571 MPI_Request *send_requests = requests;
2572 MPI_Request *recv_requests = requests + num_face_nbrs;
2573 MPI_Status *statuses = new MPI_Status[num_face_nbrs];
2574
2575 // allocate buffer and copy the vertices to be sent
2577 for (int i = 0; i < send_vertices.Size(); i++)
2578 {
2579 send_vertices[i] = vertices[send_face_nbr_vertices.GetJ()[i]];
2580 }
2581
2582 // send and receive the vertices
2583 for (int fn = 0; fn < num_face_nbrs; fn++)
2584 {
2585 int nbr_rank = GetFaceNbrRank(fn);
2586 int tag = 0;
2587
2588 MPI_Isend(send_vertices[send_face_nbr_vertices.GetI()[fn]](),
2590 MPITypeMap<real_t>::mpi_type, nbr_rank, tag, MyComm, &send_requests[fn]);
2591
2593 3*(face_nbr_vertices_offset[fn+1] -
2595 MPITypeMap<real_t>::mpi_type, nbr_rank, tag, MyComm, &recv_requests[fn]);
2596 }
2597
2598 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
2599 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2600
2601 delete [] statuses;
2602 delete [] requests;
2603 }
2604 else
2605 {
2606 ParGridFunction *pNodes = dynamic_cast<ParGridFunction *>(Nodes);
2607 MFEM_VERIFY(pNodes != NULL, "Nodes are not ParGridFunction!");
2608 pNodes->ExchangeFaceNbrData();
2609 }
2610}
2611
2613{
2614 STable3D *sfaces_tbl = new STable3D(face_nbr_vertices.Size());
2615 for (int i = 0; i < face_nbr_elements.Size(); i++)
2616 {
2617 const int *v = face_nbr_elements[i]->GetVertices();
2618 switch (face_nbr_elements[i]->GetType())
2619 {
2621 {
2622 for (int j = 0; j < 4; j++)
2623 {
2624 const int *fv = tet_t::FaceVert[j];
2625 sfaces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2626 }
2627 break;
2628 }
2629 case Element::WEDGE:
2630 {
2631 for (int j = 0; j < 2; j++)
2632 {
2633 const int *fv = pri_t::FaceVert[j];
2634 sfaces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2635 }
2636 for (int j = 2; j < 5; j++)
2637 {
2638 const int *fv = pri_t::FaceVert[j];
2639 sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2640 }
2641 break;
2642 }
2643 case Element::PYRAMID:
2644 {
2645 for (int j = 0; j < 1; j++)
2646 {
2647 const int *fv = pyr_t::FaceVert[j];
2648 sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2649 }
2650 for (int j = 1; j < 5; j++)
2651 {
2652 const int *fv = pyr_t::FaceVert[j];
2653 sfaces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2654 }
2655 break;
2656 }
2658 {
2659 // find the face by the vertices with the smallest 3 numbers
2660 // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
2661 for (int j = 0; j < 6; j++)
2662 {
2663 const int *fv = hex_t::FaceVert[j];
2664 sfaces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2665 }
2666 break;
2667 }
2668 default:
2669 MFEM_ABORT("Unexpected type of Element.");
2670 }
2671 }
2672 return sfaces_tbl;
2673}
2674
2675template <int N>
2676void
2678 const std::unique_ptr<STable3D> &faces,
2679 const std::unique_ptr<STable3D> &shared_faces,
2680 int elem, int start, int end, const int fverts[][N])
2681{
2682 for (int i = start; i < end; ++i)
2683 {
2684 // Reference face vertices.
2685 const auto fv = fverts[i];
2686 // Element specific face vertices.
2687 const Vert3 elem_fv(elem_vertices[fv[0]], elem_vertices[fv[1]],
2688 elem_vertices[fv[2]]);
2689
2690 // Check amongst the faces of elements local to this rank for this set of vertices
2691 const int lf = faces->Index(elem_fv.v[0], elem_fv.v[1], elem_fv.v[2]);
2692
2693 // If the face wasn't found amonst processor local elements, search the
2694 // ghosts for this set of vertices.
2695 const int sf = lf < 0 ? shared_faces->Index(elem_fv.v[0], elem_fv.v[1],
2696 elem_fv.v[2]) : -1;
2697 // If find local face -> use that
2698 // else if find shared face -> shift and use that
2699 // else no face found -> set to -1
2700 const int face_to_add = lf < 0 ? (sf >= 0 ? sf + NumOfFaces : -1) : lf;
2701
2702 MFEM_ASSERT(sf >= 0 ||
2703 lf >= 0, "Face must be from a local or a face neighbor element");
2704
2705 // Add this discovered face to the list of faces of this face neighbor element
2706 face_nbr_el_to_face->Push(elem, face_to_add);
2707 }
2708}
2709
2711{
2712 const auto faces = std::unique_ptr<STable3D>(GetFacesTable());
2713 const auto shared_faces = std::unique_ptr<STable3D>(GetSharedFacesTable());
2714
2715 face_nbr_el_to_face.reset(new Table(face_nbr_elements.Size(), 6));
2716
2717 Array<int> v;
2718
2719 // Helper for adding quadrilateral faces.
2720 auto add_quad_faces = [&faces, &shared_faces, &v, this]
2721 (int elem, int start, int end, const int fverts[][4])
2722 {
2723 for (int i = start; i < end; ++i)
2724 {
2725 const int * const fv = fverts[i];
2726 int k = 0;
2727 int max = v[fv[0]];
2728
2729 if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2730 if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2731 if (max < v[fv[3]]) { k = 3; }
2732
2733 int v0 = -1, v1 = -1, v2 = -1;
2734 switch (k)
2735 {
2736 case 0:
2737 v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2738 break;
2739 case 1:
2740 v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2741 break;
2742 case 2:
2743 v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2744 break;
2745 case 3:
2746 v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2747 break;
2748 }
2749 int lf = faces->Index(v0, v1, v2);
2750 if (lf < 0)
2751 {
2752 lf = shared_faces->Index(v0, v1, v2);
2753 if (lf >= 0)
2754 {
2755 lf += NumOfFaces;
2756 }
2757 }
2758 face_nbr_el_to_face->Push(elem, lf);
2759 }
2760 };
2761
2762 for (int i = 0; i < face_nbr_elements.Size(); i++)
2763 {
2764 face_nbr_elements[i]->GetVertices(v);
2765 switch (face_nbr_elements[i]->GetType())
2766 {
2768 {
2769 AddTriFaces(v, faces, shared_faces, i, 0, 4, tet_t::FaceVert);
2770 break;
2771 }
2772 case Element::WEDGE:
2773 {
2774 AddTriFaces(v, faces, shared_faces, i, 0, 2, pri_t::FaceVert);
2775 add_quad_faces(i, 2, 5, pri_t::FaceVert);
2776 break;
2777 }
2778 case Element::PYRAMID:
2779 {
2780 add_quad_faces(i, 0, 1, pyr_t::FaceVert);
2781 AddTriFaces(v, faces, shared_faces, i, 1, 5, pyr_t::FaceVert);
2782 break;
2783 }
2785 {
2786 add_quad_faces(i, 0, 6, hex_t::FaceVert);
2787 break;
2788 }
2789 default:
2790 MFEM_ABORT("Unexpected type of Element.");
2791 }
2792 }
2793 face_nbr_el_to_face->Finalize();
2794}
2795
2797{
2798 if (Conforming())
2799 {
2800 int nbr_group = face_nbr_group[fn];
2801 const int *nbs = gtopo.GetGroup(nbr_group);
2802 int nbr_lproc = (nbs[0]) ? nbs[0] : nbs[1];
2803 int nbr_rank = gtopo.GetNeighborRank(nbr_lproc);
2804 return nbr_rank;
2805 }
2806 else
2807 {
2808 // NC: simplified handling of face neighbor ranks
2809 return face_nbr_group[fn];
2810 }
2811}
2812
2813void
2815 Array<int> &orientations) const
2816{
2817 int el_nbr = i - GetNE();
2818 if (face_nbr_el_to_face != nullptr && el_nbr < face_nbr_el_to_face->Size())
2819 {
2820 face_nbr_el_to_face->GetRow(el_nbr, faces);
2821 }
2822 else
2823 {
2824 MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
2825 "face_nbr_el_to_face not generated correctly.");
2826 }
2827
2828 if (face_nbr_el_ori != nullptr && el_nbr < face_nbr_el_ori->Size())
2829 {
2830 face_nbr_el_ori->GetRow(el_nbr, orientations);
2831 }
2832 else
2833 {
2834 MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
2835 "face_nbr_el_ori not generated correctly.");
2836 }
2837}
2838
2840{
2841 const Array<int> *s2l_face;
2842 if (Dim == 1)
2843 {
2844 s2l_face = &svert_lvert;
2845 }
2846 else if (Dim == 2)
2847 {
2848 s2l_face = &sedge_ledge;
2849 }
2850 else
2851 {
2852 s2l_face = &sface_lface;
2853 }
2854
2855 Table *face_elem = new Table;
2856
2857 face_elem->MakeI(faces_info.Size());
2858
2859 for (int i = 0; i < faces_info.Size(); i++)
2860 {
2861 if (faces_info[i].Elem2No >= 0)
2862 {
2863 face_elem->AddColumnsInRow(i, 2);
2864 }
2865 else
2866 {
2867 face_elem->AddAColumnInRow(i);
2868 }
2869 }
2870 for (int i = 0; i < s2l_face->Size(); i++)
2871 {
2872 face_elem->AddAColumnInRow((*s2l_face)[i]);
2873 }
2874
2875 face_elem->MakeJ();
2876
2877 for (int i = 0; i < faces_info.Size(); i++)
2878 {
2879 face_elem->AddConnection(i, faces_info[i].Elem1No);
2880 if (faces_info[i].Elem2No >= 0)
2881 {
2882 face_elem->AddConnection(i, faces_info[i].Elem2No);
2883 }
2884 }
2885 for (int i = 0; i < s2l_face->Size(); i++)
2886 {
2887 int lface = (*s2l_face)[i];
2888 int nbr_elem_idx = -1 - faces_info[lface].Elem2No;
2889 face_elem->AddConnection(lface, NumOfElements + nbr_elem_idx);
2890 }
2891
2892 face_elem->ShiftUpI();
2893
2894 return face_elem;
2895}
2896
2904
2909 int mask) const
2910{
2911 if (FaceNo < GetNumFaces())
2912 {
2913 Mesh::GetFaceElementTransformations(FaceNo, FElTr, ElTr1, ElTr2, mask);
2914 }
2915 else
2916 {
2917 const bool fill2 = mask & 10; // Elem2 and/or Loc2
2918 GetSharedFaceTransformationsByLocalIndex(FaceNo, FElTr, ElTr1, ElTr2,
2919 fill2);
2920 }
2921}
2922
2930
2935 bool fill2) const
2936{
2937 int FaceNo = GetSharedFace(sf);
2938 GetSharedFaceTransformationsByLocalIndex(FaceNo, FElTr, ElTr1, ElTr2, fill2);
2939}
2940
2948
2950 int FaceNo, FaceElementTransformations &FElTr,
2952 bool fill2) const
2953{
2954 const FaceInfo &face_info = faces_info[FaceNo];
2955 MFEM_VERIFY(face_info.Elem2Inf >= 0, "The face must be shared.");
2956
2957 bool is_slave = Nonconforming() && IsSlaveFace(face_info);
2958 bool is_ghost = Nonconforming() && FaceNo >= GetNumFaces();
2959
2960 int mask = 0;
2961 FElTr.SetConfigurationMask(0);
2962 FElTr.Elem1 = NULL;
2963 FElTr.Elem2 = NULL;
2964
2965 int local_face =
2966 is_ghost ? nc_faces_info[face_info.NCFace].MasterFace : FaceNo;
2967 Element::Type face_type = GetFaceElementType(local_face);
2968 Geometry::Type face_geom = GetFaceGeometry(local_face);
2969
2970 // setup the transformation for the first element
2971 FElTr.Elem1No = face_info.Elem1No;
2972 GetElementTransformation(FElTr.Elem1No, &ElTr1);
2973 FElTr.Elem1 = &ElTr1;
2975
2976 // setup the transformation for the second (neighbor) element
2977 int Elem2NbrNo;
2978 if (fill2)
2979 {
2980 Elem2NbrNo = -1 - face_info.Elem2No;
2981 // Store the "shifted index" for element 2 in FElTr.Elem2No.
2982 // `Elem2NbrNo` is the index of the face neighbor (starting from 0),
2983 // and `FElTr.Elem2No` will be offset by the number of (local)
2984 // elements in the mesh.
2985 FElTr.Elem2No = NumOfElements + Elem2NbrNo;
2986 GetFaceNbrElementTransformation(Elem2NbrNo, ElTr2);
2987 FElTr.Elem2 = &ElTr2;
2989 }
2990 else
2991 {
2992 FElTr.Elem2No = -1;
2993 }
2994
2995 // setup the face transformation if the face is not a ghost
2996 if (!is_ghost)
2997 {
2998 GetFaceTransformation(FaceNo, &FElTr);
2999 // NOTE: The above call overwrites FElTr.Loc1
3001 }
3002 else
3003 {
3004 FElTr.SetGeometryType(face_geom);
3005 }
3006
3007 // setup Loc1 & Loc2
3008 int elem_type = GetElementType(face_info.Elem1No);
3009 GetLocalFaceTransformation(face_type, elem_type, FElTr.Loc1.Transf,
3010 face_info.Elem1Inf);
3012
3013 if (fill2)
3014 {
3015 elem_type = face_nbr_elements[Elem2NbrNo]->GetType();
3016 GetLocalFaceTransformation(face_type, elem_type, FElTr.Loc2.Transf,
3017 face_info.Elem2Inf);
3019 }
3020
3021 // adjust Loc1 or Loc2 of the master face if this is a slave face
3022 if (is_slave)
3023 {
3024 if (is_ghost || fill2)
3025 {
3026 // is_ghost -> modify side 1, otherwise -> modify side 2:
3027 ApplyLocalSlaveTransformation(FElTr, face_info, is_ghost);
3028 }
3029 }
3030
3031 // for ghost faces we need a special version of GetFaceTransformation
3032 if (is_ghost)
3033 {
3034 GetGhostFaceTransformation(FElTr, face_type, face_geom);
3036 }
3037
3038 FElTr.SetConfigurationMask(mask);
3039
3040 // This check can be useful for internal debugging, however it will fail on
3041 // periodic boundary faces, so we keep it disabled in general.
3042#if 0
3043#ifdef MFEM_DEBUG
3044 real_t dist = FElTr.CheckConsistency();
3045 if (dist >= 1e-12)
3046 {
3047 mfem::out << "\nInternal error: face id = " << FaceNo
3048 << ", dist = " << dist << ", rank = " << MyRank << '\n';
3049 FElTr.CheckConsistency(1); // print coordinates
3050 MFEM_ABORT("internal error");
3051 }
3052#endif
3053#endif
3054}
3055
3057 FaceElementTransformations &FElTr, Element::Type face_type,
3058 Geometry::Type face_geom) const
3059{
3060 // calculate composition of FElTr.Loc1 and FElTr.Elem1
3061 DenseMatrix &face_pm = FElTr.GetPointMat();
3062 FElTr.Reset();
3063 if (Nodes == NULL)
3064 {
3065 FElTr.Elem1->Transform(FElTr.Loc1.Transf.GetPointMat(), face_pm);
3066 FElTr.SetFE(GetTransformationFEforElementType(face_type));
3067 }
3068 else
3069 {
3070 const FiniteElement* face_el =
3071 Nodes->FESpace()->GetTraceElement(FElTr.Elem1No, face_geom);
3072 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
3073 "Mesh requires nodal Finite Element.");
3074
3075#if 0 // TODO: handle the case of non-interpolatory Nodes
3076 DenseMatrix I;
3077 face_el->Project(Transformation.GetFE(), FElTr.Loc1.Transf, I);
3078 MultABt(Transformation.GetPointMat(), I, pm_face);
3079#else
3080 IntegrationRule eir(face_el->GetDof());
3081 FElTr.Loc1.Transform(face_el->GetNodes(), eir);
3082 Nodes->GetVectorValues(*FElTr.Elem1, eir, face_pm);
3083#endif
3084 FElTr.SetFE(face_el);
3085 }
3086}
3087
3093
3095 int FaceNo, IsoparametricTransformation &ElTr) const
3096{
3097 DenseMatrix &pointmat = ElTr.GetPointMat();
3098 Element *elem = face_nbr_elements[FaceNo];
3099
3100 ElTr.Attribute = elem->GetAttribute();
3101 ElTr.ElementNo = NumOfElements + FaceNo;
3103 ElTr.mesh = this;
3104 ElTr.Reset();
3105
3106 if (Nodes == NULL)
3107 {
3108 const int nv = elem->GetNVertices();
3109 const int *v = elem->GetVertices();
3110
3111 pointmat.SetSize(spaceDim, nv);
3112 for (int k = 0; k < spaceDim; k++)
3113 {
3114 for (int j = 0; j < nv; j++)
3115 {
3116 pointmat(k, j) = face_nbr_vertices[v[j]](k);
3117 }
3118 }
3119
3121 }
3122 else
3123 {
3124 Array<int> vdofs;
3125 ParGridFunction *pNodes = dynamic_cast<ParGridFunction *>(Nodes);
3126 if (pNodes)
3127 {
3128 pNodes->ParFESpace()->GetFaceNbrElementVDofs(FaceNo, vdofs);
3129 int n = vdofs.Size()/spaceDim;
3130 pointmat.SetSize(spaceDim, n);
3131 for (int k = 0; k < spaceDim; k++)
3132 {
3133 for (int j = 0; j < n; j++)
3134 {
3135 pointmat(k,j) = (pNodes->FaceNbrData())(vdofs[n*k+j]);
3136 }
3137 }
3138
3139 ElTr.SetFE(pNodes->ParFESpace()->GetFaceNbrFE(FaceNo));
3140 }
3141 else
3142 {
3143 MFEM_ABORT("Nodes are not ParGridFunction!");
3144 }
3145 }
3146}
3147
3152
3154{
3155 if (Conforming())
3156 {
3157 switch (Dim)
3158 {
3159 case 1: return svert_lvert.Size();
3160 case 2: return sedge_ledge.Size();
3161 default: return sface_lface.Size();
3162 }
3163 }
3164 else
3165 {
3166 MFEM_ASSERT(Dim > 1, "");
3167 const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
3168 return shared.conforming.Size() + shared.slaves.Size();
3169 }
3170}
3171
3172int ParMesh::GetSharedFace(int sface) const
3173{
3174 if (Conforming())
3175 {
3176 switch (Dim)
3177 {
3178 case 1: return svert_lvert[sface];
3179 case 2: return sedge_ledge[sface];
3180 default: return sface_lface[sface];
3181 }
3182 }
3183 else
3184 {
3185 MFEM_ASSERT(Dim > 1, "");
3186 const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
3187 int csize = shared.conforming.Size();
3188 return sface < csize
3189 ? shared.conforming[sface].index
3190 : shared.slaves[sface - csize].index;
3191 }
3192}
3193
3195{
3196 const_cast<ParMesh*>(this)->ExchangeFaceNbrData();
3197 return Mesh::GetNFbyType(type);
3198}
3199
3200// shift cyclically 3 integers a, b, c, so that the smallest of
3201// order[a], order[b], order[c] is first
3202static inline
3203void Rotate3Indirect(int &a, int &b, int &c,
3204 const Array<std::int64_t> &order)
3205{
3206 if (order[a] < order[b])
3207 {
3208 if (order[a] > order[c])
3209 {
3210 ShiftRight(a, b, c);
3211 }
3212 }
3213 else
3214 {
3215 if (order[b] < order[c])
3216 {
3217 ShiftRight(c, b, a);
3218 }
3219 else
3220 {
3221 ShiftRight(a, b, c);
3222 }
3223 }
3224}
3225
3227{
3228 if (Dim != 3 || !(meshgen & 1))
3229 {
3230 return;
3231 }
3232
3233 ResetLazyData();
3234
3235 DSTable *old_v_to_v = NULL;
3236 Table *old_elem_vert = NULL;
3237
3238 if (Nodes)
3239 {
3240 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3241 }
3242
3243 // create a GroupCommunicator over shared vertices
3244 GroupCommunicator svert_comm(gtopo);
3245 GetSharedVertexCommunicator(0, svert_comm);
3246
3247 // communicate the local index of each shared vertex from the group master to
3248 // other ranks in the group
3249 Array<int> svert_master_rank(svert_lvert.Size());
3250 Array<int> svert_master_index(svert_lvert);
3251 for (int i = 0; i < group_svert.Size(); i++)
3252 {
3253 int rank = gtopo.GetGroupMasterRank(i+1);
3254 for (int j = 0; j < group_svert.RowSize(i); j++)
3255 {
3256 svert_master_rank[group_svert.GetRow(i)[j]] = rank;
3257 }
3258 }
3259 svert_comm.Bcast(svert_master_index);
3260
3261 // the pairs (master rank, master local index) define a globally consistent
3262 // vertex ordering
3263 Array<std::int64_t> glob_vert_order(vertices.Size());
3264 {
3265 Array<int> lvert_svert(vertices.Size());
3266 lvert_svert = -1;
3267 for (int i = 0; i < svert_lvert.Size(); i++)
3268 {
3269 lvert_svert[svert_lvert[i]] = i;
3270 }
3271
3272 for (int i = 0; i < vertices.Size(); i++)
3273 {
3274 int s = lvert_svert[i];
3275 if (s >= 0)
3276 {
3277 glob_vert_order[i] =
3278 (std::int64_t(svert_master_rank[s]) << 32) + svert_master_index[s];
3279 }
3280 else
3281 {
3282 glob_vert_order[i] = (std::int64_t(MyRank) << 32) + i;
3283 }
3284 }
3285 }
3286
3287 // rotate tetrahedra so that vertex zero is the lowest (global) index vertex,
3288 // vertex 1 is the second lowest (global) index and vertices 2 and 3 preserve
3289 // positive orientation of the element
3290 for (int i = 0; i < NumOfElements; i++)
3291 {
3293 {
3294 int *v = elements[i]->GetVertices();
3295
3296 Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3297
3298 if (glob_vert_order[v[0]] < glob_vert_order[v[3]])
3299 {
3300 Rotate3Indirect(v[1], v[2], v[3], glob_vert_order);
3301 }
3302 else
3303 {
3304 ShiftRight(v[0], v[1], v[3]);
3305 }
3306 }
3307 }
3308
3309 // rotate also boundary triangles
3310 for (int i = 0; i < NumOfBdrElements; i++)
3311 {
3313 {
3314 int *v = boundary[i]->GetVertices();
3315
3316 Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3317 }
3318 }
3319
3320 const bool check_consistency = true;
3321 if (check_consistency)
3322 {
3323 // create a GroupCommunicator on the shared triangles
3324 GroupCommunicator stria_comm(gtopo);
3325 GetSharedTriCommunicator(0, stria_comm);
3326
3327 Array<int> stria_flag(shared_trias.Size());
3328 for (int i = 0; i < stria_flag.Size(); i++)
3329 {
3330 const int *v = shared_trias[i].v;
3331 if (glob_vert_order[v[0]] < glob_vert_order[v[1]])
3332 {
3333 stria_flag[i] = (glob_vert_order[v[0]] < glob_vert_order[v[2]]) ? 0 : 2;
3334 }
3335 else // v[1] < v[0]
3336 {
3337 stria_flag[i] = (glob_vert_order[v[1]] < glob_vert_order[v[2]]) ? 1 : 2;
3338 }
3339 }
3340
3341 Array<int> stria_master_flag(stria_flag);
3342 stria_comm.Bcast(stria_master_flag);
3343 for (int i = 0; i < stria_flag.Size(); i++)
3344 {
3345 const int *v = shared_trias[i].v;
3346 MFEM_VERIFY(stria_flag[i] == stria_master_flag[i],
3347 "inconsistent vertex ordering found, shared triangle "
3348 << i << ": ("
3349 << v[0] << ", " << v[1] << ", " << v[2] << "), "
3350 << "local flag: " << stria_flag[i]
3351 << ", master flag: " << stria_master_flag[i]);
3352 }
3353 }
3354
3355 // rotate shared triangle faces
3356 for (int i = 0; i < shared_trias.Size(); i++)
3357 {
3358 int *v = shared_trias[i].v;
3359
3360 Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3361 }
3362
3363 // finalize
3364 if (!Nodes)
3365 {
3367 GenerateFaces();
3368 if (el_to_edge)
3369 {
3371 }
3372 }
3373 else
3374 {
3375 DoNodeReorder(old_v_to_v, old_elem_vert);
3376 delete old_elem_vert;
3377 delete old_v_to_v;
3378 }
3379
3380 // the local edge and face numbering is changed therefore we need to
3381 // update sedge_ledge and sface_lface.
3383}
3384
3385void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
3386{
3387 if (pncmesh)
3388 {
3389 MFEM_ABORT("Local and nonconforming refinements cannot be mixed.");
3390 }
3391
3393
3395
3396 if (Dim == 3)
3397 {
3398 int uniform_refinement = 0;
3399 if (type < 0)
3400 {
3401 type = -type;
3402 uniform_refinement = 1;
3403 }
3404
3405 // 1. Hash table of vertex to vertex connections corresponding to refined
3406 // edges.
3407 HashTable<Hashed2> v_to_v;
3408
3409 // 2. Do the red refinement.
3410 switch (type)
3411 {
3412 case 1:
3413 for (int i = 0; i < marked_el.Size(); i++)
3414 {
3415 Bisection(marked_el[i], v_to_v);
3416 }
3417 break;
3418 case 2:
3419 for (int i = 0; i < marked_el.Size(); i++)
3420 {
3421 Bisection(marked_el[i], v_to_v);
3422
3423 Bisection(NumOfElements - 1, v_to_v);
3424 Bisection(marked_el[i], v_to_v);
3425 }
3426 break;
3427 case 3:
3428 for (int i = 0; i < marked_el.Size(); i++)
3429 {
3430 Bisection(marked_el[i], v_to_v);
3431
3432 int j = NumOfElements - 1;
3433 Bisection(j, v_to_v);
3434 Bisection(NumOfElements - 1, v_to_v);
3435 Bisection(j, v_to_v);
3436
3437 Bisection(marked_el[i], v_to_v);
3438 Bisection(NumOfElements-1, v_to_v);
3439 Bisection(marked_el[i], v_to_v);
3440 }
3441 break;
3442 }
3443
3444 // 3. Do the green refinement (to get conforming mesh).
3445 int need_refinement;
3446 int max_faces_in_group = 0;
3447 // face_splittings identify how the shared faces have been split
3448 Array<unsigned> *face_splittings = new Array<unsigned>[GetNGroups()-1];
3449 for (int i = 0; i < GetNGroups()-1; i++)
3450 {
3451 const int faces_in_group = GroupNTriangles(i+1);
3452 face_splittings[i].Reserve(faces_in_group);
3453 if (faces_in_group > max_faces_in_group)
3454 {
3455 max_faces_in_group = faces_in_group;
3456 }
3457 }
3458 int neighbor;
3459 Array<unsigned> iBuf(max_faces_in_group);
3460
3461 MPI_Request *requests = new MPI_Request[GetNGroups()-1];
3462 MPI_Status status;
3463
3464#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3465 int ref_loops_all = 0, ref_loops_par = 0;
3466#endif
3467 do
3468 {
3469 need_refinement = 0;
3470 for (int i = 0; i < NumOfElements; i++)
3471 {
3472 if (elements[i]->NeedRefinement(v_to_v))
3473 {
3474 need_refinement = 1;
3475 Bisection(i, v_to_v);
3476 }
3477 }
3478#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3479 ref_loops_all++;
3480#endif
3481
3482 if (uniform_refinement)
3483 {
3484 continue;
3485 }
3486
3487 // if the mesh is locally conforming start making it globally
3488 // conforming
3489 if (need_refinement == 0)
3490 {
3491#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3492 ref_loops_par++;
3493#endif
3494 // MPI_Barrier(MyComm);
3495 const int tag = 293;
3496
3497 // (a) send the type of interface splitting
3498 int req_count = 0;
3499 for (int i = 0; i < GetNGroups()-1; i++)
3500 {
3501 const int *group_faces = group_stria.GetRow(i);
3502 const int faces_in_group = group_stria.RowSize(i);
3503 // it is enough to communicate through the faces
3504 if (faces_in_group == 0) { continue; }
3505
3506 face_splittings[i].SetSize(0);
3507 for (int j = 0; j < faces_in_group; j++)
3508 {
3509 GetFaceSplittings(shared_trias[group_faces[j]].v, v_to_v,
3510 face_splittings[i]);
3511 }
3512 const int *nbs = gtopo.GetGroup(i+1);
3513 neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
3514 MPI_Isend(face_splittings[i], face_splittings[i].Size(),
3515 MPI_UNSIGNED, neighbor, tag, MyComm,
3516 &requests[req_count++]);
3517 }
3518
3519 // (b) receive the type of interface splitting
3520 for (int i = 0; i < GetNGroups()-1; i++)
3521 {
3522 const int *group_faces = group_stria.GetRow(i);
3523 const int faces_in_group = group_stria.RowSize(i);
3524 if (faces_in_group == 0) { continue; }
3525
3526 const int *nbs = gtopo.GetGroup(i+1);
3527 neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
3528 MPI_Probe(neighbor, tag, MyComm, &status);
3529 int count;
3530 MPI_Get_count(&status, MPI_UNSIGNED, &count);
3531 iBuf.SetSize(count);
3532 MPI_Recv(iBuf, count, MPI_UNSIGNED, neighbor, tag, MyComm,
3533 MPI_STATUS_IGNORE);
3534
3535 for (int j = 0, pos = 0; j < faces_in_group; j++)
3536 {
3537 const int *v = shared_trias[group_faces[j]].v;
3538 need_refinement |= DecodeFaceSplittings(v_to_v, v, iBuf, pos);
3539 }
3540 }
3541
3542 int nr = need_refinement;
3543 MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
3544
3545 MPI_Waitall(req_count, requests, MPI_STATUSES_IGNORE);
3546 }
3547 }
3548 while (need_refinement == 1);
3549
3550#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3551 {
3552 int i = ref_loops_all;
3553 MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
3554 if (MyRank == 0)
3555 {
3556 mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3557 << ref_loops_all << ", ref_loops_par = " << ref_loops_par
3558 << '\n' << endl;
3559 }
3560 }
3561#endif
3562
3563 delete [] requests;
3564 iBuf.DeleteAll();
3565 delete [] face_splittings;
3566
3567 // 4. Update the boundary elements.
3568 do
3569 {
3570 need_refinement = 0;
3571 for (int i = 0; i < NumOfBdrElements; i++)
3572 {
3573 if (boundary[i]->NeedRefinement(v_to_v))
3574 {
3575 need_refinement = 1;
3576 BdrBisection(i, v_to_v);
3577 }
3578 }
3579 }
3580 while (need_refinement == 1);
3581
3582 if (NumOfBdrElements != boundary.Size())
3583 {
3584 mfem_error("ParMesh::LocalRefinement :"
3585 " (NumOfBdrElements != boundary.Size())");
3586 }
3587
3588 ResetLazyData();
3589
3590 const int old_nv = NumOfVertices;
3591 NumOfVertices = vertices.Size();
3592
3593 RefineGroups(old_nv, v_to_v);
3594
3595 // 5. Update the groups after refinement.
3596 if (el_to_face != NULL)
3597 {
3599 GenerateFaces();
3600 }
3601
3602 // 6. Update element-to-edge relations.
3603 if (el_to_edge != NULL)
3604 {
3606 }
3607 } // 'if (Dim == 3)'
3608
3609
3610 if (Dim == 2)
3611 {
3612 int uniform_refinement = 0;
3613 if (type < 0)
3614 {
3615 // type = -type; // not used
3616 uniform_refinement = 1;
3617 }
3618
3619 // 1. Get table of vertex to vertex connections.
3620 DSTable v_to_v(NumOfVertices);
3621 GetVertexToVertexTable(v_to_v);
3622
3623 // 2. Get edge to element connections in arrays edge1 and edge2
3624 int nedges = v_to_v.NumberOfEntries();
3625 int *edge1 = new int[nedges];
3626 int *edge2 = new int[nedges];
3627 int *middle = new int[nedges];
3628
3629 for (int i = 0; i < nedges; i++)
3630 {
3631 edge1[i] = edge2[i] = middle[i] = -1;
3632 }
3633
3634 for (int i = 0; i < NumOfElements; i++)
3635 {
3636 int *v = elements[i]->GetVertices();
3637 for (int j = 0; j < 3; j++)
3638 {
3639 int ind = v_to_v(v[j], v[(j+1)%3]);
3640 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
3641 }
3642 }
3643
3644 // 3. Do the red refinement.
3645 for (int i = 0; i < marked_el.Size(); i++)
3646 {
3647 RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
3648 }
3649
3650 // 4. Do the green refinement (to get conforming mesh).
3651 int need_refinement;
3652 int edges_in_group, max_edges_in_group = 0;
3653 // edge_splittings identify how the shared edges have been split
3654 int **edge_splittings = new int*[GetNGroups()-1];
3655 for (int i = 0; i < GetNGroups()-1; i++)
3656 {
3657 edges_in_group = GroupNEdges(i+1);
3658 edge_splittings[i] = new int[edges_in_group];
3659 if (edges_in_group > max_edges_in_group)
3660 {
3661 max_edges_in_group = edges_in_group;
3662 }
3663 }
3664 int neighbor, *iBuf = new int[max_edges_in_group];
3665
3666 Array<int> group_edges;
3667
3668 MPI_Request request;
3669 MPI_Status status;
3670 Vertex V;
3671 V(2) = 0.0;
3672
3673#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3674 int ref_loops_all = 0, ref_loops_par = 0;
3675#endif
3676 do
3677 {
3678 need_refinement = 0;
3679 for (int i = 0; i < nedges; i++)
3680 {
3681 if (middle[i] != -1 && edge1[i] != -1)
3682 {
3683 need_refinement = 1;
3684 GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
3685 }
3686 }
3687#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3688 ref_loops_all++;
3689#endif
3690
3691 if (uniform_refinement)
3692 {
3693 continue;
3694 }
3695
3696 // if the mesh is locally conforming start making it globally
3697 // conforming
3698 if (need_refinement == 0)
3699 {
3700#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3701 ref_loops_par++;
3702#endif
3703 // MPI_Barrier(MyComm);
3704
3705 // (a) send the type of interface splitting
3706 for (int i = 0; i < GetNGroups()-1; i++)
3707 {
3708 group_sedge.GetRow(i, group_edges);
3709 edges_in_group = group_edges.Size();
3710 // it is enough to communicate through the edges
3711 if (edges_in_group != 0)
3712 {
3713 for (int j = 0; j < edges_in_group; j++)
3714 {
3715 edge_splittings[i][j] =
3716 GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
3717 middle);
3718 }
3719 const int *nbs = gtopo.GetGroup(i+1);
3720 if (nbs[0] == 0)
3721 {
3722 neighbor = gtopo.GetNeighborRank(nbs[1]);
3723 }
3724 else
3725 {
3726 neighbor = gtopo.GetNeighborRank(nbs[0]);
3727 }
3728 MPI_Isend(edge_splittings[i], edges_in_group, MPI_INT,
3729 neighbor, 0, MyComm, &request);
3730 }
3731 }
3732
3733 // (b) receive the type of interface splitting
3734 for (int i = 0; i < GetNGroups()-1; i++)
3735 {
3736 group_sedge.GetRow(i, group_edges);
3737 edges_in_group = group_edges.Size();
3738 if (edges_in_group != 0)
3739 {
3740 const int *nbs = gtopo.GetGroup(i+1);
3741 if (nbs[0] == 0)
3742 {
3743 neighbor = gtopo.GetNeighborRank(nbs[1]);
3744 }
3745 else
3746 {
3747 neighbor = gtopo.GetNeighborRank(nbs[0]);
3748 }
3749 MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
3750 MPI_ANY_TAG, MyComm, &status);
3751
3752 for (int j = 0; j < edges_in_group; j++)
3753 {
3754 if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
3755 {
3756 int *v = shared_edges[group_edges[j]]->GetVertices();
3757 int ii = v_to_v(v[0], v[1]);
3758#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3759 if (middle[ii] != -1)
3760 {
3761 mfem_error("ParMesh::LocalRefinement (triangles) : "
3762 "Oops!");
3763 }
3764#endif
3765 need_refinement = 1;
3766 middle[ii] = NumOfVertices++;
3767 for (int c = 0; c < 2; c++)
3768 {
3769 V(c) = 0.5 * (vertices[v[0]](c) + vertices[v[1]](c));
3770 }
3771 vertices.Append(V);
3772 }
3773 }
3774 }
3775 }
3776
3777 int nr = need_refinement;
3778 MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
3779 }
3780 }
3781 while (need_refinement == 1);
3782
3783#ifdef MFEM_DEBUG_PARMESH_LOCALREF
3784 {
3785 int i = ref_loops_all;
3786 MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
3787 if (MyRank == 0)
3788 {
3789 mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3790 << ref_loops_all << ", ref_loops_par = " << ref_loops_par
3791 << '\n' << endl;
3792 }
3793 }
3794#endif
3795
3796 for (int i = 0; i < GetNGroups()-1; i++)
3797 {
3798 delete [] edge_splittings[i];
3799 }
3800 delete [] edge_splittings;
3801
3802 delete [] iBuf;
3803
3804 // 5. Update the boundary elements.
3805 int v1[2], v2[2], bisect, temp;
3806 temp = NumOfBdrElements;
3807 for (int i = 0; i < temp; i++)
3808 {
3809 int *v = boundary[i]->GetVertices();
3810 bisect = v_to_v(v[0], v[1]);
3811 if (middle[bisect] != -1)
3812 {
3813 // the element was refined (needs updating)
3814 if (boundary[i]->GetType() == Element::SEGMENT)
3815 {
3816 v1[0] = v[0]; v1[1] = middle[bisect];
3817 v2[0] = middle[bisect]; v2[1] = v[1];
3818
3819 boundary[i]->SetVertices(v1);
3820 boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
3821 }
3822 else
3823 {
3824 mfem_error("Only bisection of segment is implemented for bdr"
3825 " elem.");
3826 }
3827 }
3828 }
3829 NumOfBdrElements = boundary.Size();
3830
3831 ResetLazyData();
3832
3833 // 5a. Update the groups after refinement.
3834 RefineGroups(v_to_v, middle);
3835
3836 // 6. Free the allocated memory.
3837 delete [] edge1;
3838 delete [] edge2;
3839 delete [] middle;
3840
3841 if (el_to_edge != NULL)
3842 {
3844 GenerateFaces();
3845 }
3846 } // 'if (Dim == 2)'
3847
3848 if (Dim == 1) // --------------------------------------------------------
3849 {
3850 int cne = NumOfElements, cnv = NumOfVertices;
3851 NumOfVertices += marked_el.Size();
3852 NumOfElements += marked_el.Size();
3853 vertices.SetSize(NumOfVertices);
3854 elements.SetSize(NumOfElements);
3856
3857 for (int j = 0; j < marked_el.Size(); j++)
3858 {
3859 int i = marked_el[j];
3860 Segment *c_seg = (Segment *)elements[i];
3861 int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
3862 int new_v = cnv + j, new_e = cne + j;
3863 AverageVertices(vert, 2, new_v);
3864 elements[new_e] = new Segment(new_v, vert[1], attr);
3865 vert[1] = new_v;
3866
3869 }
3870
3871 static real_t seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
3873 UseExternalData(seg_children, 1, 2, 3);
3874
3875 GenerateFaces();
3876 } // end of 'if (Dim == 1)'
3877
3879 sequence++;
3880
3881 UpdateNodes();
3882
3883#ifdef MFEM_DEBUG
3886#endif
3887}
3888
3890 int nc_limit)
3891{
3892 if (NURBSext)
3893 {
3894 MFEM_ABORT("NURBS meshes are not supported. Please project the "
3895 "NURBS to Nodes first with SetCurvature().");
3896 }
3897
3898 if (!pncmesh)
3899 {
3900 MFEM_ABORT("Can't convert conforming ParMesh to nonconforming ParMesh "
3901 "(you need to initialize the ParMesh from a nonconforming "
3902 "serial Mesh)");
3903 }
3904
3905 ResetLazyData();
3906
3908
3909 // NOTE: no check of !refinements.Size(), in parallel we would have to reduce
3910
3911 // do the refinements
3913 pncmesh->Refine(refinements);
3914
3915 if (nc_limit > 0)
3916 {
3917 pncmesh->LimitNCLevel(nc_limit);
3918 }
3919
3920 // create a second mesh containing the finest elements from 'pncmesh'
3921 ParMesh* pmesh2 = new ParMesh(*pncmesh);
3922 pncmesh->OnMeshUpdated(pmesh2);
3923
3924 attributes.Copy(pmesh2->attributes);
3926
3927 // Copy attribute and bdr_attribute names
3930
3931 // now swap the meshes, the second mesh will become the old coarse mesh
3932 // and this mesh will be the new fine mesh
3933 Mesh::Swap(*pmesh2, false);
3934
3935 delete pmesh2; // NOTE: old face neighbors destroyed here
3936
3938
3940
3942 sequence++;
3943
3944 UpdateNodes();
3945}
3946
3948 real_t threshold, int nc_limit, int op)
3949{
3950 MFEM_VERIFY(pncmesh, "Only supported for non-conforming meshes.");
3951 MFEM_VERIFY(!NURBSext, "Derefinement of NURBS meshes is not supported. "
3952 "Project the NURBS to Nodes first.");
3953
3954 const Table &dt = pncmesh->GetDerefinementTable();
3955
3956 pncmesh->SynchronizeDerefinementData(elem_error, dt);
3957
3958 Array<int> level_ok;
3959 if (nc_limit > 0)
3960 {
3961 pncmesh->CheckDerefinementNCLevel(dt, level_ok, nc_limit);
3962 }
3963
3964 Array<int> derefs;
3965 for (int i = 0; i < dt.Size(); i++)
3966 {
3967 if (nc_limit > 0 && !level_ok[i]) { continue; }
3968
3969 real_t error =
3970 AggregateError(elem_error, dt.GetRow(i), dt.RowSize(i), op);
3971
3972 if (error < threshold) { derefs.Append(i); }
3973 }
3974
3975 long long glob_size = ReduceInt(derefs.Size());
3976 if (!glob_size) { return false; }
3977
3978 // Destroy face-neighbor data only when actually de-refining.
3980
3981 pncmesh->Derefine(derefs);
3982
3983 ParMesh* mesh2 = new ParMesh(*pncmesh);
3984 pncmesh->OnMeshUpdated(mesh2);
3985
3986 attributes.Copy(mesh2->attributes);
3988
3989 // Copy attribute and bdr_attribute names
3992
3993 Mesh::Swap(*mesh2, false);
3994 delete mesh2;
3995
3997
3999
4001 sequence++;
4002
4003 UpdateNodes();
4004
4005 return true;
4006}
4007
4008
4010{
4011 RebalanceImpl(NULL); // default SFC-based partition
4012}
4013
4014void ParMesh::Rebalance(const Array<int> &partition)
4015{
4016 RebalanceImpl(&partition);
4017}
4018
4020{
4021 if (Conforming())
4022 {
4023 MFEM_ABORT("Load balancing is currently not supported for conforming"
4024 " meshes.");
4025 }
4026
4027 if (Nodes)
4028 {
4029 // check that Nodes use a parallel FE space, so we can call UpdateNodes()
4030 MFEM_VERIFY(dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace())
4031 != NULL, "internal error");
4032 }
4033
4035
4036 pncmesh->Rebalance(partition);
4037
4038 ParMesh* pmesh2 = new ParMesh(*pncmesh);
4039 pncmesh->OnMeshUpdated(pmesh2);
4040
4041 attributes.Copy(pmesh2->attributes);
4043
4044 // Copy attribute and bdr_attribute names
4047
4048 Mesh::Swap(*pmesh2, false);
4049 delete pmesh2;
4050
4052
4054
4056 sequence++;
4057
4058 UpdateNodes();
4059}
4060
4061void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
4062{
4063 // Refine groups after LocalRefinement in 2D (triangle meshes)
4064
4065 MFEM_ASSERT(Dim == 2 && meshgen == 1, "internal error");
4066
4067 Array<int> group_verts, group_edges;
4068
4069 // To update the groups after a refinement, we observe that:
4070 // - every (new and old) vertex, edge and face belongs to exactly one group
4071 // - the refinement does not create new groups
4072 // - a new vertex appears only as the middle of a refined edge
4073 // - a face can be refined 2, 3 or 4 times producing new edges and faces
4074
4075 int *I_group_svert, *J_group_svert;
4076 int *I_group_sedge, *J_group_sedge;
4077
4078 I_group_svert = Memory<int>(GetNGroups()+1);
4079 I_group_sedge = Memory<int>(GetNGroups()+1);
4080
4081 I_group_svert[0] = I_group_svert[1] = 0;
4082 I_group_sedge[0] = I_group_sedge[1] = 0;
4083
4084 // overestimate the size of the J arrays
4085 J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4087 J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections());
4088
4089 for (int group = 0; group < GetNGroups()-1; group++)
4090 {
4091 // Get the group shared objects
4092 group_svert.GetRow(group, group_verts);
4093 group_sedge.GetRow(group, group_edges);
4094
4095 // Check which edges have been refined
4096 for (int i = 0; i < group_sedge.RowSize(group); i++)
4097 {
4098 int *v = shared_edges[group_edges[i]]->GetVertices();
4099 const int ind = middle[v_to_v(v[0], v[1])];
4100 if (ind != -1)
4101 {
4102 // add a vertex
4103 group_verts.Append(svert_lvert.Append(ind)-1);
4104 // update the edges
4105 const int attr = shared_edges[group_edges[i]]->GetAttribute();
4106 shared_edges.Append(new Segment(v[1], ind, attr));
4107 group_edges.Append(sedge_ledge.Append(-1)-1);
4108 v[1] = ind;
4109 }
4110 }
4111
4112 I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4113 I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4114
4115 int *J;
4116 J = J_group_svert+I_group_svert[group];
4117 for (int i = 0; i < group_verts.Size(); i++)
4118 {
4119 J[i] = group_verts[i];
4120 }
4121 J = J_group_sedge+I_group_sedge[group];
4122 for (int i = 0; i < group_edges.Size(); i++)
4123 {
4124 J[i] = group_edges[i];
4125 }
4126 }
4127
4129
4130 group_svert.SetIJ(I_group_svert, J_group_svert);
4131 group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4132}
4133
4134void ParMesh::RefineGroups(int old_nv, const HashTable<Hashed2> &v_to_v)
4135{
4136 // Refine groups after LocalRefinement in 3D (tetrahedral meshes)
4137
4138 MFEM_ASSERT(Dim == 3 && meshgen == 1, "internal error");
4139
4140 Array<int> group_verts, group_edges, group_trias;
4141
4142 // To update the groups after a refinement, we observe that:
4143 // - every (new and old) vertex, edge and face belongs to exactly one group
4144 // - the refinement does not create new groups
4145 // - a new vertex appears only as the middle of a refined edge
4146 // - a face can be refined multiple times producing new edges and faces
4147
4148 Array<Segment *> sedge_stack;
4149 Array<Vert3> sface_stack;
4150
4151 Array<int> I_group_svert, J_group_svert;
4152 Array<int> I_group_sedge, J_group_sedge;
4153 Array<int> I_group_stria, J_group_stria;
4154
4155 I_group_svert.SetSize(GetNGroups());
4156 I_group_sedge.SetSize(GetNGroups());
4157 I_group_stria.SetSize(GetNGroups());
4158
4159 I_group_svert[0] = 0;
4160 I_group_sedge[0] = 0;
4161 I_group_stria[0] = 0;
4162
4163 for (int group = 0; group < GetNGroups()-1; group++)
4164 {
4165 // Get the group shared objects
4166 group_svert.GetRow(group, group_verts);
4167 group_sedge.GetRow(group, group_edges);
4168 group_stria.GetRow(group, group_trias);
4169
4170 // Check which edges have been refined
4171 for (int i = 0; i < group_sedge.RowSize(group); i++)
4172 {
4173 int *v = shared_edges[group_edges[i]]->GetVertices();
4174 int ind = v_to_v.FindId(v[0], v[1]);
4175 if (ind == -1) { continue; }
4176
4177 // This shared edge is refined: walk the whole refinement tree
4178 const int attr = shared_edges[group_edges[i]]->GetAttribute();
4179 do
4180 {
4181 ind += old_nv;
4182 // Add new shared vertex
4183 group_verts.Append(svert_lvert.Append(ind)-1);
4184 // Put the right sub-edge on top of the stack
4185 sedge_stack.Append(new Segment(ind, v[1], attr));
4186 // The left sub-edge replaces the original edge
4187 v[1] = ind;
4188 ind = v_to_v.FindId(v[0], ind);
4189 }
4190 while (ind != -1);
4191 // Process all edges in the edge stack
4192 do
4193 {
4194 Segment *se = sedge_stack.Last();
4195 v = se->GetVertices();
4196 ind = v_to_v.FindId(v[0], v[1]);
4197 if (ind == -1)
4198 {
4199 // The edge 'se' is not refined
4200 sedge_stack.DeleteLast();
4201 // Add new shared edge
4202 shared_edges.Append(se);
4203 group_edges.Append(sedge_ledge.Append(-1)-1);
4204 }
4205 else
4206 {
4207 // The edge 'se' is refined
4208 ind += old_nv;
4209 // Add new shared vertex
4210 group_verts.Append(svert_lvert.Append(ind)-1);
4211 // Put the left sub-edge on top of the stack
4212 sedge_stack.Append(new Segment(v[0], ind, attr));
4213 // The right sub-edge replaces the original edge
4214 v[0] = ind;
4215 }
4216 }
4217 while (sedge_stack.Size() > 0);
4218 }
4219
4220 // Check which triangles have been refined
4221 for (int i = 0; i < group_stria.RowSize(group); i++)
4222 {
4223 int *v = shared_trias[group_trias[i]].v;
4224 int ind = v_to_v.FindId(v[0], v[1]);
4225 if (ind == -1) { continue; }
4226
4227 // This shared face is refined: walk the whole refinement tree
4228 const int edge_attr = 1;
4229 do
4230 {
4231 ind += old_nv;
4232 // Add the refinement edge to the edge stack
4233 sedge_stack.Append(new Segment(v[2], ind, edge_attr));
4234 // Put the right sub-triangle on top of the face stack
4235 sface_stack.Append(Vert3(v[1], v[2], ind));
4236 // The left sub-triangle replaces the original one
4237 v[1] = v[0]; v[0] = v[2]; v[2] = ind;
4238 ind = v_to_v.FindId(v[0], v[1]);
4239 }
4240 while (ind != -1);
4241 // Process all faces (triangles) in the face stack
4242 do
4243 {
4244 Vert3 &st = sface_stack.Last();
4245 v = st.v;
4246 ind = v_to_v.FindId(v[0], v[1]);
4247 if (ind == -1)
4248 {
4249 // The triangle 'st' is not refined
4250 // Add new shared face
4251 shared_trias.Append(st);
4252 group_trias.Append(sface_lface.Append(-1)-1);
4253 sface_stack.DeleteLast();
4254 }
4255 else
4256 {
4257 // The triangle 'st' is refined
4258 ind += old_nv;
4259 // Add the refinement edge to the edge stack
4260 sedge_stack.Append(new Segment(v[2], ind, edge_attr));
4261 // Put the left sub-triangle on top of the face stack
4262 sface_stack.Append(Vert3(v[2], v[0], ind));
4263 // Note that the above Append() may invalidate 'v'
4264 v = sface_stack[sface_stack.Size()-2].v;
4265 // The right sub-triangle replaces the original one
4266 v[0] = v[1]; v[1] = v[2]; v[2] = ind;
4267 }
4268 }
4269 while (sface_stack.Size() > 0);
4270 // Process all edges in the edge stack (same code as above)
4271 do
4272 {
4273 Segment *se = sedge_stack.Last();
4274 v = se->GetVertices();
4275 ind = v_to_v.FindId(v[0], v[1]);
4276 if (ind == -1)
4277 {
4278 // The edge 'se' is not refined
4279 sedge_stack.DeleteLast();
4280 // Add new shared edge
4281 shared_edges.Append(se);
4282 group_edges.Append(sedge_ledge.Append(-1)-1);
4283 }
4284 else
4285 {
4286 // The edge 'se' is refined
4287 ind += old_nv;
4288 // Add new shared vertex
4289 group_verts.Append(svert_lvert.Append(ind)-1);
4290 // Put the left sub-edge on top of the stack
4291 sedge_stack.Append(new Segment(v[0], ind, edge_attr));
4292 // The right sub-edge replaces the original edge
4293 v[0] = ind;
4294 }
4295 }
4296 while (sedge_stack.Size() > 0);
4297 }
4298
4299 I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4300 I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4301 I_group_stria[group+1] = I_group_stria[group] + group_trias.Size();
4302
4303 J_group_svert.Append(group_verts);
4304 J_group_sedge.Append(group_edges);
4305 J_group_stria.Append(group_trias);
4306 }
4307
4309
4310 group_svert.SetIJ(I_group_svert, J_group_svert);
4311 group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4312 group_stria.SetIJ(I_group_stria, J_group_stria);
4313 I_group_svert.LoseData(); J_group_svert.LoseData();
4314 I_group_sedge.LoseData(); J_group_sedge.LoseData();
4315 I_group_stria.LoseData(); J_group_stria.LoseData();
4316}
4317
4319{
4320 Array<int> sverts, sedges;
4321
4322 int *I_group_svert, *J_group_svert;
4323 int *I_group_sedge, *J_group_sedge;
4324
4325 I_group_svert = Memory<int>(GetNGroups());
4326 I_group_sedge = Memory<int>(GetNGroups());
4327
4328 I_group_svert[0] = 0;
4329 I_group_sedge[0] = 0;
4330
4331 // compute the size of the J arrays
4332 J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4334 J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections());
4335
4336 for (int group = 0; group < GetNGroups()-1; group++)
4337 {
4338 // Get the group shared objects
4339 group_svert.GetRow(group, sverts);
4340 group_sedge.GetRow(group, sedges);
4341
4342 // Process all the edges
4343 for (int i = 0; i < group_sedge.RowSize(group); i++)
4344 {
4345 int *v = shared_edges[sedges[i]]->GetVertices();
4346 const int ind = old_nv + sedge_ledge[sedges[i]];
4347 // add a vertex
4348 sverts.Append(svert_lvert.Append(ind)-1);
4349 // update the edges
4350 const int attr = shared_edges[sedges[i]]->GetAttribute();
4351 shared_edges.Append(new Segment(v[1], ind, attr));
4352 sedges.Append(sedge_ledge.Append(-1)-1);
4353 v[1] = ind;
4354 }
4355
4356 I_group_svert[group+1] = I_group_svert[group] + sverts.Size();
4357 I_group_sedge[group+1] = I_group_sedge[group] + sedges.Size();
4358
4359 sverts.CopyTo(J_group_svert + I_group_svert[group]);
4360 sedges.CopyTo(J_group_sedge + I_group_sedge[group]);
4361 }
4362
4364
4365 group_svert.SetIJ(I_group_svert, J_group_svert);
4366 group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4367}
4368
4369void ParMesh::UniformRefineGroups3D(int old_nv, int old_nedges,
4370 const DSTable &old_v_to_v,
4371 const STable3D &old_faces,
4372 Array<int> *f2qf)
4373{
4374 // f2qf can be NULL if all faces are quads or there are no quad faces
4375
4376 Array<int> group_verts, group_edges, group_trias, group_quads;
4377
4378 int *I_group_svert, *J_group_svert;
4379 int *I_group_sedge, *J_group_sedge;
4380 int *I_group_stria, *J_group_stria;
4381 int *I_group_squad, *J_group_squad;
4382
4383 I_group_svert = Memory<int>(GetNGroups());
4384 I_group_sedge = Memory<int>(GetNGroups());
4385 I_group_stria = Memory<int>(GetNGroups());
4386 I_group_squad = Memory<int>(GetNGroups());
4387
4388 I_group_svert[0] = 0;
4389 I_group_sedge[0] = 0;
4390 I_group_stria[0] = 0;
4391 I_group_squad[0] = 0;
4392
4393 // compute the size of the J arrays
4394 J_group_svert = Memory<int>(group_svert.Size_of_connections() +
4397 J_group_sedge = Memory<int>(2*group_sedge.Size_of_connections() +
4400 J_group_stria = Memory<int>(4*group_stria.Size_of_connections());
4401 J_group_squad = Memory<int>(4*group_squad.Size_of_connections());
4402
4403 const int oface = old_nv + old_nedges;
4404
4405 for (int group = 0; group < GetNGroups()-1; group++)
4406 {
4407 // Get the group shared objects
4408 group_svert.GetRow(group, group_verts);
4409 group_sedge.GetRow(group, group_edges);
4410 group_stria.GetRow(group, group_trias);
4411 group_squad.GetRow(group, group_quads);
4412
4413 // Process the edges that have been refined
4414 for (int i = 0; i < group_sedge.RowSize(group); i++)
4415 {
4416 int *v = shared_edges[group_edges[i]]->GetVertices();
4417 const int ind = old_nv + old_v_to_v(v[0], v[1]);
4418 // add a vertex
4419 group_verts.Append(svert_lvert.Append(ind)-1);
4420 // update the edges
4421 const int attr = shared_edges[group_edges[i]]->GetAttribute();
4422 shared_edges.Append(new Segment(v[1], ind, attr));
4423 group_edges.Append(sedge_ledge.Append(-1)-1);
4424 v[1] = ind; // v[0] remains the same
4425 }
4426
4427 // Process the triangles that have been refined
4428 for (int i = 0; i < group_stria.RowSize(group); i++)
4429 {
4430 int m[3];
4431 const int stria = group_trias[i];
4432 int *v = shared_trias[stria].v;
4433 // add the refinement edges
4434 m[0] = old_nv + old_v_to_v(v[0], v[1]);
4435 m[1] = old_nv + old_v_to_v(v[1], v[2]);
4436 m[2] = old_nv + old_v_to_v(v[2], v[0]);
4437 const int edge_attr = 1;
4438 shared_edges.Append(new Segment(m[0], m[1], edge_attr));
4439 group_edges.Append(sedge_ledge.Append(-1)-1);
4440 shared_edges.Append(new Segment(m[1], m[2], edge_attr));
4441 group_edges.Append(sedge_ledge.Append(-1)-1);
4442 shared_edges.Append(new Segment(m[0], m[2], edge_attr));
4443 group_edges.Append(sedge_ledge.Append(-1)-1);
4444 // update faces
4445 const int nst = shared_trias.Size();
4446 shared_trias.SetSize(nst+3);
4447 // The above SetSize() may invalidate 'v'
4448 v = shared_trias[stria].v;
4449 shared_trias[nst+0].Set(m[1],m[2],m[0]);
4450 shared_trias[nst+1].Set(m[0],v[1],m[1]);
4451 shared_trias[nst+2].Set(m[2],m[1],v[2]);
4452 v[1] = m[0]; v[2] = m[2]; // v[0] remains the same
4453 group_trias.Append(nst+0);
4454 group_trias.Append(nst+1);
4455 group_trias.Append(nst+2);
4456 // sface_lface is set later
4457 }
4458
4459 // Process the quads that have been refined
4460 for (int i = 0; i < group_squad.RowSize(group); i++)
4461 {
4462 int m[5];
4463 const int squad = group_quads[i];
4464 int *v = shared_quads[squad].v;
4465 const int olf = old_faces(v[0], v[1], v[2], v[3]);
4466 // f2qf can be NULL if all faces are quads
4467 m[0] = oface + (f2qf ? (*f2qf)[olf] : olf);
4468 // add a vertex
4469 group_verts.Append(svert_lvert.Append(m[0])-1);
4470 // add the refinement edges
4471 m[1] = old_nv + old_v_to_v(v[0], v[1]);
4472 m[2] = old_nv + old_v_to_v(v[1], v[2]);
4473 m[3] = old_nv + old_v_to_v(v[2], v[3]);
4474 m[4] = old_nv + old_v_to_v(v[3], v[0]);
4475 const int edge_attr = 1;
4476 shared_edges.Append(new Segment(m[1], m[0], edge_attr));
4477 group_edges.Append(sedge_ledge.Append(-1)-1);
4478 shared_edges.Append(new Segment(m[2], m[0], edge_attr));
4479 group_edges.Append(sedge_ledge.Append(-1)-1);
4480 shared_edges.Append(new Segment(m[3], m[0], edge_attr));
4481 group_edges.Append(sedge_ledge.Append(-1)-1);
4482 shared_edges.Append(new Segment(m[4], m[0], edge_attr));
4483 group_edges.Append(sedge_ledge.Append(-1)-1);
4484 // update faces
4485 const int nsq = shared_quads.Size();
4486 shared_quads.SetSize(nsq+3);
4487 // The above SetSize() may invalidate 'v'
4488 v = shared_quads[squad].v;
4489 shared_quads[nsq+0].Set(m[1],v[1],m[2],m[0]);
4490 shared_quads[nsq+1].Set(m[0],m[2],v[2],m[3]);
4491 shared_quads[nsq+2].Set(m[4],m[0],m[3],v[3]);
4492 v[1] = m[1]; v[2] = m[0]; v[3] = m[4]; // v[0] remains the same
4493 group_quads.Append(nsq+0);
4494 group_quads.Append(nsq+1);
4495 group_quads.Append(nsq+2);
4496 // sface_lface is set later
4497 }
4498
4499 I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
4500 I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
4501 I_group_stria[group+1] = I_group_stria[group] + group_trias.Size();
4502 I_group_squad[group+1] = I_group_squad[group] + group_quads.Size();
4503
4504 group_verts.CopyTo(J_group_svert + I_group_svert[group]);
4505 group_edges.CopyTo(J_group_sedge + I_group_sedge[group]);
4506 group_trias.CopyTo(J_group_stria + I_group_stria[group]);
4507 group_quads.CopyTo(J_group_squad + I_group_squad[group]);
4508 }
4509
4511
4512 group_svert.SetIJ(I_group_svert, J_group_svert);
4513 group_sedge.SetIJ(I_group_sedge, J_group_sedge);
4514 group_stria.SetIJ(I_group_stria, J_group_stria);
4515 group_squad.SetIJ(I_group_squad, J_group_squad);
4516}
4517
4519{
4521
4522 const int old_nv = NumOfVertices;
4523
4524 // call Mesh::UniformRefinement2D so that it won't update the nodes
4525 {
4526 const bool update_nodes = false;
4527 Mesh::UniformRefinement2D_base(update_nodes);
4528 }
4529
4530 // update the groups
4531 UniformRefineGroups2D(old_nv);
4532
4533 UpdateNodes();
4534
4535#ifdef MFEM_DEBUG
4536 // If there are no Nodes, the orientation is checked in the call to
4537 // UniformRefinement2D_base() above.
4538 if (Nodes) { CheckElementOrientation(false); }
4539#endif
4540}
4541
4543{
4545
4546 const int old_nv = NumOfVertices;
4547 const int old_nedges = NumOfEdges;
4548
4549 DSTable v_to_v(NumOfVertices);
4550 GetVertexToVertexTable(v_to_v);
4551 auto faces_tbl = std::unique_ptr<STable3D>(GetFacesTable());
4552
4553 // call Mesh::UniformRefinement3D_base so that it won't update the nodes
4554 Array<int> f2qf;
4555 {
4556 const bool update_nodes = false;
4557 UniformRefinement3D_base(&f2qf, &v_to_v, update_nodes);
4558 // Note: for meshes that have triangular faces, v_to_v is modified by the
4559 // above call to return different edge indices - this is used when
4560 // updating the groups. This is needed by ReorientTetMesh().
4561 }
4562
4563 // update the groups
4564 UniformRefineGroups3D(old_nv, old_nedges, v_to_v, *faces_tbl,
4565 f2qf.Size() ? &f2qf : NULL);
4566
4567 UpdateNodes();
4568}
4569
4571{
4572 if (MyRank == 0)
4573 {
4574 mfem::out << "\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
4575 }
4576}
4577
4579{
4580 if (MyRank == 0)
4581 {
4582 mfem::out << "\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
4583 }
4584}
4585
4586void ParMesh::PrintXG(std::ostream &os) const
4587{
4588 MFEM_ASSERT(Dim == spaceDim, "2D manifolds not supported");
4589 if (Dim == 3 && meshgen == 1)
4590 {
4591 int i, j, nv;
4592 const int *ind;
4593
4594 os << "NETGEN_Neutral_Format\n";
4595 // print the vertices
4596 os << NumOfVertices << '\n';
4597 for (i = 0; i < NumOfVertices; i++)
4598 {
4599 for (j = 0; j < Dim; j++)
4600 {
4601 os << " " << vertices[i](j);
4602 }
4603 os << '\n';
4604 }
4605
4606 // print the elements
4607 os << NumOfElements << '\n';
4608 for (i = 0; i < NumOfElements; i++)
4609 {
4610 nv = elements[i]->GetNVertices();
4611 ind = elements[i]->GetVertices();
4612 os << elements[i]->GetAttribute();
4613 for (j = 0; j < nv; j++)
4614 {
4615 os << " " << ind[j]+1;
4616 }
4617 os << '\n';
4618 }
4619
4620 // print the boundary + shared faces information
4621 os << NumOfBdrElements + sface_lface.Size() << '\n';
4622 // boundary
4623 for (i = 0; i < NumOfBdrElements; i++)
4624 {
4625 nv = boundary[i]->GetNVertices();
4626 ind = boundary[i]->GetVertices();
4627 os << boundary[i]->GetAttribute();
4628 for (j = 0; j < nv; j++)
4629 {
4630 os << " " << ind[j]+1;
4631 }
4632 os << '\n';
4633 }
4634 // shared faces
4635 const int sf_attr =
4636 MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
4637 for (i = 0; i < shared_trias.Size(); i++)
4638 {
4639 ind = shared_trias[i].v;
4640 os << sf_attr;
4641 for (j = 0; j < 3; j++)
4642 {
4643 os << ' ' << ind[j]+1;
4644 }
4645 os << '\n';
4646 }
4647 // There are no quad shared faces
4648 }
4649
4650 if (Dim == 3 && meshgen == 2)
4651 {
4652 int i, j, nv;
4653 const int *ind;
4654
4655 os << "TrueGrid\n"
4656 << "1 " << NumOfVertices << " " << NumOfElements
4657 << " 0 0 0 0 0 0 0\n"
4658 << "0 0 0 1 0 0 0 0 0 0 0\n"
4659 << "0 0 " << NumOfBdrElements+sface_lface.Size()
4660 << " 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
4661 << "0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
4662 << "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
4663
4664 // print the vertices
4665 for (i = 0; i < NumOfVertices; i++)
4666 {
4667 os << i+1 << " 0.0 " << vertices[i](0) << " " << vertices[i](1)
4668 << " " << vertices[i](2) << " 0.0\n";
4669 }
4670
4671 // print the elements
4672 for (i = 0; i < NumOfElements; i++)
4673 {
4674 nv = elements[i]->GetNVertices();
4675 ind = elements[i]->GetVertices();
4676 os << i+1 << " " << elements[i]->GetAttribute();
4677 for (j = 0; j < nv; j++)
4678 {
4679 os << " " << ind[j]+1;
4680 }
4681 os << '\n';
4682 }
4683
4684 // print the boundary information
4685 for (i = 0; i < NumOfBdrElements; i++)
4686 {
4687 nv = boundary[i]->GetNVertices();
4688 ind = boundary[i]->GetVertices();
4689 os << boundary[i]->GetAttribute();
4690 for (j = 0; j < nv; j++)
4691 {
4692 os << " " << ind[j]+1;
4693 }
4694 os << " 1.0 1.0 1.0 1.0\n";
4695 }
4696
4697 // print the shared faces information
4698 const int sf_attr =
4699 MyRank + 1 + (bdr_attributes.Size() > 0 ? bdr_attributes.Max() : 0);
4700 // There are no shared triangle faces
4701 for (i = 0; i < shared_quads.Size(); i++)
4702 {
4703 ind = shared_quads[i].v;
4704 os << sf_attr;
4705 for (j = 0; j < 4; j++)
4706 {
4707 os << ' ' << ind[j]+1;
4708 }
4709 os << " 1.0 1.0 1.0 1.0\n";
4710 }
4711 }
4712
4713 if (Dim == 2)
4714 {
4715 int i, j, attr;
4716 Array<int> v;
4717
4718 os << "areamesh2\n\n";
4719
4720 // print the boundary + shared edges information
4721 os << NumOfBdrElements + shared_edges.Size() << '\n';
4722 // boundary
4723 for (i = 0; i < NumOfBdrElements; i++)
4724 {
4725 attr = boundary[i]->GetAttribute();
4726 boundary[i]->GetVertices(v);
4727 os << attr << " ";
4728 for (j = 0; j < v.Size(); j++)
4729 {
4730 os << v[j] + 1 << " ";
4731 }
4732 os << '\n';
4733 }
4734 // shared edges
4735 for (i = 0; i < shared_edges.Size(); i++)
4736 {
4737 attr = shared_edges[i]->GetAttribute();
4738 shared_edges[i]->GetVertices(v);
4739 os << attr << " ";
4740 for (j = 0; j < v.Size(); j++)
4741 {
4742 os << v[j] + 1 << " ";
4743 }
4744 os << '\n';
4745 }
4746
4747 // print the elements
4748 os << NumOfElements << '\n';
4749 for (i = 0; i < NumOfElements; i++)
4750 {
4751 attr = elements[i]->GetAttribute();
4752 elements[i]->GetVertices(v);
4753
4754 os << attr << " ";
4755 if ((j = GetElementType(i)) == Element::TRIANGLE)
4756 {
4757 os << 3 << " ";
4758 }
4759 else if (j == Element::QUADRILATERAL)
4760 {
4761 os << 4 << " ";
4762 }
4763 else if (j == Element::SEGMENT)
4764 {
4765 os << 2 << " ";
4766 }
4767 for (j = 0; j < v.Size(); j++)
4768 {
4769 os << v[j] + 1 << " ";
4770 }
4771 os << '\n';
4772 }
4773
4774 // print the vertices
4775 os << NumOfVertices << '\n';
4776 for (i = 0; i < NumOfVertices; i++)
4777 {
4778 for (j = 0; j < Dim; j++)
4779 {
4780 os << vertices[i](j) << " ";
4781 }
4782 os << '\n';
4783 }
4784 }
4785 os.flush();
4786}
4787
4789{
4790 // In 2D, this is a workaround for a CPU boundary rendering artifact. We need
4791 // to skip a shared master edge if one of its slaves has the same rank.
4792
4793 const NCMesh::NCList &list = pncmesh->GetEdgeList();
4794 for (int i = master.slaves_begin; i < master.slaves_end; i++)
4795 {
4796 if (!pncmesh->IsGhost(1, list.slaves[i].index)) { return true; }
4797 }
4798 return false;
4799}
4800
4801void ParMesh::Print(std::ostream &os, const std::string &comments) const
4802{
4803 int shared_bdr_attr;
4804 Array<int> nc_shared_faces;
4805
4806 if (NURBSext)
4807 {
4808 Printer(os, "", comments); // does not print shared boundary
4809 return;
4810 }
4811
4812 const Array<int>* s2l_face;
4813 if (!pncmesh)
4814 {
4815 s2l_face = ((Dim == 1) ? &svert_lvert :
4816 ((Dim == 2) ? &sedge_ledge : &sface_lface));
4817 }
4818 else
4819 {
4820 s2l_face = &nc_shared_faces;
4821 if (Dim >= 2)
4822 {
4823 // get a list of all shared non-ghost faces
4824 const NCMesh::NCList& sfaces =
4826 const int nfaces = GetNumFaces();
4827 for (int i = 0; i < sfaces.conforming.Size(); i++)
4828 {
4829 int index = sfaces.conforming[i].index;
4830 if (index < nfaces) { nc_shared_faces.Append(index); }
4831 }
4832 for (int i = 0; i < sfaces.masters.Size(); i++)
4833 {
4834 if (Dim == 2 && WantSkipSharedMaster(sfaces.masters[i])) { continue; }
4835 int index = sfaces.masters[i].index;
4836 if (index < nfaces) { nc_shared_faces.Append(index); }
4837 }
4838 for (int i = 0; i < sfaces.slaves.Size(); i++)
4839 {
4840 int index = sfaces.slaves[i].index;
4841 if (index < nfaces) { nc_shared_faces.Append(index); }
4842 }
4843 }
4844 }
4845
4846 const bool set_names = attribute_sets.SetsExist() ||
4848
4849 os << (!set_names ? "MFEM mesh v1.0\n" : "MFEM mesh v1.3\n");
4850
4851 if (!comments.empty()) { os << '\n' << comments << '\n'; }
4852
4853 // optional
4854 os <<
4855 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
4856 "# POINT = 0\n"
4857 "# SEGMENT = 1\n"
4858 "# TRIANGLE = 2\n"
4859 "# SQUARE = 3\n"
4860 "# TETRAHEDRON = 4\n"
4861 "# CUBE = 5\n"
4862 "# PRISM = 6\n"
4863 "#\n";
4864
4865 os << "\ndimension\n" << Dim
4866 << "\n\nelements\n" << NumOfElements << '\n';
4867 for (int i = 0; i < NumOfElements; i++)
4868 {
4869 PrintElement(elements[i], os);
4870 }
4871
4872 if (set_names)
4873 {
4874 os << "\nattribute_sets\n";
4876 }
4877
4878 int num_bdr_elems = NumOfBdrElements;
4879 if (print_shared && Dim > 1)
4880 {
4881 num_bdr_elems += s2l_face->Size();
4882 }
4883 os << "\nboundary\n" << num_bdr_elems << '\n';
4884 for (int i = 0; i < NumOfBdrElements; i++)
4885 {
4886 PrintElement(boundary[i], os);
4887 }
4888
4889 if (print_shared && Dim > 1)
4890 {
4891 if (bdr_attributes.Size())
4892 {
4893 shared_bdr_attr = bdr_attributes.Max() + MyRank + 1;
4894 }
4895 else
4896 {
4897 shared_bdr_attr = MyRank + 1;
4898 }
4899 for (int i = 0; i < s2l_face->Size(); i++)
4900 {
4901 // Modify the attributes of the faces (not used otherwise?)
4902 faces[(*s2l_face)[i]]->SetAttribute(shared_bdr_attr);
4903 PrintElement(faces[(*s2l_face)[i]], os);
4904 }
4905 }
4906
4907 if (set_names)
4908 {
4909 os << "\nbdr_attribute_sets\n";
4911 }
4912
4913 os << "\nvertices\n" << NumOfVertices << '\n';
4914 if (Nodes == NULL)
4915 {
4916 os << spaceDim << '\n';
4917 for (int i = 0; i < NumOfVertices; i++)
4918 {
4919 os << vertices[i](0);
4920 for (int j = 1; j < spaceDim; j++)
4921 {
4922 os << ' ' << vertices[i](j);
4923 }
4924 os << '\n';
4925 }
4926 os.flush();
4927 }
4928 else
4929 {
4930 os << "\nnodes\n";
4931 Nodes->Save(os);
4932 }
4933
4934 if (set_names)
4935 {
4936 os << "\nmfem_mesh_end" << endl;
4937 }
4938}
4939
4940void ParMesh::Save(const std::string &fname, int precision) const
4941{
4942 ostringstream fname_with_suffix;
4943 fname_with_suffix << fname << "." << setfill('0') << setw(6) << MyRank;
4944 ofstream ofs(fname_with_suffix.str().c_str());
4945 ofs.precision(precision);
4946 Print(ofs);
4947}
4948
4949#ifdef MFEM_USE_ADIOS2
4951{
4952 Mesh::Print(os);
4953}
4954#endif
4955
4956static void dump_element(const Element* elem, Array<int> &data)
4957{
4958 data.Append(elem->GetGeometryType());
4959
4960 int nv = elem->GetNVertices();
4961 const int *v = elem->GetVertices();
4962 for (int i = 0; i < nv; i++)
4963 {
4964 data.Append(v[i]);
4965 }
4966}
4967
4968void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
4969{
4970 int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
4971 const int *v;
4972 MPI_Status status;
4973 Array<real_t> vert;
4974 Array<int> ints;
4975
4976 if (MyRank == 0)
4977 {
4978 os << "MFEM mesh v1.0\n";
4979
4980 if (!comments.empty()) { os << '\n' << comments << '\n'; }
4981
4982 // optional
4983 os <<
4984 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
4985 "# POINT = 0\n"
4986 "# SEGMENT = 1\n"
4987 "# TRIANGLE = 2\n"
4988 "# SQUARE = 3\n"
4989 "# TETRAHEDRON = 4\n"
4990 "# CUBE = 5\n"
4991 "# PRISM = 6\n"
4992 "#\n";
4993
4994 os << "\ndimension\n" << Dim;
4995 }
4996
4997 nv = NumOfElements;
4998 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
4999 if (MyRank == 0)
5000 {
5001 os << "\n\nelements\n" << ne << '\n';
5002 for (i = 0; i < NumOfElements; i++)
5003 {
5004 // processor number + 1 as attribute and geometry type
5005 os << 1 << ' ' << elements[i]->GetGeometryType();
5006 // vertices
5007 nv = elements[i]->GetNVertices();
5008 v = elements[i]->GetVertices();
5009 for (j = 0; j < nv; j++)
5010 {
5011 os << ' ' << v[j];
5012 }
5013 os << '\n';
5014 }
5015 vc = NumOfVertices;
5016 for (p = 1; p < NRanks; p++)
5017 {
5018 MPI_Recv(nv_ne, 2, MPI_INT, p, 444, MyComm, &status);
5019 ints.SetSize(ne);
5020 if (ne)
5021 {
5022 MPI_Recv(&ints[0], ne, MPI_INT, p, 445, MyComm, &status);
5023 }
5024 for (i = 0; i < ne; )
5025 {
5026 // processor number + 1 as attribute and geometry type
5027 os << p+1 << ' ' << ints[i];
5028 // vertices
5029 k = Geometries.GetVertices(ints[i++])->GetNPoints();
5030 for (j = 0; j < k; j++)
5031 {
5032 os << ' ' << vc + ints[i++];
5033 }
5034 os << '\n';
5035 }
5036 vc += nv;
5037 }
5038 }
5039 else
5040 {
5041 // for each element send its geometry type and its vertices
5042 ne = 0;
5043 for (i = 0; i < NumOfElements; i++)
5044 {
5045 ne += 1 + elements[i]->GetNVertices();
5046 }
5047 nv = NumOfVertices;
5048 MPI_Send(nv_ne, 2, MPI_INT, 0, 444, MyComm);
5049
5050 ints.Reserve(ne);
5051 ints.SetSize(0);
5052 for (i = 0; i < NumOfElements; i++)
5053 {
5054 dump_element(elements[i], ints);
5055 }
5056 MFEM_ASSERT(ints.Size() == ne, "");
5057 if (ne)
5058 {
5059 MPI_Send(&ints[0], ne, MPI_INT, 0, 445, MyComm);
5060 }
5061 }
5062
5063 // boundary + shared boundary
5064 ne = NumOfBdrElements;
5065 if (!pncmesh)
5066 {
5067 ne += GetNSharedFaces();
5068 }
5069 else if (Dim > 1)
5070 {
5071 const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
5072 ne += list.conforming.Size() + list.masters.Size() + list.slaves.Size();
5073 // In addition to the number returned by GetNSharedFaces(), include the
5074 // the master shared faces as well.
5075 }
5076 ints.Reserve(ne * (1 + 2*(Dim-1))); // just an upper bound
5077 ints.SetSize(0);
5078
5079 // for each boundary and shared boundary element send its geometry type
5080 // and its vertices
5081 ne = 0;
5082 for (i = j = 0; i < NumOfBdrElements; i++)
5083 {
5084 dump_element(boundary[i], ints); ne++;
5085 }
5086 if (!pncmesh)
5087 {
5088 switch (Dim)
5089 {
5090 case 1:
5091 for (i = 0; i < svert_lvert.Size(); i++)
5092 {
5093 ints.Append(Geometry::POINT);
5094 ints.Append(svert_lvert[i]);
5095 ne++;
5096 }
5097 break;
5098
5099 case 2:
5100 for (i = 0; i < shared_edges.Size(); i++)
5101 {
5102 dump_element(shared_edges[i], ints); ne++;
5103 }
5104 break;
5105
5106 case 3:
5107 for (i = 0; i < shared_trias.Size(); i++)
5108 {
5110 ints.Append(shared_trias[i].v, 3);
5111 ne++;
5112 }
5113 for (i = 0; i < shared_quads.Size(); i++)
5114 {
5116 ints.Append(shared_quads[i].v, 4);
5117 ne++;
5118 }
5119 break;
5120
5121 default:
5122 MFEM_ABORT("invalid dimension: " << Dim);
5123 }
5124 }
5125 else if (Dim > 1)
5126 {
5127 const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
5128 const int nfaces = GetNumFaces();
5129 for (i = 0; i < list.conforming.Size(); i++)
5130 {
5131 int index = list.conforming[i].index;
5132 if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5133 }
5134 for (i = 0; i < list.masters.Size(); i++)
5135 {
5136 int index = list.masters[i].index;
5137 if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5138 }
5139 for (i = 0; i < list.slaves.Size(); i++)
5140 {
5141 int index = list.slaves[i].index;
5142 if (index < nfaces) { dump_element(faces[index], ints); ne++; }
5143 }
5144 }
5145
5146 MPI_Reduce(&ne, &k, 1, MPI_INT, MPI_SUM, 0, MyComm);
5147 if (MyRank == 0)
5148 {
5149 os << "\nboundary\n" << k << '\n';
5150 vc = 0;
5151 for (p = 0; p < NRanks; p++)
5152 {
5153 if (p)
5154 {
5155 MPI_Recv(nv_ne, 2, MPI_INT, p, 446, MyComm, &status);
5156 ints.SetSize(ne);
5157 if (ne)
5158 {
5159 MPI_Recv(ints.GetData(), ne, MPI_INT, p, 447, MyComm, &status);
5160 }
5161 }
5162 else
5163 {
5164 ne = ints.Size();
5165 nv = NumOfVertices;
5166 }
5167 for (i = 0; i < ne; )
5168 {
5169 // processor number + 1 as bdr. attr. and bdr. geometry type
5170 os << p+1 << ' ' << ints[i];
5171 k = Geometries.NumVerts[ints[i++]];
5172 // vertices
5173 for (j = 0; j < k; j++)
5174 {
5175 os << ' ' << vc + ints[i++];
5176 }
5177 os << '\n';
5178 }
5179 vc += nv;
5180 }
5181 }
5182 else
5183 {
5184 nv = NumOfVertices;
5185 ne = ints.Size();
5186 MPI_Send(nv_ne, 2, MPI_INT, 0, 446, MyComm);
5187 if (ne)
5188 {
5189 MPI_Send(ints.GetData(), ne, MPI_INT, 0, 447, MyComm);
5190 }
5191 }
5192
5193 // vertices / nodes
5194 MPI_Reduce(&NumOfVertices, &nv, 1, MPI_INT, MPI_SUM, 0, MyComm);
5195 if (MyRank == 0)
5196 {
5197 os << "\nvertices\n" << nv << '\n';
5198 }
5199 if (Nodes == NULL)
5200 {
5201 if (MyRank == 0)
5202 {
5203 os << spaceDim << '\n';
5204 for (i = 0; i < NumOfVertices; i++)
5205 {
5206 os << vertices[i](0);
5207 for (j = 1; j < spaceDim; j++)
5208 {
5209 os << ' ' << vertices[i](j);
5210 }
5211 os << '\n';
5212 }
5213 for (p = 1; p < NRanks; p++)
5214 {
5215 MPI_Recv(&nv, 1, MPI_INT, p, 448, MyComm, &status);
5216 vert.SetSize(nv*spaceDim);
5217 if (nv)
5218 {
5219 MPI_Recv(&vert[0], nv*spaceDim, MPITypeMap<real_t>::mpi_type, p, 449, MyComm,
5220 &status);
5221 }
5222 for (i = 0; i < nv; i++)
5223 {
5224 os << vert[i*spaceDim];
5225 for (j = 1; j < spaceDim; j++)
5226 {
5227 os << ' ' << vert[i*spaceDim+j];
5228 }
5229 os << '\n';
5230 }
5231 }
5232 os.flush();
5233 }
5234 else
5235 {
5236 MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 448, MyComm);
5238 for (i = 0; i < NumOfVertices; i++)
5239 {
5240 for (j = 0; j < spaceDim; j++)
5241 {
5242 vert[i*spaceDim+j] = vertices[i](j);
5243 }
5244 }
5245 if (NumOfVertices)
5246 {
5247 MPI_Send(&vert[0], NumOfVertices*spaceDim, MPITypeMap<real_t>::mpi_type, 0, 449,
5248 MyComm);
5249 }
5250 }
5251 }
5252 else
5253 {
5254 if (MyRank == 0)
5255 {
5256 os << "\nnodes\n";
5257 }
5258 ParGridFunction *pnodes = dynamic_cast<ParGridFunction *>(Nodes);
5259 if (pnodes)
5260 {
5261 pnodes->SaveAsOne(os);
5262 }
5263 else
5264 {
5265 ParFiniteElementSpace *pfes =
5266 dynamic_cast<ParFiniteElementSpace *>(Nodes->FESpace());
5267 if (pfes)
5268 {
5269 // create a wrapper ParGridFunction
5270 ParGridFunction ParNodes(pfes, Nodes);
5271 ParNodes.SaveAsOne(os);
5272 }
5273 else
5274 {
5275 mfem_error("ParMesh::PrintAsOne : Nodes have no parallel info!");
5276 }
5277 }
5278 }
5279}
5280
5281void ParMesh::PrintAsSerial(std::ostream &os, const std::string &comments) const
5282{
5283 int save_rank = 0;
5284 Mesh serialmesh = GetSerialMesh(save_rank);
5285 if (MyRank == save_rank)
5286 {
5287 serialmesh.Printer(os, "", comments);
5288 }
5289 MPI_Barrier(MyComm);
5290}
5291
5292Mesh ParMesh::GetSerialMesh(int save_rank) const
5293{
5294 if (pncmesh || NURBSext)
5295 {
5296 MFEM_ABORT("Nonconforming meshes and NURBS meshes are not yet supported.");
5297 }
5298
5299 // Define linear H1 space for vertex numbering
5300 H1_FECollection fec_linear(1, Dim);
5301 ParMesh *pm = const_cast<ParMesh *>(this);
5302 ParFiniteElementSpace pfespace_linear(pm, &fec_linear);
5303
5304 long long ne_glob_l = GetGlobalNE(); // needs to be called by all ranks
5305 MFEM_VERIFY(int(ne_glob_l) == ne_glob_l,
5306 "overflow in the number of elements!");
5307 int ne_glob = (save_rank == MyRank) ? int(ne_glob_l) : 0;
5308
5309 long long nvertices = pfespace_linear.GetTrueVSize();
5310 long long nvertices_glob_l = 0;
5311 MPI_Reduce(&nvertices, &nvertices_glob_l, 1, MPI_LONG_LONG, MPI_SUM,
5312 save_rank, MyComm);
5313 int nvertices_glob = int(nvertices_glob_l);
5314 MFEM_VERIFY(nvertices_glob == nvertices_glob_l,
5315 "overflow in the number of vertices!");
5316
5317 long