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