MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
ncmesh.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 "mesh_headers.hpp"
14#include "../general/text.hpp"
15
16#include <string>
17#include <cmath>
18#include <map>
19
20#include "ncmesh_tables.hpp"
21
22
23namespace
24{
25/**
26 * @brief Base case of convenience variadic max function.
27 *
28 * @tparam T Base type
29 * @param arg Recursion base value
30 * @return T value to max over
31 */
32template<typename T>
33T max(T&& arg)
34{
35 return arg;
36}
37/**
38 * @brief Convenience variadic max function.
39 *
40 * @tparam T Base Type
41 * @tparam Ts Parameter pack of other types
42 * @param arg Singular argument
43 * @param args Pack of arguments
44 * @return T maximum value
45 */
46template<typename T, typename... Ts>
47T max(T arg, Ts... args)
48{
49 return std::max(std::forward<T>(arg), max(args...));
50}
51} // namespace
52namespace mfem
53{
54
55NCMesh::GeomInfo NCMesh::GI[Geometry::NumGeom];
56
58{
59 if (initialized) { return; }
60
61 auto elem = [&]()
62 {
63 switch (geom)
64 {
65 case Geometry::CUBE: return std::unique_ptr<mfem::Element>(new Hexahedron);
66 case Geometry::PRISM: return std::unique_ptr<mfem::Element>(new Wedge);
67 case Geometry::TETRAHEDRON: return std::unique_ptr<mfem::Element>
68 (new Tetrahedron);
69 case Geometry::PYRAMID: return std::unique_ptr<mfem::Element>(new Pyramid);
70 case Geometry::SQUARE: return std::unique_ptr<mfem::Element>(new Quadrilateral);
71 case Geometry::TRIANGLE: return std::unique_ptr<mfem::Element>(new Triangle);
72 case Geometry::SEGMENT: return std::unique_ptr<mfem::Element>(new Segment);
73 default: MFEM_ABORT("unsupported geometry " << geom);
74 }
75 }();
76
77 nv = elem->GetNVertices();
78 ne = elem->GetNEdges();
79 nf = elem->GetNFaces();
80 for (int i = 0; i < ne; i++)
81 {
82 for (int j = 0; j < 2; j++)
83 {
84 edges[i][j] = elem->GetEdgeVertices(i)[j];
85 }
86 }
87 for (int i = 0; i < nf; i++)
88 {
89 nfv[i] = elem->GetNFaceVertices(i);
90
91 faces[i][3] = 7; // invalid node index for 3-node faces
92 for (int j = 0; j < nfv[i]; j++)
93 {
94 faces[i][j] = elem->GetFaceVertices(i)[j];
95 }
96 }
97
98 // in 1D/2D we pretend to have faces too, so we can use NCMesh::Face::elem[2]
99 if (!nf)
100 {
101 if (ne)
102 {
103 for (int i = 0; i < ne; i++)
104 {
105 // make a degenerate face
106 faces[i][0] = faces[i][1] = edges[i][0];
107 faces[i][2] = faces[i][3] = edges[i][1];
108 nfv[i] = 2;
109 }
110 nf = ne;
111 }
112 else
113 {
114 for (int i = 0; i < nv; i++)
115 {
116 // 1D degenerate face
117 faces[i][0] = faces[i][1] = faces[i][2] = faces[i][3] = i;
118 nfv[i] = 1;
119 }
120 nf = nv;
121 }
122 }
123
124 initialized = true;
125}
126
128 : shadow(1024, 2048)
129{
130 Dim = mesh->Dimension();
131 spaceDim = mesh->SpaceDimension();
132 MyRank = 0;
133 Iso = true;
134 Legacy = false;
135
136 // create the NCMesh::Element struct for each Mesh element
137 for (int i = 0; i < mesh->GetNE(); i++)
138 {
139 const mfem::Element *elem = mesh->GetElement(i);
140
141 Geometry::Type geom = elem->GetGeometryType();
142 CheckSupportedGeom(geom);
143 GI[geom].InitGeom(geom);
144
145 // if we have pyramids we will need tets after refinement
146 if (geom == Geometry::PYRAMID)
147 {
149 }
150
151 // create NCMesh::Element for this mfem::Element
152 int root_id = AddElement(geom, elem->GetAttribute());
153 MFEM_ASSERT(root_id == i, "");
154 Element &root_elem = elements[root_id];
155
156 const int *v = elem->GetVertices();
157 for (int j = 0; j < GI[geom].nv; j++)
158 {
159 int id = v[j];
160 root_elem.node[j] = id;
161 nodes.Alloc(id, id, id);
162 // NOTE: top-level nodes are special: id == p1 == p2 == orig. vertex id
163 }
164 }
165
166 // if the user initialized any hanging nodes with Mesh::AddVertexParents,
167 // copy the hierarchy now
168 if (mesh->tmp_vertex_parents.Size())
169 {
170 for (const auto &triple : mesh->tmp_vertex_parents)
171 {
172 nodes.Reparent(triple.one, triple.two, triple.three);
173 }
174 }
175
176 // create edge nodes and faces
177 nodes.UpdateUnused();
178 for (int i = 0; i < elements.Size(); i++)
179 {
180 // increase reference count of all nodes the element is using (NOTE: this
181 // will also create and reference all edge nodes and faces)
183
184 // make links from faces back to the element
185 RegisterFaces(i);
186 }
187
188 // store boundary element attributes
189 for (int i = 0; i < mesh->GetNBE(); i++)
190 {
191 const mfem::Element *be = mesh->GetBdrElement(i);
192 const int *v = be->GetVertices();
193
194 Face* face = NULL;
195 switch (be->GetType())
196 {
198 face = faces.Find(v[0], v[1], v[2], v[3]);
199 break;
201 face = faces.Find(v[0], v[1], v[2]);
202 break;
204 face = faces.Find(v[0], v[0], v[1], v[1]);
205 break;
207 face = faces.Find(v[0], v[0], v[0], v[0]);
208 break;
209 default:
210 MFEM_ABORT("Unsupported boundary element geometry.");
211 }
212
213 MFEM_VERIFY(face, "Boundary face not found.");
214 face->attribute = be->GetAttribute();
215 }
216
217 // copy top-level vertex coordinates (leave empty if the mesh is curved)
218 if (!mesh->Nodes)
219 {
220 coordinates.SetSize(3*mesh->GetNV());
221 for (int i = 0; i < mesh->GetNV(); i++)
222 {
223 std::memcpy(&coordinates[3*i], mesh->GetVertex(i), 3*sizeof(real_t));
224 }
225 }
226
227 InitRootState(mesh->GetNE());
229
230 Update();
231}
232
234 : Dim(other.Dim)
235 , spaceDim(other.spaceDim)
236 , MyRank(other.MyRank)
237 , Iso(other.Iso)
238 , Geoms(other.Geoms)
239 , Legacy(other.Legacy)
240 , nodes(other.nodes)
241 , faces(other.faces)
242 , using_scaling(other.using_scaling)
243 , elements(other.elements)
244 , free_element_ids(other.free_element_ids)
245 , root_state(other.root_state)
246 , coordinates(other.coordinates)
247 , NEdges(other.NEdges)
248 , NFaces(other.NFaces)
249 , NGhostEdges(other.NGhostEdges)
250 , NGhostFaces(other.NGhostFaces)
251 , boundary_faces(other.boundary_faces)
252 , face_geom(other.face_geom)
253 , element_vertex(other.element_vertex)
254 , shadow(1024, 2048)
255{
256 Update();
257}
258
260{
261 Geoms = 0;
262 for (int i = 0; i < root_state.Size(); i++)
263 {
264 Geoms |= (1 << elements[i].Geom());
265 }
266}
267
269{
272
276
278}
279
281{
282#ifdef MFEM_DEBUG
283 // sign off of all faces and nodes
284 Array<int> elemFaces;
285 for (int i = 0; i < elements.Size(); i++)
286 {
287 if (elements[i].IsLeaf())
288 {
289 elemFaces.SetSize(0);
290 UnreferenceElement(i, elemFaces);
291 DeleteUnusedFaces(elemFaces);
292 }
293 }
294 // NOTE: in release mode, we just throw away all faces and nodes at once
295#endif
296}
297
298void NCMesh::Node::SetScale(real_t s, bool overwrite)
299{
300 MFEM_ASSERT(0.0 < s && s < 1.0, "Invalid scale");
301 if (!overwrite && scaleSet && std::abs((s - scale) / scale) > scaleTol)
302 {
303 MFEM_ABORT("Node scale is already set (inconsistent refinement)");
304 }
305 else if (overwrite || !scaleSet)
306 {
307 scale = s;
308 scaleSet = true;
309 }
310}
311
313{
314 MFEM_ASSERT(!vert_refc && !edge_refc, "node was not unreferenced properly, "
315 "vert_refc: " << (int) vert_refc << ", edge_refc: "
316 << (int) edge_refc);
317}
318
319void NCMesh::ReparentNode(int node, int new_p1, int new_p2, real_t scale)
320{
321 Node &nd = nodes[node];
322 nd.SetScale(GetScale(scale, new_p1 > new_p2), true);
323 int old_p1 = nd.p1, old_p2 = nd.p2;
324
325 // assign new parents
326 nodes.Reparent(node, new_p1, new_p2);
327
328 MFEM_ASSERT(shadow.FindId(old_p1, old_p2) < 0,
329 "shadow node already exists");
330
331 // store old parent pair temporarily in 'shadow'
332 int sh = shadow.GetId(old_p1, old_p2);
333 shadow[sh].vert_index = node;
334}
335
336int NCMesh::FindMidEdgeNode(int node1, int node2) const
337{
338 int mid = nodes.FindId(node1, node2);
339 if (mid < 0 && shadow.Size())
340 {
341 // if (anisotropic) refinement is underway, some nodes may temporarily be
342 // available under alternate parents (see ReparentNode)
343 mid = shadow.FindId(node1, node2);
344 if (mid >= 0)
345 {
346 mid = shadow[mid].vert_index; // index of the original node
347 }
348 }
349 return mid;
350}
351
352int NCMesh::GetMidEdgeNode(int node1, int node2)
353{
354 int mid = FindMidEdgeNode(node1, node2);
355 if (mid < 0) { mid = nodes.GetId(node1, node2); } // create if not found
356 return mid;
357}
358
359int NCMesh::GetMidFaceNode(int en1, int en2, int en3, int en4)
360{
361 // mid-face node can be created either from (en1, en3) or from (en2, en4)
362 int midf = FindMidEdgeNode(en1, en3);
363 if (midf >= 0) { return midf; }
364 return nodes.GetId(en2, en4);
365}
366
368{
369 const Element &el = elements[elem];
370 const int* node = el.node;
371 GeomInfo& gi = GI[el.Geom()];
372
373 // reference all vertices
374 for (int i = 0; i < gi.nv; i++)
375 {
376 nodes[node[i]].vert_refc++;
377 }
378
379 // reference all edges (possibly creating their nodes)
380 for (int i = 0; i < gi.ne; i++)
381 {
382 const int* ev = gi.edges[i];
383 nodes.Get(node[ev[0]], node[ev[1]])->edge_refc++;
384 }
385
386 // get all faces (possibly creating them)
387 for (int i = 0; i < gi.nf; i++)
388 {
389 const int* fv = gi.faces[i];
390 faces.GetId(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
391
392 // NOTE: face->RegisterElement called separately to avoid having to store
393 // 3 element indices temporarily in the face when refining. See also
394 // NCMesh::RegisterFaces.
395 }
396}
397
398void NCMesh::UnreferenceElement(int elem, Array<int> &elemFaces)
399{
400 Element &el = elements[elem];
401 int* node = el.node;
402 GeomInfo& gi = GI[el.Geom()];
403
404 // unreference all faces
405 for (int i = 0; i < gi.nf; i++)
406 {
407 const int* fv = gi.faces[i];
408 int face = faces.FindId(node[fv[0]], node[fv[1]],
409 node[fv[2]], node[fv[3]]);
410 MFEM_ASSERT(face >= 0, "face not found.");
411 faces[face].ForgetElement(elem);
412
413 // NOTE: faces.Delete() called later to avoid destroying and recreating
414 // faces during refinement, see NCMesh::DeleteUnusedFaces.
415 elemFaces.Append(face);
416 }
417
418 // unreference all edges (possibly destroying them)
419 for (int i = 0; i < gi.ne; i++)
420 {
421 const int* ev = gi.edges[i];
422 int enode = FindMidEdgeNode(node[ev[0]], node[ev[1]]);
423 MFEM_ASSERT(enode >= 0, "edge not found.");
424 MFEM_ASSERT(nodes.IdExists(enode), "edge does not exist.");
425 if (!nodes[enode].UnrefEdge())
426 {
427 nodes.Delete(enode);
428 }
429 }
430
431 // unreference all vertices (possibly destroying them)
432 for (int i = 0; i < gi.nv; i++)
433 {
434 if (!nodes[node[i]].UnrefVertex())
435 {
436 nodes.Delete(node[i]);
437 }
438 }
439}
440
441void NCMesh::RegisterFaces(int elem, int* fattr)
442{
443 Element &el = elements[elem];
444 GeomInfo &gi = GI[el.Geom()];
445
446 for (int i = 0; i < gi.nf; i++)
447 {
448 Face* face = GetFace(el, i);
449 MFEM_ASSERT(face, "face not found.");
450 face->RegisterElement(elem);
451 if (fattr) { face->attribute = fattr[i]; }
452 }
453}
454
456{
457 for (int i = 0; i < elemFaces.Size(); i++)
458 {
459 if (faces[elemFaces[i]].Unused())
460 {
461 faces.Delete(elemFaces[i]);
462 }
463 }
464}
465
467{
468 if (elem[0] < 0) { elem[0] = e; }
469 else if (elem[1] < 0) { elem[1] = e; }
470 else { MFEM_ABORT("can't have 3 elements in Face::elem[]."); }
471}
472
474{
475 if (elem[0] == e) { elem[0] = -1; }
476 else if (elem[1] == e) { elem[1] = -1; }
477 else { MFEM_ABORT("element " << e << " not found in Face::elem[]."); }
478}
479
481{
482 GeomInfo& gi = GI[(int) elem.geom];
483 const int* fv = gi.faces[face_no];
484 int* node = elem.node;
485 return faces.Find(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
486}
487
489{
490 if (elem[0] >= 0)
491 {
492 MFEM_ASSERT(elem[1] < 0, "not a single element face.");
493 return elem[0];
494 }
495 else
496 {
497 MFEM_ASSERT(elem[1] >= 0, "no elements in face.");
498 return elem[1];
499 }
500}
501
502
503//// Refinement ////////////////////////////////////////////////////////////////
504
505void Refinement::SetScale(const ScaledType &r)
506{
507 switch (r.first)
508 {
509 case X:
510 s[0] = r.second;
511 break;
512 case Y:
513 s[1] = r.second;
514 break;
515 case Z:
516 s[2] = r.second;
517 break;
518 case XY:
519 s[0] = r.second;
520 s[1] = r.second;
521 break;
522 case YZ:
523 s[1] = r.second;
524 s[2] = r.second;
525 break;
526 case XZ:
527 s[0] = r.second;
528 s[2] = r.second;
529 break;
530 case XYZ:
531 s[0] = r.second;
532 s[1] = r.second;
533 s[2] = r.second;
534 break;
535 default:
536 MFEM_ABORT("Unsupported refinement type.");
537 }
538}
539
541{
542 for (int i=0; i<3; ++i)
543 if (s[i] > real_t{0})
544 s[i] = scale[i];
545}
546
548 : index(index)
549{
550 for (int i=0; i<3; ++i) { s[i] = 0.0; }
551 // Default case is XYZ type with scale 0.5.
552 SetScale(ScaledType(XYZ, 0.5));
553}
554
555Refinement::Refinement(int index, const std::initializer_list<ScaledType> &refs)
556 : index(index)
557{
558 for (int i=0; i<3; ++i) { s[i] = 0.0; }
559 if (refs.size() == 0)
560 {
561 // Default case is XYZ type with scale 0.5.
562 SetScale(ScaledType(XYZ, 0.5));
563 }
564 else
565 {
566 for (const auto & ref : refs)
567 {
568 SetScale(ref);
569 }
570 }
571}
572
573Refinement::Refinement(int index, char type, real_t scale)
574 : index(index)
575{
576 for (int i=0; i<3; ++i) { s[i] = 0.0; }
577 SetScale(ScaledType(type, scale));
578}
579
581{
582 char t{0}; // Set the X, Y or Z bit
583 for (int i = 0; i < 3; ++i)
584 if (s[i] > real_t{0})
585 t |= (1 << i);
586 return t;
587}
588
589void Refinement::Set(int element, char type, real_t scale)
590{
591 index = element;
592 for (int i=0; i<3; ++i) { s[i] = 0.0; }
593 SetScale(ScaledType(type, scale));
594}
595
596void Refinement::SetType(char type, real_t scale)
597{
598 for (int i=0; i<3; ++i) { s[i] = 0.0; }
599 SetScale(ScaledType(type, scale));
600}
601
603 : geom(geom), ref_type(0), tet_type(0), flag(0), index(-1)
604 , rank(0), attribute(attr), parent(-1)
605{
606 for (int i = 0; i < MaxElemNodes; i++) { node[i] = -1; }
607 for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
608
609 // NOTE: in 2D the 8/10-element node/child arrays are not optimal, however,
610 // testing shows we would only save 17% of the total NCMesh memory if
611 // 4-element arrays were used (e.g. through templates); we thus prefer to
612 // keep the code as simple as possible.
613}
614
615int NCMesh::NewHexahedron(int n0, int n1, int n2, int n3,
616 int n4, int n5, int n6, int n7,
617 int attr,
618 int fattr0, int fattr1, int fattr2,
619 int fattr3, int fattr4, int fattr5)
620{
621 // create new element, initialize nodes
622 int new_id = AddElement(Geometry::CUBE, attr);
623 Element &el = elements[new_id];
624
625 el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
626 el.node[4] = n4, el.node[5] = n5, el.node[6] = n6, el.node[7] = n7;
627
628 // get faces and assign face attributes
630 const GeomInfo &gi_hex = GI[Geometry::CUBE];
631 for (int i = 0; i < gi_hex.nf; i++)
632 {
633 const int* fv = gi_hex.faces[i];
634 f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
635 el.node[fv[2]], el.node[fv[3]]);
636 }
637
638 f[0]->attribute = fattr0, f[1]->attribute = fattr1;
639 f[2]->attribute = fattr2, f[3]->attribute = fattr3;
640 f[4]->attribute = fattr4, f[5]->attribute = fattr5;
641
642 return new_id;
643}
644
645int NCMesh::NewWedge(int n0, int n1, int n2,
646 int n3, int n4, int n5,
647 int attr,
648 int fattr0, int fattr1,
649 int fattr2, int fattr3, int fattr4)
650{
651 // create new element, initialize nodes
652 int new_id = AddElement(Geometry::PRISM, attr);
653 Element &el = elements[new_id];
654
655 el.node[0] = n0, el.node[1] = n1, el.node[2] = n2;
656 el.node[3] = n3, el.node[4] = n4, el.node[5] = n5;
657
658 // get faces and assign face attributes
659 Face* f[5];
660 const GeomInfo &gi_wedge = GI[Geometry::PRISM];
661 for (int i = 0; i < gi_wedge.nf; i++)
662 {
663 const int* fv = gi_wedge.faces[i];
664 f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
665 el.node[fv[2]], el.node[fv[3]]);
666 }
667
668 f[0]->attribute = fattr0;
669 f[1]->attribute = fattr1;
670 f[2]->attribute = fattr2;
671 f[3]->attribute = fattr3;
672 f[4]->attribute = fattr4;
673
674 return new_id;
675}
676
677int NCMesh::NewTetrahedron(int n0, int n1, int n2, int n3, int attr,
678 int fattr0, int fattr1, int fattr2, int fattr3)
679{
680 // create new element, initialize nodes
681 int new_id = AddElement(Geometry::TETRAHEDRON, attr);
682 Element &el = elements[new_id];
683
684 el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
685
686 // get faces and assign face attributes
687 Face* f[4];
688 const GeomInfo &gi_tet = GI[Geometry::TETRAHEDRON];
689 for (int i = 0; i < gi_tet.nf; i++)
690 {
691 const int* fv = gi_tet.faces[i];
692 f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]], el.node[fv[2]]);
693 }
694
695 f[0]->attribute = fattr0;
696 f[1]->attribute = fattr1;
697 f[2]->attribute = fattr2;
698 f[3]->attribute = fattr3;
699
700 return new_id;
701}
702int NCMesh::NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr,
703 int fattr0, int fattr1, int fattr2, int fattr3,
704 int fattr4)
705{
706 // create new element, initialize nodes
707 int new_id = AddElement(Geometry::PYRAMID, attr);
708 Element &el = elements[new_id];
709
710 el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
711 el.node[4] = n4;
712
713 // get faces and assign face attributes
714 Face* f[5];
715 const GeomInfo &gi_pyr = GI[Geometry::PYRAMID];
716 for (int i = 0; i < gi_pyr.nf; i++)
717 {
718 const int* fv = gi_pyr.faces[i];
719 f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
720 el.node[fv[2]], el.node[fv[3]]);
721 }
722
723 f[0]->attribute = fattr0;
724 f[1]->attribute = fattr1;
725 f[2]->attribute = fattr2;
726 f[3]->attribute = fattr3;
727 f[4]->attribute = fattr4;
728
729 return new_id;
730}
731
732int NCMesh::NewQuadrilateral(int n0, int n1, int n2, int n3,
733 int attr,
734 int eattr0, int eattr1, int eattr2, int eattr3)
735{
736 // create new element, initialize nodes
737 int new_id = AddElement(Geometry::SQUARE, attr);
738 Element &el = elements[new_id];
739
740 el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
741
742 // get (degenerate) faces and assign face attributes
743 Face* f[4];
744 const GeomInfo &gi_quad = GI[Geometry::SQUARE];
745 for (int i = 0; i < gi_quad.nf; i++)
746 {
747 const int* fv = gi_quad.faces[i];
748 f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
749 el.node[fv[2]], el.node[fv[3]]);
750 }
751
752 f[0]->attribute = eattr0, f[1]->attribute = eattr1;
753 f[2]->attribute = eattr2, f[3]->attribute = eattr3;
754
755 return new_id;
756}
757
758int NCMesh::NewTriangle(int n0, int n1, int n2,
759 int attr, int eattr0, int eattr1, int eattr2)
760{
761 // create new element, initialize nodes
762 int new_id = AddElement(Geometry::TRIANGLE, attr);
763 Element &el = elements[new_id];
764
765 el.node[0] = n0, el.node[1] = n1, el.node[2] = n2;
766
767 // get (degenerate) faces and assign face attributes
768 Face* f[3];
769 const GeomInfo &gi_tri = GI[Geometry::TRIANGLE];
770 for (int i = 0; i < gi_tri.nf; i++)
771 {
772 const int* fv = gi_tri.faces[i];
773 f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
774 el.node[fv[2]], el.node[fv[3]]);
775 }
776
777 f[0]->attribute = eattr0;
778 f[1]->attribute = eattr1;
779 f[2]->attribute = eattr2;
780
781 return new_id;
782}
783
784int NCMesh::NewSegment(int n0, int n1, int attr, int vattr1, int vattr2)
785{
786 // create new element, initialize nodes
787 int new_id = AddElement(Geometry::SEGMENT, attr);
788 Element &el = elements[new_id];
789 el.node[0] = n0, el.node[1] = n1;
790
791 // get (degenerate) faces and assign face attributes
792 int v0 = el.node[0], v1 = el.node[1];
793 faces.Get(v0, v0, v0, v0)->attribute = vattr1;
794 faces.Get(v1, v1, v1, v1)->attribute = vattr2;
795
796 return new_id;
797}
798
799inline bool CubeFaceLeft(int node, int* n)
800{ return node == n[0] || node == n[3] || node == n[4] || node == n[7]; }
801
802inline bool CubeFaceRight(int node, int* n)
803{ return node == n[1] || node == n[2] || node == n[5] || node == n[6]; }
804
805inline bool CubeFaceFront(int node, int* n)
806{ return node == n[0] || node == n[1] || node == n[4] || node == n[5]; }
807
808inline bool CubeFaceBack(int node, int* n)
809{ return node == n[2] || node == n[3] || node == n[6] || node == n[7]; }
810
811inline bool CubeFaceBottom(int node, int* n)
812{ return node == n[0] || node == n[1] || node == n[2] || node == n[3]; }
813
814inline bool CubeFaceTop(int node, int* n)
815{ return node == n[4] || node == n[5] || node == n[6] || node == n[7]; }
816
817inline bool PrismFaceBottom(int node, int* n)
818{ return node == n[0] || node == n[1] || node == n[2]; }
819
820inline bool PrismFaceTop(int node, int* n)
821{ return node == n[3] || node == n[4] || node == n[5]; }
822
823
824void NCMesh::ForceRefinement(int vn1, int vn2, int vn3, int vn4)
825{
826 // get the element this face belongs to
827 Face* face = faces.Find(vn1, vn2, vn3, vn4);
828 if (!face) { return; }
829
830 int elem = face->GetSingleElement();
831 Element &el = elements[elem];
832 MFEM_ASSERT(!el.ref_type, "element already refined.");
833
834 int* el_nodes = el.node;
835 if (el.Geom() == Geometry::CUBE)
836 {
837 Node* node12 = nodes.Find(vn1, vn2);
838 real_t scale = 0.5;
839
840 if (node12)
841 {
842 scale = GetScale(node12->GetScale(), vn1 > vn2);
843 }
844 else
845 {
846 Node* node34 = nodes.Find(vn3, vn4);
847 MFEM_ASSERT(node34, "Scale not set in NCMesh::ForceRefinement");
848 scale = GetScale(node34->GetScale(), vn4 > vn3);
849 }
850
851 // schedule the right split depending on face orientation
852 if ((CubeFaceLeft(vn1, el_nodes) && CubeFaceRight(vn2, el_nodes)) ||
853 (CubeFaceLeft(vn2, el_nodes) && CubeFaceRight(vn1, el_nodes)))
854 {
855 // X split
856 const bool rev = CubeFaceLeft(vn2, el_nodes) &&
857 CubeFaceRight(vn1, el_nodes);
858 ref_stack.Append(Refinement(elem, 1, GetScale(scale, rev)));
859 }
860 else if ((CubeFaceFront(vn1, el_nodes) && CubeFaceBack(vn2, el_nodes)) ||
861 (CubeFaceFront(vn2, el_nodes) && CubeFaceBack(vn1, el_nodes)))
862 {
863 // Y split
864 const bool rev = CubeFaceFront(vn2, el_nodes) &&
865 CubeFaceBack(vn1, el_nodes);
866 ref_stack.Append(Refinement(elem, 2, GetScale(scale, rev)));
867 }
868 else if ((CubeFaceBottom(vn1, el_nodes) && CubeFaceTop(vn2, el_nodes)) ||
869 (CubeFaceBottom(vn2, el_nodes) && CubeFaceTop(vn1, el_nodes)))
870 {
871 // Z split
872 const bool rev = CubeFaceBottom(vn2, el_nodes) &&
873 CubeFaceTop(vn1, el_nodes);
874 ref_stack.Append(Refinement(elem, 4, GetScale(scale, rev)));
875 }
876 else
877 {
878 MFEM_ABORT("Inconsistent element/face structure.");
879 }
880 }
881 else if (el.Geom() == Geometry::PRISM)
882 {
883 if ((PrismFaceTop(vn1, el_nodes) && PrismFaceBottom(vn4, el_nodes)) ||
884 (PrismFaceTop(vn4, el_nodes) && PrismFaceBottom(vn1, el_nodes)))
885 {
886 ref_stack.Append(Refinement(elem, 3)); // XY split
887 }
888 else if ((PrismFaceTop(vn1, el_nodes) && PrismFaceBottom(vn2, el_nodes)) ||
889 (PrismFaceTop(vn2, el_nodes) && PrismFaceBottom(vn1, el_nodes)))
890 {
891 ref_stack.Append(Refinement(elem, 4)); // Z split
892 }
893 else
894 {
895 MFEM_ABORT("Inconsistent element/face structure.");
896 }
897 }
898 else
899 {
900 MFEM_ABORT("Unsupported geometry.")
901 }
902}
903
904
905void NCMesh::FindEdgeElements(int vn1, int vn2, int vn3, int vn4,
906 Array<MeshId> &elem_edge) const
907{
908 // Assuming that f = (vn1, vn2, vn3, vn4) is a quad face and e = (vn1, vn4)
909 // is its edge, this function finds the N elements sharing e, and returns the
910 // N different MeshIds of the edge (i.e., different element-local pairs
911 // describing the edge).
912
913 int ev1 = vn1, ev2 = vn4;
914
915 // follow face refinement towards 'vn1', get an existing face
916 int split, mid[5];
917 real_t scale;
918 while ((split = QuadFaceSplitType(vn1, vn2, vn3, vn4, scale, mid)) > 0)
919 {
920 if (split == 1) // vertical
921 {
922 vn2 = mid[0]; vn3 = mid[2];
923 }
924 else // horizontal
925 {
926 vn3 = mid[1]; vn4 = mid[3];
927 }
928 }
929
930 const Face *face = faces.Find(vn1, vn2, vn3, vn4);
931 MFEM_ASSERT(face != NULL, "Face not found: "
932 << vn1 << ", " << vn2 << ", " << vn3 << ", " << vn4
933 << " (edge " << ev1 << "-" << ev2 << ").");
934
935 int elem = face->GetSingleElement();
936 int local = find_node(elements[elem], vn1);
937
938 Array<int> cousins;
939 FindVertexCousins(elem, local, cousins);
940
941 elem_edge.SetSize(0);
942 for (int i = 0; i < cousins.Size(); i++)
943 {
944 local = find_element_edge(elements[cousins[i]], ev1, ev2, false);
945 if (local > 0)
946 {
947 elem_edge.Append(MeshId(-1, cousins[i], local, Geometry::SEGMENT));
948 }
949 }
950}
951
952
953void NCMesh::CheckAnisoPrism(int vn1, int vn2, int vn3, int vn4,
954 const Refinement *refs, int nref)
955{
956 MeshId buf[4];
957 Array<MeshId> eid(buf, 4);
958 FindEdgeElements(vn1, vn2, vn3, vn4, eid);
959
960 // see if there is an element that has not been force-refined yet
961 for (int i = 0, j; i < eid.Size(); i++)
962 {
963 int elem = eid[i].element;
964 for (j = 0; j < nref; j++)
965 {
966 if (refs[j].index == elem) { break; }
967 }
968 if (j == nref) // elem not found in refs[]
969 {
970 // schedule prism refinement along Z axis
971 MFEM_ASSERT(elements[elem].Geom() == Geometry::PRISM, "");
972 ref_stack.Append(Refinement(elem, 4));
973 }
974 }
975}
976
977
978void NCMesh::CheckAnisoFace(int vn1, int vn2, int vn3, int vn4,
979 int mid12, int mid34, int level)
980{
981 // When a face is getting split anisotropically (without loss of generality
982 // we assume a "vertical" split here, see picture), it is important to make
983 // sure that the mid-face vertex node (midf) has mid34 and mid12 as parents.
984 // This is necessary for the face traversal algorithm and at places like
985 // Refine() that assume the mid-edge nodes to be accessible through the right
986 // parents. However, midf may already exist under the parents mid41 and
987 // mid23. In that case we need to "reparent" midf, i.e., reinsert it to the
988 // hash-table under the correct parents. This doesn't affect other nodes as
989 // all IDs stay the same, only the face refinement "tree" is affected.
990 //
991 // vn4 mid34 vn3
992 // *------*------*
993 // | | |
994 // | |midf |
995 // mid41 *- - - *- - - * mid23
996 // | | |
997 // | | |
998 // *------*------*
999 // vn1 mid12 vn2
1000 //
1001 // This function is recursive, because the above applies to any node along
1002 // the middle vertical edge. The function calls itself again for the bottom
1003 // and upper half of the above picture.
1004
1005 const int mid23 = FindMidEdgeNode(vn2, vn3);
1006 const int mid41 = FindMidEdgeNode(vn4, vn1);
1007
1008 if (mid23 >= 0 && mid41 >= 0)
1009 {
1010 int midf = nodes.FindId(mid23, mid41);
1011 if (midf >= 0)
1012 {
1013 Node* midfNode = nodes.Find(mid23, mid41);
1014
1015 if (midfNode)
1016 {
1017 Node* node12 = nodes.Find(vn1, vn2);
1018 if (node12)
1019 {
1020 const bool rev = (vn1 < vn2) != (mid41 < mid23);
1021 midfNode->SetScale(GetScale(node12->GetScale(), rev));
1022 }
1023 else
1024 {
1025 Node* node34 = nodes.Find(vn3, vn4);
1026 if (node34)
1027 {
1028 const bool rev = (vn4 < vn3) != (mid41 < mid23);
1029 midfNode->SetScale(GetScale(node34->GetScale(), rev));
1030 }
1031 }
1032 }
1033
1034 reparents.Append(Triple<int, int, int>(midf, mid12, mid34));
1035
1036 Node* node23 = nodes.Find(vn2, vn3);
1037 reparent_scale.Append(GetScale(node23->GetScale(), vn2 > vn3));
1038
1039 int rs = ref_stack.Size();
1040
1041 CheckAnisoFace(vn1, vn2, mid23, mid41, mid12, midf, level+1);
1042 CheckAnisoFace(mid41, mid23, vn3, vn4, midf, mid34, level+1);
1043
1044 if (HavePrisms() && nodes[midf].HasEdge())
1045 {
1046 // Check if there is a prism with edge (mid23, mid41) that we may
1047 // have missed in 'CheckAnisoFace', and force-refine it if present.
1048
1049 if (ref_stack.Size() > rs)
1050 {
1051 CheckAnisoPrism(mid23, vn3, vn4, mid41,
1052 &ref_stack[rs], ref_stack.Size() - rs);
1053 }
1054 else
1055 {
1056 CheckAnisoPrism(mid23, vn3, vn4, mid41, NULL, 0);
1057 }
1058 }
1059
1060 // perform the reparents all at once at the end
1061 if (level == 0)
1062 {
1063 for (int i = 0; i < reparents.Size(); i++)
1064 {
1065 const Triple<int, int, int> &tr = reparents[i];
1066 ReparentNode(tr.one, tr.two, tr.three, reparent_scale[i]);
1067 }
1068 reparents.DeleteAll();
1070 }
1071 return;
1072 }
1073 }
1074
1075 // Also, this is the place where forced refinements begin. In the picture
1076 // above, edges mid12-midf and midf-mid34 should actually exist in the
1077 // neighboring elements, otherwise the mesh is inconsistent and needs to be
1078 // fixed. Example: suppose an element is being refined isotropically (!)
1079 // whose neighbors across some face look like this:
1080 //
1081 // *--------*--------*
1082 // | d | e |
1083 // *--------*--------*
1084 // | c |
1085 // *--------*--------*
1086 // | | |
1087 // | a | b |
1088 // | | |
1089 // *--------*--------*
1090 //
1091 // Element 'c' needs to be refined vertically for the mesh to remain valid.
1092
1093 if (level > 0)
1094 {
1095 ForceRefinement(vn1, vn2, vn3, vn4);
1096 }
1097}
1098
1099void NCMesh::CheckIsoFace(int vn1, int vn2, int vn3, int vn4,
1100 int en1, int en2, int en3, int en4, int midf)
1101{
1102 if (!Iso)
1103 {
1104 /* If anisotropic refinements are present in the mesh, we need to check
1105 isotropically split faces as well, see second comment in CheckAnisoFace
1106 above. */
1107
1108 CheckAnisoFace(vn1, vn2, en2, en4, en1, midf);
1109 CheckAnisoFace(en4, en2, vn3, vn4, midf, en3);
1110 CheckAnisoFace(vn4, vn1, en1, en3, en4, midf);
1111 CheckAnisoFace(en3, en1, vn2, vn3, midf, en2);
1112 }
1113}
1114
1115void NCMesh::SetNodeScale(int p0, int p1, real_t scale)
1116{
1117 Node* node = nodes.Find(p0, p1);
1118 if (node) { node->SetScale(GetScale(scale, p0 > p1)); }
1119}
1120
1121void NCMesh::RefineElement(int elem, char ref_type)
1122{
1123 RefineElement(Refinement(elem, ref_type));
1124}
1125
1127{
1128 const int elem = ref.index;
1129 char ref_type = ref.GetType();
1130 const real_t scale_x = ref.s[0];
1131 const real_t scale_y = ref.s[1];
1132 const real_t scale_z = ref.s[2];
1133
1134 if (!ref_type) { return; }
1135
1136 // handle elements that may have been (force-) refined already
1137 Element &el = elements[elem];
1138 if (el.ref_type)
1139 {
1140 const char remaining = ref_type & ~el.ref_type;
1141 if (remaining)
1142 {
1143 // do the remaining splits on the children
1144 for (int i = 0; i < MaxElemChildren; i++)
1145 {
1146 if (el.child[i] >= 0)
1147 {
1148 Refinement child_ref(el.child[i], remaining);
1149 child_ref.SetScaleForType(ref.s);
1150 RefineElement(child_ref);
1151 }
1152 }
1153 }
1154 return;
1155 }
1156
1157 /*mfem::out << "Refining element " << elem << " (" << el.node[0] << ", " <<
1158 el.node[1] << ", " << el.node[2] << ", " << el.node[3] << ", " <<
1159 el.node[4] << ", " << el.node[5] << ", " << el.node[6] << ", " <<
1160 el.node[7] << "), " << "ref_type " << int(ref_type) << std::endl;*/
1161
1162 int* no = el.node;
1163 int attr = el.attribute;
1164
1165 int child[MaxElemChildren];
1166 for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
1167
1168 // get parent's face attributes
1169 int fa[MaxElemFaces];
1170 GeomInfo& gi = GI[el.Geom()];
1171 for (int i = 0; i < gi.nf; i++)
1172 {
1173 const int* fv = gi.faces[i];
1174 Face* face = faces.Find(no[fv[0]], no[fv[1]], no[fv[2]], no[fv[3]]);
1175 fa[i] = face->attribute;
1176 }
1177
1178 // create child elements
1179 if (el.Geom() == Geometry::CUBE)
1180 {
1181 // Cube vertex numbering is assumed to be as follows:
1182 //
1183 // 7 6
1184 // +-----------+ Faces: 0 bottom
1185 // /| /| 1 front
1186 // 4 / | 5 / | 2 right
1187 // +-----------+ | 3 back
1188 // | | | | 4 left
1189 // | +--------|--+ 5 top
1190 // | / 3 | / 2 Z Y
1191 // |/ |/ |/
1192 // +-----------+ *--X
1193 // 0 1
1194
1195 if (ref_type == Refinement::X) // split along X axis
1196 {
1197 const int mid01 = GetMidEdgeNode(no[0], no[1]);
1198 const int mid23 = GetMidEdgeNode(no[2], no[3]);
1199 const int mid45 = GetMidEdgeNode(no[4], no[5]);
1200 const int mid67 = GetMidEdgeNode(no[6], no[7]);
1201
1202 SetNodeScale(no[0], no[1], scale_x);
1203 SetNodeScale(no[3], no[2], scale_x);
1204 SetNodeScale(no[4], no[5], scale_x);
1205 SetNodeScale(no[7], no[6], scale_x);
1206
1207 child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
1208 no[4], mid45, mid67, no[7], attr,
1209 fa[0], fa[1], -1, fa[3], fa[4], fa[5]);
1210
1211 child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
1212 mid45, no[5], no[6], mid67, attr,
1213 fa[0], fa[1], fa[2], fa[3], -1, fa[5]);
1214
1215 CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
1216 CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
1217 CheckAnisoFace(no[4], no[5], no[6], no[7], mid45, mid67);
1218 CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
1219 }
1220 else if (ref_type == Refinement::Y) // split along Y axis
1221 {
1222 const int mid12 = GetMidEdgeNode(no[1], no[2]);
1223 const int mid30 = GetMidEdgeNode(no[3], no[0]);
1224 const int mid56 = GetMidEdgeNode(no[5], no[6]);
1225 const int mid74 = GetMidEdgeNode(no[7], no[4]);
1226
1227 SetNodeScale(no[1], no[2], scale_y);
1228 SetNodeScale(no[0], no[3], scale_y);
1229 SetNodeScale(no[5], no[6], scale_y);
1230 SetNodeScale(no[4], no[7], scale_y);
1231
1232 child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
1233 no[4], no[5], mid56, mid74, attr,
1234 fa[0], fa[1], fa[2], -1, fa[4], fa[5]);
1235
1236 child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
1237 mid74, mid56, no[6], no[7], attr,
1238 fa[0], -1, fa[2], fa[3], fa[4], fa[5]);
1239
1240 CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
1241 CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
1242 CheckAnisoFace(no[5], no[6], no[7], no[4], mid56, mid74);
1243 CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
1244 }
1245 else if (ref_type == Refinement::Z) // split along Z axis
1246 {
1247 const int mid04 = GetMidEdgeNode(no[0], no[4]);
1248 const int mid15 = GetMidEdgeNode(no[1], no[5]);
1249 const int mid26 = GetMidEdgeNode(no[2], no[6]);
1250 const int mid37 = GetMidEdgeNode(no[3], no[7]);
1251
1252 SetNodeScale(no[0], no[4], scale_z);
1253 SetNodeScale(no[1], no[5], scale_z);
1254 SetNodeScale(no[2], no[6], scale_z);
1255 SetNodeScale(no[3], no[7], scale_z);
1256
1257 child[0] = NewHexahedron(no[0], no[1], no[2], no[3],
1258 mid04, mid15, mid26, mid37, attr,
1259 fa[0], fa[1], fa[2], fa[3], fa[4], -1);
1260
1261 child[1] = NewHexahedron(mid04, mid15, mid26, mid37,
1262 no[4], no[5], no[6], no[7], attr,
1263 -1, fa[1], fa[2], fa[3], fa[4], fa[5]);
1264
1265 CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
1266 CheckAnisoFace(no[5], no[1], no[2], no[6], mid15, mid26);
1267 CheckAnisoFace(no[6], no[2], no[3], no[7], mid26, mid37);
1268 CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
1269 }
1270 else if (ref_type == Refinement::XY) // XY split
1271 {
1272 const int mid01 = GetMidEdgeNode(no[0], no[1]);
1273 const int mid12 = GetMidEdgeNode(no[1], no[2]);
1274 const int mid23 = GetMidEdgeNode(no[2], no[3]);
1275 const int mid30 = GetMidEdgeNode(no[3], no[0]);
1276
1277 const int mid45 = GetMidEdgeNode(no[4], no[5]);
1278 const int mid56 = GetMidEdgeNode(no[5], no[6]);
1279 const int mid67 = GetMidEdgeNode(no[6], no[7]);
1280 const int mid74 = GetMidEdgeNode(no[7], no[4]);
1281
1282 const int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid30);
1283 const int midf5 = GetMidFaceNode(mid45, mid56, mid67, mid74);
1284
1285 SetNodeScale(no[0], no[1], scale_x);
1286 SetNodeScale(no[3], no[2], scale_x);
1287 SetNodeScale(no[4], no[5], scale_x);
1288 SetNodeScale(no[7], no[6], scale_x);
1289
1290 SetNodeScale(no[1], no[2], scale_y);
1291 SetNodeScale(no[0], no[3], scale_y);
1292 SetNodeScale(no[5], no[6], scale_y);
1293 SetNodeScale(no[4], no[7], scale_y);
1294
1295 SetNodeScale(mid30, mid12, scale_x);
1296 SetNodeScale(mid74, mid56, scale_x);
1297
1298 SetNodeScale(mid01, mid23, scale_y);
1299 SetNodeScale(mid45, mid67, scale_y);
1300
1301 child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
1302 no[4], mid45, midf5, mid74, attr,
1303 fa[0], fa[1], -1, -1, fa[4], fa[5]);
1304
1305 child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
1306 mid45, no[5], mid56, midf5, attr,
1307 fa[0], fa[1], fa[2], -1, -1, fa[5]);
1308
1309 child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
1310 midf5, mid56, no[6], mid67, attr,
1311 fa[0], -1, fa[2], fa[3], -1, fa[5]);
1312
1313 child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
1314 mid74, midf5, mid67, no[7], attr,
1315 fa[0], -1, -1, fa[3], fa[4], fa[5]);
1316
1317 CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
1318 CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
1319 CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
1320 CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
1321
1322 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1323 CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1324 }
1325 else if (ref_type == Refinement::XZ) // XZ split
1326 {
1327 const int mid01 = GetMidEdgeNode(no[0], no[1]);
1328 const int mid23 = GetMidEdgeNode(no[2], no[3]);
1329 const int mid45 = GetMidEdgeNode(no[4], no[5]);
1330 const int mid67 = GetMidEdgeNode(no[6], no[7]);
1331
1332 const int mid04 = GetMidEdgeNode(no[0], no[4]);
1333 const int mid15 = GetMidEdgeNode(no[1], no[5]);
1334 const int mid26 = GetMidEdgeNode(no[2], no[6]);
1335 const int mid37 = GetMidEdgeNode(no[3], no[7]);
1336
1337 const int midf1 = GetMidFaceNode(mid01, mid15, mid45, mid04);
1338 const int midf3 = GetMidFaceNode(mid23, mid37, mid67, mid26);
1339
1340 SetNodeScale(no[0], no[1], scale_x);
1341 SetNodeScale(no[3], no[2], scale_x);
1342 SetNodeScale(no[4], no[5], scale_x);
1343 SetNodeScale(no[7], no[6], scale_x);
1344
1345 SetNodeScale(no[0], no[4], scale_z);
1346 SetNodeScale(no[1], no[5], scale_z);
1347 SetNodeScale(no[2], no[6], scale_z);
1348 SetNodeScale(no[3], no[7], scale_z);
1349
1350 SetNodeScale(mid04, mid15, scale_x);
1351 SetNodeScale(mid37, mid26, scale_x);
1352
1353 SetNodeScale(mid01, mid45, scale_z);
1354 SetNodeScale(mid23, mid67, scale_z);
1355
1356 child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
1357 mid04, midf1, midf3, mid37, attr,
1358 fa[0], fa[1], -1, fa[3], fa[4], -1);
1359
1360 child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
1361 midf1, mid15, mid26, midf3, attr,
1362 fa[0], fa[1], fa[2], fa[3], -1, -1);
1363
1364 child[2] = NewHexahedron(midf1, mid15, mid26, midf3,
1365 mid45, no[5], no[6], mid67, attr,
1366 -1, fa[1], fa[2], fa[3], -1, fa[5]);
1367
1368 child[3] = NewHexahedron(mid04, midf1, midf3, mid37,
1369 no[4], mid45, mid67, no[7], attr,
1370 -1, fa[1], -1, fa[3], fa[4], fa[5]);
1371
1372 CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
1373 CheckAnisoFace(no[2], no[6], no[5], no[1], mid26, mid15);
1374 CheckAnisoFace(no[6], no[7], no[4], no[5], mid67, mid45);
1375 CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
1376
1377 CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1378 CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1379 }
1380 else if (ref_type == Refinement::YZ) // YZ split
1381 {
1382 const int mid12 = GetMidEdgeNode(no[1], no[2]);
1383 const int mid30 = GetMidEdgeNode(no[3], no[0]);
1384 const int mid56 = GetMidEdgeNode(no[5], no[6]);
1385 const int mid74 = GetMidEdgeNode(no[7], no[4]);
1386
1387 const int mid04 = GetMidEdgeNode(no[0], no[4]);
1388 const int mid15 = GetMidEdgeNode(no[1], no[5]);
1389 const int mid26 = GetMidEdgeNode(no[2], no[6]);
1390 const int mid37 = GetMidEdgeNode(no[3], no[7]);
1391
1392 const int midf2 = GetMidFaceNode(mid12, mid26, mid56, mid15);
1393 const int midf4 = GetMidFaceNode(mid30, mid04, mid74, mid37);
1394
1395 SetNodeScale(no[1], no[2], scale_y);
1396 SetNodeScale(no[0], no[3], scale_y);
1397 SetNodeScale(no[5], no[6], scale_y);
1398 SetNodeScale(no[4], no[7], scale_y);
1399
1400 SetNodeScale(no[0], no[4], scale_z);
1401 SetNodeScale(no[1], no[5], scale_z);
1402 SetNodeScale(no[2], no[6], scale_z);
1403 SetNodeScale(no[3], no[7], scale_z);
1404
1405 SetNodeScale(mid15, mid26, scale_y);
1406 SetNodeScale(mid04, mid37, scale_y);
1407
1408 SetNodeScale(mid12, mid56, scale_z);
1409 SetNodeScale(mid30, mid74, scale_z);
1410
1411 child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
1412 mid04, mid15, midf2, midf4, attr,
1413 fa[0], fa[1], fa[2], -1, fa[4], -1);
1414
1415 child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
1416 midf4, midf2, mid26, mid37, attr,
1417 fa[0], -1, fa[2], fa[3], fa[4], -1);
1418
1419 child[2] = NewHexahedron(mid04, mid15, midf2, midf4,
1420 no[4], no[5], mid56, mid74, attr,
1421 -1, fa[1], fa[2], -1, fa[4], fa[5]);
1422
1423 child[3] = NewHexahedron(midf4, midf2, mid26, mid37,
1424 mid74, mid56, no[6], no[7], attr,
1425 -1, -1, fa[2], fa[3], fa[4], fa[5]);
1426
1427 CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
1428 CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
1429 CheckAnisoFace(no[3], no[7], no[6], no[2], mid37, mid26);
1430 CheckAnisoFace(no[7], no[4], no[5], no[6], mid74, mid56);
1431
1432 CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1433 CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1434 }
1435 else if (ref_type == Refinement::XYZ) // full isotropic refinement
1436 {
1437 const int mid01 = GetMidEdgeNode(no[0], no[1]);
1438 const int mid12 = GetMidEdgeNode(no[1], no[2]);
1439 const int mid23 = GetMidEdgeNode(no[2], no[3]);
1440 const int mid30 = GetMidEdgeNode(no[3], no[0]);
1441
1442 const int mid45 = GetMidEdgeNode(no[4], no[5]);
1443 const int mid56 = GetMidEdgeNode(no[5], no[6]);
1444 const int mid67 = GetMidEdgeNode(no[6], no[7]);
1445 const int mid74 = GetMidEdgeNode(no[7], no[4]);
1446
1447 const int mid04 = GetMidEdgeNode(no[0], no[4]);
1448 const int mid15 = GetMidEdgeNode(no[1], no[5]);
1449 const int mid26 = GetMidEdgeNode(no[2], no[6]);
1450 const int mid37 = GetMidEdgeNode(no[3], no[7]);
1451
1452 const int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid30);
1453 const int midf1 = GetMidFaceNode(mid01, mid15, mid45, mid04);
1454 const int midf2 = GetMidFaceNode(mid12, mid26, mid56, mid15);
1455 const int midf3 = GetMidFaceNode(mid23, mid37, mid67, mid26);
1456 const int midf4 = GetMidFaceNode(mid30, mid04, mid74, mid37);
1457 const int midf5 = GetMidFaceNode(mid45, mid56, mid67, mid74);
1458
1459 const int midel = GetMidEdgeNode(midf1, midf3);
1460
1461 SetNodeScale(midf1, midf3, scale_y);
1462
1463 SetNodeScale(no[0], no[1], scale_x);
1464 SetNodeScale(no[3], no[2], scale_x);
1465 SetNodeScale(no[1], no[2], scale_y);
1466 SetNodeScale(no[0], no[3], scale_y);
1467
1468 SetNodeScale(no[4], no[5], scale_x);
1469 SetNodeScale(no[7], no[6], scale_x);
1470 SetNodeScale(no[5], no[6], scale_y);
1471 SetNodeScale(no[4], no[7], scale_y);
1472
1473 SetNodeScale(no[0], no[4], scale_z);
1474 SetNodeScale(no[1], no[5], scale_z);
1475 SetNodeScale(no[2], no[6], scale_z);
1476 SetNodeScale(no[3], no[7], scale_z);
1477
1478 SetNodeScale(mid01, mid23, scale_y);
1479 SetNodeScale(mid30, mid12, scale_x);
1480
1481 SetNodeScale(mid04, mid15, scale_x);
1482 SetNodeScale(mid01, mid45, scale_z);
1483
1484 SetNodeScale(mid15, mid26, scale_y);
1485 SetNodeScale(mid12, mid56, scale_z);
1486
1487 SetNodeScale(mid23, mid67, scale_z);
1488 SetNodeScale(mid37, mid26, scale_x);
1489
1490 SetNodeScale(mid30, mid74, scale_z);
1491 SetNodeScale(mid04, mid37, scale_y);
1492
1493 SetNodeScale(mid74, mid56, scale_x);
1494 SetNodeScale(mid45, mid67, scale_y);
1495
1496 child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
1497 mid04, midf1, midel, midf4, attr,
1498 fa[0], fa[1], -1, -1, fa[4], -1);
1499
1500 child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
1501 midf1, mid15, midf2, midel, attr,
1502 fa[0], fa[1], fa[2], -1, -1, -1);
1503
1504 child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
1505 midel, midf2, mid26, midf3, attr,
1506 fa[0], -1, fa[2], fa[3], -1, -1);
1507
1508 child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
1509 midf4, midel, midf3, mid37, attr,
1510 fa[0], -1, -1, fa[3], fa[4], -1);
1511
1512 child[4] = NewHexahedron(mid04, midf1, midel, midf4,
1513 no[4], mid45, midf5, mid74, attr,
1514 -1, fa[1], -1, -1, fa[4], fa[5]);
1515
1516 child[5] = NewHexahedron(midf1, mid15, midf2, midel,
1517 mid45, no[5], mid56, midf5, attr,
1518 -1, fa[1], fa[2], -1, -1, fa[5]);
1519
1520 child[6] = NewHexahedron(midel, midf2, mid26, midf3,
1521 midf5, mid56, no[6], mid67, attr,
1522 -1, -1, fa[2], fa[3], -1, fa[5]);
1523
1524 child[7] = NewHexahedron(midf4, midel, midf3, mid37,
1525 mid74, midf5, mid67, no[7], attr,
1526 -1, -1, -1, fa[3], fa[4], fa[5]);
1527
1528 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1529 CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1530 CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1531 CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1532 CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1533 CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1534 }
1535 else
1536 {
1537 MFEM_ABORT("invalid refinement type.");
1538 }
1539
1540 if (ref_type != Refinement::XYZ) { Iso = false; }
1541 }
1542 else if (el.Geom() == Geometry::PRISM)
1543 {
1544 // Wedge vertex numbering:
1545 //
1546 // 5
1547 // _+_
1548 // _/ | \_ Faces: 0 bottom
1549 // 3 / | \ 4 1 top
1550 // +---------+ 2 front
1551 // | | | 3 right (1 2 5 4)
1552 // | _+_ | 4 left (2 0 3 5)
1553 // | _/ 2 \_ | Z Y
1554 // |/ \| | /
1555 // +---------+ *--X
1556 // 0 1
1557
1558 if (ref_type < 4) // XY refinement (split in 4 wedges)
1559 {
1560 ref_type = Refinement::XY; // for consistency
1561
1562 int mid01 = GetMidEdgeNode(no[0], no[1]);
1563 int mid12 = GetMidEdgeNode(no[1], no[2]);
1564 int mid20 = GetMidEdgeNode(no[2], no[0]);
1565
1566 int mid34 = GetMidEdgeNode(no[3], no[4]);
1567 int mid45 = GetMidEdgeNode(no[4], no[5]);
1568 int mid53 = GetMidEdgeNode(no[5], no[3]);
1569
1570 child[0] = NewWedge(no[0], mid01, mid20,
1571 no[3], mid34, mid53, attr,
1572 fa[0], fa[1], fa[2], -1, fa[4]);
1573
1574 child[1] = NewWedge(mid01, no[1], mid12,
1575 mid34, no[4], mid45, attr,
1576 fa[0], fa[1], fa[2], fa[3], -1);
1577
1578 child[2] = NewWedge(mid20, mid12, no[2],
1579 mid53, mid45, no[5], attr,
1580 fa[0], fa[1], -1, fa[3], fa[4]);
1581
1582 child[3] = NewWedge(mid12, mid20, mid01,
1583 mid45, mid53, mid34, attr,
1584 fa[0], fa[1], -1, -1, -1);
1585
1586 CheckAnisoFace(no[0], no[1], no[4], no[3], mid01, mid34);
1587 CheckAnisoFace(no[1], no[2], no[5], no[4], mid12, mid45);
1588 CheckAnisoFace(no[2], no[0], no[3], no[5], mid20, mid53);
1589 }
1590 else if (ref_type == Refinement::Z) // Z refinement only (split in 2 wedges)
1591 {
1592 int mid03 = GetMidEdgeNode(no[0], no[3]);
1593 int mid14 = GetMidEdgeNode(no[1], no[4]);
1594 int mid25 = GetMidEdgeNode(no[2], no[5]);
1595
1596 child[0] = NewWedge(no[0], no[1], no[2],
1597 mid03, mid14, mid25, attr,
1598 fa[0], -1, fa[2], fa[3], fa[4]);
1599
1600 child[1] = NewWedge(mid03, mid14, mid25,
1601 no[3], no[4], no[5], attr,
1602 -1, fa[1], fa[2], fa[3], fa[4]);
1603
1604 CheckAnisoFace(no[3], no[0], no[1], no[4], mid03, mid14);
1605 CheckAnisoFace(no[4], no[1], no[2], no[5], mid14, mid25);
1606 CheckAnisoFace(no[5], no[2], no[0], no[3], mid25, mid03);
1607 }
1608 else // ref_type > 4, full isotropic refinement (split in 8 wedges)
1609 {
1610 ref_type = Refinement::XYZ; // for consistency
1611
1612 int mid01 = GetMidEdgeNode(no[0], no[1]);
1613 int mid12 = GetMidEdgeNode(no[1], no[2]);
1614 int mid20 = GetMidEdgeNode(no[2], no[0]);
1615
1616 int mid34 = GetMidEdgeNode(no[3], no[4]);
1617 int mid45 = GetMidEdgeNode(no[4], no[5]);
1618 int mid53 = GetMidEdgeNode(no[5], no[3]);
1619
1620 int mid03 = GetMidEdgeNode(no[0], no[3]);
1621 int mid14 = GetMidEdgeNode(no[1], no[4]);
1622 int mid25 = GetMidEdgeNode(no[2], no[5]);
1623
1624 int midf2 = GetMidFaceNode(mid01, mid14, mid34, mid03);
1625 int midf3 = GetMidFaceNode(mid12, mid25, mid45, mid14);
1626 int midf4 = GetMidFaceNode(mid20, mid03, mid53, mid25);
1627
1628 child[0] = NewWedge(no[0], mid01, mid20,
1629 mid03, midf2, midf4, attr,
1630 fa[0], -1, fa[2], -1, fa[4]);
1631
1632 child[1] = NewWedge(mid01, no[1], mid12,
1633 midf2, mid14, midf3, attr,
1634 fa[0], -1, fa[2], fa[3], -1);
1635
1636 child[2] = NewWedge(mid20, mid12, no[2],
1637 midf4, midf3, mid25, attr,
1638 fa[0], -1, -1, fa[3], fa[4]);
1639
1640 child[3] = NewWedge(mid12, mid20, mid01,
1641 midf3, midf4, midf2, attr,
1642 fa[0], -1, -1, -1, -1);
1643
1644 child[4] = NewWedge(mid03, midf2, midf4,
1645 no[3], mid34, mid53, attr,
1646 -1, fa[1], fa[2], -1, fa[4]);
1647
1648 child[5] = NewWedge(midf2, mid14, midf3,
1649 mid34, no[4], mid45, attr,
1650 -1, fa[1], fa[2], fa[3], -1);
1651
1652 child[6] = NewWedge(midf4, midf3, mid25,
1653 mid53, mid45, no[5], attr,
1654 -1, fa[1], -1, fa[3], fa[4]);
1655
1656 child[7] = NewWedge(midf3, midf4, midf2,
1657 mid45, mid53, mid34, attr,
1658 -1, fa[1], -1, -1, -1);
1659
1660 CheckIsoFace(no[0], no[1], no[4], no[3], mid01, mid14, mid34, mid03, midf2);
1661 CheckIsoFace(no[1], no[2], no[5], no[4], mid12, mid25, mid45, mid14, midf3);
1662 CheckIsoFace(no[2], no[0], no[3], no[5], mid20, mid03, mid53, mid25, midf4);
1663 }
1664
1665 if (ref_type != Refinement::XYZ) { Iso = false; }
1666 }
1667 else if (el.Geom() == Geometry::TETRAHEDRON)
1668 {
1669 // Tetrahedron vertex numbering:
1670 //
1671 // 3
1672 // + Faces: 0 back (1, 2, 3)
1673 // |\\_ 1 left (0, 3, 2)
1674 // || \_ 2 front (0, 1, 3)
1675 // | \ \_ 3 bottom (0, 1, 2)
1676 // | +__ \_
1677 // | /2 \__ \_ Z Y
1678 // |/ \__\ | /
1679 // +------------+ *--X
1680 // 0 1
1681
1682 ref_type = Refinement::XYZ; // for consistency
1683
1684 int mid01 = GetMidEdgeNode(no[0], no[1]);
1685 int mid12 = GetMidEdgeNode(no[1], no[2]);
1686 int mid02 = GetMidEdgeNode(no[2], no[0]);
1687
1688 int mid03 = GetMidEdgeNode(no[0], no[3]);
1689 int mid13 = GetMidEdgeNode(no[1], no[3]);
1690 int mid23 = GetMidEdgeNode(no[2], no[3]);
1691
1692 child[0] = NewTetrahedron(no[0], mid01, mid02, mid03, attr,
1693 -1, fa[1], fa[2], fa[3]);
1694
1695 child[1] = NewTetrahedron(mid01, no[1], mid12, mid13, attr,
1696 fa[0], -1, fa[2], fa[3]);
1697
1698 child[2] = NewTetrahedron(mid02, mid12, no[2], mid23, attr,
1699 fa[0], fa[1], -1, fa[3]);
1700
1701 child[3] = NewTetrahedron(mid03, mid13, mid23, no[3], attr,
1702 fa[0], fa[1], fa[2], -1);
1703
1704 // There are three ways to split the inner octahedron. A good strategy is
1705 // to use the shortest diagonal. At the moment we don't have the geometric
1706 // information in this class to determine which diagonal is the shortest,
1707 // but it seems that with reasonable shapes of the coarse tets and MFEM's
1708 // default tet orientation, always using tet_type == 0 produces stable
1709 // refinements. Types 1 and 2 are unused for now.
1710 el.tet_type = 0;
1711
1712 if (el.tet_type == 0) // shortest diagonal mid01--mid23
1713 {
1714 child[4] = NewTetrahedron(mid01, mid23, mid02, mid03, attr,
1715 fa[1], -1, -1, -1);
1716
1717 child[5] = NewTetrahedron(mid01, mid23, mid03, mid13, attr,
1718 -1, fa[2], -1, -1);
1719
1720 child[6] = NewTetrahedron(mid01, mid23, mid13, mid12, attr,
1721 fa[0], -1, -1, -1);
1722
1723 child[7] = NewTetrahedron(mid01, mid23, mid12, mid02, attr,
1724 -1, fa[3], -1, -1);
1725 }
1726 else if (el.tet_type == 1) // shortest diagonal mid12--mid03
1727 {
1728 child[4] = NewTetrahedron(mid03, mid01, mid02, mid12, attr,
1729 fa[3], -1, -1, -1);
1730
1731 child[5] = NewTetrahedron(mid03, mid02, mid23, mid12, attr,
1732 -1, -1, -1, fa[1]);
1733
1734 child[6] = NewTetrahedron(mid03, mid23, mid13, mid12, attr,
1735 fa[0], -1, -1, -1);
1736
1737 child[7] = NewTetrahedron(mid03, mid13, mid01, mid12, attr,
1738 -1, -1, -1, fa[2]);
1739 }
1740 else // el.tet_type == 2, shortest diagonal mid02--mid13
1741 {
1742 child[4] = NewTetrahedron(mid02, mid01, mid13, mid03, attr,
1743 fa[2], -1, -1, -1);
1744
1745 child[5] = NewTetrahedron(mid02, mid03, mid13, mid23, attr,
1746 -1, -1, fa[1], -1);
1747
1748 child[6] = NewTetrahedron(mid02, mid23, mid13, mid12, attr,
1749 fa[0], -1, -1, -1);
1750
1751 child[7] = NewTetrahedron(mid02, mid12, mid13, mid01, attr,
1752 -1, -1, fa[3], -1);
1753 }
1754 }
1755 else if (el.Geom() == Geometry::PYRAMID)
1756 {
1757 // Pyramid vertex numbering:
1758 //
1759 // 4
1760 // + \_ Faces: 0 bottom (3,2,1,0)
1761 // |\\_ \_ 1 front (0, 1, 4)
1762 // || \_ \__ 2 right (1, 2, 4)
1763 // | \ \_ \__ 3 back (2, 3, 4)
1764 // | +____\_ ____\ 4 left (3, 0, 4)
1765 // | /3 \_ 2 Z Y
1766 // |/ \ / | /
1767 // +------------+ *--X
1768 // 0 1
1769
1770 ref_type = Refinement::XYZ; // for consistency
1771
1772 int mid01 = GetMidEdgeNode(no[0], no[1]);
1773 int mid12 = GetMidEdgeNode(no[1], no[2]);
1774 int mid23 = GetMidEdgeNode(no[2], no[3]);
1775 int mid03 = GetMidEdgeNode(no[0], no[3]);
1776 int mid04 = GetMidEdgeNode(no[0], no[4]);
1777 int mid14 = GetMidEdgeNode(no[1], no[4]);
1778 int mid24 = GetMidEdgeNode(no[2], no[4]);
1779 int mid34 = GetMidEdgeNode(no[3], no[4]);
1780 int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid03);
1781
1782 child[0] = NewPyramid(no[0], mid01, midf0, mid03, mid04,
1783 attr, fa[0], fa[1], -1, -1, fa[4]);
1784
1785 child[1] = NewPyramid(mid01, no[1], mid12, midf0, mid14,
1786 attr, fa[0], fa[1], fa[2], -1, -1);
1787
1788 child[2] = NewPyramid(midf0, mid12, no[2], mid23, mid24,
1789 attr, fa[0], -1, fa[2], fa[3], -1);
1790
1791 child[3] = NewPyramid(mid03, midf0, mid23, no[3], mid34,
1792 attr, fa[0], -1, -1, fa[3], fa[4]);
1793
1794 child[4] = NewPyramid(mid24, mid14, mid04, mid34, midf0,
1795 attr, -1, -1, -1, -1, -1);
1796
1797 child[5] = NewPyramid(mid04, mid14, mid24, mid34, no[4],
1798 attr, -1, fa[1], fa[2], fa[3], fa[4]);
1799
1800 child[6] = NewTetrahedron(mid01, midf0, mid04, mid14,
1801 attr, -1, -1, -1, fa[1]);
1802
1803 child[7] = NewTetrahedron(midf0, mid14, mid12, mid24,
1804 attr, -1, -1, fa[2], -1);
1805
1806 child[8] = NewTetrahedron(midf0, mid23, mid34, mid24,
1807 attr, -1, -1, fa[3], -1);
1808
1809 child[9] = NewTetrahedron(mid03, mid04, midf0, mid34,
1810 attr, -1, fa[4], -1, -1);
1811
1812 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid03, midf0);
1813 }
1814 else if (el.Geom() == Geometry::SQUARE)
1815 {
1816 ref_type &= 0x3; // ignore Z bit
1817
1818 if (ref_type == Refinement::X) // X split
1819 {
1820 const int mid01 = nodes.GetId(no[0], no[1]);
1821 const int mid23 = nodes.GetId(no[2], no[3]);
1822
1823 SetNodeScale(no[0], no[1], scale_x);
1824 SetNodeScale(no[3], no[2], scale_x);
1825
1826 child[0] = NewQuadrilateral(no[0], mid01, mid23, no[3],
1827 attr, fa[0], -1, fa[2], fa[3]);
1828
1829 child[1] = NewQuadrilateral(mid01, no[1], no[2], mid23,
1830 attr, fa[0], fa[1], fa[2], -1);
1831 }
1832 else if (ref_type == Refinement::Y) // Y split
1833 {
1834 const int mid12 = nodes.GetId(no[1], no[2]);
1835 const int mid30 = nodes.GetId(no[3], no[0]);
1836
1837 SetNodeScale(no[1], no[2], scale_y);
1838 SetNodeScale(no[0], no[3], scale_y);
1839
1840 child[0] = NewQuadrilateral(no[0], no[1], mid12, mid30,
1841 attr, fa[0], fa[1], -1, fa[3]);
1842
1843 child[1] = NewQuadrilateral(mid30, mid12, no[2], no[3],
1844 attr, -1, fa[1], fa[2], fa[3]);
1845 }
1846 else if (ref_type == Refinement::XY) // iso split
1847 {
1848 const int mid01 = nodes.GetId(no[0], no[1]);
1849 const int mid12 = nodes.GetId(no[1], no[2]);
1850 const int mid23 = nodes.GetId(no[2], no[3]);
1851 const int mid30 = nodes.GetId(no[3], no[0]);
1852
1853 const int midel = nodes.GetId(mid01, mid23);
1854
1855 SetNodeScale(no[0], no[1], scale_x);
1856 SetNodeScale(no[3], no[2], scale_x);
1857
1858 SetNodeScale(no[1], no[2], scale_y);
1859 SetNodeScale(no[0], no[3], scale_y);
1860
1861 SetNodeScale(mid01, mid23, scale_y);
1862
1863 child[0] = NewQuadrilateral(no[0], mid01, midel, mid30,
1864 attr, fa[0], -1, -1, fa[3]);
1865
1866 child[1] = NewQuadrilateral(mid01, no[1], mid12, midel,
1867 attr, fa[0], fa[1], -1, -1);
1868
1869 child[2] = NewQuadrilateral(midel, mid12, no[2], mid23,
1870 attr, -1, fa[1], fa[2], -1);
1871
1872 child[3] = NewQuadrilateral(mid30, midel, mid23, no[3],
1873 attr, -1, -1, fa[2], fa[3]);
1874 }
1875 else
1876 {
1877 MFEM_ABORT("Invalid refinement type.");
1878 }
1879
1880 if (ref_type != Refinement::XY) { Iso = false; }
1881 }
1882 else if (el.Geom() == Geometry::TRIANGLE)
1883 {
1884 ref_type = Refinement::XY; // for consistency
1885
1886 // isotropic split - the only ref_type available for triangles
1887 int mid01 = nodes.GetId(no[0], no[1]);
1888 int mid12 = nodes.GetId(no[1], no[2]);
1889 int mid20 = nodes.GetId(no[2], no[0]);
1890
1891 child[0] = NewTriangle(no[0], mid01, mid20, attr, fa[0], -1, fa[2]);
1892 child[1] = NewTriangle(mid01, no[1], mid12, attr, fa[0], fa[1], -1);
1893 child[2] = NewTriangle(mid20, mid12, no[2], attr, -1, fa[1], fa[2]);
1894 child[3] = NewTriangle(mid12, mid20, mid01, attr, -1, -1, -1);
1895 }
1896 else if (el.Geom() == Geometry::SEGMENT)
1897 {
1898 ref_type = Refinement::X; // for consistency
1899
1900 int mid = nodes.GetId(no[0], no[1]);
1901 child[0] = NewSegment(no[0], mid, attr, fa[0], -1);
1902 child[1] = NewSegment(mid, no[1], attr, -1, fa[1]);
1903 }
1904 else
1905 {
1906 MFEM_ABORT("Unsupported element geometry.");
1907 }
1908
1909 // start using the nodes of the children, create edges & faces
1910 for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1911 {
1912 ReferenceElement(child[i]);
1913 }
1914
1915 int buf[MaxElemFaces];
1916 Array<int> parentFaces(buf, MaxElemFaces);
1917 parentFaces.SetSize(0);
1918
1919 // sign off of all nodes of the parent, clean up unused nodes, but keep faces
1920 UnreferenceElement(elem, parentFaces);
1921
1922 // register the children in their faces
1923 for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1924 {
1925 RegisterFaces(child[i]);
1926 }
1927
1928 // clean up parent faces, if unused
1929 DeleteUnusedFaces(parentFaces);
1930
1931 // make the children inherit our rank; set the parent element
1932 for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1933 {
1934 Element &ch = elements[child[i]];
1935 ch.rank = el.rank;
1936 ch.parent = elem;
1937 }
1938
1939 // finish the refinement
1940 el.ref_type = ref_type;
1941 std::memcpy(el.child, child, sizeof(el.child));
1942}
1943
1944
1945void NCMesh::Refine(const Array<Refinement>& refinements)
1946{
1947 // push all refinements on the stack in reverse order
1948 ref_stack.Reserve(refinements.Size());
1949 for (int i = refinements.Size()-1; i >= 0; i--)
1950 {
1951 Refinement ref = refinements[i]; // Copy
1952 ref.index = leaf_elements[ref.index];
1953 ref_stack.Append(ref);
1954 }
1955
1956 // keep refining as long as the stack contains something
1957 int nforced = 0;
1958 while (ref_stack.Size())
1959 {
1960 Refinement ref = ref_stack.Last();
1961 ref_stack.DeleteLast();
1962
1963 int size = ref_stack.Size();
1964 RefineElement(ref);
1965 nforced += ref_stack.Size() - size;
1966 }
1967
1968 /* TODO: the current algorithm of forced refinements is not optimal. As
1969 forced refinements spread through the mesh, some may not be necessary in
1970 the end, since the affected elements may still be scheduled for refinement
1971 that could stop the propagation. We should introduce the member
1972 Element::ref_pending that would show the intended refinement in the batch.
1973 A forced refinement would be combined with ref_pending to (possibly) stop
1974 the propagation earlier.
1975
1976 Update: what about a FIFO instead of ref_stack? */
1977
1978#if defined(MFEM_DEBUG) && !defined(MFEM_USE_MPI)
1979 mfem::out << "Refined " << refinements.Size() << " + " << nforced
1980 << " elements" << std::endl;
1981#else
1982 MFEM_CONTRACT_VAR(nforced);
1983#endif
1984
1985 ref_stack.DeleteAll();
1986 shadow.DeleteAll();
1987
1988 Update();
1989}
1990
1991
1992//// Derefinement //////////////////////////////////////////////////////////////
1993
1995{
1996 if (!el.ref_type) { return el.node[index]; }
1997
1998 // need to retrieve node from a child element (there is always a child that
1999 // inherited the parent's corner under the same index)
2000 int ch;
2001 switch (el.Geom())
2002 {
2003 case Geometry::CUBE:
2004 ch = el.child[hex_deref_table[el.ref_type - 1][index]];
2005 break;
2006
2007 case Geometry::PRISM:
2008 ch = prism_deref_table[el.ref_type - 1][index];
2009 MFEM_ASSERT(ch != -1, "");
2010 ch = el.child[ch];
2011 break;
2012
2013 case Geometry::PYRAMID:
2014 ch = pyramid_deref_table[el.ref_type - 1][index];
2015 MFEM_ASSERT(ch != -1, "");
2016 ch = el.child[ch];
2017 break;
2018
2019 case Geometry::SQUARE:
2020 ch = el.child[quad_deref_table[el.ref_type - 1][index]];
2021 break;
2022
2024 case Geometry::TRIANGLE:
2025 ch = el.child[index];
2026 break;
2027
2028 default:
2029 ch = 0; // suppress compiler warning
2030 MFEM_ABORT("Unsupported element geometry.");
2031 }
2032 return RetrieveNode(elements[ch], index);
2033}
2034
2035
2037{
2038 Element &el = elements[elem];
2039 if (!el.ref_type) { return; }
2040
2041 int child[MaxElemChildren];
2042 std::memcpy(child, el.child, sizeof(child));
2043
2044 // first make sure that all children are leaves, derefine them if not
2045 for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
2046 {
2047 if (elements[child[i]].ref_type)
2048 {
2049 DerefineElement(child[i]);
2050 }
2051 }
2052
2053 int faces_attribute[MaxElemFaces];
2054 int ref_type_key = el.ref_type - 1;
2055
2056 for (int i = 0; i < MaxElemNodes; i++) { el.node[i] = -1; }
2057
2058 // retrieve original corner nodes and face attributes from the children
2059 if (el.Geom() == Geometry::CUBE)
2060 {
2061 // Sets corner nodes from childs
2062 constexpr int nb_cube_childs = 8;
2063 for (int i = 0; i < nb_cube_childs; i++)
2064 {
2065 const int child_local_index = hex_deref_table[ref_type_key][i];
2066 const int child_global_index = child[child_local_index];
2067 Element &ch = elements[child_global_index];
2068 el.node[i] = ch.node[i];
2069 }
2070 // Sets faces attributes from childs' faces
2071 constexpr int nb_cube_faces = 6;
2072 for (int i = 0; i < nb_cube_faces; i++)
2073 {
2074 const int child_local_index = hex_deref_table[ref_type_key]
2075 [i + nb_cube_childs];
2076 const int child_global_index = child[child_local_index];
2077 Element &ch = elements[child_global_index];
2078 const int* fv = GI[el.Geom()].faces[i];
2079 faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
2080 ch.node[fv[2]], ch.node[fv[3]])
2081 ->attribute;
2082 }
2083 }
2084 else if (el.Geom() == Geometry::PRISM)
2085 {
2086 MFEM_ASSERT(prism_deref_table[ref_type_key][0] != -1,
2087 "invalid prism refinement");
2088 constexpr int nb_prism_childs = 6;
2089 for (int i = 0; i < nb_prism_childs; i++)
2090 {
2091 const int child_local_index = prism_deref_table[ref_type_key][i];
2092 const int child_global_index = child[child_local_index];
2093 Element &ch = elements[child_global_index];
2094 el.node[i] = ch.node[i];
2095 }
2096 el.node[6] = el.node[7] = -1;
2097
2098 constexpr int nb_prism_faces = 5;
2099 for (int i = 0; i < nb_prism_faces; i++)
2100 {
2101 const int child_local_index = prism_deref_table[ref_type_key]
2102 [i + nb_prism_childs];
2103 const int child_global_index = child[child_local_index];
2104 Element &ch = elements[child_global_index];
2105 const int* fv = GI[el.Geom()].faces[i];
2106 faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
2107 ch.node[fv[2]], ch.node[fv[3]])
2108 ->attribute;
2109 }
2110 }
2111 else if (el.Geom() == Geometry::PYRAMID)
2112 {
2113 MFEM_ASSERT(pyramid_deref_table[ref_type_key][0] != -1,
2114 "invalid pyramid refinement");
2115 constexpr int nb_pyramid_childs = 5;
2116 for (int i = 0; i < nb_pyramid_childs; i++)
2117 {
2118 const int child_local_index = pyramid_deref_table[ref_type_key][i];
2119 const int child_global_index = child[child_local_index];
2120 Element &ch = elements[child_global_index];
2121 el.node[i] = ch.node[i];
2122 }
2123 el.node[5] = el.node[6] = el.node[7] = -1;
2124
2125
2126 constexpr int nb_pyramid_faces = 5;
2127 for (int i = 0; i < nb_pyramid_faces; i++)
2128 {
2129 const int child_local_index = pyramid_deref_table[ref_type_key]
2130 [i + nb_pyramid_childs];
2131 const int child_global_index = child[child_local_index];
2132 Element &ch = elements[child_global_index];
2133 const int* fv = GI[el.Geom()].faces[i];
2134 faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
2135 ch.node[fv[2]], ch.node[fv[3]])
2136 ->attribute;
2137 }
2138 }
2139 else if (el.Geom() == Geometry::TETRAHEDRON)
2140 {
2141 for (int i = 0; i < 4; i++)
2142 {
2143 Element& ch1 = elements[child[i]];
2144 Element& ch2 = elements[child[(i+1) & 0x3]];
2145 el.node[i] = ch1.node[i];
2146 const int* fv = GI[el.Geom()].faces[i];
2147 faces_attribute[i] = faces.Find(ch2.node[fv[0]], ch2.node[fv[1]],
2148 ch2.node[fv[2]], ch2.node[fv[3]])
2149 ->attribute;
2150 }
2151 }
2152 else if (el.Geom() == Geometry::SQUARE)
2153 {
2154 constexpr int nb_square_childs = 4;
2155 for (int i = 0; i < nb_square_childs; i++)
2156 {
2157 const int child_local_index = quad_deref_table[ref_type_key][i];
2158 const int child_global_index = child[child_local_index];
2159 Element &ch = elements[child_global_index];
2160 el.node[i] = ch.node[i];
2161 }
2162 constexpr int nb_square_faces = 4;
2163 for (int i = 0; i < nb_square_faces; i++)
2164 {
2165 const int child_local_index = quad_deref_table[ref_type_key]
2166 [i + nb_square_childs];
2167 const int child_global_index = child[child_local_index];
2168 Element &ch = elements[child_global_index];
2169 const int* fv = GI[el.Geom()].faces[i];
2170 faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
2171 ch.node[fv[2]], ch.node[fv[3]])
2172 ->attribute;
2173 }
2174 }
2175 else if (el.Geom() == Geometry::TRIANGLE)
2176 {
2177 constexpr int nb_triangle_childs = 3;
2178 for (int i = 0; i < nb_triangle_childs; i++)
2179 {
2180 Element& ch = elements[child[i]];
2181 el.node[i] = ch.node[i];
2182 const int* fv = GI[el.Geom()].faces[i];
2183 faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
2184 ch.node[fv[2]], ch.node[fv[3]])
2185 ->attribute;
2186 }
2187 }
2188 else if (el.Geom() == Geometry::SEGMENT)
2189 {
2190 constexpr int nb_segment_childs = 2;
2191 for (int i = 0; i < nb_segment_childs; i++)
2192 {
2193 int ni = elements[child[i]].node[i];
2194 el.node[i] = ni;
2195 faces_attribute[i] = faces.Find(ni, ni, ni, ni)->attribute;
2196 }
2197 }
2198 else
2199 {
2200 MFEM_ABORT("Unsupported element geometry.");
2201 }
2202
2203 // sign in to all nodes
2204 ReferenceElement(elem);
2205
2207 Array<int> childFaces(buf, MaxElemChildren*MaxElemFaces);
2208 childFaces.SetSize(0);
2209
2210 // delete children, determine rank
2211 el.rank = std::numeric_limits<int>::max();
2212 for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
2213 {
2214 el.rank = std::min(el.rank, elements[child[i]].rank);
2215 UnreferenceElement(child[i], childFaces);
2216 FreeElement(child[i]);
2217 }
2218
2219 RegisterFaces(elem, faces_attribute);
2220
2221 // delete unused faces
2222 childFaces.Sort();
2223 childFaces.Unique();
2224 DeleteUnusedFaces(childFaces);
2225
2226 el.ref_type = 0;
2227}
2228
2229
2231{
2232 Element &el = elements[elem];
2233 if (!el.ref_type) { return; }
2234
2235 int total = 0, ref = 0, ghost = 0;
2236 for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
2237 {
2238 total++;
2239 Element &ch = elements[el.child[i]];
2240 if (ch.ref_type) { ref++; break; }
2241 if (IsGhost(ch)) { ghost++; }
2242 }
2243
2244 if (!ref && ghost < total)
2245 {
2246 // can be derefined, add to list
2247 int next_row = list.Size() ? (list.Last().from + 1) : 0;
2248 for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
2249 {
2250 Element &ch = elements[el.child[i]];
2251 list.Append(Connection(next_row, ch.index));
2252 }
2253 }
2254 else
2255 {
2256 for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
2257 {
2258 CollectDerefinements(el.child[i], list);
2259 }
2260 }
2261}
2262
2264{
2265 Array<Connection> list;
2266 list.Reserve(leaf_elements.Size());
2267
2268 for (int i = 0; i < root_state.Size(); i++)
2269 {
2270 CollectDerefinements(i, list);
2271 }
2272
2273 int size = list.Size() ? (list.Last().from + 1) : 0;
2274 derefinements.MakeFromList(size, list);
2275 return derefinements;
2276}
2277
2279 Array<int> &level_ok, int max_nc_level)
2280{
2281 level_ok.SetSize(deref_table.Size());
2282 for (int i = 0; i < deref_table.Size(); i++)
2283 {
2284 const int* fine = deref_table.GetRow(i), size = deref_table.RowSize(i);
2285 Element &parent = elements[elements[leaf_elements[fine[0]]].parent];
2286
2287 int ok = 1;
2288 for (int j = 0; j < size; j++)
2289 {
2290 int splits[3];
2291 CountSplits(leaf_elements[fine[j]], splits);
2292
2293 for (int k = 0; k < Dim; k++)
2294 {
2295 if ((parent.ref_type & (1 << k)) &&
2296 splits[k] >= max_nc_level)
2297 {
2298 ok = 0; break;
2299 }
2300 }
2301 if (!ok) { break; }
2302 }
2303 level_ok[i] = ok;
2304 }
2305}
2306
2307void NCMesh::Derefine(const Array<int> &derefs)
2308{
2309 MFEM_VERIFY(Dim < 3 || Iso,
2310 "derefinement of 3D anisotropic meshes not implemented yet.");
2311
2313
2314 Array<int> fine_coarse;
2315 leaf_elements.Copy(fine_coarse);
2316
2317 // perform the derefinements
2318 for (int i = 0; i < derefs.Size(); i++)
2319 {
2320 int row = derefs[i];
2321 MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
2322 "invalid derefinement number.");
2323
2324 const int* fine = derefinements.GetRow(row);
2325 int parent = elements[leaf_elements[fine[0]]].parent;
2326
2327 // record the relation of the fine elements to their parent
2328 SetDerefMatrixCodes(parent, fine_coarse);
2329
2330 DerefineElement(parent);
2331 }
2332
2333 // update leaf_elements, Element::index etc.
2334 Update();
2335
2336 // link old fine elements to the new coarse elements
2337 for (int i = 0; i < fine_coarse.Size(); i++)
2338 {
2339 transforms.embeddings[i].parent = elements[fine_coarse[i]].index;
2340 }
2341}
2342
2344{
2345 int nfine = leaf_elements.Size();
2346
2347 // this will tell GetDerefinementTransforms that transforms are not finished
2348 transforms.Clear();
2349
2350 transforms.embeddings.SetSize(nfine);
2351 for (int i = 0; i < nfine; i++)
2352 {
2354 emb.parent = -1;
2355 emb.matrix = 0;
2356 Element &el = elements[leaf_elements[i]];
2357 emb.geom = el.Geom();
2358 emb.ghost = IsGhost(el);
2359 }
2360}
2361
2362void NCMesh::SetDerefMatrixCodes(int parent, Array<int> &fine_coarse)
2363{
2364 // encode the ref_type and child number for GetDerefinementTransforms()
2365 Element &prn = elements[parent];
2366 for (int i = 0; i < MaxElemChildren && prn.child[i] >= 0; i++)
2367 {
2368 Element &ch = elements[prn.child[i]];
2369 if (ch.index >= 0)
2370 {
2371 int code = (prn.ref_type << 4) | i;
2372 transforms.embeddings[ch.index].matrix = code;
2373 fine_coarse[ch.index] = parent;
2374 }
2375 }
2376}
2377
2378
2379//// Mesh Interface ////////////////////////////////////////////////////////////
2380
2381void NCMesh::CollectLeafElements(int elem, int state, Array<int> &ghosts,
2382 int &counter)
2383{
2384 Element &el = elements[elem];
2385 if (!el.ref_type)
2386 {
2387 if (el.rank >= 0) // skip elements beyond the ghost layer in parallel
2388 {
2389 if (!IsGhost(el))
2390 {
2391 leaf_elements.Append(elem);
2392 }
2393 else
2394 {
2395 // in parallel (or in serial loading a parallel file), collect
2396 // elements of neighboring ranks in a separate array
2397 ghosts.Append(elem);
2398 }
2399
2400 // assign the SFC index (temporarily, will be replaced by Mesh index)
2401 el.index = counter++;
2402 }
2403 else
2404 {
2405 // elements beyond the ghost layer are invalid and don't appear in
2406 // 'leaf_elements' (also for performance reasons)
2407 el.index = -1;
2408 }
2409 }
2410 else // Refined element
2411 {
2412 // in non-leaf elements, the 'rank' and 'index' members have no meaning
2413 el.rank = -1;
2414 el.index = -1;
2415
2416 // recurse to subtrees; try to order leaf elements along a space-filling
2417 // curve by changing the order the children are visited at each level
2418 if (el.Geom() == Geometry::SQUARE && el.ref_type == Refinement::XY)
2419 {
2420 for (int i = 0; i < 4; i++)
2421 {
2422 int ch = quad_hilbert_child_order[state][i];
2423 int st = quad_hilbert_child_state[state][i];
2424 CollectLeafElements(el.child[ch], st, ghosts, counter);
2425 }
2426 }
2427 else if (el.Geom() == Geometry::CUBE && el.ref_type == Refinement::XYZ)
2428 {
2429 for (int i = 0; i < 8; i++)
2430 {
2431 int ch = hex_hilbert_child_order[state][i];
2432 int st = hex_hilbert_child_state[state][i];
2433 CollectLeafElements(el.child[ch], st, ghosts, counter);
2434 }
2435 }
2436 else // no space filling curve tables yet for remaining cases
2437 {
2438 for (int i = 0; i < MaxElemChildren; i++)
2439 {
2440 if (el.child[i] >= 0)
2441 {
2442 CollectLeafElements(el.child[i], state, ghosts, counter);
2443 }
2444 }
2445 }
2446 }
2447}
2448
2450{
2451 Array<int> ghosts;
2452
2453 // collect leaf elements in leaf_elements and ghosts elements in ghosts from
2454 // all roots
2456 for (int i = 0, counter = 0; i < root_state.Size(); i++)
2457 {
2458 CollectLeafElements(i, root_state[i], ghosts, counter);
2459 }
2460
2462 NGhostElements = ghosts.Size();
2463
2464 // append ghost elements at the end of 'leaf_element' (if any) and assign the
2465 // final (Mesh) indices of leaves
2466 leaf_elements.Append(ghosts);
2468 for (int i = 0; i < leaf_elements.Size(); i++)
2469 {
2470 Element &el = elements[leaf_elements[i]];
2471 leaf_sfc_index[i] = el.index;
2472 el.index = i;
2473 }
2474}
2475
2477{
2478#ifndef MFEM_NCMESH_OLD_VERTEX_ORDERING
2479 // This method assigns indices to vertices (Node::vert_index) that will be
2480 // seen by the Mesh class and the rest of MFEM. We must be careful to:
2481 //
2482 // 1. Stay compatible with the conforming code, which expects top-level
2483 // (original) vertices to be indexed first, otherwise GridFunctions
2484 // defined on a conforming mesh would no longer be valid when the mesh
2485 // is converted to an NC mesh.
2486 //
2487 // 2. Make sure serial NCMesh is compatible with the parallel ParNCMesh, so
2488 // it is possible to read parallel partial solutions in serial code
2489 // (e.g., serial GLVis). This means handling ghost elements, if present.
2490 //
2491 // 3. Assign vertices in a globally consistent order for parallel meshes:
2492 // if two vertices i,j are shared by two ranks r1,r2, and i<j on r1,
2493 // then i<j on r2 as well. This is true for top-level vertices but also
2494 // for the remaining shared vertices thanks to the globally consistent
2495 // SFC ordering of the leaf elements. This property reduces
2496 // communication and simplifies ParNCMesh.
2497
2498 // STEP 1: begin by splitting vertices into 4 classes:
2499 // - local top-level vertices (code -1)
2500 // - local non-top level vertices (code -2)
2501 // - ghost (non-local) vertices (code -3)
2502 // - vertices beyond the ghost layer (code -4)
2503
2504 for (auto & node : nodes)
2505 {
2506 node.vert_index = -4; // assume beyond ghost layer
2507 }
2508
2509 for (int i = 0; i < leaf_elements.Size(); i++)
2510 {
2511 Element &el = elements[leaf_elements[i]];
2512 for (int j = 0; j < GI[el.Geom()].nv; j++)
2513 {
2514 Node &nd = nodes[el.node[j]];
2515 if (el.rank == MyRank)
2516 {
2517 if (nd.p1 == nd.p2) // local top-level vertex
2518 {
2519 if (nd.vert_index < -1) { nd.vert_index = -1; }
2520 }
2521 else // local non-top-level vertex
2522 {
2523 if (nd.vert_index < -2) { nd.vert_index = -2; }
2524 }
2525 }
2526 else // ghost vertex
2527 {
2528 if (nd.vert_index < -3) { nd.vert_index = -3; }
2529 }
2530 }
2531 }
2532
2533 // STEP 2: assign indices of top-level local vertices, in original order
2534 NVertices = 0;
2535 for (auto &node : nodes)
2536 {
2537 if (node.vert_index == -1)
2538 {
2539 node.vert_index = NVertices++;
2540 }
2541 }
2542
2543 // STEP 3: go over all elements (local and ghost) in SFC order and assign
2544 // remaining local vertices in that order.
2545 Array<int> sfc_order(leaf_elements.Size());
2546 for (int i = 0; i < sfc_order.Size(); i++)
2547 {
2548 sfc_order[leaf_sfc_index[i]] = leaf_elements[i];
2549 }
2550
2551 for (int i = 0; i < sfc_order.Size(); i++)
2552 {
2553 const Element &el = elements[sfc_order[i]];
2554 for (int j = 0; j < GI[el.Geom()].nv; j++)
2555 {
2556 Node &nd = nodes[el.node[j]];
2557 if (nd.vert_index == -2) { nd.vert_index = NVertices++; }
2558 }
2559 }
2560
2561 // STEP 4: create the mapping from Mesh vertex index to NCMesh node index
2563 for (auto node = nodes.begin(); node != nodes.end(); ++node)
2564 {
2565 if (node->HasVertex() && node->vert_index >= 0)
2566 {
2567 MFEM_ASSERT(node->vert_index < vertex_nodeId.Size(), "");
2568 vertex_nodeId[node->vert_index] = node.index();
2569 }
2570 }
2571
2572 // STEP 5: assign remaining ghost vertices, ignore vertices beyond the ghost
2573 // layer
2574 NGhostVertices = 0;
2575 for (int i = 0; i < sfc_order.Size(); i++)
2576 {
2577 const Element &el = elements[sfc_order[i]];
2578 for (int j = 0; j < GI[el.Geom()].nv; j++)
2579 {
2580 Node &nd = nodes[el.node[j]];
2581 if (nd.vert_index == -3)
2582 {
2584 }
2585 }
2586 }
2587
2588#else // old ordering for debugging/testing only
2589 bool parallel = false;
2590#ifdef MFEM_USE_MPI
2591 if (dynamic_cast<ParNCMesh*>(this)) { parallel = true; }
2592#endif
2593
2594 if (!parallel)
2595 {
2596 NVertices = 0;
2597 for (auto node = nodes.begin(); node != nodes.end(); ++node)
2598 {
2599 if (node->HasVertex()) { node->vert_index = NVertices++; }
2600 }
2601
2603
2604 NVertices = 0;
2605 for (auto node = nodes.begin(); node != nodes.end(); ++node)
2606 {
2607 if (node->HasVertex()) { vertex_nodeId[NVertices++] = node.index(); }
2608 }
2609 }
2610 else
2611 {
2612 for (auto node = nodes.begin(); node != nodes.end(); ++node)
2613 {
2614 if (node->HasVertex()) { node->vert_index = -1; }
2615 }
2616
2617 NVertices = 0;
2618 for (int i = 0; i < leaf_elements.Size(); i++)
2619 {
2620 Element &el = elements[leaf_elements[i]];
2621 if (el.rank == MyRank)
2622 {
2623 for (int j = 0; j < GI[el.Geom()].nv; j++)
2624 {
2625 int &vindex = nodes[el.node[j]].vert_index;
2626 if (vindex < 0) { vindex = NVertices++; }
2627 }
2628 }
2629 }
2630
2632 for (auto &node : nodes)
2633 {
2634 if (node.HasVertex() && node.vert_index >= 0)
2635 {
2636 vertex_nodeId[node.vert_index] = node.index();
2637 }
2638 }
2639
2640 NGhostVertices = 0;
2641 for (auto &node : nodes)
2642 {
2643 if (node.HasVertex() && node.vert_index < 0)
2644 {
2645 node.vert_index = NVertices + (NGhostVertices++);
2646 }
2647 }
2648 }
2649#endif
2650}
2651
2652void NCMesh::InitRootState(int root_count)
2653{
2654 root_state.SetSize(root_count);
2655 root_state = 0;
2656
2657 if (elements.Size() == 0) { return; }
2658
2659 char* node_order;
2660 int nch;
2661
2662 switch (elements[0].Geom()) // TODO: mixed meshes
2663 {
2664 case Geometry::SQUARE:
2665 nch = 4;
2666 node_order = (char*) quad_hilbert_child_order;
2667 break;
2668
2669 case Geometry::CUBE:
2670 nch = 8;
2671 node_order = (char*) hex_hilbert_child_order;
2672 break;
2673
2674 default:
2675 return; // do nothing, all states stay zero
2676 }
2677
2678 int entry_node = -2;
2679
2680 // process the root element sequence
2681 for (int i = 0; i < root_count; i++)
2682 {
2683 Element &el = elements[i];
2684
2685 int v_in = FindNodeExt(el, entry_node, false);
2686 if (v_in < 0) { v_in = 0; }
2687
2688 // determine which nodes are shared with the next element
2689 bool shared[MaxElemNodes];
2690 for (int ni = 0; ni < MaxElemNodes; ++ni) { shared[ni] = 0; }
2691 if (i+1 < root_count)
2692 {
2693 Element &next = elements[i+1];
2694 for (int j = 0; j < nch; j++)
2695 {
2696 int node = FindNodeExt(el, RetrieveNode(next, j), false);
2697 if (node >= 0) { shared[node] = true; }
2698 }
2699 }
2700
2701 // select orientation that starts in v_in and exits in shared node
2702 int state = Dim*v_in;
2703 for (int j = 0; j < Dim; j++)
2704 {
2705 if (shared[(int) node_order[nch*(state + j) + nch-1]])
2706 {
2707 state += j;
2708 break;
2709 }
2710 }
2711
2712 root_state[i] = state;
2713
2714 entry_node = RetrieveNode(el, node_order[nch*state + nch-1]);
2715 }
2716}
2717
2719{
2720 switch (geom)
2721 {
2722 case Geometry::CUBE: return new mfem::Hexahedron;
2723 case Geometry::PRISM: return new mfem::Wedge;
2724 case Geometry::PYRAMID: return new mfem::Pyramid;
2726 case Geometry::SQUARE: return new mfem::Quadrilateral;
2727 case Geometry::TRIANGLE: return new mfem::Triangle;
2728 }
2729 MFEM_ABORT("invalid geometry");
2730 return NULL;
2731}
2732
2733const real_t* NCMesh::CalcVertexPos(int node) const
2734{
2735 const Node &nd = nodes[node];
2736 if (nd.p1 == nd.p2) // top-level vertex
2737 {
2738 return &coordinates[3*nd.p1];
2739 }
2740
2741 TmpVertex &tv = tmp_vertex[node];
2742 if (tv.valid) { return tv.pos; }
2743
2744 MFEM_VERIFY(tv.visited == false, "cyclic vertex dependencies.");
2745 tv.visited = true;
2746
2747 const real_t* pos1 = CalcVertexPos(nd.p1);
2748 const real_t* pos2 = CalcVertexPos(nd.p2);
2749
2750 for (int i = 0; i < 3; i++)
2751 {
2752 tv.pos[i] = ((1.0 - nd.GetScale()) * pos1[i]) + (nd.GetScale() * pos2[i]);
2753 }
2754 tv.valid = true;
2755 return tv.pos;
2756}
2757
2759{
2760 mesh.vertices.SetSize(vertex_nodeId.Size());
2761 if (coordinates.Size())
2762 {
2763 // calculate vertex positions from stored top-level vertex coordinates
2764 tmp_vertex = new TmpVertex[nodes.NumIds()];
2765 for (int i = 0; i < mesh.vertices.Size(); i++)
2766 {
2767 mesh.vertices[i].SetCoords(spaceDim, CalcVertexPos(vertex_nodeId[i]));
2768 }
2769 delete [] tmp_vertex;
2770 }
2771 // NOTE: if the mesh is curved ('coordinates' is empty), mesh.vertices are
2772 // left uninitialized here; they will be initialized later by the Mesh from
2773 // Nodes -- here we just make sure mesh.vertices has the correct size.
2774
2775 for (auto &elem : mesh.elements)
2776 {
2777 mesh.FreeElement(elem);
2778 }
2779 mesh.elements.SetSize(0);
2780
2781 for (auto &elem : mesh.boundary)
2782 {
2783 mesh.FreeElement(elem);
2784 }
2785 mesh.boundary.SetSize(0);
2786
2787 // Save off boundary face vertices to make boundary elements later.
2788 std::map<int, mfem::Array<int>> unique_boundary_faces;
2789
2790 // create an mfem::Element for each leaf Element
2791 for (int i = 0; i < NElements; i++)
2792 {
2793 const Element &nc_elem = elements[leaf_elements[i]];
2794
2795 const int* node = nc_elem.node;
2796 GeomInfo& gi = GI[(int) nc_elem.geom];
2797
2798 mfem::Element* elem = mesh.NewElement(nc_elem.geom);
2799 mesh.elements.Append(elem);
2800
2801 elem->SetAttribute(nc_elem.attribute);
2802 for (int j = 0; j < gi.nv; j++)
2803 {
2804 elem->GetVertices()[j] = nodes[node[j]].vert_index;
2805 }
2806
2807 // Loop over faces and collect those marked as boundaries
2808 for (int k = 0; k < gi.nf; ++k)
2809 {
2810 const int nfv = gi.nfv[k];
2811 const int * const fv = gi.faces[k];
2812 const auto id = faces.FindId(node[fv[0]], node[fv[1]], node[fv[2]],
2813 node[fv[3]]);
2814 if (id >= 0 && faces[id].Boundary())
2815 {
2816 const auto &face = faces[id];
2817 if (face.elem[0] >= 0 && face.elem[1] >= 0 &&
2818 nc_elem.rank != std::min(elements[face.elem[0]].rank,
2819 elements[face.elem[1]].rank))
2820 {
2821 // This is a conformal internal face, but this element is not the
2822 // lowest ranking attached processor, thus not the owner of the
2823 // face. Consequently, we do not add this face to avoid double
2824 // counting.
2825 continue;
2826 }
2827
2828 // Add in all boundary faces that are actual boundaries or not
2829 // masters of another face. The fv[2] in the edge split is on
2830 // purpose. A point cannot have a split level, thus do not check for
2831 // master/slave relation.
2832 if ((nfv == 4 &&
2833 !QuadFaceIsMaster(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]))
2834 || (nfv == 3 && !TriFaceIsMaster(node[fv[0]], node[fv[1]], node[fv[2]]))
2835 || (nfv == 2 && EdgeSplitLevel(node[fv[0]], node[fv[2]]) == 0) || (nfv == 1))
2836 {
2837 // This face has no split faces below, it is conformal or a
2838 // slave.
2839 unique_boundary_faces[id].SetSize(nfv);
2840 for (int v = 0; v < nfv; ++v)
2841 {
2842 // Using a map overwrites if a face is visited twice. The
2843 // nfv==2 is necessary because faces of 2D are storing the
2844 // second index in the 2 slot, not the 1 slot.
2845 unique_boundary_faces[id][v] = nodes[node[fv[(nfv==2) ? 2*v : v]]].vert_index;
2846 }
2847 }
2848 }
2849 }
2850 }
2851
2852 auto geom_from_nfv = [](int nfv)
2853 {
2854 switch (nfv)
2855 {
2856 case 1: return Geometry::POINT;
2857 case 2: return Geometry::SEGMENT;
2858 case 3: return Geometry::TRIANGLE;
2859 case 4: return Geometry::SQUARE;
2860 }
2861 return Geometry::INVALID;
2862 };
2863
2864 for (const auto &fv : unique_boundary_faces)
2865 {
2866 const auto f = fv.first;
2867 const auto &v = fv.second;
2868 const auto &face = faces.At(f);
2869
2870 auto geom = geom_from_nfv(v.Size());
2871
2872 MFEM_ASSERT(geom != Geometry::INVALID,
2873 "nfv: " << v.Size() <<
2874 " does not match a valid face geometry: Quad, Tri, Segment, Point");
2875
2876 // Add a new boundary element, with matching attribute and vertices
2877 mesh.boundary.Append(mesh.NewElement(geom));
2878 auto * const be = mesh.boundary.Last();
2879 be->SetAttribute(face.attribute);
2880 be->SetVertices(v);
2881 }
2882}
2883
2884
2886{
2887 //// PART 1: pull indices of regular edges/faces from the Mesh
2888
2889 NEdges = mesh->GetNEdges();
2890 NFaces = mesh->GetNumFaces();
2891 if (Dim < 2) { NFaces = 0; }
2892 // clear Node::edge_index and Face::index
2893 for (auto &node : nodes)
2894 {
2895 if (node.HasEdge()) { node.edge_index = -1; }
2896 }
2897 for (auto &face : faces)
2898 {
2899 face.index = -1;
2900 }
2901
2902 // get edge enumeration from the Mesh
2903 Table *edge_vertex = mesh->GetEdgeVertexTable();
2904 for (int i = 0; i < edge_vertex->Size(); i++)
2905 {
2906 const int *ev = edge_vertex->GetRow(i);
2907 Node* node = nodes.Find(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
2908 MFEM_ASSERT(node && node->HasEdge(),
2909 "edge (" << ev[0] << "," << ev[1] << ") not found, "
2910 "node = " << node << " node->HasEdge() "
2911 << (node != nullptr ? node->HasEdge() : false));
2912 node->edge_index = i;
2913 }
2914
2915 // get face enumeration from the Mesh, initialize 'face_geom'
2917 for (int i = 0; i < NFaces; i++)
2918 {
2919 const int* fv = mesh->GetFace(i)->GetVertices();
2920 const int nfv = mesh->GetFace(i)->GetNVertices();
2921
2922 Face* face;
2923 if (Dim == 3)
2924 {
2925 if (nfv == 4)
2926 {
2928 face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2929 vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
2930 }
2931 else
2932 {
2933 MFEM_ASSERT(nfv == 3, "");
2935 face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2936 vertex_nodeId[fv[2]]);
2937 }
2938 }
2939 else
2940 {
2941 MFEM_ASSERT(nfv == 2, "");
2943 int n0 = vertex_nodeId[fv[0]], n1 = vertex_nodeId[fv[1]];
2944 face = faces.Find(n0, n0, n1, n1); // look up degenerate face
2945
2946#ifdef MFEM_DEBUG
2947 // (non-ghost) edge and face numbers must match in 2D
2948 const int *ev = edge_vertex->GetRow(i);
2949 MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
2950 (ev[1] == fv[0] && ev[0] == fv[1]), "");
2951#endif
2952 }
2953
2954 MFEM_VERIFY(face, "face not found.");
2955 face->index = i;
2956 }
2957
2958 //// PART 2: assign indices of ghost edges/faces, if any
2959
2960 // count ghost edges and assign their indices
2961 NGhostEdges = 0;
2962 for (auto &node : nodes)
2963 {
2964 if (node.HasEdge() && node.edge_index < 0)
2965 {
2966 node.edge_index = NEdges + (NGhostEdges++);
2967 }
2968 }
2969
2970 // count ghost faces
2971 NGhostFaces = 0;
2972 for (auto &face : faces)
2973 {
2974 if (face.index < 0) { NGhostFaces++; }
2975 }
2976
2977 if (Dim == 2)
2978 {
2979 // in 2D we have fake faces because of DG
2980 MFEM_ASSERT(NFaces == NEdges, "");
2981 MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
2982 }
2983
2984 // resize face_geom (default_geom is for slave faces beyond the ghost layer)
2985 Geometry::Type default_geom = Geometry::SQUARE;
2986 face_geom.SetSize(NFaces + NGhostFaces, default_geom);
2987
2988 // update 'face_geom' for ghost faces, assign ghost face indices
2989 int nghosts = 0;
2990 for (int i = 0; i < NGhostElements; i++)
2991 {
2992 Element &el = elements[leaf_elements[NElements + i]]; // ghost element
2993 GeomInfo &gi = GI[el.Geom()];
2994
2995 for (int j = 0; j < gi.nf; j++)
2996 {
2997 const int *fv = gi.faces[j];
2998 const int fid = faces.FindId(el.node[fv[0]], el.node[fv[1]],
2999 el.node[fv[2]], el.node[fv[3]]);
3000 MFEM_ASSERT(fid >= 0, "face not found!");
3001 auto &face = faces[fid];
3002
3003 if (face.index < 0)
3004 {
3005 face.index = NFaces + (nghosts++);
3006 // store the face geometry
3007 static const Geometry::Type types[5] =
3008 {
3011 };
3012 face_geom[face.index] = types[gi.nfv[j]];
3013 }
3014 }
3015 }
3016
3017 // assign valid indices also to faces beyond the ghost layer
3018 for (auto &face : faces)
3019 {
3020 if (face.index < 0) { face.index = NFaces + (nghosts++); }
3021 }
3022 MFEM_ASSERT(nghosts == NGhostFaces, "");
3023}
3024
3025
3026//// Face/edge lists ///////////////////////////////////////////////////////////
3027
3028int NCMesh::QuadFaceSplitType(int v1, int v2, int v3, int v4, real_t & s,
3029 int mid[5]) const
3030{
3031 MFEM_ASSERT(Dim >= 3, "");
3032
3033 // find edge nodes
3034 const int e1 = FindMidEdgeNode(v1, v2);
3035 const int e2 = FindMidEdgeNode(v2, v3);
3036 const int e3 = (e1 >= 0 &&
3037 nodes[e1].HasVertex()) ? FindMidEdgeNode(v3, v4) : -1;
3038 const int e4 = (e2 >= 0 &&
3039 nodes[e2].HasVertex()) ? FindMidEdgeNode(v4, v1) : -1;
3040
3041 // optional: return the mid-edge nodes if requested
3042 if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
3043
3044 // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
3045 int midf1 = -1, midf2 = -1;
3046 if (e1 >= 0 && e3 >= 0) { midf1 = FindMidEdgeNode(e1, e3); }
3047 if (e2 >= 0 && e4 >= 0) { midf2 = FindMidEdgeNode(e2, e4); }
3048
3049 // get proper node if shadow node exists
3050 if (midf1 >= 0 && midf1 == midf2)
3051 {
3052 const Node &nd = nodes[midf1];
3053 if (nd.p1 != e1 && nd.p2 != e1) { midf1 = -1; }
3054 if (nd.p1 != e2 && nd.p2 != e2) { midf2 = -1; }
3055 }
3056
3057 // only one way to access the mid-face node must always exist
3058 MFEM_ASSERT(!(midf1 >= 0 && midf2 >= 0), "incorrectly split face!");
3059
3060 if (midf1 < 0 && midf2 < 0) // face not split
3061 {
3062 if (mid) { mid[4] = -1; }
3063 s = 0.0;
3064 return 0;
3065 }
3066 else if (midf1 >= 0) // face split "vertically"
3067 {
3068 if (mid) { mid[4] = midf1; }
3069 s = nodes[e1].GetScale();
3070 if (v1 > v2) { s = 1.0 - s; }
3071 return 1;
3072 }
3073 else // face split "horizontally"
3074 {
3075 if (mid) { mid[4] = midf2; }
3076 s = nodes[e2].GetScale();
3077 if (v2 > v3) { s = 1.0 - s; }
3078 return 2;
3079 }
3080}
3081
3082bool NCMesh::TriFaceSplit(int v1, int v2, int v3, int mid[3]) const
3083{
3084 int e1 = nodes.FindId(v1, v2);
3085 if (e1 < 0 || !nodes[e1].HasVertex()) { return false; }
3086
3087 int e2 = nodes.FindId(v2, v3);
3088 if (e2 < 0 || !nodes[e2].HasVertex()) { return false; }
3089
3090 int e3 = nodes.FindId(v3, v1);
3091 if (e3 < 0 || !nodes[e3].HasVertex()) { return false; }
3092
3093 if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3; }
3094
3095 // This is necessary but not sufficient to determine if a face has been
3096 // split. All edges might have been split due to edge attached faces being
3097 // refined. Need to check for existence of face made up of midpoints.
3098 return true;
3099}
3100
3101bool contains_node(const std::array<int, 4> &nodes, int n)
3102{
3103 return std::find(nodes.begin(), nodes.end(), n) != nodes.end();
3104};
3105
3106int NCMesh::ParentFaceNodes(std::array<int, 4> &face_nodes) const
3107{
3108 const bool is_tri = face_nodes[3] == -1;
3109 const bool is_segment = (face_nodes[0] == face_nodes[1] &&
3110 face_nodes[2] == face_nodes[3]);
3111 const bool is_quad = *std::min_element(face_nodes.begin(),
3112 face_nodes.end()) >= 0;
3113
3114 MFEM_ASSERT((is_tri && !is_segment && !is_quad)
3115 || (!is_tri && is_segment && !is_quad) || (!is_tri && !is_segment &&
3116 is_quad), "Inconsistent node geometry");
3117
3118 bool all_nodes_root = true;
3119 for (auto x : face_nodes)
3120 {
3121 all_nodes_root = all_nodes_root && (x < 0 || (nodes[x].p1 == nodes[x].p2));
3122 }
3123 // This face is a root face -> nothing to do.
3124 if (all_nodes_root) { return -1; }
3125
3126 int child = -1; // The index into parent.child that this face corresponds to.
3127 auto parent_nodes = face_nodes;
3128 if (is_quad)
3129 {
3130 // Logic for coarsening anisotropic faces is more complex, needs
3131 // identification and handling of multiple "crux" points. Will require
3132 // inspection of edge nodes.
3133 MFEM_VERIFY(Iso,
3134 "ParentFaceNodes does not support anisotropic refinement yet!");
3135
3136 // Finds the first node whose parents aren't in the face_nodes. This is
3137 // also the index of the child location in the parent face. Treated
3138 // separately as ultimately multiple crux will need to be handled for
3139 // anisotropic faces.
3140 const auto crux = [&]()
3141 {
3142 for (int i = 0; i < static_cast<int>(face_nodes.size()); i++)
3143 {
3144 if ((!contains_node(face_nodes, nodes[face_nodes[i]].p1)
3145 && !contains_node(face_nodes, nodes[face_nodes[i]].p2))
3146 || (nodes[face_nodes[i]].p1 == nodes[face_nodes[i]].p2) /* top level node */)
3147 {
3148 return i;
3149 }
3150 }
3151 return -1;
3152 }();
3153 MFEM_ASSERT(crux != -1, "A root face should have been returned early");
3154
3155 // Loop over nodes, starting from diagonal to child, wrapping and skipping
3156 // child. This will visit the node opposite child twice, thereby
3157 // coarsening to the diagonally opposite. NOTE: This assumes that the
3158 // nodes for a square are numbered (0 -> 1 -> 2 -> 3 -> 0).
3159 for (int i = 0; i < static_cast<int>(face_nodes.size()) + 1; i++)
3160 {
3161 int ind = (crux + i + 2) %
3162 4; // Start and end with coarsening of the diagonally opposite
3163 if (ind == crux) { continue; }
3164 auto &x = parent_nodes[ind];
3165
3166 // Check against parent_nodes rather than face_nodes so on second lap
3167 // the node opposite crux will coarsen again to the diagonally across
3168 // in the parent face. A top level node has p1 == p2, thus these
3169 // modifications do nothing.
3170 if (contains_node(parent_nodes, nodes[x].p1))
3171 {
3172 MFEM_ASSERT(nodes[x].p2 == nodes[x].p1 ||
3173 !contains_node(parent_nodes, nodes[x].p2), "!");
3174 x = nodes[x].p2;
3175 }
3176 else if (contains_node(parent_nodes, nodes[x].p2))
3177 {
3178 MFEM_ASSERT(nodes[x].p2 == nodes[x].p1 ||
3179 !contains_node(parent_nodes, nodes[x].p1), "!");
3180 x = nodes[x].p1;
3181 }
3182 else { /* do nothing */ }
3183 }
3184 }
3185 else if (is_tri)
3186 {
3187 for (int i = 0; i < 3; i++)
3188 {
3189 auto x = face_nodes[i];
3190 if (x == -1) { continue; }
3191 if (contains_node(face_nodes, nodes[x].p1))
3192 {
3193 MFEM_ASSERT(nodes[x].p2 == nodes[x].p1 ||
3194 !contains_node(face_nodes, nodes[x].p2), "!");
3195 parent_nodes[i] = nodes[x].p2;
3196 }
3197 else if (contains_node(face_nodes, nodes[x].p2))
3198 {
3199 MFEM_ASSERT(nodes[x].p2 == nodes[x].p1 ||
3200 !contains_node(face_nodes, nodes[x].p1), "!");
3201 parent_nodes[i] = nodes[x].p1;
3202 }
3203 else { /* do nothing */ }
3204 }
3205
3206 if (std::equal(face_nodes.begin(), face_nodes.end(), parent_nodes.begin()))
3207 {
3208 // Having excluded root faces, this must be an interior face. We need
3209 // to handle the special case of the interior face of the parent face.
3210 std::array<std::array<int, 2>, 6> parent_pairs;
3211 for (std::size_t i = 0; i < face_nodes.size() - 1; i++)
3212 {
3213 parent_pairs[i][0] = nodes[face_nodes[i]].p1;
3214 parent_pairs[i][1] = nodes[face_nodes[i]].p2;
3215 }
3216 // Each node gets mapped to the common node from its parents and the
3217 // predecessor node's parents.
3218 for (int i = 0; i < 3; i++)
3219 {
3220 // Parenting convention here assumes parent face has the SAME
3221 // orientation as the original. This is true on exterior boundaries,
3222 // but for an interior boundary the master face will have an
3223 // opposing orientation. TODO: Possibly fix for interior boundaries.
3224 const auto &prev = parent_pairs[(i - 1 + 3) % 3]; // (0 -> 2, 1 -> 0, 2 -> 1)
3225 const auto &next = parent_pairs[(i + 1 + 3) % 3]; // (0 -> 1, 1 -> 2, 2 -> 0)
3226 for (auto x : next)
3227 {
3228 if (std::find(prev.begin(), prev.end(), x) != prev.end()) { parent_nodes[i] = x; }
3229 }
3230 }
3231 child = 3; // The interior face is the final child.
3232 }
3233 }
3234 else if (is_segment)
3235 {
3236 // Given this isn't a root face, one node must be the parent of the other.
3237 if (face_nodes[0] == nodes[face_nodes[1]].p1)
3238 {
3239 face_nodes[1] = nodes[face_nodes[1]].p2;
3240 }
3241 else if (face_nodes[0] == nodes[face_nodes[1]].p2)
3242 {
3243 face_nodes[1] = nodes[face_nodes[1]].p1;
3244 }
3245 else if (face_nodes[1] == nodes[face_nodes[0]].p1)
3246 {
3247 face_nodes[0] = nodes[face_nodes[0]].p2;
3248 }
3249 else if (face_nodes[1] == nodes[face_nodes[0]].p2)
3250 {
3251 face_nodes[0] = nodes[face_nodes[0]].p1;
3252 }
3253 else
3254 {
3255 MFEM_ABORT("Internal logic error!");
3256 }
3257 }
3258 else
3259 {
3260 MFEM_ABORT("Unrecognized face geometry!");
3261 }
3262 for (int i = 0; i < 4 && face_nodes[i] >= 0; i++)
3263 {
3264 if (face_nodes[i] == parent_nodes[i])
3265 {
3266 MFEM_ASSERT(child == -1,
3267 "This face cannot be more than one child of the parent face!");
3268 child = i;
3269 }
3270 }
3271 MFEM_ASSERT(child != -1, "Root elements must have exited early!");
3272 std::swap(face_nodes, parent_nodes);
3273 return child;
3274}
3275
3276int NCMesh::find_node(const Element &el, int node)
3277{
3278 for (int i = 0; i < MaxElemNodes; i++)
3279 {
3280 if (el.node[i] == node) { return i; }
3281 }
3282 MFEM_ABORT("Node not found.");
3283 return -1;
3284}
3285
3286int NCMesh::FindNodeExt(const Element &el, int node, bool abort)
3287{
3288 for (int i = 0; i < GI[el.Geom()].nv; i++)
3289 {
3290 if (RetrieveNode(el, i) == node) { return i; }
3291 }
3292 if (abort) { MFEM_ABORT("Node not found."); }
3293 return -1;
3294}
3295
3296int NCMesh::find_element_edge(const Element &el, int vn0, int vn1, bool abort)
3297{
3298 MFEM_ASSERT(!el.ref_type, "");
3299
3300 GeomInfo &gi = GI[el.Geom()];
3301 for (int i = 0; i < gi.ne; i++)
3302 {
3303 const int* ev = gi.edges[i];
3304 int n0 = el.node[ev[0]];
3305 int n1 = el.node[ev[1]];
3306 if ((n0 == vn0 && n1 == vn1) ||
3307 (n0 == vn1 && n1 == vn0)) { return i; }
3308 }
3309
3310 if (abort) { MFEM_ABORT("Edge (" << vn0 << ", " << vn1 << ") not found"); }
3311 return -1;
3312}
3313
3314int NCMesh::find_local_face(int geom, int a, int b, int c)
3315{
3316 GeomInfo &gi = GI[geom];
3317 for (int i = 0; i < gi.nf; i++)
3318 {
3319 const int* fv = gi.faces[i];
3320 if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
3321 (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
3322 (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
3323 {
3324 return i;
3325 }
3326 }
3327 MFEM_ABORT("Face not found.");
3328 return -1;
3329}
3330
3331namespace
3332{
3333template <typename T> struct IntHash;
3334template <> struct IntHash<float>
3335{
3336 using int_type = uint32_t;
3337 static constexpr int_type initial_value = 0xc4a016dd; // random value;
3338};
3339template <> struct IntHash<double>
3340{
3341 using int_type = uint64_t;
3342 static constexpr int_type initial_value = 0xf9ca9ba106acbba9; // random value
3343};
3344}
3345
3346/// Hash function for a PointMatrix, used in MatrixMap::map.
3347struct PointMatrixHash
3348{
3349 std::size_t operator()(const NCMesh::PointMatrix &pm) const
3350 {
3351 // This is a variation on "Hashing an array of floats" from here:
3352 // https://cs.stackexchange.com/questions/37952
3353
3354 // Make sure (at compile time) that the types have compatible sizes
3355 static_assert(sizeof(IntHash<real_t>::int_type) == sizeof(real_t), "");
3356 // Suppress maybe unused warnings
3357 MFEM_CONTRACT_VAR(IntHash<float>::initial_value);
3358 MFEM_CONTRACT_VAR(IntHash<double>::initial_value);
3359
3360 auto int_bit_cast = [](real_t val)
3361 {
3362 // std::memcpy is the proper way of doing type punning, see e.g.
3363 // https://gist.github.com/shafik/848ae25ee209f698763cffee272a58f8
3364 IntHash<real_t>::int_type int_val;
3365 std::memcpy(&int_val, &val, sizeof(real_t));
3366 return int_val;
3367 };
3368
3369 IntHash<real_t>::int_type hash = IntHash<real_t>::initial_value;
3370
3371 for (int i = 0; i < pm.np; i++)
3372 {
3373 for (int j = 0; j < pm.points[i].dim; j++)
3374 {
3375 // mix the doubles by adding their binary representations many times
3376 // over (note: 31 is 11111 in binary)
3377 real_t coord = pm.points[i].coord[j];
3378 hash = 31*hash + int_bit_cast(coord);
3379 }
3380 }
3381 return hash; // return the lowest bits of the huge sum
3382 }
3383};
3384
3385/** Helper container to keep track of point matrices encountered during
3386 * face/edge traversal and to assign unique indices to them.
3387 */
3388struct MatrixMap
3389{
3390 int GetIndex(const NCMesh::PointMatrix &pm)
3391 {
3392 int &index = map[pm];
3393 if (!index) { index = static_cast<int>(map.size()); }
3394 return index - 1;
3395 }
3396
3397 void ExportMatrices(Array<DenseMatrix*> &point_matrices) const
3398 {
3399 point_matrices.SetSize(static_cast<int>(map.size()));
3400 for (const auto &pair : map)
3401 {
3402 DenseMatrix* mat = new DenseMatrix();
3403 pair.first.GetMatrix(*mat);
3404 point_matrices[pair.second - 1] = mat;
3405 }
3406 }
3407
3408 void DumpBucketSizes() const
3409 {
3410 for (unsigned i = 0; i < map.bucket_count(); i++)
3411 {
3412 mfem::out << map.bucket_size(i) << " ";
3413 }
3414 }
3415
3416private:
3417 std::unordered_map<NCMesh::PointMatrix, int, PointMatrixHash> map;
3418};
3419
3420
3421int NCMesh::ReorderFacePointMat(int v0, int v1, int v2, int v3,
3422 int elem, const PointMatrix &pm,
3423 PointMatrix &reordered) const
3424{
3425 const Element &el = elements[elem];
3426 int master[4] =
3427 {
3428 find_node(el, v0), find_node(el, v1), find_node(el, v2),
3429 (v3 >= 0) ? find_node(el, v3) : -1
3430 };
3431 int nfv = (v3 >= 0) ? 4 : 3;
3432
3433 int local = find_local_face(el.Geom(), master[0], master[1], master[2]);
3434 const int* fv = GI[el.Geom()].faces[local];
3435
3436 reordered.np = pm.np;
3437 for (int i = 0, j; i < nfv; i++)
3438 {
3439 for (j = 0; j < nfv; j++)
3440 {
3441 if (fv[i] == master[j])
3442 {
3443 reordered.points[i] = pm.points[j];
3444 break;
3445 }
3446 }
3447 MFEM_ASSERT(j != nfv, "node not found.");
3448 }
3449 return local;
3450}
3451
3452void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
3453 const PointMatrix& pm, int level,
3454 Face* eface[4], MatrixMap &matrix_map)
3455{
3456 if (level > 0)
3457 {
3458 // check if we made it to a face that is not split further
3459 Face* fa = faces.Find(vn0, vn1, vn2, vn3);
3460 if (fa)
3461 {
3462 // we have a slave face, add it to the list
3463 int elem = fa->GetSingleElement();
3464 face_list.slaves.Append(
3465 Slave(fa->index, elem, -1, Geometry::SQUARE));
3466 Slave &sl = face_list.slaves.Last();
3467
3468 // reorder the point matrix according to slave face orientation
3469 PointMatrix pm_r;
3470 sl.local = ReorderFacePointMat(vn0, vn1, vn2, vn3, elem, pm, pm_r);
3471 sl.matrix = matrix_map.GetIndex(pm_r);
3472
3473 eface[0] = eface[2] = fa;
3474 eface[1] = eface[3] = fa;
3475
3476 return;
3477 }
3478 }
3479
3480 // we need to recurse deeper
3481 int mid[5];
3482 real_t scale;
3483 const int split = QuadFaceSplitType(vn0, vn1, vn2, vn3, scale, mid);
3484
3485 Face *ef[2][4];
3486 if (split == 1) // "X" split face
3487 {
3488 Point pmid0(pm(0), pm(1), scale), pmid2(pm(2), pm(3), 1.0 - scale);
3489
3490 TraverseQuadFace(vn0, mid[0], mid[2], vn3,
3491 PointMatrix(pm(0), pmid0, pmid2, pm(3)),
3492 level+1, ef[0], matrix_map);
3493
3494 TraverseQuadFace(mid[0], vn1, vn2, mid[2],
3495 PointMatrix(pmid0, pm(1), pm(2), pmid2),
3496 level+1, ef[1], matrix_map);
3497
3498 eface[1] = ef[1][1];
3499 eface[3] = ef[0][3];
3500 eface[0] = eface[2] = NULL;
3501 }
3502 else if (split == 2) // "Y" split face
3503 {
3504 Point pmid1(pm(1), pm(2), scale), pmid3(pm(3), pm(0), 1.0 - scale);
3505
3506 TraverseQuadFace(vn0, vn1, mid[1], mid[3],
3507 PointMatrix(pm(0), pm(1), pmid1, pmid3),
3508 level+1, ef[0], matrix_map);
3509
3510 TraverseQuadFace(mid[3], mid[1], vn2, vn3,
3511 PointMatrix(pmid3, pmid1, pm(2), pm(3)),
3512 level+1, ef[1], matrix_map);
3513
3514 eface[0] = ef[0][0];
3515 eface[2] = ef[1][2];
3516 eface[1] = eface[3] = NULL;
3517 }
3518
3519 // check for a prism edge constrained by the master face
3520 if (HavePrisms() && mid[4] >= 0)
3521 {
3522 Node& enode = nodes[mid[4]];
3523 if (enode.HasEdge())
3524 {
3525 // process the edge only if it's not shared by slave faces within this
3526 // master face (i.e. the edge is "hidden")
3527 const int fi[3][2] = {{0, 0}, {1, 3}, {2, 0}};
3528 if (!ef[0][fi[split][0]] && !ef[1][fi[split][1]])
3529 {
3530 MFEM_ASSERT(enode.edge_refc == 1, "");
3531
3532 MeshId buf[4];
3533 Array<MeshId> eid(buf, 4);
3534
3535 (split == 1) ? FindEdgeElements(mid[0], vn1, vn2, mid[2], eid)
3536 /* */ : FindEdgeElements(mid[3], vn0, vn1, mid[1], eid);
3537
3538 MFEM_ASSERT(eid.Size() > 0, "edge prism not found");
3539 MFEM_ASSERT(eid.Size() < 2, "non-unique edge prism");
3540
3541 // create a slave face record with a degenerate point matrix
3542 face_list.slaves.Append(
3543 Slave(-1 - enode.edge_index,
3544 eid[0].element, eid[0].local, Geometry::SQUARE));
3545 Slave &sl = face_list.slaves.Last();
3546
3547 if (split == 1)
3548 {
3549 Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
3550 int v1 = nodes[mid[0]].vert_index;
3551 int v2 = nodes[mid[2]].vert_index;
3552 sl.matrix =
3553 matrix_map.GetIndex(
3554 (v1 < v2) ? PointMatrix(mid0, mid2, mid2, mid0) :
3555 /* */ PointMatrix(mid2, mid0, mid0, mid2));
3556 }
3557 else
3558 {
3559 Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
3560 int v1 = nodes[mid[1]].vert_index;
3561 int v2 = nodes[mid[3]].vert_index;
3562 sl.matrix =
3563 matrix_map.GetIndex(
3564 (v1 < v2) ? PointMatrix(mid1, mid3, mid3, mid1) :
3565 /* */ PointMatrix(mid3, mid1, mid1, mid3));
3566 }
3567 }
3568 }
3569 }
3570}
3571
3572void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
3573 MatrixMap &matrix_map)
3574{
3575 int mid = nodes.FindId(vn0, vn1);
3576 if (mid < 0) { return; }
3577
3578 const Node &nd = nodes[mid];
3579 if (nd.HasEdge())
3580 {
3581 // check if the edge is already a master in 'edge_list'
3582 const auto eid_and_type = edge_list.GetMeshIdAndType(nd.edge_index);
3583 if (eid_and_type.type == NCList::MeshIdType::MASTER
3584 || eid_and_type.type == NCList::MeshIdType::CONFORMING)
3585 {
3586 // in this case we need to add an edge-face constraint, because the
3587 // non-slave edge is really a (face-)slave itself.
3588 const MeshId &eid = *eid_and_type.id;
3589 face_list.slaves.Append(
3590 Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
3591
3592 int v0index = nodes[vn0].vert_index;
3593 int v1index = nodes[vn1].vert_index;
3594
3595 face_list.slaves.Last().matrix =
3596 matrix_map.GetIndex((v0index < v1index) ? PointMatrix(p0, p1, p0)
3597 /* */ : PointMatrix(p1, p0, p1));
3598
3599 return; // no need to continue deeper
3600 }
3601 }
3602
3603 // recurse deeper
3604 Point pmid(p0, p1);
3605 TraverseTetEdge(vn0, mid, p0, pmid, matrix_map);
3606 TraverseTetEdge(mid, vn1, pmid, p1, matrix_map);
3607}
3608
3610 int vn2,
3611 const PointMatrix& pm, int level,
3612 MatrixMap &matrix_map)
3613{
3614 if (level > 0)
3615 {
3616 // check if we made it to a face that is not split further
3617 Face* fa = faces.Find(vn0, vn1, vn2);
3618 if (fa)
3619 {
3620 // we have a slave face, add it to the list
3621 int elem = fa->GetSingleElement();
3622 face_list.slaves.Append(
3623 Slave(fa->index, elem, -1, Geometry::TRIANGLE));
3624 Slave &sl = face_list.slaves.Last();
3625
3626 // reorder the point matrix according to slave face orientation
3627 PointMatrix pm_r;
3628 sl.local = ReorderFacePointMat(vn0, vn1, vn2, -1, elem, pm, pm_r);
3629 sl.matrix = matrix_map.GetIndex(pm_r);
3630
3631 return {true, elements[elem].rank != MyRank};
3632 }
3633 }
3634
3635 int mid[3];
3636 if (TriFaceSplit(vn0, vn1, vn2, mid))
3637 {
3638 Point pmid0(pm(0), pm(1)), pmid1(pm(1), pm(2)), pmid2(pm(2), pm(0));
3640
3641 b[0] = TraverseTriFace(vn0, mid[0], mid[2],
3642 PointMatrix(pm(0), pmid0, pmid2),
3643 level+1, matrix_map);
3644
3645 b[1] = TraverseTriFace(mid[0], vn1, mid[1],
3646 PointMatrix(pmid0, pm(1), pmid1),
3647 level+1, matrix_map);
3648
3649 b[2] = TraverseTriFace(mid[2], mid[1], vn2,
3650 PointMatrix(pmid2, pmid1, pm(2)),
3651 level+1, matrix_map);
3652
3653 b[3] = TraverseTriFace(mid[1], mid[2], mid[0],
3654 PointMatrix(pmid1, pmid2, pmid0),
3655 level+1, matrix_map);
3656
3657 // Traverse possible tet edges constrained by the master face. This needs
3658 // to occur if none of these first NC level faces are split further, OR if
3659 // they are on different processors. The different processor constraint is
3660 // needed in the case of local elements constrained by this face via the
3661 // edge alone. Cannot know this a priori, so just constrain any edge
3662 // attached to two neighbors.
3663 if (HaveTets() && (!b[3].unsplit || b[3].ghost_neighbor))
3664 {
3665 // If the faces have no further splits, so would not be captured by
3666 // normal face relations, add possible edge constraints.
3667 if (!b[1].unsplit || b[1].ghost_neighbor) { TraverseTetEdge(mid[0],mid[1], pmid0,pmid1, matrix_map); }
3668 if (!b[2].unsplit || b[2].ghost_neighbor) { TraverseTetEdge(mid[1],mid[2], pmid1,pmid2, matrix_map); }
3669 if (!b[0].unsplit || b[0].ghost_neighbor) { TraverseTetEdge(mid[2],mid[0], pmid2,pmid0, matrix_map); }
3670 }
3671 }
3672 return {false, false};
3673}
3674
3676{
3677 face_list.Clear();
3678 if (Dim < 3) { return; }
3679
3680 if (HaveTets()) { GetEdgeList(); } // needed by TraverseTetEdge()
3681
3683
3684 Array<char> processed_faces(faces.NumIds());
3685 processed_faces = 0;
3686
3687 MatrixMap matrix_maps[Geometry::NumGeom];
3688
3689 // visit faces of leaf elements
3690 for (int i = 0; i < leaf_elements.Size(); i++)
3691 {
3692 int elem = leaf_elements[i];
3693 Element &el = elements[elem];
3694 MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3695
3696 GeomInfo& gi = GI[el.Geom()];
3697 for (int j = 0; j < gi.nf; j++)
3698 {
3699 // get nodes for this face
3700 int node[4];
3701 for (int k = 0; k < 4; k++)
3702 {
3703 node[k] = el.node[gi.faces[j][k]];
3704 }
3705
3706 int face = faces.FindId(node[0], node[1], node[2], node[3]);
3707 MFEM_ASSERT(face >= 0, "face not found!");
3708
3709 // tell ParNCMesh about the face
3710 ElementSharesFace(elem, j, face);
3711
3712 // have we already processed this face? skip if yes
3713 if (processed_faces[face]) { continue; }
3714 processed_faces[face] = 1;
3715
3716 int fgeom = (node[3] >= 0) ? Geometry::SQUARE : Geometry::TRIANGLE;
3717
3718 Face &fa = faces[face];
3719 bool is_master = false;
3720 if (fa.elem[0] >= 0 && fa.elem[1] >= 0)
3721 {
3722 // this is a conforming face, add it to the list
3723 face_list.conforming.Append(MeshId(fa.index, elem, j, fgeom));
3724 }
3725 else
3726 {
3727 // this is either a master face or a slave face, but we can't tell
3728 // until we traverse the face refinement 'tree'...
3729 int sb = face_list.slaves.Size();
3730 if (fgeom == Geometry::SQUARE)
3731 {
3732 Face* dummy[4];
3733 TraverseQuadFace(node[0], node[1], node[2], node[3],
3734 pm_quad_identity, 0, dummy, matrix_maps[fgeom]);
3735 }
3736 else
3737 {
3738 TraverseTriFace(node[0], node[1], node[2],
3739 pm_tri_identity, 0, matrix_maps[fgeom]);
3740 }
3741
3742 int se = face_list.slaves.Size();
3743 if (sb < se)
3744 {
3745 // found slaves, so this is a master face; add it to the list
3746 is_master = true;
3747 face_list.masters.Append(
3748 Master(fa.index, elem, j, fgeom, sb, se));
3749
3750 // also, set the master index for the slaves
3751 for (int ii = sb; ii < se; ii++)
3752 {
3753 face_list.slaves[ii].master = fa.index;
3754 }
3755 }
3756 }
3757
3758 // To support internal boundaries can only insert non-master faces.
3759 if (fa.Boundary() && !is_master) { boundary_faces.Append(face); }
3760 }
3761 }
3762
3763 // export unique point matrices
3764 for (int i = 0; i < Geometry::NumGeom; i++)
3765 {
3766 matrix_maps[i].ExportMatrices(face_list.point_matrices[i]);
3767 }
3768}
3769
3770void NCMesh::TraverseEdge(int vn0, int vn1, real_t t0, real_t t1, int flags,
3771 int level, MatrixMap &matrix_map)
3772{
3773 int mid = nodes.FindId(vn0, vn1);
3774 if (mid < 0) { return; }
3775
3776 Node &nd = nodes[mid];
3777 if (nd.HasEdge() && level > 0)
3778 {
3779 // we have a slave edge, add it to the list
3780 edge_list.slaves.Append(Slave(nd.edge_index, -1, -1, Geometry::SEGMENT));
3781
3782 Slave &sl = edge_list.slaves.Last();
3783 sl.matrix = matrix_map.GetIndex(PointMatrix(Point(t0), Point(t1)));
3784
3785 // handle slave edge orientation
3786 sl.edge_flags = flags;
3787 int v0index = nodes[vn0].vert_index;
3788 int v1index = nodes[vn1].vert_index;
3789 if (v0index > v1index) { sl.edge_flags |= 2; }
3790 }
3791
3792 // recurse deeper
3793 const real_t scale = GetScale(nd.GetScale(), vn0 > vn1);
3794
3795 const real_t tmid = ((1.0 - scale) * t0) + (scale * t1);
3796 TraverseEdge(vn0, mid, t0, tmid, flags, level+1, matrix_map);
3797 TraverseEdge(mid, vn1, tmid, t1, flags, level+1, matrix_map);
3798}
3799
3801{
3802 edge_list.Clear();
3803 if (Dim < 3) { boundary_faces.SetSize(0); }
3804
3805 Array<char> processed_edges(nodes.NumIds());
3806 processed_edges = 0;
3807
3808 Array<int> edge_element(nodes.NumIds());
3809 Array<signed char> edge_local(nodes.NumIds());
3810 edge_local = -1;
3811
3812 MatrixMap matrix_map;
3813
3814 // visit edges of leaf elements
3815 for (int i = 0; i < leaf_elements.Size(); i++)
3816 {
3817 int elem = leaf_elements[i];
3818 Element &el = elements[elem];
3819 MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3820
3821 GeomInfo& gi = GI[el.Geom()];
3822 for (int j = 0; j < gi.ne; j++)
3823 {
3824 // get nodes for this edge
3825 const int* ev = gi.edges[j];
3826 int node[2] = { el.node[ev[0]], el.node[ev[1]] };
3827
3828 int enode = nodes.FindId(node[0], node[1]);
3829 MFEM_ASSERT(enode >= 0, "edge node not found!");
3830
3831 Node &nd = nodes[enode];
3832 MFEM_ASSERT(nd.HasEdge(), "edge not found!");
3833
3834 // tell ParNCMesh about the edge
3835 ElementSharesEdge(elem, j, enode);
3836
3837 // store element/local for later
3838 edge_element[nd.edge_index] = elem;
3839 edge_local[nd.edge_index] = j;
3840
3841 // skip slave edges here, they will be reached from their masters
3842 if (GetEdgeMaster(enode) >= 0)
3843 {
3844 // (2D only, store internal boundary faces)
3845 if (Dim <= 2)
3846 {
3847 int face = faces.FindId(node[0], node[0], node[1], node[1]);
3848 MFEM_ASSERT(face >= 0, "face not found!");
3849 if (faces[face].Boundary()) { boundary_faces.Append(face); }
3850 }
3851 continue;
3852 }
3853
3854 // have we already processed this edge? skip if yes
3855 if (processed_edges[enode]) { continue; }
3856 processed_edges[enode] = 1;
3857
3858 // prepare edge interval for slave traversal, handle orientation
3859 real_t t0 = 0.0, t1 = 1.0;
3860 int v0index = nodes[node[0]].vert_index;
3861 int v1index = nodes[node[1]].vert_index;
3862 int flags = (v0index > v1index) ? 1 : 0;
3863
3864 // try traversing the edge to find slave edges
3865 int sb = edge_list.slaves.Size();
3866 TraverseEdge(node[0], node[1], t0, t1, flags, 0, matrix_map);
3867
3868 int se = edge_list.slaves.Size();
3869 if (sb < se)
3870 {
3871 // found slaves, this is a master face; add it to the list
3872 edge_list.masters.Append(
3873 Master(nd.edge_index, elem, j, Geometry::SEGMENT, sb, se));
3874
3875 // also, set the master index for the slaves
3876 for (int ii = sb; ii < se; ii++)
3877 {
3878 edge_list.slaves[ii].master = nd.edge_index;
3879 }
3880 }
3881 else
3882 {
3883 // no slaves, this is a conforming edge
3884 edge_list.conforming.Append(MeshId(nd.edge_index, elem, j));
3885 // (2D only, store boundary faces)
3886 if (Dim <= 2)
3887 {
3888 int face = faces.FindId(node[0], node[0], node[1], node[1]);
3889 MFEM_ASSERT(face >= 0, "face not found!");
3890 if (faces[face].Boundary()) { boundary_faces.Append(face); }
3891 }
3892 }
3893 }
3894 }
3895
3896 // fix up slave edge element/local
3897 for (int i = 0; i < edge_list.slaves.Size(); i++)
3898 {
3899 Slave &sl = edge_list.slaves[i];
3900 int local = edge_local[sl.index];
3901 if (local >= 0)
3902 {
3903 sl.local = local;
3904 sl.element = edge_element[sl.index];
3905 }
3906 }
3907
3908 // export unique point matrices
3909 matrix_map.ExportMatrices(edge_list.point_matrices[Geometry::SEGMENT]);
3910}
3911
3913{
3914 int total = NVertices + NGhostVertices;
3915
3917 vertex_list.conforming.Reserve(total);
3918
3919 Array<char> processed_vertices(total);
3920 processed_vertices = 0;
3921
3922 // analogously to above, visit vertices of leaf elements
3923 for (int i = 0; i < leaf_elements.Size(); i++)
3924 {
3925 int elem = leaf_elements[i];
3926 Element &el = elements[elem];
3927
3928 for (int j = 0; j < GI[el.Geom()].nv; j++)
3929 {
3930 int node = el.node[j];
3931 Node &nd = nodes[node];
3932
3933 int index = nd.vert_index;
3934 if (index >= 0)
3935 {
3936 ElementSharesVertex(elem, j, node);
3937
3938 if (processed_vertices[index]) { continue; }
3939 processed_vertices[index] = 1;
3940
3941 vertex_list.conforming.Append(MeshId(index, elem, j));
3942 }
3943 }
3944 }
3945}
3946
3948 DenseMatrix &oriented_matrix) const
3949{
3950 oriented_matrix = *(point_matrices[slave.Geom()][slave.matrix]);
3951
3952 if (slave.edge_flags)
3953 {
3954 MFEM_ASSERT(oriented_matrix.Height() == 1 &&
3955 oriented_matrix.Width() == 2, "not an edge point matrix");
3956
3957 if (slave.edge_flags & 1) // master inverted
3958 {
3959 oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
3960 oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
3961 }
3962 if (slave.edge_flags & 2) // slave inverted
3963 {
3964 std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
3965 }
3966 }
3967}
3968
3970{
3971 conforming.DeleteAll();
3972 masters.DeleteAll();
3973 slaves.DeleteAll();
3974
3975 for (int i = 0; i < Geometry::NumGeom; i++)
3976 {
3977 for (int j = 0; j < point_matrices[i].Size(); j++)
3978 {
3979 delete point_matrices[i][j];
3980 }
3981 point_matrices[i].DeleteAll();
3982 }
3983
3984 inv_index.clear();
3985}
3986
3989{
3990 BuildIndex();
3991 const auto it = inv_index.find(index);
3992 auto ft = it != inv_index.end() ? it->second.first : MeshIdType::UNRECOGNIZED;
3993 switch (ft)
3994 {
3995 case MeshIdType::CONFORMING:
3996 return {&conforming[it->second.second], it->second.first};
3997 case MeshIdType::MASTER:
3998 return {&masters[it->second.second], it->second.first};
3999 case MeshIdType::SLAVE:
4000 return {&slaves[it->second.second], it->second.first};
4001 case MeshIdType::UNRECOGNIZED:
4002 default:
4003 return {nullptr, MeshIdType::UNRECOGNIZED};
4004 }
4005}
4006
4009{
4010 BuildIndex();
4011 auto it = inv_index.find(index);
4012 return (it != inv_index.end()) ? it->second.first : MeshIdType::UNRECOGNIZED;
4013}
4014
4015bool
4017{
4018 return GetMeshIdType(index) == ft;
4019}
4020
4021void
4022NCMesh::NCList::BuildIndex() const
4023{
4024 if (inv_index.size() == 0)
4025 {
4026 auto index_compare = [](const MeshId &a, const MeshId &b) { return a.index < b.index; };
4027 auto max_conforming = std::max_element(conforming.begin(), conforming.end(),
4028 index_compare);
4029 auto max_master = std::max_element(masters.begin(), masters.end(),
4030 index_compare);
4031 auto max_slave = std::max_element(slaves.begin(), slaves.end(), index_compare);
4032
4033 int max_conforming_index = max_conforming != nullptr ? max_conforming->index :
4034 -1;
4035 int max_master_index = max_master != nullptr ? max_master->index : -1;
4036 int max_slave_index = max_slave != nullptr ? max_slave->index : -1;
4037
4038 inv_index.reserve(max(max_conforming_index, max_master_index, max_slave_index,
4039 0));
4040 for (int i = 0; i < conforming.Size(); i++)
4041 {
4042 inv_index.emplace(conforming[i].index, std::make_pair(MeshIdType::CONFORMING,
4043 i));
4044 }
4045 for (int i = 0; i < masters.Size(); i++)
4046 {
4047 inv_index.emplace(masters[i].index, std::make_pair(MeshIdType::MASTER, i));
4048 }
4049 for (int i = 0; i < slaves.Size(); i++)
4050 {
4051 inv_index.emplace(slaves[i].index, std::make_pair(MeshIdType::SLAVE, i));
4052 }
4053 }
4054}
4055
4056//// Neighbors /////////////////////////////////////////////////////////////////
4057
4058void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
4059{
4060 int mid = nodes.FindId(v0, v1);
4061 if (mid >= 0 && nodes[mid].HasVertex())
4062 {
4063 indices.Append(mid);
4064
4065 CollectEdgeVertices(v0, mid, indices);
4066 CollectEdgeVertices(mid, v1, indices);
4067 }
4068}
4069
4070void NCMesh::CollectTriFaceVertices(int v0, int v1, int v2, Array<int> &indices)
4071{
4072 int mid[3];
4073 if (TriFaceSplit(v0, v1, v2, mid))
4074 {
4075 for (int i = 0; i < 3; i++)
4076 {
4077 indices.Append(mid[i]);
4078 }
4079
4080 CollectTriFaceVertices(v0, mid[0], mid[2], indices);
4081 CollectTriFaceVertices(mid[0], v1, mid[1], indices);
4082 CollectTriFaceVertices(mid[2], mid[1], v2, indices);
4083 CollectTriFaceVertices(mid[0], mid[1], mid[2], indices);
4084
4085 if (HaveTets()) // possible edge-face contact
4086 {
4087 CollectEdgeVertices(mid[0], mid[1], indices);
4088 CollectEdgeVertices(mid[1], mid[2], indices);
4089 CollectEdgeVertices(mid[2], mid[0], indices);
4090 }
4091 }
4092}
4093
4094void NCMesh::CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
4095 Array<int> &indices)
4096{
4097 int mid[5];
4098 real_t scale;
4099 switch (QuadFaceSplitType(v0, v1, v2, v3, scale, mid))
4100 {
4101 case 1:
4102 indices.Append(mid[0]);
4103 indices.Append(mid[2]);
4104
4105 CollectQuadFaceVertices(v0, mid[0], mid[2], v3, indices);
4106 CollectQuadFaceVertices(mid[0], v1, v2, mid[2], indices);
4107
4108 if (HavePrisms()) // possible edge-face contact
4109 {
4110 CollectEdgeVertices(mid[0], mid[2], indices);
4111 }
4112 break;
4113
4114 case 2:
4115 indices.Append(mid[1]);
4116 indices.Append(mid[3]);
4117
4118 CollectQuadFaceVertices(v0, v1, mid[1], mid[3], indices);
4119 CollectQuadFaceVertices(mid[3], mid[1], v2, v3, indices);
4120
4121 if (HavePrisms()) // possible edge-face contact
4122 {
4123 CollectEdgeVertices(mid[1], mid[3], indices);
4124 }
4125 break;
4126 }
4127}
4128
4130{
4131 int nrows = leaf_elements.Size();
4132 int* I = Memory<int>(nrows + 1);
4133 int** JJ = new int*[nrows];
4134
4135 Array<int> indices;
4136 indices.Reserve(128);
4137
4138 // collect vertices coinciding with each element, including hanging vertices
4139 for (int i = 0; i < leaf_elements.Size(); i++)
4140 {
4141 int elem = leaf_elements[i];
4142 Element &el = elements[elem];
4143 MFEM_ASSERT(!el.ref_type, "not a leaf element.");
4144
4145 GeomInfo& gi = GI[el.Geom()];
4146 int* node = el.node;
4147
4148 indices.SetSize(0);
4149 for (int j = 0; j < gi.ne; j++)
4150 {
4151 const int* ev = gi.edges[j];
4152 CollectEdgeVertices(node[ev[0]], node[ev[1]], indices);
4153 }
4154
4155 if (Dim >= 3)
4156 {
4157 for (int j = 0; j < gi.nf; j++)
4158 {
4159 const int* fv = gi.faces[j];
4160 if (gi.nfv[j] == 4)
4161 {
4162 CollectQuadFaceVertices(node[fv[0]], node[fv[1]],
4163 node[fv[2]], node[fv[3]], indices);
4164 }
4165 else
4166 {
4167 CollectTriFaceVertices(node[fv[0]], node[fv[1]], node[fv[2]],
4168 indices);
4169 }
4170 }
4171 }
4172
4173 // temporarily store one row of the table
4174 indices.Sort();
4175 indices.Unique();
4176 int size = indices.Size();
4177 I[i] = size;
4178 JJ[i] = new int[size];
4179 std::memcpy(JJ[i], indices.GetData(), size * sizeof(int));
4180 }
4181
4182 // finalize the I array of the table
4183 int nnz = 0;
4184 for (int i = 0; i < nrows; i++)
4185 {
4186 int cnt = I[i];
4187 I[i] = nnz;
4188 nnz += cnt;
4189 }
4190 I[nrows] = nnz;
4191
4192 // copy the temporarily stored rows into one J array
4193 int *J = Memory<int>(nnz);
4194 nnz = 0;
4195 for (int i = 0; i < nrows; i++)
4196 {
4197 int cnt = I[i+1] - I[i];
4198 std::memcpy(J+nnz, JJ[i], cnt * sizeof(int));
4199 delete [] JJ[i];
4200 nnz += cnt;
4201 }
4202
4203 element_vertex.SetIJ(I, J, nrows);
4204
4205 delete [] JJ;
4206}
4207
4208
4210 Array<int> *neighbors,
4211 Array<char> *neighbor_set)
4212{
4213 // If A is the element-to-vertex table (see 'element_vertex') listing all
4214 // vertices touching each element, including hanging vertices, then A*A^T is
4215 // the element-to-neighbor table. Multiplying the element set with A*A^T
4216 // gives the neighbor set. To save memory, this function only computes the
4217 // action of A*A^T, the product itself is not stored anywhere.
4218
4219 // Optimization: the 'element_vertex' table does not store the obvious corner
4220 // nodes in it. The table is therefore empty for conforming meshes.
4221
4223
4224 int nleaves = leaf_elements.Size();
4225 MFEM_VERIFY(elem_set.Size() == nleaves, "");
4226 MFEM_ASSERT(element_vertex.Size() == nleaves, "");
4227
4228 // step 1: vertices = A^T * elem_set, i.e, find all vertices touching the
4229 // element set
4230
4231 Array<char> vmark(nodes.NumIds());
4232 vmark = 0;
4233
4234 for (int i = 0; i < nleaves; i++)
4235 {
4236 if (elem_set[i])
4237 {
4238 int *v = element_vertex.GetRow(i);
4239 int nv = element_vertex.RowSize(i);
4240 for (int j = 0; j < nv; j++)
4241 {
4242 vmark[v[j]] = 1;
4243 }
4244
4245 Element &el = elements[leaf_elements[i]];
4246 nv = GI[el.Geom()].nv;
4247 for (int j = 0; j < nv; j++)
4248 {
4249 vmark[el.node[j]] = 1;
4250 }
4251 }
4252 }
4253
4254 // step 2: neighbors = A * vertices, i.e., find all elements coinciding with
4255 // vertices from step 1; NOTE: in the result we don't include elements from
4256 // the original set
4257
4258 if (neighbor_set)
4259 {
4260 neighbor_set->SetSize(nleaves);
4261 *neighbor_set = 0;
4262 }
4263
4264 for (int i = 0; i < nleaves; i++)
4265 {
4266 if (!elem_set[i])
4267 {
4268 bool hit = false;
4269
4270 int *v = element_vertex.GetRow(i);
4271 int nv = element_vertex.RowSize(i);
4272 for (int j = 0; j < nv; j++)
4273 {
4274 if (vmark[v[j]]) { hit = true; break; }
4275 }
4276
4277 if (!hit)
4278 {
4279 Element &el = elements[leaf_elements[i]];
4280 nv = GI[el.Geom()].nv;
4281 for (int j = 0; j < nv; j++)
4282 {
4283 if (vmark[el.node[j]]) { hit = true; break; }
4284 }
4285 }
4286
4287 if (hit)
4288 {
4289 if (neighbors) { neighbors->Append(leaf_elements[i]); }
4290 if (neighbor_set) { (*neighbor_set)[i] = 1; }
4291 }
4292 }
4293 }
4294}
4295
4296static bool sorted_lists_intersect(const int* a, const int* b, int na, int nb)
4297{
4298 // pointers to "end" sentinel, not last entry. Not for dereferencing.
4299 const int * const a_end = a + na;
4300 const int * const b_end = b + nb;
4301 while (a != a_end && b != b_end)
4302 {
4303 if (*a < *b)
4304 {
4305 ++a;
4306 }
4307 else if (*b < *a)
4308 {
4309 ++b;
4310 }
4311 else
4312 {
4313 return true; // neither *a < *b nor *b < *a thus a == b
4314 }
4315 }
4316 return false; // no common element found
4317}
4318
4319
4320void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
4321 const Array<int> *search_set)
4322{
4323 // TODO future: this function is inefficient. For a single element, an octree
4324 // neighbor search algorithm would be better. However, the octree neighbor
4325 // algorithm is hard to get right in the multi-octree case due to the
4326 // different orientations of the octrees (i.e., the root elements).
4327
4329
4330 // sorted list of all vertex nodes touching 'elem'
4331 Array<int> vert;
4332 vert.Reserve(128);
4333
4334 // support for non-leaf 'elem', collect vertices of all children
4335 Array<int> stack;
4336 stack.Reserve(64);
4337 stack.Append(elem);
4338
4339 while (stack.Size())
4340 {
4341 Element &el = elements[stack.Last()];
4342 stack.DeleteLast();
4343
4344 if (!el.ref_type)
4345 {
4346 int *v = element_vertex.GetRow(el.index);
4347 int nv = element_vertex.RowSize(el.index);
4348 for (int i = 0; i < nv; i++)
4349 {
4350 vert.Append(v[i]);
4351 }
4352
4353 nv = GI[el.Geom()].nv;
4354 for (int i = 0; i < nv; i++)
4355 {
4356 vert.Append(el.node[i]);
4357 }
4358 }
4359 else
4360 {
4361 for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
4362 {
4363 stack.Append(el.child[i]);
4364 }
4365 }
4366 }
4367
4368 vert.Sort();
4369 vert.Unique();
4370
4371 int *v1 = vert.GetData();
4372 int nv1 = vert.Size();
4373
4374 if (!search_set) { search_set = &leaf_elements; }
4375
4376 // test *all* potential neighbors from the search set
4377 for (int i = 0; i < search_set->Size(); i++)
4378 {
4379 int testme = (*search_set)[i];
4380 if (testme != elem)
4381 {
4382 Element &el = elements[testme];
4383 int *v2 = element_vertex.GetRow(el.index);
4384 int nv2 = element_vertex.RowSize(el.index);
4385
4386 bool hit = sorted_lists_intersect(v1, v2, nv1, nv2);
4387
4388 if (!hit)
4389 {
4390 int nv = GI[el.Geom()].nv;
4391 for (int j = 0; j < nv; j++)
4392 {
4393 hit = sorted_lists_intersect(&el.node[j], v1, 1, nv1);
4394 if (hit) { break; }
4395 }
4396 }
4397
4398 if (hit) { neighbors.Append(testme); }
4399 }
4400 }
4401}
4402
4404 Array<int> &expanded,
4405 const Array<int> *search_set)
4406{
4408
4409 Array<char> vmark(nodes.NumIds());
4410 vmark = 0;
4411
4412 for (int i = 0; i < elems.Size(); i++)
4413 {
4414 Element &el = elements[elems[i]];
4415
4416 int *v = element_vertex.GetRow(el.index);
4417 int nv = element_vertex.RowSize(el.index);
4418 for (int j = 0; j < nv; j++)
4419 {
4420 vmark[v[j]] = 1;
4421 }
4422
4423 nv = GI[el.Geom()].nv;
4424 for (int j = 0; j < nv; j++)
4425 {
4426 vmark[el.node[j]] = 1;
4427 }
4428 }
4429
4430 if (!search_set)
4431 {
4432 search_set = &leaf_elements;
4433 }
4434
4435 expanded.SetSize(0);
4436 for (int i = 0; i < search_set->Size(); i++)
4437 {
4438 int testme = (*search_set)[i];
4439 Element &el = elements[testme];
4440 bool hit = false;
4441
4442 int *v = element_vertex.GetRow(el.index);
4443 int nv = element_vertex.RowSize(el.index);
4444 for (int j = 0; j < nv; j++)
4445 {
4446 if (vmark[v[j]]) { hit = true; break; }
4447 }
4448
4449 if (!hit)
4450 {
4451 nv = GI[el.Geom()].nv;
4452 for (int j = 0; j < nv; j++)
4453 {
4454 if (vmark[el.node[j]]) { hit = true; break; }
4455 }
4456 }
4457
4458 if (hit) { expanded.Append(testme); }
4459 }
4460}
4461
4462int NCMesh::GetVertexRootCoord(int elem, RefCoord coord[3]) const
4463{
4464 while (1)
4465 {
4466 const Element &el = elements[elem];
4467 if (el.parent < 0) { return elem; }
4468
4469 const Element &pa = elements[el.parent];
4470 MFEM_ASSERT(pa.ref_type, "internal error");
4471
4472 int ch = 0;
4473 while (ch < MaxElemChildren && pa.child[ch] != elem) { ch++; }
4474 MFEM_ASSERT(ch < MaxElemChildren, "internal error");
4475
4476 MFEM_ASSERT(geom_parent[el.Geom()], "unsupported geometry");
4477 const RefTrf &tr = geom_parent[el.Geom()][(int) pa.ref_type][ch];
4478 tr.Apply(coord, coord);
4479
4480 elem = el.parent;
4481 }
4482}
4483
4484static bool RefPointInside(Geometry::Type geom, const RefCoord pt[3])
4485{
4486 switch (geom)
4487 {
4488 case Geometry::SQUARE:
4489 return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
4490 (pt[1] >= 0) && (pt[1] <= T_ONE);
4491
4492 case Geometry::CUBE:
4493 return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
4494 (pt[1] >= 0) && (pt[1] <= T_ONE) &&
4495 (pt[2] >= 0) && (pt[2] <= T_ONE);
4496
4497 case Geometry::TRIANGLE:
4498 return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE);
4499
4500 case Geometry::PRISM:
4501 return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE) &&
4502 (pt[2] >= 0) && (pt[2] <= T_ONE);
4503
4504 case Geometry::PYRAMID:
4505 return (pt[0] >= 0) && (pt[1] >= 0) && (pt[2] >= 0.0) &&
4506 (pt[0] + pt[2] <= T_ONE) && (pt[1] + pt[2] <= T_ONE) &&
4507 (pt[2] <= T_ONE);
4508
4509 default:
4510 MFEM_ABORT("unsupported geometry");
4511 return false;
4512 }
4513}
4514
4515void NCMesh::CollectIncidentElements(int elem, const RefCoord coord[3],
4516 Array<int> &list) const
4517{
4518 const Element &el = elements[elem];
4519 if (!el.ref_type)
4520 {
4521 list.Append(elem);
4522 return;
4523 }
4524
4525 RefCoord tcoord[3];
4526 for (int ch = 0; ch < MaxElemChildren && el.child[ch] >= 0; ch++)
4527 {
4528 const RefTrf &tr = geom_child[el.Geom()][(int) el.ref_type][ch];
4529 tr.Apply(coord, tcoord);
4530
4531 if (RefPointInside(el.Geom(), tcoord))
4532 {
4533 CollectIncidentElements(el.child[ch], tcoord, list);
4534 }
4535 }
4536}
4537
4538void NCMesh::FindVertexCousins(int elem, int local, Array<int> &cousins) const
4539{
4540 const Element &el = elements[elem];
4541
4542 RefCoord coord[3];
4543 MFEM_ASSERT(geom_corners[el.Geom()], "unsupported geometry");
4544 std::memcpy(coord, geom_corners[el.Geom()][local], sizeof(coord));
4545
4546 int root = GetVertexRootCoord(elem, coord);
4547
4548 cousins.SetSize(0);
4549 CollectIncidentElements(root, coord, cousins);
4550}
4551
4552
4553//// Coarse/fine transformations ///////////////////////////////////////////////
4554
4556{
4557 MFEM_ASSERT(np == pm.np, "");
4558 for (int i = 0; i < np; i++)
4559 {
4560 MFEM_ASSERT(points[i].dim == pm.points[i].dim, "");
4561 for (int j = 0; j < points[i].dim; j++)
4562 {
4563 if (points[i].coord[j] != pm.points[i].coord[j]) { return false; }
4564 }
4565 }
4566 return true;
4567}
4568
4570{
4571 point_matrix.SetSize(points[0].dim, np);
4572 for (int i = 0; i < np; i++)
4573 {
4574 for (int j = 0; j < points[0].dim; j++)
4575 {
4576 point_matrix(j, i) = points[i].coord[j];
4577 }
4578 }
4579}
4580
4582 Point(0), Point(1)
4583);
4585 Point(0, 0), Point(1, 0), Point(0, 1)
4586);
4588 Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)
4589);
4591 Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)
4592);
4594 Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0),
4595 Point(0, 0, 1), Point(1, 0, 1), Point(0, 1, 1)
4596);
4598 Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0),
4599 Point(0, 1, 0), Point(0, 0, 1)
4600);
4602 Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
4603 Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
4604);
4605
4607{
4608 switch (geom)
4609 {
4616 case Geometry::CUBE: return pm_hex_identity;
4617 default:
4618 MFEM_ABORT("unsupported geometry " << geom);
4619 return pm_tri_identity;
4620 }
4621}
4622
4623void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
4624 DenseMatrix& matrix) const
4625{
4626 PointMatrix pm = GetGeomIdentity(geom);
4627
4628 while (*ref_path)
4629 {
4630 int ref_type = *ref_path++;
4631 int child = *ref_path++;
4632
4633 // TODO: do this with the new child transform tables
4634
4635 if (geom == Geometry::CUBE)
4636 {
4637 if (ref_type == 1) // split along X axis
4638 {
4639 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4640 Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
4641
4642 if (child == 0)
4643 {
4644 pm = PointMatrix(pm(0), mid01, mid23, pm(3),
4645 pm(4), mid45, mid67, pm(7));
4646 }
4647 else if (child == 1)
4648 {
4649 pm = PointMatrix(mid01, pm(1), pm(2), mid23,
4650 mid45, pm(5), pm(6), mid67);
4651 }
4652 }
4653 else if (ref_type == 2) // split along Y axis
4654 {
4655 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4656 Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4657
4658 if (child == 0)
4659 {
4660 pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4661 pm(4), pm(5), mid56, mid74);
4662 }
4663 else if (child == 1)
4664 {
4665 pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4666 mid74, mid56, pm(6), pm(7));
4667 }
4668 }
4669 else if (ref_type == 4) // split along Z axis
4670 {
4671 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4672 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4673
4674 if (child == 0)
4675 {
4676 pm = PointMatrix(pm(0), pm(1), pm(2), pm(3),
4677 mid04, mid15, mid26, mid37);
4678 }
4679 else if (child == 1)
4680 {
4681 pm = PointMatrix(mid04, mid15, mid26, mid37,
4682 pm(4), pm(5), pm(6), pm(7));
4683 }
4684 }
4685 else if (ref_type == 3) // XY split
4686 {
4687 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4688 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4689 Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4690 Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4691
4692 Point midf0(mid23, mid12, mid01, mid30);
4693 Point midf5(mid45, mid56, mid67, mid74);
4694
4695 if (child == 0)
4696 {
4697 pm = PointMatrix(pm(0), mid01, midf0, mid30,
4698 pm(4), mid45, midf5, mid74);
4699 }
4700 else if (child == 1)
4701 {
4702 pm = PointMatrix(mid01, pm(1), mid12, midf0,
4703 mid45, pm(5), mid56, midf5);
4704 }
4705 else if (child == 2)
4706 {
4707 pm = PointMatrix(midf0, mid12, pm(2), mid23,
4708 midf5, mid56, pm(6), mid67);
4709 }
4710 else if (child == 3)
4711 {
4712 pm = PointMatrix(mid30, midf0, mid23, pm(3),
4713 mid74, midf5, mid67, pm(7));
4714 }
4715 }
4716 else if (ref_type == 5) // XZ split
4717 {
4718 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4719 Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
4720 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4721 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4722
4723 Point midf1(mid01, mid15, mid45, mid04);
4724 Point midf3(mid23, mid37, mid67, mid26);
4725
4726 if (child == 0)
4727 {
4728 pm = PointMatrix(pm(0), mid01, mid23, pm(3),
4729 mid04, midf1, midf3, mid37);
4730 }
4731 else if (child == 1)
4732 {
4733 pm = PointMatrix(mid01, pm(1), pm(2), mid23,
4734 midf1, mid15, mid26, midf3);
4735 }
4736 else if (child == 2)
4737 {
4738 pm = PointMatrix(midf1, mid15, mid26, midf3,
4739 mid45, pm(5), pm(6), mid67);
4740 }
4741 else if (child == 3)
4742 {
4743 pm = PointMatrix(mid04, midf1, midf3, mid37,
4744 pm(4), mid45, mid67, pm(7));
4745 }
4746 }
4747 else if (ref_type == 6) // YZ split
4748 {
4749 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4750 Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4751 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4752 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4753
4754 Point midf2(mid12, mid26, mid56, mid15);
4755 Point midf4(mid30, mid04, mid74, mid37);
4756
4757 if (child == 0)
4758 {
4759 pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4760 mid04, mid15, midf2, midf4);
4761 }
4762 else if (child == 1)
4763 {
4764 pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4765 midf4, midf2, mid26, mid37);
4766 }
4767 else if (child == 2)
4768 {
4769 pm = PointMatrix(mid04, mid15, midf2, midf4,
4770 pm(4), pm(5), mid56, mid74);
4771 }
4772 else if (child == 3)
4773 {
4774 pm = PointMatrix(midf4, midf2, mid26, mid37,
4775 mid74, mid56, pm(6), pm(7));
4776 }
4777 }
4778 else if (ref_type == 7) // full isotropic refinement
4779 {
4780 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4781 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4782 Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4783 Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4784 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4785 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4786
4787 Point midf0(mid23, mid12, mid01, mid30);
4788 Point midf1(mid01, mid15, mid45, mid04);
4789 Point midf2(mid12, mid26, mid56, mid15);
4790 Point midf3(mid23, mid37, mid67, mid26);
4791 Point midf4(mid30, mid04, mid74, mid37);
4792 Point midf5(mid45, mid56, mid67, mid74);
4793
4794 Point midel(midf1, midf3);
4795
4796 if (child == 0)
4797 {
4798 pm = PointMatrix(pm(0), mid01, midf0, mid30,
4799 mid04, midf1, midel, midf4);
4800 }
4801 else if (child == 1)
4802 {
4803 pm = PointMatrix(mid01, pm(1), mid12, midf0,
4804 midf1, mid15, midf2, midel);
4805 }
4806 else if (child == 2)
4807 {
4808 pm = PointMatrix(midf0, mid12, pm(2), mid23,
4809 midel, midf2, mid26, midf3);
4810 }
4811 else if (child == 3)
4812 {
4813 pm = PointMatrix(mid30, midf0, mid23, pm(3),
4814 midf4, midel, midf3, mid37);
4815 }
4816 else if (child == 4)
4817 {
4818 pm = PointMatrix(mid04, midf1, midel, midf4,
4819 pm(4), mid45, midf5, mid74);
4820 }
4821 else if (child == 5)
4822 {
4823 pm = PointMatrix(midf1, mid15, midf2, midel,
4824 mid45, pm(5), mid56, midf5);
4825 }
4826 else if (child == 6)
4827 {
4828 pm = PointMatrix(midel, midf2, mid26, midf3,
4829 midf5, mid56, pm(6), mid67);
4830 }
4831 else if (child == 7)
4832 {
4833 pm = PointMatrix(midf4, midel, midf3, mid37,
4834 mid74, midf5, mid67, pm(7));
4835 }
4836 }
4837 }
4838 else if (geom == Geometry::PRISM)
4839 {
4840 if (ref_type < 4) // XY split
4841 {
4842 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4843 Point mid20(pm(2), pm(0)), mid34(pm(3), pm(4));
4844 Point mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4845
4846 if (child == 0)
4847 {
4848 pm = PointMatrix(pm(0), mid01, mid20, pm(3), mid34, mid53);
4849 }
4850 else if (child == 1)
4851 {
4852 pm = PointMatrix(mid01, pm(1), mid12, mid34, pm(4), mid45);
4853 }
4854 else if (child == 2)
4855 {
4856 pm = PointMatrix(mid20, mid12, pm(2), mid53, mid45, pm(5));
4857 }
4858 else if (child == 3)
4859 {
4860 pm = PointMatrix(mid12, mid20, mid01, mid45, mid53, mid34);
4861 }
4862 }
4863 else if (ref_type == 4) // Z split
4864 {
4865 Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4866
4867 if (child == 0)
4868 {
4869 pm = PointMatrix(pm(0), pm(1), pm(2), mid03, mid14, mid25);
4870 }
4871 else if (child == 1)
4872 {
4873 pm = PointMatrix(mid03, mid14, mid25, pm(3), pm(4), pm(5));
4874 }
4875 }
4876 else // ref_type > 4, iso split
4877 {
4878 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4879 Point mid34(pm(3), pm(4)), mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4880 Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4881
4882 Point midf2(mid01, mid14, mid34, mid03);
4883 Point midf3(mid12, mid25, mid45, mid14);
4884 Point midf4(mid20, mid03, mid53, mid25);
4885
4886 if (child == 0)
4887 {
4888 pm = PointMatrix(pm(0), mid01, mid20, mid03, midf2, midf4);
4889 }
4890 else if (child == 1)
4891 {
4892 pm = PointMatrix(mid01, pm(1), mid12, midf2, mid14, midf3);
4893 }
4894 else if (child == 2)
4895 {
4896 pm = PointMatrix(mid20, mid12, pm(2), midf4, midf3, mid25);
4897 }
4898 else if (child == 3)
4899 {
4900 pm = PointMatrix(mid12, mid20, mid01, midf3, midf4, midf2);
4901 }
4902 else if (child == 4)
4903 {
4904 pm = PointMatrix(mid03, midf2, midf4, pm(3), mid34, mid53);
4905 }
4906 else if (child == 5)
4907 {
4908 pm = PointMatrix(midf2, mid14, midf3, mid34, pm(4), mid45);
4909 }
4910 else if (child == 6)
4911 {
4912 pm = PointMatrix(midf4, midf3, mid25, mid53, mid45, pm(5));
4913 }
4914 else if (child == 7)
4915 {
4916 pm = PointMatrix(midf3, midf4, midf2, mid45, mid53, mid34);
4917 }
4918 }
4919 }
4920 else if (geom == Geometry::PYRAMID)
4921 {
4922 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4923 Point mid03(pm(0), pm(3)), mid12(pm(1), pm(2));
4924 Point mid04(pm(0), pm(4)), mid14(pm(1), pm(4));
4925 Point mid24(pm(2), pm(4)), mid34(pm(3), pm(4));
4926 Point midf0(mid23, mid12, mid01, mid03);
4927
4928 if (child == 0) // Pyramid
4929 {
4930 pm = PointMatrix(pm(0), mid01, midf0, mid03, mid04);
4931 }
4932 else if (child == 1) // Pyramid
4933 {
4934 pm = PointMatrix(mid01, pm(1), mid12, midf0, mid14);
4935 }
4936 else if (child == 2) // Pyramid
4937 {
4938 pm = PointMatrix(midf0, mid12, pm(2), mid23, mid24);
4939 }
4940 else if (child == 3) // Pyramid
4941 {
4942 pm = PointMatrix(mid03, midf0, mid23, pm(3), mid34);
4943 }
4944 else if (child == 4) // Pyramid
4945 {
4946 pm = PointMatrix(mid24, mid14, mid04, mid34, midf0);
4947 }
4948 else if (child == 5) // Pyramid
4949 {
4950 pm = PointMatrix(mid04, mid14, mid24, mid34, pm(4));
4951 }
4952 else if (child == 6) // Tet
4953 {
4954 pm = PointMatrix(mid01, midf0, mid04, mid14);
4955 }
4956 else if (child == 7) // Tet
4957 {
4958 pm = PointMatrix(midf0, mid14, mid12, mid24);
4959 }
4960 else if (child == 8) // Tet
4961 {
4962 pm = PointMatrix(midf0, mid23, mid34, mid24);
4963 }
4964 else if (child == 9) // Tet
4965 {
4966 pm = PointMatrix(mid03, mid04, midf0, mid34);
4967 }
4968 }
4969 else if (geom == Geometry::TETRAHEDRON)
4970 {
4971 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid02(pm(2), pm(0));
4972 Point mid03(pm(0), pm(3)), mid13(pm(1), pm(3)), mid23(pm(2), pm(3));
4973
4974 if (child == 0)
4975 {
4976 pm = PointMatrix(pm(0), mid01, mid02, mid03);
4977 }
4978 else if (child == 1)
4979 {
4980 pm = PointMatrix(mid01, pm(1), mid12, mid13);
4981 }
4982 else if (child == 2)
4983 {
4984 pm = PointMatrix(mid02, mid12, pm(2), mid23);
4985 }
4986 else if (child == 3)
4987 {
4988 pm = PointMatrix(mid03, mid13, mid23, pm(3));
4989 }
4990 else if (child == 4)
4991 {
4992 pm = PointMatrix(mid01, mid23, mid02, mid03);
4993 }
4994 else if (child == 5)
4995 {
4996 pm = PointMatrix(mid01, mid23, mid03, mid13);
4997 }
4998 else if (child == 6)
4999 {
5000 pm = PointMatrix(mid01, mid23, mid13, mid12);
5001 }
5002 else if (child == 7)
5003 {
5004 pm = PointMatrix(mid01, mid23, mid12, mid02);
5005 }
5006 }
5007 else if (geom == Geometry::SQUARE)
5008 {
5009 if (ref_type == 1) // X split
5010 {
5011 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
5012
5013 if (child == 0)
5014 {
5015 pm = PointMatrix(pm(0), mid01, mid23, pm(3));
5016 }
5017 else if (child == 1)
5018 {
5019 pm = PointMatrix(mid01, pm(1), pm(2), mid23);
5020 }
5021 }
5022 else if (ref_type == 2) // Y split
5023 {
5024 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
5025
5026 if (child == 0)
5027 {
5028 pm = PointMatrix(pm(0), pm(1), mid12, mid30);
5029 }
5030 else if (child == 1)
5031 {
5032 pm = PointMatrix(mid30, mid12, pm(2), pm(3));
5033 }
5034 }
5035 else if (ref_type == 3) // iso split
5036 {
5037 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
5038 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
5039 Point midel(mid01, mid23);
5040
5041 if (child == 0)
5042 {
5043 pm = PointMatrix(pm(0), mid01, midel, mid30);
5044 }
5045 else if (child == 1)
5046 {
5047 pm = PointMatrix(mid01, pm(1), mid12, midel);
5048 }
5049 else if (child == 2)
5050 {
5051 pm = PointMatrix(midel, mid12, pm(2), mid23);
5052 }
5053 else if (child == 3)
5054 {
5055 pm = PointMatrix(mid30, midel, mid23, pm(3));
5056 }
5057 }
5058 }
5059 else if (geom == Geometry::TRIANGLE)
5060 {
5061 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
5062
5063 if (child == 0)
5064 {
5065 pm = PointMatrix(pm(0), mid01, mid20);
5066 }
5067 else if (child == 1)
5068 {
5069 pm = PointMatrix(mid01, pm(1), mid12);
5070 }
5071 else if (child == 2)
5072 {
5073 pm = PointMatrix(mid20, mid12, pm(2));
5074 }
5075 else if (child == 3)
5076 {
5077 pm = PointMatrix(mid12, mid20, mid01);
5078 }
5079 }
5080 else if (geom == Geometry::SEGMENT)
5081 {
5082 Point mid01(pm(0), pm(1));
5083
5084 if (child == 0)
5085 {
5086 pm = PointMatrix(pm(0), mid01);
5087 }
5088 else if (child == 1)
5089 {
5090 pm = PointMatrix(mid01, pm(1));
5091 }
5092 }
5093 }
5094
5095 // write the points to the matrix
5096 for (int i = 0; i < pm.np; i++)
5097 {
5098 for (int j = 0; j < pm(i).dim; j++)
5099 {
5100 matrix(j, i) = pm(i).coord[j];
5101 }
5102 }
5103}
5104
5106{
5109
5110 for (int i = 0; i < leaf_elements.Size(); i++)
5111 {
5112 int elem = leaf_elements[i];
5113 if (!IsGhost(elements[elem])) { coarse_elements.Append(elem); }
5114 }
5115
5116 transforms.embeddings.DeleteAll();
5117}
5118
5119void NCMesh::TraverseRefinements(int elem, int coarse_index,
5120 std::string &ref_path, RefPathMap &map) const
5121{
5122 const Element &el = elements[elem];
5123 if (!el.ref_type)
5124 {
5125 int &matrix = map[ref_path];
5126 if (!matrix) { matrix = static_cast<int>(map.size()); }
5127
5129 emb.parent = coarse_index;
5130 emb.matrix = matrix - 1;
5131 emb.geom = el.Geom();
5132 emb.ghost = IsGhost(el);
5133 }
5134 else
5135 {
5136 MFEM_ASSERT(el.tet_type == 0, "not implemented");
5137
5138 ref_path.push_back(el.ref_type);
5139 ref_path.push_back(0);
5140
5141 for (int i = 0; i < MaxElemChildren; i++)
5142 {
5143 if (el.child[i] >= 0)
5144 {
5145 ref_path[ref_path.length()-1] = i;
5146 TraverseRefinements(el.child[i], coarse_index, ref_path, map);
5147 }
5148 }
5149 ref_path.resize(ref_path.length()-2);
5150 }
5151}
5152
5154{
5155 MFEM_VERIFY(coarse_elements.Size() || !leaf_elements.Size(),
5156 "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
5157 " and Refine().");
5158
5159 if (!transforms.embeddings.Size())
5160 {
5161 transforms.Clear();
5163
5164 std::string ref_path;
5165 ref_path.reserve(100);
5166
5167 RefPathMap path_map[Geometry::NumGeom];
5168 for (int g = 0; g < Geometry::NumGeom; g++)
5169 {
5170 path_map[g][ref_path] = 1; // empty path == identity
5171 }
5172
5173 int used_geoms = 0;
5174 for (int i = 0; i < coarse_elements.Size(); i++)
5175 {
5176 int geom = elements[coarse_elements[i]].geom;
5177 TraverseRefinements(coarse_elements[i], i, ref_path, path_map[geom]);
5178 used_geoms |= (1 << geom);
5179 }
5180
5181 for (int g = 0; g < Geometry::NumGeom; g++)
5182 {
5183 if (used_geoms & (1 << g))
5184 {
5186 const PointMatrix &identity = GetGeomIdentity(geom);
5187
5189 .SetSize(Dim, identity.np, static_cast<int>(path_map[g].size()));
5190
5191 // calculate the point matrices
5192 RefPathMap::iterator it;
5193 for (it = path_map[g].begin(); it != path_map[g].end(); ++it)
5194 {
5195 GetPointMatrix(geom, it->first.c_str(),
5196 transforms.point_matrices[g](it->second-1));
5197 }
5198 }
5199 }
5200 }
5201 return transforms;
5202}
5203
5205{
5206 MFEM_VERIFY(transforms.embeddings.Size() || !leaf_elements.Size(),
5207 "GetDerefinementTransforms() must be preceded by Derefine().");
5208
5210 {
5211 std::map<int, int> mat_no[Geometry::NumGeom];
5212 for (int g = 0; g < Geometry::NumGeom; g++)
5213 {
5214 mat_no[g][0] = 1; // 0 == identity
5215 }
5216
5217 // assign numbers to the different matrices used
5218 for (int i = 0; i < transforms.embeddings.Size(); i++)
5219 {
5221 int code = emb.matrix; // see SetDerefMatrixCodes()
5222 if (code)
5223 {
5224 int &matrix = mat_no[emb.geom][code];
5225 if (!matrix) { matrix = static_cast<int>(mat_no[emb.geom].size()); }
5226
5227 emb.matrix = matrix - 1;
5228 }
5229 }
5230
5231 for (int g = 0; g < Geometry::NumGeom; g++)
5232 {
5233 if (Geoms & (1 << g))
5234 {
5236 const PointMatrix &identity = GetGeomIdentity(geom);
5237
5239 .SetSize(Dim, identity.np, static_cast<int>(mat_no[geom].size()));
5240
5241 // calculate point matrices
5242 for (auto it = mat_no[geom].begin(); it != mat_no[geom].end(); ++it)
5243 {
5244 char path[3] = { 0, 0, 0 };
5245
5246 int code = it->first;
5247 if (code)
5248 {
5249 path[0] = code >> 4; // ref_type (see SetDerefMatrixCodes())
5250 path[1] = code & 0xf; // child
5251 }
5252
5253 GetPointMatrix(geom, path,
5254 transforms.point_matrices[geom](it->second-1));
5255 }
5256 }
5257 }
5258 }
5259 return transforms;
5260}
5261
5263 bool want_ghosts) const
5264{
5265 Array<Connection> conn;
5266 conn.Reserve(embeddings.Size());
5267
5268 int max_parent = -1;
5269 for (int i = 0; i < embeddings.Size(); i++)
5270 {
5271 const Embedding &emb = embeddings[i];
5272 if ((emb.parent >= 0) &&
5273 (!emb.ghost || want_ghosts))
5274 {
5275 conn.Append(Connection(emb.parent, i));
5276 max_parent = std::max(emb.parent, max_parent);
5277 }
5278 }
5279
5280 conn.Sort(); // NOTE: unique is not necessary
5281 coarse_to_fine.MakeFromList(max_parent+1, conn);
5282}
5283
5289
5291{
5292 for (int i = 0; i < Geometry::NumGeom; i++)
5293 {
5294 point_matrices[i].SetSize(0, 0, 0);
5295 }
5296 embeddings.DeleteAll();
5297}
5298
5300{
5301 // return true if point matrices are present for any geometry
5302 for (int i = 0; i < Geometry::NumGeom; i++)
5303 {
5304 if (point_matrices[i].SizeK()) { return true; }
5305 }
5306 return false;
5307}
5308
5310{
5311 for (int g = 0; g < Geometry::NumGeom; ++g)
5312 {
5313 a.point_matrices[g].Swap(b.point_matrices[g]);
5314 }
5315 Swap(a.embeddings, b.embeddings);
5316}
5317
5318
5319//// SFC Ordering //////////////////////////////////////////////////////////////
5320
5321static int sgn(int x)
5322{
5323 return (x < 0) ? -1 : (x > 0) ? 1 : 0;
5324}
5325
5326static void HilbertSfc2D(int x, int y, int ax, int ay, int bx, int by,
5327 Array<int> &coords)
5328{
5329 int w = std::abs(ax + ay);
5330 int h = std::abs(bx + by);
5331
5332 int dax = sgn(ax), day = sgn(ay); // unit major direction ("right")
5333 int dbx = sgn(bx), dby = sgn(by); // unit orthogonal direction ("up")
5334
5335 if (h == 1) // trivial row fill
5336 {
5337 for (int i = 0; i < w; i++, x += dax, y += day)
5338 {
5339 coords.Append(x);
5340 coords.Append(y);
5341 }
5342 return;
5343 }
5344 if (w == 1) // trivial column fill
5345 {
5346 for (int i = 0; i < h; i++, x += dbx, y += dby)
5347 {
5348 coords.Append(x);
5349 coords.Append(y);
5350 }
5351 return;
5352 }
5353
5354 int ax2 = ax/2, ay2 = ay/2;
5355 int bx2 = bx/2, by2 = by/2;
5356
5357 int w2 = std::abs(ax2 + ay2);
5358 int h2 = std::abs(bx2 + by2);
5359
5360 if (2*w > 3*h) // long case: split in two parts only
5361 {
5362 if ((w2 & 0x1) && (w > 2))
5363 {
5364 ax2 += dax, ay2 += day; // prefer even steps
5365 }
5366
5367 HilbertSfc2D(x, y, ax2, ay2, bx, by, coords);
5368 HilbertSfc2D(x+ax2, y+ay2, ax-ax2, ay-ay2, bx, by, coords);
5369 }
5370 else // standard case: one step up, one long horizontal step, one step down
5371 {
5372 if ((h2 & 0x1) && (h > 2))
5373 {
5374 bx2 += dbx, by2 += dby; // prefer even steps
5375 }
5376
5377 HilbertSfc2D(x, y, bx2, by2, ax2, ay2, coords);
5378 HilbertSfc2D(x+bx2, y+by2, ax, ay, bx-bx2, by-by2, coords);
5379 HilbertSfc2D(x+(ax-dax)+(bx2-dbx), y+(ay-day)+(by2-dby),
5380 -bx2, -by2, -(ax-ax2), -(ay-ay2), coords);
5381 }
5382}
5383
5384static void HilbertSfc3D(int x, int y, int z,
5385 int ax, int ay, int az,
5386 int bx, int by, int bz,
5387 int cx, int cy, int cz,
5388 Array<int> &coords)
5389{
5390 int w = std::abs(ax + ay + az);
5391 int h = std::abs(bx + by + bz);
5392 int d = std::abs(cx + cy + cz);
5393
5394 int dax = sgn(ax), day = sgn(ay), daz = sgn(az); // unit major dir ("right")
5395 int dbx = sgn(bx), dby = sgn(by), dbz = sgn(bz); // unit ortho dir ("forward")
5396 int dcx = sgn(cx), dcy = sgn(cy), dcz = sgn(cz); // unit ortho dir ("up")
5397
5398 // trivial row/column fills
5399 if (h == 1 && d == 1)
5400 {
5401 for (int i = 0; i < w; i++, x += dax, y += day, z += daz)
5402 {
5403 coords.Append(x);
5404 coords.Append(y);
5405 coords.Append(z);
5406 }
5407 return;
5408 }
5409 if (w == 1 && d == 1)
5410 {
5411 for (int i = 0; i < h; i++, x += dbx, y += dby, z += dbz)
5412 {
5413 coords.Append(x);
5414 coords.Append(y);
5415 coords.Append(z);
5416 }
5417 return;
5418 }
5419 if (w == 1 && h == 1)
5420 {
5421 for (int i = 0; i < d; i++, x += dcx, y += dcy, z += dcz)
5422 {
5423 coords.Append(x);
5424 coords.Append(y);
5425 coords.Append(z);
5426 }
5427 return;
5428 }
5429
5430 int ax2 = ax/2, ay2 = ay/2, az2 = az/2;
5431 int bx2 = bx/2, by2 = by/2, bz2 = bz/2;
5432 int cx2 = cx/2, cy2 = cy/2, cz2 = cz/2;
5433
5434 int w2 = std::abs(ax2 + ay2 + az2);
5435 int h2 = std::abs(bx2 + by2 + bz2);
5436 int d2 = std::abs(cx2 + cy2 + cz2);
5437
5438 // prefer even steps
5439 if ((w2 & 0x1) && (w > 2))
5440 {
5441 ax2 += dax, ay2 += day, az2 += daz;
5442 }
5443 if ((h2 & 0x1) && (h > 2))
5444 {
5445 bx2 += dbx, by2 += dby, bz2 += dbz;
5446 }
5447 if ((d2 & 0x1) && (d > 2))
5448 {
5449 cx2 += dcx, cy2 += dcy, cz2 += dcz;
5450 }
5451
5452 // wide case, split in w only
5453 if ((2*w > 3*h) && (2*w > 3*d))
5454 {
5455 HilbertSfc3D(x, y, z,
5456 ax2, ay2, az2,
5457 bx, by, bz,
5458 cx, cy, cz, coords);
5459
5460 HilbertSfc3D(x+ax2, y+ay2, z+az2,
5461 ax-ax2, ay-ay2, az-az2,
5462 bx, by, bz,
5463 cx, cy, cz, coords);
5464 }
5465 // do not split in d
5466 else if (3*h > 4*d)
5467 {
5468 HilbertSfc3D(x, y, z,
5469 bx2, by2, bz2,
5470 cx, cy, cz,
5471 ax2, ay2, az2, coords);
5472
5473 HilbertSfc3D(x+bx2, y+by2, z+bz2,
5474 ax, ay, az,
5475 bx-bx2, by-by2, bz-bz2,
5476 cx, cy, cz, coords);
5477
5478 HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
5479 y+(ay-day)+(by2-dby),
5480 z+(az-daz)+(bz2-dbz),
5481 -bx2, -by2, -bz2,
5482 cx, cy, cz,
5483 -(ax-ax2), -(ay-ay2), -(az-az2), coords);
5484 }
5485 // do not split in h
5486 else if (3*d > 4*h)
5487 {
5488 HilbertSfc3D(x, y, z,
5489 cx2, cy2, cz2,
5490 ax2, ay2, az2,
5491 bx, by, bz, coords);
5492
5493 HilbertSfc3D(x+cx2, y+cy2, z+cz2,
5494 ax, ay, az,
5495 bx, by, bz,
5496 cx-cx2, cy-cy2, cz-cz2, coords);
5497
5498 HilbertSfc3D(x+(ax-dax)+(cx2-dcx),
5499 y+(ay-day)+(cy2-dcy),
5500 z+(az-daz)+(cz2-dcz),
5501 -cx2, -cy2, -cz2,
5502 -(ax-ax2), -(ay-ay2), -(az-az2),
5503 bx, by, bz, coords);
5504 }
5505 // regular case, split in all w/h/d
5506 else
5507 {
5508 HilbertSfc3D(x, y, z,
5509 bx2, by2, bz2,
5510 cx2, cy2, cz2,
5511 ax2, ay2, az2, coords);
5512
5513 HilbertSfc3D(x+bx2, y+by2, z+bz2,
5514 cx, cy, cz,
5515 ax2, ay2, az2,
5516 bx-bx2, by-by2, bz-bz2, coords);
5517
5518 HilbertSfc3D(x+(bx2-dbx)+(cx-dcx),
5519 y+(by2-dby)+(cy-dcy),
5520 z+(bz2-dbz)+(cz-dcz),
5521 ax, ay, az,
5522 -bx2, -by2, -bz2,
5523 -(cx-cx2), -(cy-cy2), -(cz-cz2), coords);
5524
5525 HilbertSfc3D(x+(ax-dax)+bx2+(cx-dcx),
5526 y+(ay-day)+by2+(cy-dcy),
5527 z+(az-daz)+bz2+(cz-dcz),
5528 -cx, -cy, -cz,
5529 -(ax-ax2), -(ay-ay2), -(az-az2),
5530 bx-bx2, by-by2, bz-bz2, coords);
5531
5532 HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
5533 y+(ay-day)+(by2-dby),
5534 z+(az-daz)+(bz2-dbz),
5535 -bx2, -by2, -bz2,
5536 cx2, cy2, cz2,
5537 -(ax-ax2), -(ay-ay2), -(az-az2), coords);
5538 }
5539}
5540
5541void NCMesh::GridSfcOrdering2D(int width, int height, Array<int> &coords)
5542{
5543 coords.SetSize(0);
5544 coords.Reserve(2*width*height);
5545
5546 if (width >= height)
5547 {
5548 HilbertSfc2D(0, 0, width, 0, 0, height, coords);
5549 }
5550 else
5551 {
5552 HilbertSfc2D(0, 0, 0, height, width, 0, coords);
5553 }
5554}
5555
5556void NCMesh::GridSfcOrdering3D(int width, int height, int depth,
5557 Array<int> &coords)
5558{
5559 coords.SetSize(0);
5560 coords.Reserve(3*width*height*depth);
5561
5562 if (width >= height && width >= depth)
5563 {
5564 HilbertSfc3D(0, 0, 0,
5565 width, 0, 0,
5566 0, height, 0,
5567 0, 0, depth, coords);
5568 }
5569 else if (height >= width && height >= depth)
5570 {
5571 HilbertSfc3D(0, 0, 0,
5572 0, height, 0,
5573 width, 0, 0,
5574 0, 0, depth, coords);
5575 }
5576 else // depth >= width && depth >= height
5577 {
5578 HilbertSfc3D(0, 0, 0,
5579 0, 0, depth,
5580 width, 0, 0,
5581 0, height, 0, coords);
5582 }
5583}
5584
5585
5586//// Utility ///////////////////////////////////////////////////////////////////
5587
5588void NCMesh::GetEdgeVertices(const MeshId &edge_id, int vert_index[2],
5589 bool oriented) const
5590{
5591 const Element &el = elements[edge_id.element];
5592 const GeomInfo& gi = GI[el.Geom()];
5593 const int* ev = gi.edges[(int) edge_id.local];
5594
5595 int n0 = el.node[ev[0]], n1 = el.node[ev[1]];
5596 if (n0 > n1) { std::swap(n0, n1); }
5597
5598 vert_index[0] = nodes[n0].vert_index;
5599 vert_index[1] = nodes[n1].vert_index;
5600
5601 if (oriented && vert_index[0] > vert_index[1])
5602 {
5603 std::swap(vert_index[0], vert_index[1]);
5604 }
5605}
5606
5608{
5609 const Element &el = elements[edge_id.element];
5610 const GeomInfo& gi = GI[el.Geom()];
5611 const int* ev = gi.edges[(int) edge_id.local];
5612
5613 int v0 = nodes[el.node[ev[0]]].vert_index;
5614 int v1 = nodes[el.node[ev[1]]].vert_index;
5615
5616 return ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1])) ? -1 : 1;
5617}
5618
5620 int vert_index[4], int edge_index[4],
5621 int edge_orientation[4]) const
5622{
5623 MFEM_ASSERT(Dim >= 3, "");
5624
5625 const Element &el = elements[face_id.element];
5626 const GeomInfo& gi = GI[el.Geom()];
5627
5628 const int *fv = gi.faces[(int) face_id.local];
5629 const int nfv = gi.nfv[(int) face_id.local];
5630
5631 vert_index[3] = edge_index[3] = -1;
5632 edge_orientation[3] = 0;
5633
5634 for (int i = 0; i < nfv; i++)
5635 {
5636 vert_index[i] = nodes[el.node[fv[i]]].vert_index;
5637 }
5638
5639 for (int i = 0; i < nfv; i++)
5640 {
5641 int j = i+1;
5642 if (j >= nfv) { j = 0; }
5643
5644 int n1 = el.node[fv[i]];
5645 int n2 = el.node[fv[j]];
5646
5647 const Node* en = nodes.Find(n1, n2);
5648 MFEM_ASSERT(en != NULL, "edge not found.");
5649
5650 edge_index[i] = en->edge_index;
5651 edge_orientation[i] = (vert_index[i] < vert_index[j]) ? 1 : -1;
5652 }
5653
5654 return nfv;
5655}
5656
5657int NCMesh::GetEdgeMaster(int node) const
5658{
5659 MFEM_ASSERT(node >= 0, "edge node not found.");
5660 const Node &nd = nodes[node];
5661
5662 int p1 = nd.p1, p2 = nd.p2;
5663 MFEM_ASSERT(p1 != p2, "invalid edge node.");
5664
5665 const Node &n1 = nodes[p1], &n2 = nodes[p2];
5666
5667 int n1p1 = n1.p1, n1p2 = n1.p2;
5668 int n2p1 = n2.p1, n2p2 = n2.p2;
5669
5670 if ((n2p1 != n2p2) && (p1 == n2p1 || p1 == n2p2))
5671 {
5672 // n1 is parent of n2: (n1)--(nd)--(n2)------(*)
5673 if (n2.HasEdge()) { return p2; }
5674 else { return GetEdgeMaster(p2); }
5675 }
5676
5677 if ((n1p1 != n1p2) && (p2 == n1p1 || p2 == n1p2))
5678 {
5679 // n2 is parent of n1: (n2)--(nd)--(n1)------(*)
5680 if (n1.HasEdge()) { return p1; }
5681 else { return GetEdgeMaster(p1); }
5682 }
5683
5684 return -1;
5685}
5686
5687int NCMesh::GetEdgeMaster(int v1, int v2) const
5688{
5689 int node = nodes.FindId(vertex_nodeId[v1], vertex_nodeId[v2]);
5690 MFEM_ASSERT(node >= 0 && nodes[node].HasEdge(), "(v1, v2) is not an edge.");
5691
5692 int master = GetEdgeMaster(node);
5693 return (master >= 0) ? nodes[master].edge_index : -1;
5694}
5695
5697{
5698 int elem = leaf_elements[i];
5699 int depth = 0, parent;
5700 while ((parent = elements[elem].parent) != -1)
5701 {
5702 elem = parent;
5703 depth++;
5704 }
5705 return depth;
5706}
5707
5709{
5710 int elem = leaf_elements[i];
5711 int parent, reduction = 1;
5712 while ((parent = elements[elem].parent) != -1)
5713 {
5714 if (elements[parent].ref_type & 1) { reduction *= 2; }
5715 if (elements[parent].ref_type & 2) { reduction *= 2; }
5716 if (elements[parent].ref_type & 4) { reduction *= 2; }
5717 elem = parent;
5718 }
5719 return reduction;
5720}
5721
5723 Array<int> &face_indices,
5724 Array<int> &face_attribs) const
5725{
5726 const Element &el = elements[leaf_elements[leaf_elem]];
5727 const GeomInfo& gi = GI[el.Geom()];
5728
5729 face_indices.SetSize(gi.nf);
5730 face_attribs.SetSize(gi.nf);
5731
5732 for (int i = 0; i < gi.nf; i++)
5733 {
5734 const int* fv = gi.faces[i];
5735 const Face *face = faces.Find(el.node[fv[0]], el.node[fv[1]],
5736 el.node[fv[2]], el.node[fv[3]]);
5737 MFEM_ASSERT(face, "face not found");
5738 face_indices[i] = face->index;
5739 face_attribs[i] = face->attribute;
5740 }
5741}
5742void NCMesh::FindFaceNodes(int face, int node[4]) const
5743{
5744 auto tmp = FindFaceNodes(face);
5745 std::copy(tmp.begin(), tmp.end(), node);
5746}
5747
5748std::array<int, 4> NCMesh::FindFaceNodes(int face) const
5749{
5750 return FindFaceNodes(faces[face]);
5751}
5752
5753std::array<int, 4> NCMesh::FindFaceNodes(const Face &fa) const
5754{
5755 // Obtain face nodes from one of its elements (note that face->p1, p2, p3
5756 // cannot be used directly since they are not in order and p4 is missing).
5757 int elem = fa.elem[0];
5758 if (elem < 0) { elem = fa.elem[1]; }
5759 MFEM_ASSERT(elem >= 0, "Face has no elements?");
5760
5761 const Element &el = elements[elem];
5762 int f = find_local_face(el.Geom(),
5763 find_node(el, fa.p1),
5764 find_node(el, fa.p2),
5765 find_node(el, fa.p3));
5766
5767 const int* fv = GI[el.Geom()].faces[f];
5768 std::array<int, 4> node;
5769 for (int i = 0; i < 4; i++)
5770 {
5771 node[i] = el.node[fv[i]];
5772 }
5773 return node;
5774}
5775
5776void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
5777 Array<int> &bdr_vertices, Array<int> &bdr_edges,
5778 Array<int> &bdr_faces)
5779{
5780 bdr_vertices.SetSize(0);
5781 bdr_edges.SetSize(0);
5782 bdr_faces.SetSize(0);
5783
5784 if (Dim == 3)
5785 {
5786 GetFaceList(); // make sure 'boundary_faces' is up to date
5787 for (int f : boundary_faces)
5788 {
5789 if (bdr_attr_is_ess[faces[f].attribute - 1])
5790 {
5791 auto node = FindFaceNodes(f);
5792 int nfv = (node[3] < 0) ? 3 : 4;
5793
5794 for (int j = 0; j < nfv; j++)
5795 {
5796 bdr_vertices.Append(nodes[node[j]].vert_index);
5797
5798 int enode = nodes.FindId(node[j], node[(j+1) % nfv]);
5799 MFEM_ASSERT(enode >= 0 && nodes[enode].HasEdge(), "Edge not found.");
5800 bdr_edges.Append(nodes[enode].edge_index);
5801
5802 while ((enode = GetEdgeMaster(enode)) >= 0)
5803 {
5804 // append master edges that may not be accessible from any
5805 // boundary element, this happens in 3D in re-entrant corners
5806 bdr_edges.Append(nodes[enode].edge_index);
5807 }
5808 }
5809
5810 // If the face is a slave face, collect its non-ghost master face
5811 const Face &face = faces[f];
5812
5813 const auto id_and_type = GetFaceList().GetMeshIdAndType(face.index);
5814 if (id_and_type.type == NCList::MeshIdType::SLAVE)
5815 {
5816 // A slave face must mark its masters
5817 const auto &slave_face_id = static_cast<const Slave&>(*id_and_type.id);
5818 bdr_faces.Append(slave_face_id.master);
5819 }
5820 }
5821 }
5822 }
5823 else if (Dim == 2)
5824 {
5825 GetFaceList();
5826 GetEdgeList(); // make sure 'boundary_faces' is up to date
5827
5828 for (int f : boundary_faces)
5829 {
5830 Face &face = faces[f];
5831 if (bdr_attr_is_ess[face.attribute - 1])
5832 {
5833 bdr_vertices.Append(nodes[face.p1].vert_index);
5834 bdr_vertices.Append(nodes[face.p3].vert_index);
5835 }
5836
5837 const auto id_and_type = GetEdgeList().GetMeshIdAndType(face.index);
5838 if (id_and_type.type == NCList::MeshIdType::SLAVE)
5839 {
5840 // A slave face must mark its masters
5841 const auto &slave_edge_id = static_cast<const Slave&>(*id_and_type.id);
5842 bdr_edges.Append(slave_edge_id.master);
5843 }
5844 }
5845 }
5846
5847 // Filter, sort and unique an array, so it contains only local unique values.
5848 auto FilterSortUnique = [](Array<int> &v, int N)
5849 {
5850 // Perform the O(N) filter before the O(NlogN) sort. begin -> it is only
5851 // entries < N.
5852 auto it = std::remove_if(v.begin(), v.end(), [N](int i) { return i >= N; });
5853 std::sort(v.begin(), it);
5854 v.SetSize(std::distance(v.begin(), std::unique(v.begin(), it)));
5855 };
5856
5857 FilterSortUnique(bdr_vertices, NVertices);
5858 FilterSortUnique(bdr_edges, NEdges);
5859 FilterSortUnique(bdr_faces, NFaces);
5860}
5861
5862int NCMesh::EdgeSplitLevel(int vn1, int vn2) const
5863{
5864 int mid = nodes.FindId(vn1, vn2);
5865 if (mid < 0 || !nodes[mid].HasVertex()) { return 0; }
5866 return 1 + std::max(EdgeSplitLevel(vn1, mid), EdgeSplitLevel(mid, vn2));
5867}
5868
5869int NCMesh::TriFaceSplitLevel(int vn1, int vn2, int vn3) const
5870{
5871 int mid[3];
5872 if (TriFaceSplit(vn1, vn2, vn3, mid) &&
5873 faces.FindId(vn1, vn2, vn3) < 0)
5874 {
5875 return 1 + max(TriFaceSplitLevel(vn1, mid[0], mid[2]),
5876 TriFaceSplitLevel(mid[0], vn2, mid[1]),
5877 TriFaceSplitLevel(mid[2], mid[1], vn3),
5878 TriFaceSplitLevel(mid[0], mid[1], mid[2]));
5879 }
5880
5881 return 0; // not split
5882}
5883
5884void NCMesh::QuadFaceSplitLevel(int vn1, int vn2, int vn3, int vn4,
5885 int& h_level, int& v_level) const
5886{
5887 int hl1, hl2, vl1, vl2;
5888 int mid[5];
5889 real_t scale;
5890
5891 switch (QuadFaceSplitType(vn1, vn2, vn3, vn4, scale, mid))
5892 {
5893 case 0: // not split
5894 h_level = v_level = 0;
5895 break;
5896
5897 case 1: // vertical
5898 QuadFaceSplitLevel(vn1, mid[0], mid[2], vn4, hl1, vl1);
5899 QuadFaceSplitLevel(mid[0], vn2, vn3, mid[2], hl2, vl2);
5900 h_level = std::max(hl1, hl2);
5901 v_level = std::max(vl1, vl2) + 1;
5902 break;
5903
5904 default: // horizontal
5905 QuadFaceSplitLevel(vn1, vn2, mid[1], mid[3], hl1, vl1);
5906 QuadFaceSplitLevel(mid[3], mid[1], vn3, vn4, hl2, vl2);
5907 h_level = std::max(hl1, hl2) + 1;
5908 v_level = std::max(vl1, vl2);
5909 }
5910}
5911
5912int NCMesh::QuadFaceSplitLevel(int vn1, int vn2, int vn3, int vn4) const
5913{
5914 int h_level, v_level;
5915 QuadFaceSplitLevel(vn1, vn2, vn3, vn4, h_level, v_level);
5916 return h_level + v_level;
5917}
5918
5919void NCMesh::CountSplits(int elem, int splits[3]) const
5920{
5921 const Element &el = elements[elem];
5922 const int* node = el.node;
5923 GeomInfo& gi = GI[el.Geom()];
5924
5925 int elevel[MaxElemEdges];
5926 for (int i = 0; i < gi.ne; i++)
5927 {
5928 const int* ev = gi.edges[i];
5929 elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
5930 }
5931
5932 int flevel[MaxElemFaces][2];
5933 if (Dim >= 3)
5934 {
5935 for (int i = 0; i < gi.nf; i++)
5936 {
5937 const int* fv = gi.faces[i];
5938 if (gi.nfv[i] == 4)
5939 {
5940 QuadFaceSplitLevel(node[fv[0]], node[fv[1]],
5941 node[fv[2]], node[fv[3]],
5942 flevel[i][1], flevel[i][0]);
5943 }
5944 else
5945 {
5946 flevel[i][1] = 0;
5947 flevel[i][0] =
5948 TriFaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]]);
5949 }
5950 }
5951 }
5952
5953 if (el.Geom() == Geometry::CUBE)
5954 {
5955 splits[0] = max(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
5956 elevel[0], elevel[2], elevel[4], elevel[6]);
5957
5958 splits[1] = max(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
5959 elevel[1], elevel[3], elevel[5], elevel[7]);
5960
5961 splits[2] = max(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
5962 elevel[8], elevel[9], elevel[10], elevel[11]);
5963 }
5964 else if (el.Geom() == Geometry::PRISM)
5965 {
5966 splits[0] = splits[1] = max(flevel[0][0], flevel[1][0], 0,
5967 flevel[2][0], flevel[3][0], flevel[4][0],
5968 elevel[0], elevel[1], elevel[2],
5969 elevel[3], elevel[4], elevel[5]);
5970
5971 splits[2] = max(flevel[2][1], flevel[3][1], flevel[4][1],
5972 elevel[6], elevel[7], elevel[8]);
5973 }
5974 else if (el.Geom() == Geometry::PYRAMID)
5975 {
5976 splits[0] = max(flevel[0][0], flevel[1][0], 0,
5977 flevel[2][0], flevel[3][0], flevel[4][0],
5978 elevel[0], elevel[1], elevel[2],
5979 elevel[3], elevel[4], elevel[5],
5980 elevel[6], elevel[7]);
5981
5982 splits[1] = splits[0];
5983 splits[2] = splits[0];
5984 }
5985 else if (el.Geom() == Geometry::TETRAHEDRON)
5986 {
5987 splits[0] = max(flevel[0][0], flevel[1][0], flevel[2][0], flevel[3][0],
5988 elevel[0], elevel[1], elevel[2], elevel[3], elevel[4], elevel[5]);
5989
5990 splits[1] = splits[0];
5991 splits[2] = splits[0];
5992 }
5993 else if (el.Geom() == Geometry::SQUARE)
5994 {
5995 splits[0] = max(elevel[0], elevel[2]);
5996 splits[1] = max(elevel[1], elevel[3]);
5997 }
5998 else if (el.Geom() == Geometry::TRIANGLE)
5999 {
6000 splits[0] = max(elevel[0], elevel[1], elevel[2]);
6001 splits[1] = splits[0];
6002 }
6003 else
6004 {
6005 MFEM_ABORT("Unsupported element geometry.");
6006 }
6007}
6008
6009void NCMesh::GetLimitRefinements(Array<Refinement> &refinements, int max_level)
6010{
6011 for (int i = 0; i < leaf_elements.Size(); i++)
6012 {
6013 if (IsGhost(elements[leaf_elements[i]])) { break; } // TODO: NElements
6014
6015 int splits[3];
6016 CountSplits(leaf_elements[i], splits);
6017
6018 char ref_type = 0;
6019 for (int k = 0; k < Dim; k++)
6020 {
6021 if (splits[k] > max_level)
6022 {
6023 ref_type |= (1 << k);
6024 }
6025 }
6026
6027 if (ref_type)
6028 {
6029 if (Iso)
6030 {
6031 // iso meshes should only be modified by iso refinements
6032 ref_type = 7;
6033 }
6034 refinements.Append(Refinement(i, ref_type));
6035 }
6036 }
6037}
6038
6039void NCMesh::LimitNCLevel(int max_nc_level)
6040{
6041 MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
6042
6043 while (1)
6044 {
6045 Array<Refinement> refinements;
6046 GetLimitRefinements(refinements, max_nc_level);
6047 if (!refinements.Size()) { break; }
6048 Refine(refinements);
6049 }
6050}
6051
6052
6053//// I/O ///////////////////////////////////////////////////////////////////////
6054
6055int NCMesh::PrintVertexParents(std::ostream *os) const
6056{
6057 if (!os)
6058 {
6059 // count vertex nodes with parents
6060 int nv = 0;
6061 for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
6062 {
6063 if (node->HasVertex() && node->p1 != node->p2) { nv++; }
6064 }
6065 return nv;
6066 }
6067 else
6068 {
6069 // print the relations
6070 bool uniform_scaling = true; // Check whether all nodes have scale 0.5
6071 for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
6072 {
6073 if (node->HasVertex() && node->p1 != node->p2)
6074 {
6075 MFEM_ASSERT(nodes[node->p1].HasVertex(), "");
6076 MFEM_ASSERT(nodes[node->p2].HasVertex(), "");
6077
6078 (*os) << node.index() << " " << node->p1 << " " << node->p2;
6079 if (node->GetScale() != 0.5) { uniform_scaling = false; }
6080 if (using_scaling) { (*os) << " " << node->GetScale(); }
6081 (*os) << "\n";
6082 }
6083 }
6084 MFEM_VERIFY(using_scaling || uniform_scaling, "NCMesh has nonuniform "
6085 "scaling. Call Mesh::SetScaledNCMesh first.");
6086 return 0;
6087 }
6088}
6089
6090void NCMesh::LoadVertexParents(std::istream &input)
6091{
6092 int nv;
6093 input >> nv;
6094 while (nv--)
6095 {
6096 int id, p1, p2;
6097 real_t s{0.5};
6098 input >> id >> p1 >> p2;
6099 if (using_scaling) { input >> s; }
6100 MFEM_VERIFY(input, "problem reading vertex parents.");
6101
6102 MFEM_VERIFY(nodes.IdExists(id), "vertex " << id << " not found.");
6103 MFEM_VERIFY(nodes.IdExists(p1), "parent " << p1 << " not found.");
6104 MFEM_VERIFY(nodes.IdExists(p2), "parent " << p2 << " not found.");
6105
6106 int check = nodes.FindId(p1, p2);
6107 MFEM_VERIFY(check < 0, "parents (" << p1 << ", " << p2 << ") already "
6108 "assigned to node " << check);
6109
6110 // assign new parents for the node
6111 nodes.Reparent(id, p1, p2);
6112
6113 nodes[id].SetScale(s);
6114 }
6115}
6116
6117int NCMesh::PrintBoundary(std::ostream *os) const
6118{
6119 static const int nfv2geom[5] =
6120 {
6123 };
6124 int deg = (Dim == 2) ? 2 : 1; // for degenerate faces in 2D
6125
6126 int count = 0;
6127 for (int i = 0; i < elements.Size(); i++)
6128 {
6129 const Element &el = elements[i];
6130 if (!el.IsLeaf()) { continue; }
6131
6132 GeomInfo& gi = GI[el.Geom()];
6133 for (int k = 0; k < gi.nf; k++)
6134 {
6135 const int* fv = gi.faces[k];
6136 const int nfv = gi.nfv[k];
6137 const Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
6138 el.node[fv[2]], el.node[fv[3]]);
6139 MFEM_ASSERT(face != NULL, "face not found");
6140 if (face->Boundary())
6141 {
6142 if (!os) { count++; continue; }
6143
6144 (*os) << face->attribute << " " << nfv2geom[nfv];
6145 for (int j = 0; j < nfv; j++)
6146 {
6147 (*os) << " " << el.node[fv[j*deg]];
6148 }
6149 (*os) << "\n";
6150 }
6151 }
6152 }
6153 return count;
6154}
6155
6156void NCMesh::LoadBoundary(std::istream &input)
6157{
6158 int nb, attr, geom;
6159 input >> nb;
6160 for (int i = 0; i < nb; i++)
6161 {
6162 input >> attr >> geom;
6163
6164 int v1, v2, v3, v4;
6165 if (geom == Geometry::SQUARE)
6166 {
6167 input >> v1 >> v2 >> v3 >> v4;
6168 Face* face = faces.Get(v1, v2, v3, v4);
6169 face->attribute = attr;
6170 }
6171 else if (geom == Geometry::TRIANGLE)
6172 {
6173 input >> v1 >> v2 >> v3;
6174 Face* face = faces.Get(v1, v2, v3);
6175 face->attribute = attr;
6176 }
6177 else if (geom == Geometry::SEGMENT)
6178 {
6179 input >> v1 >> v2;
6180 Face* face = faces.Get(v1, v1, v2, v2);
6181 face->attribute = attr;
6182 }
6183 else if (geom == Geometry::POINT)
6184 {
6185 input >> v1;
6186 Face* face = faces.Get(v1, v1, v1, v1);
6187 face->attribute = attr;
6188 }
6189 else
6190 {
6191 MFEM_ABORT("unsupported boundary element geometry: " << geom);
6192 }
6193 }
6194}
6195
6196void NCMesh::PrintCoordinates(std::ostream &os) const
6197{
6198 int nv = coordinates.Size()/3;
6199 os << nv << "\n";
6200 if (!nv) { return; }
6201
6202 os << spaceDim << "\n";
6203 for (int i = 0; i < nv; i++)
6204 {
6205 os << coordinates[3*i];
6206 for (int j = 1; j < spaceDim; j++)
6207 {
6208 os << " " << coordinates[3*i + j];
6209 }
6210 os << "\n";
6211 }
6212}
6213
6214void NCMesh::LoadCoordinates(std::istream &input)
6215{
6216 int nv;
6217 input >> nv;
6218 if (!nv) { return; }
6219
6220 input >> spaceDim;
6221
6222 coordinates.SetSize(nv * 3);
6223 coordinates = 0.0;
6224
6225 for (int i = 0; i < nv; i++)
6226 {
6227 for (int j = 0; j < spaceDim; j++)
6228 {
6229 input >> coordinates[3*i + j];
6230 MFEM_VERIFY(input.good(), "unexpected EOF");
6231 }
6232 }
6233}
6234
6236{
6237 for (int i = 0; i < root_state.Size(); i++)
6238 {
6239 if (root_state[i]) { return false; }
6240 }
6241 return true;
6242}
6243
6244void NCMesh::Print(std::ostream &os, const std::string &comments) const
6245{
6246 if (using_scaling)
6247 {
6248 os << "MFEM NC mesh v1.1\n\n";
6249 }
6250 else
6251 {
6252 os << "MFEM NC mesh v1.0\n\n";
6253 }
6254
6255 if (!comments.empty()) { os << comments << "\n\n"; }
6256
6257 os <<
6258 "# NCMesh supported geometry types:\n"
6259 "# SEGMENT = 1\n"
6260 "# TRIANGLE = 2\n"
6261 "# SQUARE = 3\n"
6262 "# TETRAHEDRON = 4\n"
6263 "# CUBE = 5\n"
6264 "# PRISM = 6\n"
6265 "# PYRAMID = 7\n";
6266
6267 os << "\ndimension\n" << Dim << "\n";
6268
6269#ifndef MFEM_USE_MPI
6270 if (MyRank != 0) // don't print this section in serial: default rank is 0
6271#endif
6272 {
6273 os << "\nrank\n" << MyRank << "\n";
6274 }
6275
6276 os << "\n# rank attr geom ref_type nodes/children";
6277 os << "\nelements\n" << elements.Size() << "\n";
6278
6279 for (int i = 0; i < elements.Size(); i++)
6280 {
6281 const Element &el = elements[i];
6282 os << el.rank << " " << el.attribute << " ";
6283 if (el.parent == -2) { os << "-1\n"; continue; } // unused element
6284
6285 os << int(el.geom) << " " << int(el.ref_type);
6286 for (int j = 0; j < MaxElemNodes && el.node[j] >= 0; j++)
6287 {
6288 os << " " << el.node[j];
6289 }
6290 os << "\n";
6291 }
6292
6293 int nb = PrintBoundary(NULL);
6294 if (nb)
6295 {
6296 os << "\n# attr geom nodes";
6297 os << "\nboundary\n" << nb << "\n";
6298
6299 PrintBoundary(&os);
6300 }
6301
6302 int nvp = PrintVertexParents(NULL);
6303 if (nvp)
6304 {
6305 os << "\n# vert_id p1 p2";
6306 os << "\nvertex_parents\n" << nvp << "\n";
6307
6308 PrintVertexParents(&os);
6309 }
6310
6311 if (!ZeroRootStates()) // root_state section is optional
6312 {
6313 os << "\n# root element orientation";
6314 os << "\nroot_state\n" << root_state.Size() << "\n";
6315
6316 for (int i = 0; i < root_state.Size(); i++)
6317 {
6318 os << root_state[i] << "\n";
6319 }
6320 }
6321
6322 if (coordinates.Size())
6323 {
6324 os << "\n# top-level node coordinates";
6325 os << "\ncoordinates\n";
6326
6327 PrintCoordinates(os);
6328 }
6329 else
6330 {
6331 // 'nodes' will be printed one level up by Mesh::Printer()
6332 }
6333}
6334
6336{
6337 // initialize Element::parent
6338 for (int i = 0; i < elements.Size(); i++)
6339 {
6340 Element &el = elements[i];
6341 if (el.ref_type)
6342 {
6343 for (int j = 0; j < MaxElemChildren && el.child[j] >= 0; j++)
6344 {
6345 int child = el.child[j];
6346 MFEM_VERIFY(child < elements.Size(), "invalid mesh file: "
6347 "element " << i << " references invalid child " << child);
6348 elements[child].parent = i;
6349 }
6350 }
6351 }
6352
6353 // count the root elements
6354 int nroots = 0;
6355 for (const auto &e : elements)
6356 if (e.parent == -1)
6357 {
6358 ++nroots;
6359 }
6360 MFEM_VERIFY(nroots > 0 ||
6361 elements.Size() == 0,
6362 "invalid mesh file: no root elements in non-empty mesh found.");
6363
6364
6365 // check that only the first 'nroot' elements are roots (have no parent)
6366 for (int i = nroots; i < elements.Size(); i++)
6367 {
6368 MFEM_VERIFY(elements[i].parent != -1,
6369 "invalid mesh file: only the first M elements can be roots. "
6370 "Found element " << i << " with no parent.");
6371 }
6372
6373 // default root state is zero
6374 root_state.SetSize(nroots);
6375 root_state = 0;
6376}
6377
6379{
6380 int ntop = 0;
6381 for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
6382 {
6383 if (node->p1 == node->p2) { ntop = node.index() + 1; }
6384 }
6385 return ntop;
6386}
6387
6388NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
6389 : spaceDim(0), MyRank(0), Iso(true), Legacy(false),
6390 using_scaling(version == 11)
6391{
6392 is_nc = 1;
6393 if (version == 1) // old MFEM mesh v1.1 format
6394 {
6395 LoadLegacyFormat(input, curved, is_nc);
6396 Legacy = true;
6397 return;
6398 }
6399
6400 MFEM_ASSERT(version == 10 || version == 11, "");
6401 std::string ident;
6402 int count;
6403
6404 // Skip the version string
6405 skip_comment_lines(input, 'M');
6406
6407 // load dimension
6408 skip_comment_lines(input, '#');
6409 input >> ident;
6410 MFEM_VERIFY(ident == "dimension", "Invalid mesh file: " << ident);
6411 input >> Dim;
6412
6413 // load rank, if present
6414 skip_comment_lines(input, '#');
6415 input >> ident;
6416 if (ident == "rank")
6417 {
6418 input >> MyRank;
6419 MFEM_VERIFY(MyRank >= 0, "Invalid rank");
6420
6421 skip_comment_lines(input, '#');
6422 input >> ident;
6423 }
6424
6425 // load file SFC version, if present (for future changes to SFC ordering)
6426 if (ident == "sfc_version")
6427 {
6428 int sfc_version; // TODO future: store as class member
6429 input >> sfc_version;
6430 MFEM_VERIFY(sfc_version == 0,
6431 "Unsupported mesh file SFC version (" << sfc_version << "). "
6432 "Please update MFEM.");
6433
6434 skip_comment_lines(input, '#');
6435 input >> ident;
6436 }
6437
6438 // load elements
6439 MFEM_VERIFY(ident == "elements", "Invalid mesh file: " << ident);
6440 input >> count;
6441 for (int i = 0; i < count; i++)
6442 {
6443 int rank, attr, geom, ref_type;
6444 input >> rank >> attr >> geom;
6445
6446 Geometry::Type type = Geometry::Type(geom);
6447 elements.Append(Element(type, attr));
6448
6449 MFEM_ASSERT(elements.Size() == i+1, "");
6450 Element &el = elements[i];
6451 el.rank = rank;
6452
6453 if (geom >= 0)
6454 {
6455 CheckSupportedGeom(type);
6456 GI[geom].InitGeom(type);
6457
6458 input >> ref_type;
6459 MFEM_VERIFY(ref_type >= 0 && ref_type < 8, "");
6460 el.ref_type = ref_type;
6461
6462 if (ref_type) // refined element
6463 {
6464 for (int j = 0; j < ref_type_num_children[ref_type]; j++)
6465 {
6466 input >> el.child[j];
6467 }
6468 if (Dim == 3 && ref_type != 7) { Iso = false; }
6469 }
6470 else // leaf element
6471 {
6472 for (int j = 0; j < GI[geom].nv; j++)
6473 {
6474 int id;
6475 input >> id;
6476 el.node[j] = id;
6477 nodes.Alloc(id, id, id);
6478 // NOTE: nodes that won't get parents assigned will stay hashed
6479 // with p1 == p2 == id (top-level nodes)
6480 }
6481 }
6482 }
6483 else
6484 {
6485 el.parent = -2; // mark as unused
6487 }
6488 }
6489
6491 InitGeomFlags();
6492
6493 skip_comment_lines(input, '#');
6494 input >> ident;
6495
6496 // load boundary
6497 if (ident == "boundary")
6498 {
6499 LoadBoundary(input);
6500
6501 skip_comment_lines(input, '#');
6502 input >> ident;
6503 }
6504
6505 // load vertex hierarchy
6506 if (ident == "vertex_parents")
6507 {
6508 LoadVertexParents(input);
6509
6510 skip_comment_lines(input, '#');
6511 input >> ident;
6512 }
6513
6514 // load root states
6515 if (ident == "root_state")
6516 {
6517 input >> count;
6518 MFEM_VERIFY(count <= root_state.Size(), "Too many root states");
6519 for (int i = 0; i < count; i++)
6520 {
6521 input >> root_state[i];
6522 }
6523
6524 skip_comment_lines(input, '#');
6525 input >> ident;
6526 }
6527
6528 // load coordinates or nodes
6529 if (ident == "coordinates")
6530 {
6531 LoadCoordinates(input);
6532
6533 MFEM_VERIFY(coordinates.Size() >= 3*CountTopLevelNodes(),
6534 "Invalid mesh file: not all top-level nodes are covered by "
6535 "the 'coordinates' section of the mesh file: " << coordinates.Size() << ' ' <<
6536 3*CountTopLevelNodes());
6537 curved = 0;
6538 }
6539 else if (ident == "nodes")
6540 {
6541 coordinates.SetSize(0); // this means the mesh is curved
6542
6543 // prepare to read the nodes
6544 input >> std::ws;
6545 curved = 1;
6546 }
6547 else
6548 {
6549 MFEM_ABORT("Invalid mesh file: either 'coordinates' or "
6550 "'nodes' must be present");
6551 }
6552
6553 // create edge nodes and faces
6554 nodes.UpdateUnused();
6555 for (int i = 0; i < elements.Size(); i++)
6556 {
6557 if (elements[i].IsLeaf())
6558 {
6560 RegisterFaces(i);
6561 }
6562 }
6563
6564 Update();
6565}
6566
6568 const BlockArray<Element> &tmp_elements)
6569{
6570 Element &el = elements[elem];
6571 if (el.ref_type)
6572 {
6573 for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
6574 {
6575 int old_id = el.child[i];
6576 // here we know 'free_element_ids' is empty
6577 int new_id = elements.Append(tmp_elements[old_id]);
6578 el.child[i] = new_id;
6579 elements[new_id].parent = elem;
6580 CopyElements(new_id, tmp_elements);
6581 }
6582 }
6583}
6584
6585void NCMesh::LoadCoarseElements(std::istream &input)
6586{
6587 int ne;
6588 input >> ne;
6589
6590 bool iso = true;
6591
6592 // load the coarse elements
6593 while (ne--)
6594 {
6595 int ref_type;
6596 input >> ref_type;
6597
6598 int elem = AddElement(Geometry::INVALID, 0);
6599 Element &el = elements[elem];
6600 el.ref_type = ref_type;
6601
6602 if (Dim == 3 && ref_type != 7) { iso = false; }
6603
6604 // load child IDs and make parent-child links
6605 int nch = ref_type_num_children[ref_type];
6606 for (int i = 0, id; i < nch; i++)
6607 {
6608 input >> id;
6609 MFEM_VERIFY(id >= 0, "");
6610 MFEM_VERIFY(id < elements.Size(),
6611 "coarse element cannot be referenced before it is "
6612 "defined (id=" << id << ").");
6613
6614 Element &child = elements[id];
6615 MFEM_VERIFY(child.parent == -1,
6616 "element " << id << " cannot have two parents.");
6617
6618 el.child[i] = id;
6619 child.parent = elem;
6620
6621 if (!i) // copy geom and attribute from first child
6622 {
6623 el.geom = child.geom;
6624 el.attribute = child.attribute;
6625 }
6626 }
6627 }
6628
6629 // prepare for reordering the elements
6630 BlockArray<Element> tmp_elements;
6631 elements.Swap(tmp_elements);
6632
6633 // copy roots, they need to be at the beginning of 'elements'
6634 int root_count = 0;
6635 for (auto el = tmp_elements.begin(); el != tmp_elements.end(); ++el)
6636 {
6637 if (el->parent == -1)
6638 {
6639 elements.Append(*el); // same as AddElement()
6640 root_count++;
6641 }
6642 }
6643
6644 // copy the rest of the hierarchy
6645 for (int i = 0; i < root_count; i++)
6646 {
6647 CopyElements(i, tmp_elements);
6648 }
6649
6650 // set the Iso flag (must be false if there are 3D aniso refinements)
6651 Iso = iso;
6652
6653 InitRootState(root_count);
6654}
6655
6656void NCMesh::LoadLegacyFormat(std::istream &input, int &curved, int &is_nc)
6657{
6658 MFEM_ASSERT(elements.Size() == 0, "");
6659 MFEM_ASSERT(nodes.Size() == 0, "");
6660 MFEM_ASSERT(free_element_ids.Size() == 0, "");
6661
6662 std::string ident;
6663 int count, attr, geom;
6664
6665 // load dimension
6666 skip_comment_lines(input, '#');
6667 input >> ident;
6668 MFEM_VERIFY(ident == "dimension", "invalid mesh file");
6669 input >> Dim;
6670
6671 // load elements
6672 skip_comment_lines(input, '#');
6673 input >> ident;
6674 MFEM_VERIFY(ident == "elements", "invalid mesh file");
6675
6676 input >> count;
6677 for (int i = 0; i < count; i++)
6678 {
6679 input >> attr >> geom;
6680
6681 Geometry::Type type = Geometry::Type(geom);
6682 CheckSupportedGeom(type);
6683 GI[geom].InitGeom(type);
6684
6685 int eid = AddElement(type, attr);
6686 MFEM_ASSERT(eid == i, "");
6687
6688 Element &el = elements[eid];
6689 for (int j = 0; j < GI[geom].nv; j++)
6690 {
6691 int id;
6692 input >> id;
6693 el.node[j] = id;
6694 nodes.Alloc(id, id, id); // see comment in NCMesh::NCMesh
6695 }
6696 el.index = i; // needed for file leaf order below
6697 }
6698
6699 // load boundary
6700 skip_comment_lines(input, '#');
6701 input >> ident;
6702 MFEM_VERIFY(ident == "boundary", "invalid mesh file");
6703
6704 LoadBoundary(input);
6705
6706 // load vertex hierarchy
6707 skip_comment_lines(input, '#');
6708 input >> ident;
6709 if (ident == "vertex_parents")
6710 {
6711 LoadVertexParents(input);
6712 is_nc = 1;
6713
6714 skip_comment_lines(input, '#');
6715 input >> ident;
6716 }
6717 else
6718 {
6719 // no "vertex_parents" section: this file needs to be treated as a
6720 // conforming mesh for complete backward compatibility with MFEM 4.2
6721 is_nc = 0;
6722 }
6723
6724 // load element hierarchy
6725 if (ident == "coarse_elements")
6726 {
6727 LoadCoarseElements(input);
6728
6729 skip_comment_lines(input, '#');
6730 input >> ident;
6731 }
6732 else
6733 {
6734 // no element hierarchy -> all elements are roots
6735 InitRootState(elements.Size());
6736 }
6737 InitGeomFlags();
6738
6739 // load vertices
6740 MFEM_VERIFY(ident == "vertices", "invalid mesh file");
6741 int nvert;
6742 input >> nvert;
6743 input >> std::ws >> ident;
6744 if (ident != "nodes")
6745 {
6746 spaceDim = atoi(ident.c_str());
6747
6748 coordinates.SetSize(3*nvert);
6749 coordinates = 0.0;
6750
6751 for (int i = 0; i < nvert; i++)
6752 {
6753 for (int j = 0; j < spaceDim; j++)
6754 {
6755 input >> coordinates[3*i + j];
6756 MFEM_VERIFY(input.good(), "unexpected EOF");
6757 }
6758 }
6759
6760 // truncate extra coordinates (legacy vertices section is longer)
6761 int ntop = CountTopLevelNodes();
6762 if (3*ntop < coordinates.Size())
6763 {
6764 coordinates.SetSize(3*ntop);
6765 }
6766 }
6767 else
6768 {
6770
6771 // prepare to read the nodes
6772 input >> std::ws;
6773 curved = 1;
6774 }
6775
6776 // create edge nodes and faces
6777 nodes.UpdateUnused();
6778 int leaf_count = 0;
6779 for (int i = 0; i < elements.Size(); i++)
6780 {
6781 if (elements[i].IsLeaf())
6782 {
6784 RegisterFaces(i);
6785 leaf_count++;
6786 }
6787 }
6788
6789 // v1.1 honors file leaf order on load, prepare legacy 'leaf_elements'
6790 Array<int> file_leaf_elements(leaf_count);
6791 file_leaf_elements = -1;
6792 for (int i = 0; i < elements.Size(); i++)
6793 {
6794 if (elements[i].IsLeaf())
6795 {
6796 file_leaf_elements[elements[i].index] = i;
6797 }
6798 }
6799 MFEM_ASSERT(file_leaf_elements.Min() >= 0, "");
6800
6801 Update();
6802
6803 // force file leaf order
6804 Swap(leaf_elements, file_leaf_elements);
6805
6806 // make sure Mesh::NVertices is equal to "nvert" from the file (in case of
6807 // unused vertices), see also GetMeshComponents
6808 if (nvert > vertex_nodeId.Size())
6809 {
6810 vertex_nodeId.SetSize(nvert, -1);
6811 }
6812}
6813
6815{
6816 order.SetSize(NVertices);
6817 order = -1;
6818
6819 int count = 0;
6820 for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
6821 {
6822 if (node->HasVertex())
6823 {
6824 MFEM_ASSERT(node.index() >= 0, "");
6825 MFEM_ASSERT(node.index() < order.Size(), "");
6826 MFEM_ASSERT(order[node.index()] == -1, "");
6827
6828 order[node.index()] = node->vert_index;
6829 count++;
6830 }
6831 }
6832 MFEM_ASSERT(count == order.Size(), "");
6833 MFEM_CONTRACT_VAR(count);
6834}
6835
6836
6837////////////////////////////////////////////////////////////////////////////////
6838
6840{
6842 face_list.Clear();
6843 edge_list.Clear();
6844
6847
6849
6850 // TODO future: consider trimming unused blocks at the end of 'elements' and
6851 // maybe also of 'nodes' and 'faces'.
6852}
6853
6855{
6856 int pm_size = 0;
6857 for (int i = 0; i < Geometry::NumGeom; i++)
6858 {
6859 for (int j = 0; j < point_matrices[i].Size(); i++)
6860 {
6861 pm_size += static_cast<int>(point_matrices[i][j]->MemoryUsage());
6862 }
6863 pm_size += static_cast<int>(point_matrices[i].MemoryUsage());
6864 }
6865
6866 return conforming.MemoryUsage() +
6867 masters.MemoryUsage() +
6868 slaves.MemoryUsage() +
6869 pm_size;
6870}
6871
6873{
6874 long mem = embeddings.MemoryUsage();
6875 for (int i = 0; i < Geometry::NumGeom; i++)
6876 {
6877 mem += point_matrices[i].MemoryUsage();
6878 }
6879 return mem;
6880}
6881
6883{
6884 return nodes.MemoryUsage() +
6885 faces.MemoryUsage() +
6886 elements.MemoryUsage() +
6898 ref_stack.MemoryUsage() +
6902 sizeof(*this);
6903}
6904
6906{
6907 nodes.PrintMemoryDetail(); mfem::out << " nodes\n";
6908 faces.PrintMemoryDetail(); mfem::out << " faces\n";
6909
6910 mfem::out << elements.MemoryUsage() << " elements\n"
6911 << free_element_ids.MemoryUsage() << " free_element_ids\n"
6912 << root_state.MemoryUsage() << " root_state\n"
6913 << coordinates.MemoryUsage() << " top_vertex_pos\n"
6914 << leaf_elements.MemoryUsage() << " leaf_elements\n"
6915 << leaf_sfc_index.MemoryUsage() << " leaf_sfc_index\n"
6916 << vertex_nodeId.MemoryUsage() << " vertex_nodeId\n"
6917 << face_list.MemoryUsage() << " face_list\n"
6918 << edge_list.MemoryUsage() << " edge_list\n"
6919 << vertex_list.MemoryUsage() << " vertex_list\n"
6920 << boundary_faces.MemoryUsage() << " boundary_faces\n"
6921 << element_vertex.MemoryUsage() << " element_vertex\n"
6922 << ref_stack.MemoryUsage() << " ref_stack\n"
6923 << derefinements.MemoryUsage() << " derefinements\n"
6924 << transforms.MemoryUsage() << " transforms\n"
6925 << coarse_elements.MemoryUsage() << " coarse_elements\n"
6926 << sizeof(*this) << " NCMesh"
6927 << std::endl;
6928
6929 return elements.Size() - free_element_ids.Size();
6930}
6931
6932#ifdef MFEM_DEBUG
6933void NCMesh::DebugLeafOrder(std::ostream &os) const
6934{
6935 tmp_vertex = new TmpVertex[nodes.NumIds()];
6936 for (int i = 0; i < leaf_elements.Size(); i++)
6937 {
6938 const Element* elem = &elements[leaf_elements[i]];
6939 for (int j = 0; j < Dim; j++)
6940 {
6941 real_t sum = 0.0;
6942 int count = 0;
6943 for (int k = 0; k < MaxElemNodes; k++)
6944 {
6945 if (elem->node[k] >= 0)
6946 {
6947 sum += CalcVertexPos(elem->node[k])[j];
6948 count++;
6949 }
6950 }
6951 os << sum / count << " ";
6952 }
6953 os << "\n";
6954 }
6955 delete [] tmp_vertex;
6956}
6957
6958void NCMesh::DebugDump(std::ostream &os) const
6959{
6960 // dump nodes
6961 tmp_vertex = new TmpVertex[nodes.NumIds()];
6962 os << nodes.Size() << "\n";
6963 for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
6964 {
6965 const real_t *pos = CalcVertexPos(node.index());
6966 os << node.index() << " "
6967 << pos[0] << " " << pos[1] << " " << pos[2] << " "
6968 << node->p1 << " " << node->p2 << " "
6969 << node->vert_index << " " << node->edge_index << " "
6970 << 0 << "\n";
6971 }
6972 delete [] tmp_vertex;
6973 os << "\n";
6974
6975 // dump elements
6976 int nleaves = 0;
6977 for (int i = 0; i < elements.Size(); i++)
6978 {
6979 if (elements[i].IsLeaf()) { nleaves++; }
6980 }
6981 os << nleaves << "\n";
6982 for (int i = 0; i < elements.Size(); i++)
6983 {
6984 const Element &el = elements[i];
6985 if (el.IsLeaf())
6986 {
6987 const GeomInfo& gi = GI[el.Geom()];
6988 os << gi.nv << " ";
6989 for (int j = 0; j < gi.nv; j++)
6990 {
6991 os << el.node[j] << " ";
6992 }
6993 os << el.attribute << " " << el.rank << " " << i << "\n";
6994 }
6995 }
6996 os << "\n";
6997
6998 // dump faces
6999 os << faces.Size() << "\n";
7000 for (const auto &face : faces)
7001 {
7002 int elem = face.elem[0];
7003 if (elem < 0) { elem = face.elem[1]; }
7004 MFEM_ASSERT(elem >= 0, "");
7005 const Element &el = elements[elem];
7006
7007 int lf = find_local_face(el.Geom(),
7008 find_node(el, face.p1),
7009 find_node(el, face.p2),
7010 find_node(el, face.p3));
7011
7012 const int* fv = GI[el.Geom()].faces[lf];
7013 const int nfv = GI[el.Geom()].nfv[lf];
7014
7015 os << nfv;
7016 for (int i = 0; i < nfv; i++)
7017 {
7018 os << " " << el.node[fv[i]];
7019 }
7020 //os << " # face " << face.index() << ", index " << face.index << "\n";
7021 os << "\n";
7022 }
7023}
7024#endif
7025
7026} // namespace mfem
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition array.hpp:278
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition array.hpp:165
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
Definition array.cpp:85
int Size() const
Return the logical size of the array.
Definition array.hpp:147
void DeleteAll()
Delete the whole array.
Definition array.hpp:925
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
T * GetData()
Returns the data.
Definition array.hpp:121
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition array.hpp:935
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
Definition array.hpp:286
T * end()
STL-like end. Returns pointer after the last element of the array.
Definition array.hpp:325
T * begin()
STL-like begin. Returns pointer to the first element of the array.
Definition array.hpp:322
std::size_t MemoryUsage() const
Returns the number of bytes allocated for the array including any reserve.
Definition array.hpp:334
iterator begin()
Definition array.hpp:636
void Swap(BlockArray< T > &other)
Definition array.hpp:1123
iterator end()
Definition array.hpp:637
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Abstract data type element.
Definition element.hpp:29
Geometry::Type GetGeometryType() const
Definition element.hpp:55
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Definition element.hpp:61
virtual Type GetType() const =0
Returns element's type.
int GetAttribute() const
Return element's attribute.
Definition element.hpp:58
static const int NumGeom
Definition geom.hpp:46
Data type hexahedron element.
Class used by MFEM to store pointers to host and/or device memory.
Mesh data type.
Definition mesh.hpp:64
Array< Vertex > vertices
Definition mesh.hpp:105
Element * NewElement(int geom)
Definition mesh.cpp:4672
int GetNEdges() const
Return the number of edges.
Definition mesh.hpp:1288
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6531
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1339
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
const Element * GetFace(int i) const
Return pointer to the i'th face element object.
Definition mesh.hpp:1366
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1354
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition mesh.hpp:1279
GridFunction * Nodes
Definition mesh.hpp:260
Array< Triple< int, int, int > > tmp_vertex_parents
Definition mesh.hpp:274
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1285
void FreeElement(Element *E)
Definition mesh.cpp:13381
Array< Element * > boundary
Definition mesh.hpp:106
Table * GetEdgeVertexTable() const
Definition mesh.cpp:7388
Array< Element * > elements
Definition mesh.hpp:100
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1321
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
Definition ncmesh.hpp:140
int GetEdgeMaster(int v1, int v2) const
Definition ncmesh.cpp:5687
static constexpr int MaxElemEdges
Number of edges an element can have.
Definition ncmesh.hpp:490
static PointMatrix pm_tet_identity
Definition ncmesh.hpp:1179
void OnMeshUpdated(Mesh *mesh)
Definition ncmesh.cpp:2885
void GetElementFacesAttributes(int i, Array< int > &faces, Array< int > &fattr) const
Return the faces and face attributes of leaf element i.
Definition ncmesh.cpp:5722
static GeomInfo GI[Geometry::NumGeom]
Definition ncmesh.hpp:1340
virtual void Update()
Definition ncmesh.cpp:268
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
Definition ncmesh.cpp:4320
void FreeElement(int id)
Definition ncmesh.hpp:826
virtual void Trim()
Save memory by releasing all non-essential and cached data.
Definition ncmesh.cpp:6839
virtual void ElementSharesVertex(int elem, int local, int vnode)
Definition ncmesh.hpp:1014
HashTable< Node > shadow
temporary storage for reparented nodes
Definition ncmesh.hpp:792
void SetDerefMatrixCodes(int parent, Array< int > &fine_coarse)
Definition ncmesh.cpp:2362
void CopyElements(int elem, const BlockArray< Element > &tmp_elements)
Definition ncmesh.cpp:6567
Array< Refinement > ref_stack
stack of scheduled refinements (temporary)
Definition ncmesh.hpp:791
std::int64_t RefCoord
Definition ncmesh.hpp:486
const Face & GetFace(int i) const
Access a Face.
Definition ncmesh.hpp:652
mfem::Element * NewMeshElement(int geom) const
Definition ncmesh.cpp:2718
void SetNodeScale(int p0, int p1, real_t scale)
Helper function to set scale for a node with parents p0, p1.
Definition ncmesh.cpp:1115
int NewTriangle(int n0, int n1, int n2, int attr, int eattr0, int eattr1, int eattr2)
Definition ncmesh.cpp:758
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
Definition ncmesh.cpp:2758
NCMesh()=default
int NGhostElements
Definition ncmesh.hpp:721
void LoadCoarseElements(std::istream &input)
Load the element refinement hierarchy from a legacy mesh file.
Definition ncmesh.cpp:6585
void NeighborExpand(const Array< int > &elems, Array< int > &expanded, const Array< int > *search_set=NULL)
Definition ncmesh.cpp:4403
HashTable< Node > nodes
Definition ncmesh.hpp:619
int PrintVertexParents(std::ostream *out) const
Print the "vertex_parents" section of the mesh file.
Definition ncmesh.cpp:6055
static PointMatrix pm_prism_identity
Definition ncmesh.hpp:1180
int NewSegment(int n0, int n1, int attr, int vattr1, int vattr2)
Definition ncmesh.cpp:784
static int find_node(const Element &el, int node)
Definition ncmesh.cpp:3276
void DeleteUnusedFaces(const Array< int > &elemFaces)
Definition ncmesh.cpp:455
const CoarseFineTransformations & GetRefinementTransforms() const
Definition ncmesh.cpp:5153
void BuildElementToVertexTable()
Definition ncmesh.cpp:4129
void TraverseEdge(int vn0, int vn1, real_t t0, real_t t1, int flags, int level, MatrixMap &matrix_map)
Definition ncmesh.cpp:3770
bool using_scaling
Definition ncmesh.hpp:622
static int find_element_edge(const Element &el, int vn0, int vn1, bool abort=true)
Definition ncmesh.cpp:3296
Array< int > free_element_ids
Definition ncmesh.hpp:625
void LoadBoundary(std::istream &input)
Load the "boundary" section of the mesh file.
Definition ncmesh.cpp:6156
virtual void ElementSharesFace(int elem, int local, int face)
Definition ncmesh.hpp:1012
int PrintMemoryDetail() const
Definition ncmesh.cpp:6905
Array< real_t > reparent_scale
scale associated with reparents (tmp)
Definition ncmesh.hpp:794
static constexpr int MaxElemNodes
Number of nodes an element can have.
Definition ncmesh.hpp:488
int NewHexahedron(int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int attr, int fattr0, int fattr1, int fattr2, int fattr3, int fattr4, int fattr5)
Definition ncmesh.cpp:615
HashTable< Face > faces
Definition ncmesh.hpp:620
bool HaveTets() const
Return true if the mesh contains tetrahedral elements.
Definition ncmesh.hpp:784
void GetEdgeVertices(const MeshId &edge_id, int vert_index[2], bool oriented=true) const
Return Mesh vertex indices of an edge identified by 'edge_id'.
Definition ncmesh.cpp:5588
static PointMatrix pm_seg_identity
Definition ncmesh.hpp:1176
bool QuadFaceIsMaster(int n1, int n2, int n3, int n4) const
Determine if a Quad face is a master face.
Definition ncmesh.hpp:930
Array< int > boundary_faces
subset of all faces, set by BuildFaceList
Definition ncmesh.hpp:731
friend struct MatrixMap
Definition ncmesh.hpp:1349
std::array< int, 4 > FindFaceNodes(int face) const
Method for finding the nodes associated to a face.
Definition ncmesh.cpp:5748
static PointMatrix pm_quad_identity
Definition ncmesh.hpp:1178
void ReferenceElement(int elem)
Add references to all nodes, edges and faces of the element.
Definition ncmesh.cpp:367
void UpdateElementToVertexTable()
Definition ncmesh.hpp:1051
BlockArray< Element > elements
Definition ncmesh.hpp:624
Array< int > coarse_elements
state of leaf_elements before Refine(), set by MarkCoarseLevel()
Definition ncmesh.hpp:1198
const CoarseFineTransformations & GetDerefinementTransforms() const
Definition ncmesh.cpp:5204
Table derefinements
possible derefinements, see GetDerefinementTable
Definition ncmesh.hpp:796
int FindNodeExt(const Element &el, int node, bool abort=true)
Extended version of find_node: works if 'el' is refined.
Definition ncmesh.cpp:3286
int ReorderFacePointMat(int v0, int v1, int v2, int v3, int elem, const PointMatrix &pm, PointMatrix &reordered) const
Definition ncmesh.cpp:3421
int FindMidEdgeNode(int node1, int node2) const
Definition ncmesh.cpp:336
friend struct PointMatrixHash
Definition ncmesh.hpp:1350
Array< char > face_geom
face geometry by face index, set by OnMeshUpdated
Definition ncmesh.hpp:732
int GetEdgeNCOrientation(const MeshId &edge_id) const
Definition ncmesh.cpp:5607
Table element_vertex
leaf-element to vertex table, see FindSetNeighbors
Definition ncmesh.hpp:734
Array< Triple< int, int, int > > reparents
scheduled node reparents (tmp)
Definition ncmesh.hpp:793
virtual void BuildFaceList()
Definition ncmesh.cpp:3675
void CheckAnisoPrism(int vn1, int vn2, int vn3, int vn4, const Refinement *refs, int nref)
Definition ncmesh.cpp:953
void LoadVertexParents(std::istream &input)
Load the vertex parent hierarchy from a mesh file.
Definition ncmesh.cpp:6090
TmpVertex * tmp_vertex
Definition ncmesh.hpp:1211
static void GridSfcOrdering3D(int width, int height, int depth, Array< int > &coords)
Definition ncmesh.cpp:5556
Array< int > leaf_elements
finest elements, in Mesh ordering (+ ghosts)
Definition ncmesh.hpp:723
const real_t * CalcVertexPos(int node) const
Definition ncmesh.cpp:2733
void CollectLeafElements(int elem, int state, Array< int > &ghosts, int &counter)
Definition ncmesh.cpp:2381
int CountTopLevelNodes() const
Return the index of the last top-level node plus one.
Definition ncmesh.cpp:6378
void QuadFaceSplitLevel(int vn1, int vn2, int vn3, int vn4, int &h_level, int &v_level) const
Computes the number of horizontal and vertical splits of this quad that have occurred in the NCMesh....
Definition ncmesh.cpp:5884
int Geoms
bit mask of element geometries present, see InitGeomFlags()
Definition ncmesh.hpp:527
void InitDerefTransforms()
Definition ncmesh.cpp:2343
int GetElementSizeReduction(int i) const
Definition ncmesh.cpp:5708
virtual void LimitNCLevel(int max_nc_level)
Definition ncmesh.cpp:6039
int NGhostVertices
Definition ncmesh.hpp:721
bool ZeroRootStates() const
Return true if all root_states are zero.
Definition ncmesh.cpp:6235
Array< int > vertex_nodeId
vertex-index to node-id map, see UpdateVertices
Definition ncmesh.hpp:725
int ParentFaceNodes(std::array< int, 4 > &nodes) const
Given a set of nodes defining a face, traverse the nodes structure to find the nodes that make up the...
Definition ncmesh.cpp:3106
int RetrieveNode(const Element &el, int index)
Return el.node[index] correctly, even if the element is refined.
Definition ncmesh.cpp:1994
void FindEdgeElements(int vn1, int vn2, int vn3, int vn4, Array< MeshId > &prisms) const
Definition ncmesh.cpp:905
virtual ~NCMesh()
Definition ncmesh.cpp:280
void CollectDerefinements(int elem, Array< Connection > &list)
Definition ncmesh.cpp:2230
virtual void ElementSharesEdge(int elem, int local, int enode)
Definition ncmesh.hpp:1013
int GetFaceVerticesEdges(const MeshId &face_id, int vert_index[4], int edge_index[4], int edge_orientation[4]) const
Definition ncmesh.cpp:5619
void RefineElement(int elem, char ref_type)
Refine the element elem with the refinement type ref_type.
Definition ncmesh.cpp:1121
void InitRootState(int root_count)
Definition ncmesh.cpp:2652
int NewQuadrilateral(int n0, int n1, int n2, int n3, int attr, int eattr0, int eattr1, int eattr2, int eattr3)
Definition ncmesh.cpp:732
Array< int > leaf_sfc_index
natural tree ordering of leaf elements
Definition ncmesh.hpp:724
void ForceRefinement(int vn1, int vn2, int vn3, int vn4)
Definition ncmesh.cpp:824
void CollectEdgeVertices(int v0, int v1, Array< int > &indices)
Definition ncmesh.cpp:4058
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
Definition ncmesh.hpp:728
static PointMatrix pm_pyramid_identity
Definition ncmesh.hpp:1181
Array< int > root_state
Definition ncmesh.hpp:701
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Definition ncmesh.hpp:319
real_t GetScale(real_t s, bool reverse) const
Return directed scale in (0,1).
Definition ncmesh.hpp:1221
void CollectTriFaceVertices(int v0, int v1, int v2, Array< int > &indices)
Definition ncmesh.cpp:4070
void InitRootElements()
Count root elements and initialize root_state.
Definition ncmesh.cpp:6335
int GetMidEdgeNode(int node1, int node2)
Definition ncmesh.cpp:352
virtual void BuildEdgeList()
Definition ncmesh.cpp:3800
void DebugDump(std::ostream &out) const
Definition ncmesh.cpp:6958
void DebugLeafOrder(std::ostream &out) const
Definition ncmesh.cpp:6933
std::map< std::string, int > RefPathMap
Definition ncmesh.hpp:1189
int GetElementDepth(int i) const
Return the distance of leaf i from the root.
Definition ncmesh.cpp:5696
void TraverseRefinements(int elem, int coarse_index, std::string &ref_path, RefPathMap &map) const
Definition ncmesh.cpp:5119
void CheckIsoFace(int vn1, int vn2, int vn3, int vn4, int en1, int en2, int en3, int en4, int midf)
Definition ncmesh.cpp:1099
void GetLimitRefinements(Array< Refinement > &refinements, int max_level)
Definition ncmesh.cpp:6009
int NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr, int fattr0, int fattr1, int fattr2, int fattr3, int fattr4)
Definition ncmesh.cpp:702
virtual void Derefine(const Array< int > &derefs)
Definition ncmesh.cpp:2307
bool Iso
true if the mesh only contains isotropic refinements
Definition ncmesh.hpp:526
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges, Array< int > &bdr_faces)
Get a list of vertices (2D/3D), edges (3D) and faces (3D) that coincide with boundary elements with t...
Definition ncmesh.cpp:5776
bool TriFaceSplit(int n1, int n2, int n3, int mid[3]=NULL) const
Given a tri face defined by three vertices, establish whether the edges that make up this face have b...
Definition ncmesh.cpp:3082
bool Legacy
true if the mesh was loaded from the legacy v1.1 format
Definition ncmesh.hpp:528
virtual void BuildVertexList()
Definition ncmesh.cpp:3912
static const PointMatrix & GetGeomIdentity(Geometry::Type geom)
Definition ncmesh.cpp:4606
static constexpr int MaxElemFaces
Number of faces an element can have.
Definition ncmesh.hpp:492
void CheckAnisoFace(int vn1, int vn2, int vn3, int vn4, int mid12, int mid34, int level=0)
Definition ncmesh.cpp:978
int NewWedge(int n0, int n1, int n2, int n3, int n4, int n5, int attr, int fattr0, int fattr1, int fattr2, int fattr3, int fattr4)
Definition ncmesh.cpp:645
int GetVertexRootCoord(int elem, RefCoord coord[3]) const
Definition ncmesh.cpp:4462
static void CheckSupportedGeom(Geometry::Type geom)
Definition ncmesh.hpp:1283
Array< real_t > coordinates
Definition ncmesh.hpp:705
void FindVertexCousins(int elem, int local, Array< int > &cousins) const
Definition ncmesh.cpp:4538
bool IsGhost(const Element &el) const
Return true if the Element el is a ghost element.
Definition ncmesh.hpp:787
static constexpr int MaxElemChildren
Number of children an element can have.
Definition ncmesh.hpp:494
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition ncmesh.hpp:326
bool TriFaceIsMaster(int n1, int n2, int n3) const
Determine if a Triangle face is a master face.
Definition ncmesh.hpp:910
int MyRank
used in parallel, or when loading a parallel file in serial
Definition ncmesh.hpp:525
void Print(std::ostream &out, const std::string &comments="") const
Definition ncmesh.cpp:6244
int AddElement(const Element &el)
Add an Element el to the NCMesh, optimized to reuse freed elements.
Definition ncmesh.hpp:812
void CountSplits(int elem, int splits[3]) const
Definition ncmesh.cpp:5919
int spaceDim
dimensions of the elements and the vertex coordinates
Definition ncmesh.hpp:524
void UnreferenceElement(int elem, Array< int > &elemFaces)
Definition ncmesh.cpp:398
void LegacyToNewVertexOrdering(Array< int > &order) const
I/O: Return a map from old (v1.1) vertex indices to new vertex indices.
Definition ncmesh.cpp:6814
int NGhostEdges
Definition ncmesh.hpp:721
void CollectIncidentElements(int elem, const RefCoord coord[3], Array< int > &list) const
Definition ncmesh.cpp:4515
void TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1, MatrixMap &matrix_map)
Definition ncmesh.cpp:3572
int QuadFaceSplitType(int n1, int n2, int n3, int n4, real_t &s, int mid[5]=NULL) const
Given a quad face defined by four vertices, establish which edges of this face have been split,...
Definition ncmesh.cpp:3028
void RegisterFaces(int elem, int *fattr=NULL)
Definition ncmesh.cpp:441
int PrintBoundary(std::ostream *out) const
Definition ncmesh.cpp:6117
CoarseFineTransformations transforms
storage for data returned by Get[De]RefinementTransforms()
Definition ncmesh.hpp:1195
void MarkCoarseLevel()
Definition ncmesh.cpp:5105
int TriFaceSplitLevel(int vn1, int vn2, int vn3) const
Return the number of splits of this triangle that have occurred in the NCMesh. If zero,...
Definition ncmesh.cpp:5869
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
Definition ncmesh.cpp:2278
void UpdateVertices()
This method assigns indices to vertices (Node::vert_index) that will be seen by the Mesh class and th...
Definition ncmesh.cpp:2476
NCList vertex_list
lazy-initialized list of vertices, see GetVertexList
Definition ncmesh.hpp:729
void InitGeomFlags()
Definition ncmesh.cpp:259
int NewTetrahedron(int n0, int n1, int n2, int n3, int attr, int fattr0, int fattr1, int fattr2, int fattr3)
Definition ncmesh.cpp:677
void LoadCoordinates(std::istream &input)
Load the "coordinates" section of the mesh file.
Definition ncmesh.cpp:6214
void GetPointMatrix(Geometry::Type geom, const char *ref_path, DenseMatrix &matrix) const
Definition ncmesh.cpp:4623
void TraverseQuadFace(int vn0, int vn1, int vn2, int vn3, const PointMatrix &pm, int level, Face *eface[4], MatrixMap &matrix_map)
Definition ncmesh.cpp:3452
const Table & GetDerefinementTable()
Definition ncmesh.cpp:2263
void ClearTransforms()
Free all internal data created by the above three functions.
Definition ncmesh.cpp:5284
void FindSetNeighbors(const Array< char > &elem_set, Array< int > *neighbors, Array< char > *neighbor_set=NULL)
Definition ncmesh.cpp:4209
void CollectQuadFaceVertices(int v0, int v1, int v2, int v3, Array< int > &indices)
Definition ncmesh.cpp:4094
void PrintCoordinates(std::ostream &out) const
Print the "coordinates" section of the mesh file.
Definition ncmesh.cpp:6196
TriFaceTraverseResults TraverseTriFace(int vn0, int vn1, int vn2, const PointMatrix &pm, int level, MatrixMap &matrix_map)
Definition ncmesh.cpp:3609
int GetMidFaceNode(int en1, int en2, int en3, int en4)
Definition ncmesh.cpp:359
bool HavePrisms() const
Return true if the mesh contains prism elements.
Definition ncmesh.hpp:778
long MemoryUsage() const
Return total number of bytes allocated.
Definition ncmesh.cpp:6882
void DerefineElement(int elem)
Derefine the element elem, does nothing on leaf elements.
Definition ncmesh.cpp:2036
NCList face_list
lazy-initialized list of faces, see GetFaceList
Definition ncmesh.hpp:727
virtual void Refine(const Array< Refinement > &refinements)
Definition ncmesh.cpp:1945
void ReparentNode(int node, int new_p1, int new_p2, real_t scale)
Definition ncmesh.cpp:319
static PointMatrix pm_tri_identity
Definition ncmesh.hpp:1177
int NGhostFaces
Definition ncmesh.hpp:721
static PointMatrix pm_hex_identity
Definition ncmesh.hpp:1182
void UpdateLeafElements()
Update the leaf elements indices in leaf_elements.
Definition ncmesh.cpp:2449
void LoadLegacyFormat(std::istream &input, int &curved, int &is_nc)
Load the deprecated MFEM mesh v1.1 format for backward compatibility.
Definition ncmesh.cpp:6656
int EdgeSplitLevel(int vn1, int vn2) const
Return the number of splits of this edge that have occurred in the NCMesh. If zero,...
Definition ncmesh.cpp:5862
static void GridSfcOrdering2D(int width, int height, Array< int > &coords)
Definition ncmesh.cpp:5541
static int find_local_face(int geom, int a, int b, int c)
Definition ncmesh.cpp:3314
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
A parallel extension of the NCMesh class.
Definition pncmesh.hpp:63
Data type Pyramid element.
Definition pyramid.hpp:23
Data type quadrilateral element.
Data type line segment element.
Definition segment.hpp:23
int RowSize(int i) const
Definition table.hpp:116
void Clear()
Definition table.cpp:453
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition table.cpp:259
void SetIJ(int *newI, int *newJ, int newsize=-1)
Replace the I and J arrays with the given newI and newJ arrays.
Definition table.cpp:279
int Size() const
Returns the number of TYPE I elements.
Definition table.hpp:100
std::size_t MemoryUsage() const
Definition table.cpp:474
void MakeFromList(int nrows, const Array< Connection > &list)
Definition table.cpp:352
Data type tetrahedron element.
Data type triangle element.
Definition triangle.hpp:24
A triple of objects.
Data type Wedge element.
Definition wedge.hpp:23
int dim
Definition ex24.cpp:53
int index(int i, int j, int nx, int ny)
Definition life.cpp:236
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
bool PrismFaceBottom(int node, int *n)
Definition ncmesh.cpp:817
void Swap(Array< T > &, Array< T > &)
Definition array.hpp:672
bool PrismFaceTop(int node, int *n)
Definition ncmesh.cpp:820
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
bool CubeFaceBottom(int node, int *n)
Definition ncmesh.cpp:811
bool contains_node(const std::array< int, 4 > &nodes, int n)
Definition ncmesh.cpp:3101
bool CubeFaceBack(int node, int *n)
Definition ncmesh.cpp:808
bool CubeFaceLeft(int node, int *n)
Definition ncmesh.cpp:799
bool CubeFaceRight(int node, int *n)
Definition ncmesh.cpp:802
bool CubeFaceFront(int node, int *n)
Definition ncmesh.cpp:805
NCMesh::RefCoord RefCoord
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
Definition text.hpp:31
bool CubeFaceTop(int node, int *n)
Definition ncmesh.cpp:814
Defines the coarse-fine transformations of all fine elements.
Definition ncmesh.hpp:90
Array< Embedding > embeddings
Fine element positions in their parents.
Definition ncmesh.hpp:92
void MakeCoarseToFineTable(Table &coarse_to_fine, bool want_ghosts=false) const
Definition ncmesh.cpp:5262
DenseTensor point_matrices[Geometry::NumGeom]
Definition ncmesh.hpp:96
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition table.hpp:28
Defines the position of a fine element within a coarse element.
Definition ncmesh.hpp:69
unsigned geom
Definition ncmesh.hpp:77
unsigned ghost
For internal use: 0 if regular fine element, 1 if parallel ghost element.
Definition ncmesh.hpp:81
int parent
Coarse Element index in the coarse mesh.
Definition ncmesh.hpp:71
unsigned matrix
Definition ncmesh.hpp:78
int rank
processor number (ParNCMesh), -1 if undefined/unknown
Definition ncmesh.hpp:602
bool IsLeaf() const
Definition ncmesh.hpp:613
int child[MaxElemChildren]
2-10 children (if ref_type != 0)
Definition ncmesh.hpp:607
char tet_type
tetrahedron split type, currently always 0
Definition ncmesh.hpp:599
Element(Geometry::Type geom, int attr)
Definition ncmesh.cpp:602
int node[MaxElemNodes]
element corners (if ref_type == 0)
Definition ncmesh.hpp:606
char ref_type
bit mask of X,Y,Z refinements (bits 0,1,2 respectively)
Definition ncmesh.hpp:598
char geom
Geometry::Type of the element (char for storage only)
Definition ncmesh.hpp:597
int index
element number in the Mesh, -1 if refined
Definition ncmesh.hpp:601
int parent
parent element, -1 if this is a root element, -2 if free'd
Definition ncmesh.hpp:609
Geometry::Type Geom() const
Definition ncmesh.hpp:612
void ForgetElement(int e)
Definition ncmesh.cpp:473
int elem[2]
up to 2 elements sharing the face
Definition ncmesh.hpp:576
void RegisterElement(int e)
Definition ncmesh.cpp:466
bool Boundary() const
Definition ncmesh.hpp:580
int index
face number in the Mesh
Definition ncmesh.hpp:575
int GetSingleElement() const
Return one of elem[0] or elem[1] and make sure the other is -1.
Definition ncmesh.cpp:488
int attribute
boundary element attribute, -1 if internal face
Definition ncmesh.hpp:574
This holds in one place the constants about the geometries we support.
Definition ncmesh.hpp:1328
int nfv[MaxElemFaces]
Definition ncmesh.hpp:1332
int faces[MaxElemFaces][4]
Definition ncmesh.hpp:1331
int edges[MaxElemEdges][2]
Definition ncmesh.hpp:1330
void InitGeom(Geometry::Type geom)
Definition ncmesh.cpp:57
Identifies a vertex/edge/face in both Mesh and NCMesh.
Definition ncmesh.hpp:210
int element
NCMesh::Element containing this vertex/edge/face.
Definition ncmesh.hpp:212
int index
Mesh number.
Definition ncmesh.hpp:211
Geometry::Type Geom() const
Definition ncmesh.hpp:216
signed char local
local number within 'element'
Definition ncmesh.hpp:213
MeshIdType GetMeshIdType(int index) const
Return a face type for a given nc index.
Definition ncmesh.cpp:4008
Array< MeshId > conforming
All MeshIds corresponding to conformal faces.
Definition ncmesh.hpp:252
long MemoryUsage() const
Definition ncmesh.cpp:6854
Array< Slave > slaves
All MeshIds corresponding to slave faces.
Definition ncmesh.hpp:254
bool CheckMeshIdType(int index, MeshIdType type) const
Given an index, check if this is a certain face type.
Definition ncmesh.cpp:4016
void Clear()
Erase the contents of the conforming, master and slave arrays.
Definition ncmesh.cpp:3969
void OrientedPointMatrix(const Slave &slave, DenseMatrix &oriented_matrix) const
Definition ncmesh.cpp:3947
Array< Master > masters
All MeshIds corresponding to master faces.
Definition ncmesh.hpp:253
Array< DenseMatrix * > point_matrices[Geometry::NumGeom]
List of unique point matrices for each slave geometry.
Definition ncmesh.hpp:257
MeshIdAndType GetMeshIdAndType(int index) const
Return a mesh id and type for a given nc index.
Definition ncmesh.cpp:3988
bool HasEdge() const
Definition ncmesh.hpp:549
real_t GetScale() const
Definition ncmesh.hpp:555
void SetScale(real_t s, bool overwrite=false)
Definition ncmesh.cpp:298
The PointMatrix stores the coordinates of the slave face using the master face coordinate as referenc...
Definition ncmesh.hpp:1130
void GetMatrix(DenseMatrix &point_matrix) const
Definition ncmesh.cpp:4569
Point points[MaxElemNodes]
Definition ncmesh.hpp:1132
bool operator==(const PointMatrix &pm) const
Definition ncmesh.cpp:4555
Nonconforming edge/face within a bigger edge/face.
Definition ncmesh.hpp:237
unsigned edge_flags
orientation flags, see OrientedPointMatrix
Definition ncmesh.hpp:240
unsigned matrix
index into NCList::point_matrices[geom]
Definition ncmesh.hpp:239
Refinement()=default
Refinement scale in each dimension.
real_t s[3]
Definition ncmesh.hpp:45
void SetScaleForType(const real_t *scale)
Set the scale in the directions for the currently set type.
Definition ncmesh.cpp:540
std::pair< char, real_t > ScaledType
Definition ncmesh.hpp:44
int index
Mesh element number.
Definition ncmesh.hpp:42
void Set(int element, char type, real_t scale=0.5)
Set the element, type, and scale.
Definition ncmesh.cpp:589
char GetType() const
Return the type as char.
Definition ncmesh.cpp:580
void SetType(char type, real_t scale=0.5)
Set the type and scale, assuming the element is already set.
Definition ncmesh.cpp:596
std::array< int, NCMesh::MaxFaceNodes > nodes