MFEM  v4.5.2
Finite element discretization library
ncmesh.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "mesh_headers.hpp"
13 #include "../general/sort_pairs.hpp"
14 #include "../general/text.hpp"
15 
16 #include <string>
17 #include <cmath>
18 #include <map>
19 
20 #include "ncmesh_tables.hpp"
21 
22 namespace mfem
23 {
24 
25 NCMesh::GeomInfo NCMesh::GI[Geometry::NumGeom];
26 
28 {
29  if (initialized) { return; }
30 
31  mfem::Element *elem = NULL;
32  switch (geom)
33  {
34  case Geometry::CUBE: elem = new Hexahedron; break;
35  case Geometry::PRISM: elem = new Wedge; break;
36  case Geometry::TETRAHEDRON: elem = new Tetrahedron; break;
37  case Geometry::PYRAMID: elem = new Pyramid; break;
38  case Geometry::SQUARE: elem = new Quadrilateral; break;
39  case Geometry::TRIANGLE: elem = new Triangle; break;
40  case Geometry::SEGMENT: elem = new Segment; break;
41  default: MFEM_ABORT("unsupported geometry " << geom);
42  }
43 
44  nv = elem->GetNVertices();
45  ne = elem->GetNEdges();
46  nf = elem->GetNFaces();
47 
48  for (int i = 0; i < ne; i++)
49  {
50  for (int j = 0; j < 2; j++)
51  {
52  edges[i][j] = elem->GetEdgeVertices(i)[j];
53  }
54  }
55  for (int i = 0; i < nf; i++)
56  {
57  nfv[i] = elem->GetNFaceVertices(i);
58 
59  faces[i][3] = 7; // invalid node index for 3-node faces
60  for (int j = 0; j < nfv[i]; j++)
61  {
62  faces[i][j] = elem->GetFaceVertices(i)[j];
63  }
64  }
65 
66  // in 1D/2D we pretend to have faces too, so we can use NCMesh::Face::elem[2]
67  if (!nf)
68  {
69  if (ne)
70  {
71  for (int i = 0; i < ne; i++)
72  {
73  // make a degenerate face
74  faces[i][0] = faces[i][1] = edges[i][0];
75  faces[i][2] = faces[i][3] = edges[i][1];
76  nfv[i] = 2;
77  }
78  nf = ne;
79  }
80  else
81  {
82  for (int i = 0; i < nv; i++)
83  {
84  // 1D degenerate face
85  faces[i][0] = faces[i][1] = faces[i][2] = faces[i][3] = i;
86  nfv[i] = 1;
87  }
88  nf = nv;
89  }
90  }
91 
92  delete elem;
93  initialized = true;
94 }
95 
96 static void CheckSupportedGeom(Geometry::Type geom)
97 {
98  MFEM_VERIFY(geom == Geometry::SEGMENT ||
99  geom == Geometry::TRIANGLE || geom == Geometry::SQUARE ||
100  geom == Geometry::CUBE || geom == Geometry::PRISM ||
101  geom == Geometry::PYRAMID || geom == Geometry::TETRAHEDRON,
102  "Element type " << geom << " is not supported by NCMesh.");
103 }
104 
105 NCMesh::NCMesh(const Mesh *mesh)
106  : shadow(1024, 2048)
107 {
108  Dim = mesh->Dimension();
109  spaceDim = mesh->SpaceDimension();
110  MyRank = 0;
111  Iso = true;
112  Legacy = false;
113 
114  // create the NCMesh::Element struct for each Mesh element
115  for (int i = 0; i < mesh->GetNE(); i++)
116  {
117  const mfem::Element *elem = mesh->GetElement(i);
118 
119  Geometry::Type geom = elem->GetGeometryType();
120  CheckSupportedGeom(geom);
121  GI[geom].InitGeom(geom);
122 
123  // if we have pyramids we will need tets after refinement
124  if (geom == Geometry::PYRAMID)
125  {
127  }
128 
129  // create NCMesh::Element for this mfem::Element
130  int root_id = AddElement(Element(geom, elem->GetAttribute()));
131  MFEM_ASSERT(root_id == i, "");
132  Element &root_elem = elements[root_id];
133 
134  const int *v = elem->GetVertices();
135  for (int j = 0; j < GI[geom].nv; j++)
136  {
137  int id = v[j];
138  root_elem.node[j] = id;
139  nodes.Alloc(id, id, id);
140  // NOTE: top-level nodes are special: id == p1 == p2 == orig. vertex id
141  }
142  }
143 
144  // if the user initialized any hanging nodes with Mesh::AddVertexParents,
145  // copy the hierarchy now
146  if (mesh->tmp_vertex_parents.Size())
147  {
148  for (const auto &triple : mesh->tmp_vertex_parents)
149  {
150  nodes.Reparent(triple.one, triple.two, triple.three);
151  }
152  }
153 
154  // create edge nodes and faces
155  nodes.UpdateUnused();
156  for (int i = 0; i < elements.Size(); i++)
157  {
158  // increase reference count of all nodes the element is using
159  // (NOTE: this will also create and reference all edge nodes and faces)
160  ReferenceElement(i);
161 
162  // make links from faces back to the element
163  RegisterFaces(i);
164  }
165 
166  // store boundary element attributes
167  for (int i = 0; i < mesh->GetNBE(); i++)
168  {
169  const mfem::Element *be = mesh->GetBdrElement(i);
170  const int *v = be->GetVertices();
171 
172  Face* face = NULL;
173  switch (be->GetType())
174  {
176  face = faces.Find(v[0], v[1], v[2], v[3]);
177  break;
179  face = faces.Find(v[0], v[1], v[2]);
180  break;
182  face = faces.Find(v[0], v[0], v[1], v[1]);
183  break;
185  face = faces.Find(v[0], v[0], v[0], v[0]);
186  break;
187  default:
188  MFEM_ABORT("Unsupported boundary element geometry.");
189  }
190 
191  MFEM_VERIFY(face, "Boundary face not found.");
192  face->attribute = be->GetAttribute();
193  }
194 
195  // copy top-level vertex coordinates (leave empty if the mesh is curved)
196  if (!mesh->Nodes)
197  {
198  coordinates.SetSize(3*mesh->GetNV());
199  for (int i = 0; i < mesh->GetNV(); i++)
200  {
201  std::memcpy(&coordinates[3*i], mesh->GetVertex(i), 3*sizeof(double));
202  }
203  }
204 
205  InitRootState(mesh->GetNE());
206  InitGeomFlags();
207 
208  Update();
209 }
210 
211 NCMesh::NCMesh(const NCMesh &other)
212  : Dim(other.Dim)
213  , spaceDim(other.spaceDim)
214  , MyRank(other.MyRank)
215  , Iso(other.Iso)
216  , Geoms(other.Geoms)
217  , Legacy(other.Legacy)
218  , nodes(other.nodes)
219  , faces(other.faces)
220  , elements(other.elements)
221  , shadow(1024, 2048)
222 {
224  other.root_state.Copy(root_state);
226  Update();
227 }
228 
230 {
231  Geoms = 0;
232  for (int i = 0; i < root_state.Size(); i++)
233  {
234  Geoms |= (1 << elements[i].Geom());
235  }
236 }
237 
239 {
241  UpdateVertices();
242 
243  vertex_list.Clear();
244  face_list.Clear();
245  edge_list.Clear();
246 
248 }
249 
251 {
252 #ifdef MFEM_DEBUG
253  // sign off of all faces and nodes
254  Array<int> elemFaces;
255  for (int i = 0; i < elements.Size(); i++)
256  {
257  if (elements[i].IsLeaf())
258  {
259  elemFaces.SetSize(0);
260  UnreferenceElement(i, elemFaces);
261  DeleteUnusedFaces(elemFaces);
262  }
263  }
264  // NOTE: in release mode, we just throw away all faces and nodes at once
265 #endif
266 }
267 
269 {
270  MFEM_ASSERT(!vert_refc && !edge_refc, "node was not unreferenced properly, "
271  "vert_refc: " << (int) vert_refc << ", edge_refc: "
272  << (int) edge_refc);
273 }
274 
275 void NCMesh::ReparentNode(int node, int new_p1, int new_p2)
276 {
277  Node &nd = nodes[node];
278  int old_p1 = nd.p1, old_p2 = nd.p2;
279 
280  // assign new parents
281  nodes.Reparent(node, new_p1, new_p2);
282 
283  MFEM_ASSERT(shadow.FindId(old_p1, old_p2) < 0,
284  "shadow node already exists");
285 
286  // store old parent pair temporarily in 'shadow'
287  int sh = shadow.GetId(old_p1, old_p2);
288  shadow[sh].vert_index = node;
289 }
290 
291 int NCMesh::FindMidEdgeNode(int node1, int node2) const
292 {
293  int mid = nodes.FindId(node1, node2);
294  if (mid < 0 && shadow.Size())
295  {
296  // if (anisotropic) refinement is underway, some nodes may temporarily
297  // be available under alternate parents (see ReparentNode)
298  mid = shadow.FindId(node1, node2);
299  if (mid >= 0)
300  {
301  mid = shadow[mid].vert_index; // index of the original node
302  }
303  }
304  return mid;
305 }
306 
307 int NCMesh::GetMidEdgeNode(int node1, int node2)
308 {
309  int mid = FindMidEdgeNode(node1, node2);
310  if (mid < 0) { mid = nodes.GetId(node1, node2); } // create if not found
311  return mid;
312 }
313 
314 int NCMesh::GetMidFaceNode(int en1, int en2, int en3, int en4)
315 {
316  // mid-face node can be created either from (en1, en3) or from (en2, en4)
317  int midf = FindMidEdgeNode(en1, en3);
318  if (midf >= 0) { return midf; }
319  return nodes.GetId(en2, en4);
320 }
321 
323 {
324  Element &el = elements[elem];
325  int* node = el.node;
326  GeomInfo& gi = GI[el.Geom()];
327 
328  // reference all vertices
329  for (int i = 0; i < gi.nv; i++)
330  {
331  nodes[node[i]].vert_refc++;
332  }
333 
334  // reference all edges (possibly creating their nodes)
335  for (int i = 0; i < gi.ne; i++)
336  {
337  const int* ev = gi.edges[i];
338  nodes.Get(node[ev[0]], node[ev[1]])->edge_refc++;
339  }
340 
341  // get all faces (possibly creating them)
342  for (int i = 0; i < gi.nf; i++)
343  {
344  const int* fv = gi.faces[i];
345  faces.GetId(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
346 
347  // NOTE: face->RegisterElement called separately to avoid having
348  // to store 3 element indices temporarily in the face when refining.
349  // See also NCMesh::RegisterFaces.
350  }
351 }
352 
353 void NCMesh::UnreferenceElement(int elem, Array<int> &elemFaces)
354 {
355  Element &el = elements[elem];
356  int* node = el.node;
357  GeomInfo& gi = GI[el.Geom()];
358 
359  // unreference all faces
360  for (int i = 0; i < gi.nf; i++)
361  {
362  const int* fv = gi.faces[i];
363  int face = faces.FindId(node[fv[0]], node[fv[1]],
364  node[fv[2]], node[fv[3]]);
365  MFEM_ASSERT(face >= 0, "face not found.");
366  faces[face].ForgetElement(elem);
367 
368  // NOTE: faces.Delete() called later to avoid destroying and
369  // recreating faces during refinement, see NCMesh::DeleteUnusedFaces.
370  elemFaces.Append(face);
371  }
372 
373  // unreference all edges (possibly destroying them)
374  for (int i = 0; i < gi.ne; i++)
375  {
376  const int* ev = gi.edges[i];
377  int enode = FindMidEdgeNode(node[ev[0]], node[ev[1]]);
378  MFEM_ASSERT(enode >= 0, "edge not found.");
379  MFEM_ASSERT(nodes.IdExists(enode), "edge does not exist.");
380  if (!nodes[enode].UnrefEdge())
381  {
382  nodes.Delete(enode);
383  }
384  }
385 
386  // unreference all vertices (possibly destroying them)
387  for (int i = 0; i < gi.nv; i++)
388  {
389  if (!nodes[node[i]].UnrefVertex())
390  {
391  nodes.Delete(node[i]);
392  }
393  }
394 }
395 
396 void NCMesh::RegisterFaces(int elem, int* fattr)
397 {
398  Element &el = elements[elem];
399  GeomInfo &gi = GI[el.Geom()];
400 
401  for (int i = 0; i < gi.nf; i++)
402  {
403  Face* face = GetFace(el, i);
404  MFEM_ASSERT(face, "face not found.");
405  face->RegisterElement(elem);
406  if (fattr) { face->attribute = fattr[i]; }
407  }
408 }
409 
410 void NCMesh::DeleteUnusedFaces(const Array<int> &elemFaces)
411 {
412  for (int i = 0; i < elemFaces.Size(); i++)
413  {
414  if (faces[elemFaces[i]].Unused())
415  {
416  faces.Delete(elemFaces[i]);
417  }
418  }
419 }
420 
422 {
423  if (elem[0] < 0) { elem[0] = e; }
424  else if (elem[1] < 0) { elem[1] = e; }
425  else { MFEM_ABORT("can't have 3 elements in Face::elem[]."); }
426 }
427 
429 {
430  if (elem[0] == e) { elem[0] = -1; }
431  else if (elem[1] == e) { elem[1] = -1; }
432  else { MFEM_ABORT("element " << e << " not found in Face::elem[]."); }
433 }
434 
436 {
437  GeomInfo& gi = GI[(int) elem.geom];
438  const int* fv = gi.faces[face_no];
439  int* node = elem.node;
440  return faces.Find(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
441 }
442 
444 {
445  if (elem[0] >= 0)
446  {
447  MFEM_ASSERT(elem[1] < 0, "not a single element face.");
448  return elem[0];
449  }
450  else
451  {
452  MFEM_ASSERT(elem[1] >= 0, "no elements in face.");
453  return elem[1];
454  }
455 }
456 
457 
458 //// Refinement ////////////////////////////////////////////////////////////////
459 
461  : geom(geom), ref_type(0), tet_type(0), flag(0), index(-1)
462  , rank(0), attribute(attr), parent(-1)
463 {
464  for (int i = 0; i < MaxElemNodes; i++) { node[i] = -1; }
465  for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
466 
467  // NOTE: in 2D the 8/10-element node/child arrays are not optimal, however,
468  // testing shows we would only save 17% of the total NCMesh memory if
469  // 4-element arrays were used (e.g. through templates); we thus prefer to
470  // keep the code as simple as possible.
471 }
472 
473 int NCMesh::NewHexahedron(int n0, int n1, int n2, int n3,
474  int n4, int n5, int n6, int n7,
475  int attr,
476  int fattr0, int fattr1, int fattr2,
477  int fattr3, int fattr4, int fattr5)
478 {
479  // create new element, initialize nodes
480  int new_id = AddElement(Element(Geometry::CUBE, attr));
481  Element &el = elements[new_id];
482 
483  el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
484  el.node[4] = n4, el.node[5] = n5, el.node[6] = n6, el.node[7] = n7;
485 
486  // get faces and assign face attributes
487  Face* f[MaxElemFaces];
488  const GeomInfo &gi_hex = GI[Geometry::CUBE];
489  for (int i = 0; i < gi_hex.nf; i++)
490  {
491  const int* fv = gi_hex.faces[i];
492  f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
493  el.node[fv[2]], el.node[fv[3]]);
494  }
495 
496  f[0]->attribute = fattr0, f[1]->attribute = fattr1;
497  f[2]->attribute = fattr2, f[3]->attribute = fattr3;
498  f[4]->attribute = fattr4, f[5]->attribute = fattr5;
499 
500  return new_id;
501 }
502 
503 int NCMesh::NewWedge(int n0, int n1, int n2,
504  int n3, int n4, int n5,
505  int attr,
506  int fattr0, int fattr1,
507  int fattr2, int fattr3, int fattr4)
508 {
509  // create new element, initialize nodes
510  int new_id = AddElement(Element(Geometry::PRISM, attr));
511  Element &el = elements[new_id];
512 
513  el.node[0] = n0, el.node[1] = n1, el.node[2] = n2;
514  el.node[3] = n3, el.node[4] = n4, el.node[5] = n5;
515 
516  // get faces and assign face attributes
517  Face* f[5];
518  const GeomInfo &gi_wedge = GI[Geometry::PRISM];
519  for (int i = 0; i < gi_wedge.nf; i++)
520  {
521  const int* fv = gi_wedge.faces[i];
522  f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
523  el.node[fv[2]], el.node[fv[3]]);
524  }
525 
526  f[0]->attribute = fattr0;
527  f[1]->attribute = fattr1;
528  f[2]->attribute = fattr2;
529  f[3]->attribute = fattr3;
530  f[4]->attribute = fattr4;
531 
532  return new_id;
533 }
534 
535 int NCMesh::NewTetrahedron(int n0, int n1, int n2, int n3, int attr,
536  int fattr0, int fattr1, int fattr2, int fattr3)
537 {
538  // create new element, initialize nodes
539  int new_id = AddElement(Element(Geometry::TETRAHEDRON, attr));
540  Element &el = elements[new_id];
541 
542  el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
543 
544  // get faces and assign face attributes
545  Face* f[4];
546  const GeomInfo &gi_tet = GI[Geometry::TETRAHEDRON];
547  for (int i = 0; i < gi_tet.nf; i++)
548  {
549  const int* fv = gi_tet.faces[i];
550  f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]], el.node[fv[2]]);
551  }
552 
553  f[0]->attribute = fattr0;
554  f[1]->attribute = fattr1;
555  f[2]->attribute = fattr2;
556  f[3]->attribute = fattr3;
557 
558  return new_id;
559 }
560 int NCMesh::NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr,
561  int fattr0, int fattr1, int fattr2, int fattr3,
562  int fattr4)
563 {
564  // create new element, initialize nodes
565  int new_id = AddElement(Element(Geometry::PYRAMID, attr));
566  Element &el = elements[new_id];
567 
568  el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
569  el.node[4] = n4;
570 
571  // get faces and assign face attributes
572  Face* f[5];
573  const GeomInfo &gi_pyr = GI[Geometry::PYRAMID];
574  for (int i = 0; i < gi_pyr.nf; i++)
575  {
576  const int* fv = gi_pyr.faces[i];
577  f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
578  el.node[fv[2]], el.node[fv[3]]);
579  }
580 
581  f[0]->attribute = fattr0;
582  f[1]->attribute = fattr1;
583  f[2]->attribute = fattr2;
584  f[3]->attribute = fattr3;
585  f[4]->attribute = fattr4;
586 
587  return new_id;
588 }
589 
590 int NCMesh::NewQuadrilateral(int n0, int n1, int n2, int n3,
591  int attr,
592  int eattr0, int eattr1, int eattr2, int eattr3)
593 {
594  // create new element, initialize nodes
595  int new_id = AddElement(Element(Geometry::SQUARE, attr));
596  Element &el = elements[new_id];
597 
598  el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
599 
600  // get (degenerate) faces and assign face attributes
601  Face* f[4];
602  const GeomInfo &gi_quad = GI[Geometry::SQUARE];
603  for (int i = 0; i < gi_quad.nf; i++)
604  {
605  const int* fv = gi_quad.faces[i];
606  f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
607  el.node[fv[2]], el.node[fv[3]]);
608  }
609 
610  f[0]->attribute = eattr0, f[1]->attribute = eattr1;
611  f[2]->attribute = eattr2, f[3]->attribute = eattr3;
612 
613  return new_id;
614 }
615 
616 int NCMesh::NewTriangle(int n0, int n1, int n2,
617  int attr, int eattr0, int eattr1, int eattr2)
618 {
619  // create new element, initialize nodes
620  int new_id = AddElement(Element(Geometry::TRIANGLE, attr));
621  Element &el = elements[new_id];
622 
623  el.node[0] = n0, el.node[1] = n1, el.node[2] = n2;
624 
625  // get (degenerate) faces and assign face attributes
626  Face* f[3];
627  const GeomInfo &gi_tri = GI[Geometry::TRIANGLE];
628  for (int i = 0; i < gi_tri.nf; i++)
629  {
630  const int* fv = gi_tri.faces[i];
631  f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
632  el.node[fv[2]], el.node[fv[3]]);
633  }
634 
635  f[0]->attribute = eattr0;
636  f[1]->attribute = eattr1;
637  f[2]->attribute = eattr2;
638 
639  return new_id;
640 }
641 
642 int NCMesh::NewSegment(int n0, int n1, int attr, int vattr1, int vattr2)
643 {
644  // create new element, initialize nodes
645  int new_id = AddElement(Element(Geometry::SEGMENT, attr));
646  Element &el = elements[new_id];
647  el.node[0] = n0, el.node[1] = n1;
648 
649  // get (degenerate) faces and assign face attributes
650  int v0 = el.node[0], v1 = el.node[1];
651  faces.Get(v0, v0, v0, v0)->attribute = vattr1;
652  faces.Get(v1, v1, v1, v1)->attribute = vattr2;
653 
654  return new_id;
655 }
656 
657 inline bool CubeFaceLeft(int node, int* n)
658 { return node == n[0] || node == n[3] || node == n[4] || node == n[7]; }
659 
660 inline bool CubeFaceRight(int node, int* n)
661 { return node == n[1] || node == n[2] || node == n[5] || node == n[6]; }
662 
663 inline bool CubeFaceFront(int node, int* n)
664 { return node == n[0] || node == n[1] || node == n[4] || node == n[5]; }
665 
666 inline bool CubeFaceBack(int node, int* n)
667 { return node == n[2] || node == n[3] || node == n[6] || node == n[7]; }
668 
669 inline bool CubeFaceBottom(int node, int* n)
670 { return node == n[0] || node == n[1] || node == n[2] || node == n[3]; }
671 
672 inline bool CubeFaceTop(int node, int* n)
673 { return node == n[4] || node == n[5] || node == n[6] || node == n[7]; }
674 
675 inline bool PrismFaceBottom(int node, int* n)
676 { return node == n[0] || node == n[1] || node == n[2]; }
677 
678 inline bool PrismFaceTop(int node, int* n)
679 { return node == n[3] || node == n[4] || node == n[5]; }
680 
681 
682 void NCMesh::ForceRefinement(int vn1, int vn2, int vn3, int vn4)
683 {
684  // get the element this face belongs to
685  Face* face = faces.Find(vn1, vn2, vn3, vn4);
686  if (!face) { return; }
687 
688  int elem = face->GetSingleElement();
689  Element &el = elements[elem];
690  MFEM_ASSERT(!el.ref_type, "element already refined.");
691 
692  int* el_nodes = el.node;
693  if (el.Geom() == Geometry::CUBE)
694  {
695  // schedule the right split depending on face orientation
696  if ((CubeFaceLeft(vn1, el_nodes) && CubeFaceRight(vn2, el_nodes)) ||
697  (CubeFaceLeft(vn2, el_nodes) && CubeFaceRight(vn1, el_nodes)))
698  {
699  ref_stack.Append(Refinement(elem, 1)); // X split
700  }
701  else if ((CubeFaceFront(vn1, el_nodes) && CubeFaceBack(vn2, el_nodes)) ||
702  (CubeFaceFront(vn2, el_nodes) && CubeFaceBack(vn1, el_nodes)))
703  {
704  ref_stack.Append(Refinement(elem, 2)); // Y split
705  }
706  else if ((CubeFaceBottom(vn1, el_nodes) && CubeFaceTop(vn2, el_nodes)) ||
707  (CubeFaceBottom(vn2, el_nodes) && CubeFaceTop(vn1, el_nodes)))
708  {
709  ref_stack.Append(Refinement(elem, 4)); // Z split
710  }
711  else
712  {
713  MFEM_ABORT("Inconsistent element/face structure.");
714  }
715  }
716  else if (el.Geom() == Geometry::PRISM)
717  {
718  if ((PrismFaceTop(vn1, el_nodes) && PrismFaceBottom(vn4, el_nodes)) ||
719  (PrismFaceTop(vn4, el_nodes) && PrismFaceBottom(vn1, el_nodes)))
720  {
721  ref_stack.Append(Refinement(elem, 3)); // XY split
722  }
723  else if ((PrismFaceTop(vn1, el_nodes) && PrismFaceBottom(vn2, el_nodes)) ||
724  (PrismFaceTop(vn2, el_nodes) && PrismFaceBottom(vn1, el_nodes)))
725  {
726  ref_stack.Append(Refinement(elem, 4)); // Z split
727  }
728  else
729  {
730  MFEM_ABORT("Inconsistent element/face structure.");
731  }
732  }
733  else
734  {
735  MFEM_ABORT("Unsupported geometry.")
736  }
737 }
738 
739 
740 void NCMesh::FindEdgeElements(int vn1, int vn2, int vn3, int vn4,
741  Array<MeshId> &elem_edge) const
742 {
743  // Assuming that f = (vn1, vn2, vn3, vn4) is a quad face and
744  // e = (vn1, vn4) is its edge, this function finds the N elements
745  // sharing e, and returns the N different MeshIds of the edge (i.e.,
746  // different element-local pairs describing the edge).
747 
748  int ev1 = vn1, ev2 = vn4;
749 
750  // follow face refinement towards 'vn1', get an existing face
751  int split, mid[5];
752  while ((split = QuadFaceSplitType(vn1, vn2, vn3, vn4, mid)) > 0)
753  {
754  if (split == 1) // vertical
755  {
756  vn2 = mid[0]; vn3 = mid[2];
757  }
758  else // horizontal
759  {
760  vn3 = mid[1]; vn4 = mid[3];
761  }
762  }
763 
764  const Face *face = faces.Find(vn1, vn2, vn3, vn4);
765  MFEM_ASSERT(face != NULL, "Face not found: "
766  << vn1 << ", " << vn2 << ", " << vn3 << ", " << vn4
767  << " (edge " << ev1 << "-" << ev2 << ").");
768 
769  int elem = face->GetSingleElement();
770  int local = find_node(elements[elem], vn1);
771 
772  Array<int> cousins;
773  FindVertexCousins(elem, local, cousins);
774 
775  elem_edge.SetSize(0);
776  for (int i = 0; i < cousins.Size(); i++)
777  {
778  local = find_element_edge(elements[cousins[i]], ev1, ev2, false);
779  if (local > 0)
780  {
781  elem_edge.Append(MeshId(-1, cousins[i], local, Geometry::SEGMENT));
782  }
783  }
784 }
785 
786 
787 void NCMesh::CheckAnisoPrism(int vn1, int vn2, int vn3, int vn4,
788  const Refinement *refs, int nref)
789 {
790  MeshId buf[4];
791  Array<MeshId> eid(buf, 4);
792  FindEdgeElements(vn1, vn2, vn3, vn4, eid);
793 
794  // see if there is an element that has not been force-refined yet
795  for (int i = 0, j; i < eid.Size(); i++)
796  {
797  int elem = eid[i].element;
798  for (j = 0; j < nref; j++)
799  {
800  if (refs[j].index == elem) { break; }
801  }
802  if (j == nref) // elem not found in refs[]
803  {
804  // schedule prism refinement along Z axis
805  MFEM_ASSERT(elements[elem].Geom() == Geometry::PRISM, "");
806  ref_stack.Append(Refinement(elem, 4));
807  }
808  }
809 }
810 
811 
812 void NCMesh::CheckAnisoFace(int vn1, int vn2, int vn3, int vn4,
813  int mid12, int mid34, int level)
814 {
815  // When a face is getting split anisotropically (without loss of generality
816  // we assume a "vertical" split here, see picture), it is important to make
817  // sure that the mid-face vertex node (midf) has mid34 and mid12 as parents.
818  // This is necessary for the face traversal algorithm and at places like
819  // Refine() that assume the mid-edge nodes to be accessible through the right
820  // parents. However, midf may already exist under the parents mid41 and
821  // mid23. In that case we need to "reparent" midf, i.e., reinsert it to the
822  // hash-table under the correct parents. This doesn't affect other nodes as
823  // all IDs stay the same, only the face refinement "tree" is affected.
824  //
825  // vn4 mid34 vn3
826  // *------*------*
827  // | | |
828  // | |midf |
829  // mid41 *- - - *- - - * mid23
830  // | | |
831  // | | |
832  // *------*------*
833  // vn1 mid12 vn2
834  //
835  // This function is recursive, because the above applies to any node along
836  // the middle vertical edge. The function calls itself again for the bottom
837  // and upper half of the above picture.
838 
839  int mid23 = FindMidEdgeNode(vn2, vn3);
840  int mid41 = FindMidEdgeNode(vn4, vn1);
841  if (mid23 >= 0 && mid41 >= 0)
842  {
843  int midf = nodes.FindId(mid23, mid41);
844  if (midf >= 0)
845  {
846  reparents.Append(Triple<int, int, int>(midf, mid12, mid34));
847 
848  int rs = ref_stack.Size();
849 
850  CheckAnisoFace(vn1, vn2, mid23, mid41, mid12, midf, level+1);
851  CheckAnisoFace(mid41, mid23, vn3, vn4, midf, mid34, level+1);
852 
853  if (HavePrisms() && nodes[midf].HasEdge())
854  {
855  // Check if there is a prism with edge (mid23, mid41) that we may
856  // have missed in 'CheckAnisoFace', and force-refine it if present.
857 
858  if (ref_stack.Size() > rs)
859  {
860  CheckAnisoPrism(mid23, vn3, vn4, mid41,
861  &ref_stack[rs], ref_stack.Size() - rs);
862  }
863  else
864  {
865  CheckAnisoPrism(mid23, vn3, vn4, mid41, NULL, 0);
866  }
867  }
868 
869  // perform the reparents all at once at the end
870  if (level == 0)
871  {
872  for (int i = 0; i < reparents.Size(); i++)
873  {
874  const Triple<int, int, int> &tr = reparents[i];
875  ReparentNode(tr.one, tr.two, tr.three);
876  }
877  reparents.DeleteAll();
878  }
879  return;
880  }
881  }
882 
883  // Also, this is the place where forced refinements begin. In the picture
884  // above, edges mid12-midf and midf-mid34 should actually exist in the
885  // neighboring elements, otherwise the mesh is inconsistent and needs to be
886  // fixed. Example: suppose an element is being refined isotropically (!)
887  // whose neighbors across some face look like this:
888  //
889  // *--------*--------*
890  // | d | e |
891  // *--------*--------*
892  // | c |
893  // *--------*--------*
894  // | | |
895  // | a | b |
896  // | | |
897  // *--------*--------*
898  //
899  // Element 'c' needs to be refined vertically for the mesh to remain valid.
900 
901  if (level > 0)
902  {
903  ForceRefinement(vn1, vn2, vn3, vn4);
904  }
905 }
906 
907 void NCMesh::CheckIsoFace(int vn1, int vn2, int vn3, int vn4,
908  int en1, int en2, int en3, int en4, int midf)
909 {
910  if (!Iso)
911  {
912  /* If anisotropic refinements are present in the mesh, we need to check
913  isotropically split faces as well, see second comment in
914  CheckAnisoFace above. */
915 
916  CheckAnisoFace(vn1, vn2, en2, en4, en1, midf);
917  CheckAnisoFace(en4, en2, vn3, vn4, midf, en3);
918  CheckAnisoFace(vn4, vn1, en1, en3, en4, midf);
919  CheckAnisoFace(en3, en1, vn2, vn3, midf, en2);
920  }
921 }
922 
923 
924 void NCMesh::RefineElement(int elem, char ref_type)
925 {
926  if (!ref_type) { return; }
927 
928  // handle elements that may have been (force-) refined already
929  Element &el = elements[elem];
930  if (el.ref_type)
931  {
932  char remaining = ref_type & ~el.ref_type;
933 
934  // do the remaining splits on the children
935  for (int i = 0; i < MaxElemChildren; i++)
936  {
937  if (el.child[i] >= 0) { RefineElement(el.child[i], remaining); }
938  }
939  return;
940  }
941 
942  /*mfem::out << "Refining element " << elem << " ("
943  << el.node[0] << ", " << el.node[1] << ", "
944  << el.node[2] << ", " << el.node[3] << ", "
945  << el.node[4] << ", " << el.node[5] << ", "
946  << el.node[6] << ", " << el.node[7] << "), "
947  << "ref_type " << int(ref_type) << std::endl;*/
948 
949  int* no = el.node;
950  int attr = el.attribute;
951 
952  int child[MaxElemChildren];
953  for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
954 
955  // get parent's face attributes
956  int fa[MaxElemFaces];
957  GeomInfo& gi = GI[el.Geom()];
958  for (int i = 0; i < gi.nf; i++)
959  {
960  const int* fv = gi.faces[i];
961  Face* face = faces.Find(no[fv[0]], no[fv[1]], no[fv[2]], no[fv[3]]);
962  fa[i] = face->attribute;
963  }
964 
965  // create child elements
966  if (el.Geom() == Geometry::CUBE)
967  {
968  // Vertex numbering is assumed to be as follows:
969  //
970  // 7 6
971  // +-----------+ Faces: 0 bottom
972  // /| /| 1 front
973  // 4 / | 5 / | 2 right
974  // +-----------+ | 3 back
975  // | | | | 4 left
976  // | +--------|--+ 5 top
977  // | / 3 | / 2 Z Y
978  // |/ |/ |/
979  // +-----------+ *--X
980  // 0 1
981 
982  if (ref_type == Refinement::X) // split along X axis
983  {
984  int mid01 = GetMidEdgeNode(no[0], no[1]);
985  int mid23 = GetMidEdgeNode(no[2], no[3]);
986  int mid45 = GetMidEdgeNode(no[4], no[5]);
987  int mid67 = GetMidEdgeNode(no[6], no[7]);
988 
989  child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
990  no[4], mid45, mid67, no[7], attr,
991  fa[0], fa[1], -1, fa[3], fa[4], fa[5]);
992 
993  child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
994  mid45, no[5], no[6], mid67, attr,
995  fa[0], fa[1], fa[2], fa[3], -1, fa[5]);
996 
997  CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
998  CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
999  CheckAnisoFace(no[4], no[5], no[6], no[7], mid45, mid67);
1000  CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
1001  }
1002  else if (ref_type == Refinement::Y) // split along Y axis
1003  {
1004  int mid12 = GetMidEdgeNode(no[1], no[2]);
1005  int mid30 = GetMidEdgeNode(no[3], no[0]);
1006  int mid56 = GetMidEdgeNode(no[5], no[6]);
1007  int mid74 = GetMidEdgeNode(no[7], no[4]);
1008 
1009  child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
1010  no[4], no[5], mid56, mid74, attr,
1011  fa[0], fa[1], fa[2], -1, fa[4], fa[5]);
1012 
1013  child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
1014  mid74, mid56, no[6], no[7], attr,
1015  fa[0], -1, fa[2], fa[3], fa[4], fa[5]);
1016 
1017  CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
1018  CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
1019  CheckAnisoFace(no[5], no[6], no[7], no[4], mid56, mid74);
1020  CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
1021  }
1022  else if (ref_type == Refinement::Z) // split along Z axis
1023  {
1024  int mid04 = GetMidEdgeNode(no[0], no[4]);
1025  int mid15 = GetMidEdgeNode(no[1], no[5]);
1026  int mid26 = GetMidEdgeNode(no[2], no[6]);
1027  int mid37 = GetMidEdgeNode(no[3], no[7]);
1028 
1029  child[0] = NewHexahedron(no[0], no[1], no[2], no[3],
1030  mid04, mid15, mid26, mid37, attr,
1031  fa[0], fa[1], fa[2], fa[3], fa[4], -1);
1032 
1033  child[1] = NewHexahedron(mid04, mid15, mid26, mid37,
1034  no[4], no[5], no[6], no[7], attr,
1035  -1, fa[1], fa[2], fa[3], fa[4], fa[5]);
1036 
1037  CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
1038  CheckAnisoFace(no[5], no[1], no[2], no[6], mid15, mid26);
1039  CheckAnisoFace(no[6], no[2], no[3], no[7], mid26, mid37);
1040  CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
1041  }
1042  else if (ref_type == Refinement::XY) // XY split
1043  {
1044  int mid01 = GetMidEdgeNode(no[0], no[1]);
1045  int mid12 = GetMidEdgeNode(no[1], no[2]);
1046  int mid23 = GetMidEdgeNode(no[2], no[3]);
1047  int mid30 = GetMidEdgeNode(no[3], no[0]);
1048 
1049  int mid45 = GetMidEdgeNode(no[4], no[5]);
1050  int mid56 = GetMidEdgeNode(no[5], no[6]);
1051  int mid67 = GetMidEdgeNode(no[6], no[7]);
1052  int mid74 = GetMidEdgeNode(no[7], no[4]);
1053 
1054  int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid30);
1055  int midf5 = GetMidFaceNode(mid45, mid56, mid67, mid74);
1056 
1057  child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
1058  no[4], mid45, midf5, mid74, attr,
1059  fa[0], fa[1], -1, -1, fa[4], fa[5]);
1060 
1061  child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
1062  mid45, no[5], mid56, midf5, attr,
1063  fa[0], fa[1], fa[2], -1, -1, fa[5]);
1064 
1065  child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
1066  midf5, mid56, no[6], mid67, attr,
1067  fa[0], -1, fa[2], fa[3], -1, fa[5]);
1068 
1069  child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
1070  mid74, midf5, mid67, no[7], attr,
1071  fa[0], -1, -1, fa[3], fa[4], fa[5]);
1072 
1073  CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
1074  CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
1075  CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
1076  CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
1077 
1078  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1079  CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1080  }
1081  else if (ref_type == Refinement::XZ) // XZ split
1082  {
1083  int mid01 = GetMidEdgeNode(no[0], no[1]);
1084  int mid23 = GetMidEdgeNode(no[2], no[3]);
1085  int mid45 = GetMidEdgeNode(no[4], no[5]);
1086  int mid67 = GetMidEdgeNode(no[6], no[7]);
1087 
1088  int mid04 = GetMidEdgeNode(no[0], no[4]);
1089  int mid15 = GetMidEdgeNode(no[1], no[5]);
1090  int mid26 = GetMidEdgeNode(no[2], no[6]);
1091  int mid37 = GetMidEdgeNode(no[3], no[7]);
1092 
1093  int midf1 = GetMidFaceNode(mid01, mid15, mid45, mid04);
1094  int midf3 = GetMidFaceNode(mid23, mid37, mid67, mid26);
1095 
1096  child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
1097  mid04, midf1, midf3, mid37, attr,
1098  fa[0], fa[1], -1, fa[3], fa[4], -1);
1099 
1100  child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
1101  midf1, mid15, mid26, midf3, attr,
1102  fa[0], fa[1], fa[2], fa[3], -1, -1);
1103 
1104  child[2] = NewHexahedron(midf1, mid15, mid26, midf3,
1105  mid45, no[5], no[6], mid67, attr,
1106  -1, fa[1], fa[2], fa[3], -1, fa[5]);
1107 
1108  child[3] = NewHexahedron(mid04, midf1, midf3, mid37,
1109  no[4], mid45, mid67, no[7], attr,
1110  -1, fa[1], -1, fa[3], fa[4], fa[5]);
1111 
1112  CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
1113  CheckAnisoFace(no[2], no[6], no[5], no[1], mid26, mid15);
1114  CheckAnisoFace(no[6], no[7], no[4], no[5], mid67, mid45);
1115  CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
1116 
1117  CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1118  CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1119  }
1120  else if (ref_type == Refinement::YZ) // YZ split
1121  {
1122  int mid12 = GetMidEdgeNode(no[1], no[2]);
1123  int mid30 = GetMidEdgeNode(no[3], no[0]);
1124  int mid56 = GetMidEdgeNode(no[5], no[6]);
1125  int mid74 = GetMidEdgeNode(no[7], no[4]);
1126 
1127  int mid04 = GetMidEdgeNode(no[0], no[4]);
1128  int mid15 = GetMidEdgeNode(no[1], no[5]);
1129  int mid26 = GetMidEdgeNode(no[2], no[6]);
1130  int mid37 = GetMidEdgeNode(no[3], no[7]);
1131 
1132  int midf2 = GetMidFaceNode(mid12, mid26, mid56, mid15);
1133  int midf4 = GetMidFaceNode(mid30, mid04, mid74, mid37);
1134 
1135  child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
1136  mid04, mid15, midf2, midf4, attr,
1137  fa[0], fa[1], fa[2], -1, fa[4], -1);
1138 
1139  child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
1140  midf4, midf2, mid26, mid37, attr,
1141  fa[0], -1, fa[2], fa[3], fa[4], -1);
1142 
1143  child[2] = NewHexahedron(mid04, mid15, midf2, midf4,
1144  no[4], no[5], mid56, mid74, attr,
1145  -1, fa[1], fa[2], -1, fa[4], fa[5]);
1146 
1147  child[3] = NewHexahedron(midf4, midf2, mid26, mid37,
1148  mid74, mid56, no[6], no[7], attr,
1149  -1, -1, fa[2], fa[3], fa[4], fa[5]);
1150 
1151  CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
1152  CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
1153  CheckAnisoFace(no[3], no[7], no[6], no[2], mid37, mid26);
1154  CheckAnisoFace(no[7], no[4], no[5], no[6], mid74, mid56);
1155 
1156  CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1157  CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1158  }
1159  else if (ref_type == Refinement::XYZ) // full isotropic refinement
1160  {
1161  int mid01 = GetMidEdgeNode(no[0], no[1]);
1162  int mid12 = GetMidEdgeNode(no[1], no[2]);
1163  int mid23 = GetMidEdgeNode(no[2], no[3]);
1164  int mid30 = GetMidEdgeNode(no[3], no[0]);
1165 
1166  int mid45 = GetMidEdgeNode(no[4], no[5]);
1167  int mid56 = GetMidEdgeNode(no[5], no[6]);
1168  int mid67 = GetMidEdgeNode(no[6], no[7]);
1169  int mid74 = GetMidEdgeNode(no[7], no[4]);
1170 
1171  int mid04 = GetMidEdgeNode(no[0], no[4]);
1172  int mid15 = GetMidEdgeNode(no[1], no[5]);
1173  int mid26 = GetMidEdgeNode(no[2], no[6]);
1174  int mid37 = GetMidEdgeNode(no[3], no[7]);
1175 
1176  int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid30);
1177  int midf1 = GetMidFaceNode(mid01, mid15, mid45, mid04);
1178  int midf2 = GetMidFaceNode(mid12, mid26, mid56, mid15);
1179  int midf3 = GetMidFaceNode(mid23, mid37, mid67, mid26);
1180  int midf4 = GetMidFaceNode(mid30, mid04, mid74, mid37);
1181  int midf5 = GetMidFaceNode(mid45, mid56, mid67, mid74);
1182 
1183  int midel = GetMidEdgeNode(midf1, midf3);
1184 
1185  child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
1186  mid04, midf1, midel, midf4, attr,
1187  fa[0], fa[1], -1, -1, fa[4], -1);
1188 
1189  child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
1190  midf1, mid15, midf2, midel, attr,
1191  fa[0], fa[1], fa[2], -1, -1, -1);
1192 
1193  child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
1194  midel, midf2, mid26, midf3, attr,
1195  fa[0], -1, fa[2], fa[3], -1, -1);
1196 
1197  child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
1198  midf4, midel, midf3, mid37, attr,
1199  fa[0], -1, -1, fa[3], fa[4], -1);
1200 
1201  child[4] = NewHexahedron(mid04, midf1, midel, midf4,
1202  no[4], mid45, midf5, mid74, attr,
1203  -1, fa[1], -1, -1, fa[4], fa[5]);
1204 
1205  child[5] = NewHexahedron(midf1, mid15, midf2, midel,
1206  mid45, no[5], mid56, midf5, attr,
1207  -1, fa[1], fa[2], -1, -1, fa[5]);
1208 
1209  child[6] = NewHexahedron(midel, midf2, mid26, midf3,
1210  midf5, mid56, no[6], mid67, attr,
1211  -1, -1, fa[2], fa[3], -1, fa[5]);
1212 
1213  child[7] = NewHexahedron(midf4, midel, midf3, mid37,
1214  mid74, midf5, mid67, no[7], attr,
1215  -1, -1, -1, fa[3], fa[4], fa[5]);
1216 
1217  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1218  CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1219  CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1220  CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1221  CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1222  CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1223  }
1224  else
1225  {
1226  MFEM_ABORT("invalid refinement type.");
1227  }
1228 
1229  if (ref_type != Refinement::XYZ) { Iso = false; }
1230  }
1231  else if (el.Geom() == Geometry::PRISM)
1232  {
1233  // Wedge vertex numbering:
1234  //
1235  // 5
1236  // _+_
1237  // _/ | \_ Faces: 0 bottom
1238  // 3 / | \ 4 1 top
1239  // +---------+ 2 front
1240  // | | | 3 right (1 2 5 4)
1241  // | _+_ | 4 left (2 0 3 5)
1242  // | _/ 2 \_ | Z Y
1243  // |/ \| | /
1244  // +---------+ *--X
1245  // 0 1
1246 
1247  if (ref_type < 4) // XY refinement (split in 4 wedges)
1248  {
1249  ref_type = Refinement::XY; // for consistence
1250 
1251  int mid01 = GetMidEdgeNode(no[0], no[1]);
1252  int mid12 = GetMidEdgeNode(no[1], no[2]);
1253  int mid20 = GetMidEdgeNode(no[2], no[0]);
1254 
1255  int mid34 = GetMidEdgeNode(no[3], no[4]);
1256  int mid45 = GetMidEdgeNode(no[4], no[5]);
1257  int mid53 = GetMidEdgeNode(no[5], no[3]);
1258 
1259  child[0] = NewWedge(no[0], mid01, mid20,
1260  no[3], mid34, mid53, attr,
1261  fa[0], fa[1], fa[2], -1, fa[4]);
1262 
1263  child[1] = NewWedge(mid01, no[1], mid12,
1264  mid34, no[4], mid45, attr,
1265  fa[0], fa[1], fa[2], fa[3], -1);
1266 
1267  child[2] = NewWedge(mid20, mid12, no[2],
1268  mid53, mid45, no[5], attr,
1269  fa[0], fa[1], -1, fa[3], fa[4]);
1270 
1271  child[3] = NewWedge(mid12, mid20, mid01,
1272  mid45, mid53, mid34, attr,
1273  fa[0], fa[1], -1, -1, -1);
1274 
1275  CheckAnisoFace(no[0], no[1], no[4], no[3], mid01, mid34);
1276  CheckAnisoFace(no[1], no[2], no[5], no[4], mid12, mid45);
1277  CheckAnisoFace(no[2], no[0], no[3], no[5], mid20, mid53);
1278  }
1279  else if (ref_type == Refinement::Z) // Z refinement only (split in 2 wedges)
1280  {
1281  int mid03 = GetMidEdgeNode(no[0], no[3]);
1282  int mid14 = GetMidEdgeNode(no[1], no[4]);
1283  int mid25 = GetMidEdgeNode(no[2], no[5]);
1284 
1285  child[0] = NewWedge(no[0], no[1], no[2],
1286  mid03, mid14, mid25, attr,
1287  fa[0], -1, fa[2], fa[3], fa[4]);
1288 
1289  child[1] = NewWedge(mid03, mid14, mid25,
1290  no[3], no[4], no[5], attr,
1291  -1, fa[1], fa[2], fa[3], fa[4]);
1292 
1293  CheckAnisoFace(no[3], no[0], no[1], no[4], mid03, mid14);
1294  CheckAnisoFace(no[4], no[1], no[2], no[5], mid14, mid25);
1295  CheckAnisoFace(no[5], no[2], no[0], no[3], mid25, mid03);
1296  }
1297  else // ref_type > 4, full isotropic refinement (split in 8 wedges)
1298  {
1299  ref_type = Refinement::XYZ; // for consistence
1300 
1301  int mid01 = GetMidEdgeNode(no[0], no[1]);
1302  int mid12 = GetMidEdgeNode(no[1], no[2]);
1303  int mid20 = GetMidEdgeNode(no[2], no[0]);
1304 
1305  int mid34 = GetMidEdgeNode(no[3], no[4]);
1306  int mid45 = GetMidEdgeNode(no[4], no[5]);
1307  int mid53 = GetMidEdgeNode(no[5], no[3]);
1308 
1309  int mid03 = GetMidEdgeNode(no[0], no[3]);
1310  int mid14 = GetMidEdgeNode(no[1], no[4]);
1311  int mid25 = GetMidEdgeNode(no[2], no[5]);
1312 
1313  int midf2 = GetMidFaceNode(mid01, mid14, mid34, mid03);
1314  int midf3 = GetMidFaceNode(mid12, mid25, mid45, mid14);
1315  int midf4 = GetMidFaceNode(mid20, mid03, mid53, mid25);
1316 
1317  child[0] = NewWedge(no[0], mid01, mid20,
1318  mid03, midf2, midf4, attr,
1319  fa[0], -1, fa[2], -1, fa[4]);
1320 
1321  child[1] = NewWedge(mid01, no[1], mid12,
1322  midf2, mid14, midf3, attr,
1323  fa[0], -1, fa[2], fa[3], -1);
1324 
1325  child[2] = NewWedge(mid20, mid12, no[2],
1326  midf4, midf3, mid25, attr,
1327  fa[0], -1, -1, fa[3], fa[4]);
1328 
1329  child[3] = NewWedge(mid12, mid20, mid01,
1330  midf3, midf4, midf2, attr,
1331  fa[0], -1, -1, -1, -1);
1332 
1333  child[4] = NewWedge(mid03, midf2, midf4,
1334  no[3], mid34, mid53, attr,
1335  -1, fa[1], fa[2], -1, fa[4]);
1336 
1337  child[5] = NewWedge(midf2, mid14, midf3,
1338  mid34, no[4], mid45, attr,
1339  -1, fa[1], fa[2], fa[3], -1);
1340 
1341  child[6] = NewWedge(midf4, midf3, mid25,
1342  mid53, mid45, no[5], attr,
1343  -1, fa[1], -1, fa[3], fa[4]);
1344 
1345  child[7] = NewWedge(midf3, midf4, midf2,
1346  mid45, mid53, mid34, attr,
1347  -1, fa[1], -1, -1, -1);
1348 
1349  CheckIsoFace(no[0], no[1], no[4], no[3], mid01, mid14, mid34, mid03, midf2);
1350  CheckIsoFace(no[1], no[2], no[5], no[4], mid12, mid25, mid45, mid14, midf3);
1351  CheckIsoFace(no[2], no[0], no[3], no[5], mid20, mid03, mid53, mid25, midf4);
1352  }
1353 
1354  if (ref_type != Refinement::XYZ) { Iso = false; }
1355  }
1356  else if (el.Geom() == Geometry::TETRAHEDRON)
1357  {
1358  // Tetrahedron vertex numbering:
1359  //
1360  // 3
1361  // + Faces: 0 back (1, 2, 3)
1362  // |\\_ 1 left (0, 3, 2)
1363  // || \_ 2 front (0, 1, 3)
1364  // | \ \_ 3 bottom (0, 1, 2)
1365  // | +__ \_
1366  // | /2 \__ \_ Z Y
1367  // |/ \__\ | /
1368  // +------------+ *--X
1369  // 0 1
1370 
1371  ref_type = Refinement::XYZ; // for consistence
1372 
1373  int mid01 = GetMidEdgeNode(no[0], no[1]);
1374  int mid12 = GetMidEdgeNode(no[1], no[2]);
1375  int mid02 = GetMidEdgeNode(no[2], no[0]);
1376 
1377  int mid03 = GetMidEdgeNode(no[0], no[3]);
1378  int mid13 = GetMidEdgeNode(no[1], no[3]);
1379  int mid23 = GetMidEdgeNode(no[2], no[3]);
1380 
1381  child[0] = NewTetrahedron(no[0], mid01, mid02, mid03, attr,
1382  -1, fa[1], fa[2], fa[3]);
1383 
1384  child[1] = NewTetrahedron(mid01, no[1], mid12, mid13, attr,
1385  fa[0], -1, fa[2], fa[3]);
1386 
1387  child[2] = NewTetrahedron(mid02, mid12, no[2], mid23, attr,
1388  fa[0], fa[1], -1, fa[3]);
1389 
1390  child[3] = NewTetrahedron(mid03, mid13, mid23, no[3], attr,
1391  fa[0], fa[1], fa[2], -1);
1392 
1393  // There are three ways to split the inner octahedron. A good strategy is
1394  // to use the shortest diagonal. At the moment we don't have the geometric
1395  // information in this class to determine which diagonal is the shortest,
1396  // but it seems that with reasonable shapes of the coarse tets and MFEM's
1397  // default tet orientation, always using tet_type == 0 produces stable
1398  // refinements. Types 1 and 2 are unused for now.
1399  el.tet_type = 0;
1400 
1401  if (el.tet_type == 0) // shortest diagonal mid01--mid23
1402  {
1403  child[4] = NewTetrahedron(mid01, mid23, mid02, mid03, attr,
1404  fa[1], -1, -1, -1);
1405 
1406  child[5] = NewTetrahedron(mid01, mid23, mid03, mid13, attr,
1407  -1, fa[2], -1, -1);
1408 
1409  child[6] = NewTetrahedron(mid01, mid23, mid13, mid12, attr,
1410  fa[0], -1, -1, -1);
1411 
1412  child[7] = NewTetrahedron(mid01, mid23, mid12, mid02, attr,
1413  -1, fa[3], -1, -1);
1414  }
1415  else if (el.tet_type == 1) // shortest diagonal mid12--mid03
1416  {
1417  child[4] = NewTetrahedron(mid03, mid01, mid02, mid12, attr,
1418  fa[3], -1, -1, -1);
1419 
1420  child[5] = NewTetrahedron(mid03, mid02, mid23, mid12, attr,
1421  -1, -1, -1, fa[1]);
1422 
1423  child[6] = NewTetrahedron(mid03, mid23, mid13, mid12, attr,
1424  fa[0], -1, -1, -1);
1425 
1426  child[7] = NewTetrahedron(mid03, mid13, mid01, mid12, attr,
1427  -1, -1, -1, fa[2]);
1428  }
1429  else // el.tet_type == 2, shortest diagonal mid02--mid13
1430  {
1431  child[4] = NewTetrahedron(mid02, mid01, mid13, mid03, attr,
1432  fa[2], -1, -1, -1);
1433 
1434  child[5] = NewTetrahedron(mid02, mid03, mid13, mid23, attr,
1435  -1, -1, fa[1], -1);
1436 
1437  child[6] = NewTetrahedron(mid02, mid23, mid13, mid12, attr,
1438  fa[0], -1, -1, -1);
1439 
1440  child[7] = NewTetrahedron(mid02, mid12, mid13, mid01, attr,
1441  -1, -1, fa[3], -1);
1442  }
1443  }
1444  else if (el.Geom() == Geometry::PYRAMID)
1445  {
1446  // Pyramid vertex numbering:
1447  //
1448  // 4
1449  // + \_ Faces: 0 bottom (3,2,1,0)
1450  // |\\_ \_ 1 front (0, 1, 4)
1451  // || \_ \__ 2 right (1, 2, 4)
1452  // | \ \_ \__ 3 back (2, 3, 4)
1453  // | +____\_ ____\ 4 left (3, 0, 4)
1454  // | /3 \_ 2 Z Y
1455  // |/ \ / | /
1456  // +------------+ *--X
1457  // 0 1
1458 
1459  ref_type = Refinement::XYZ; // for consistence
1460 
1461  int mid01 = GetMidEdgeNode(no[0], no[1]);
1462  int mid12 = GetMidEdgeNode(no[1], no[2]);
1463  int mid23 = GetMidEdgeNode(no[2], no[3]);
1464  int mid03 = GetMidEdgeNode(no[0], no[3]);
1465  int mid04 = GetMidEdgeNode(no[0], no[4]);
1466  int mid14 = GetMidEdgeNode(no[1], no[4]);
1467  int mid24 = GetMidEdgeNode(no[2], no[4]);
1468  int mid34 = GetMidEdgeNode(no[3], no[4]);
1469  int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid03);
1470 
1471  child[0] = NewPyramid(no[0], mid01, midf0, mid03, mid04,
1472  attr, fa[0], fa[1], -1, -1, fa[4]);
1473 
1474  child[1] = NewPyramid(mid01, no[1], mid12, midf0, mid14,
1475  attr, fa[0], fa[1], fa[2], -1, -1);
1476 
1477  child[2] = NewPyramid(midf0, mid12, no[2], mid23, mid24,
1478  attr, fa[0], -1, fa[2], fa[3], -1);
1479 
1480  child[3] = NewPyramid(mid03, midf0, mid23, no[3], mid34,
1481  attr, fa[0], -1, -1, fa[3], fa[4]);
1482 
1483  child[4] = NewPyramid(mid24, mid14, mid04, mid34, midf0,
1484  attr, -1, -1, -1, -1, -1);
1485 
1486  child[5] = NewPyramid(mid04, mid14, mid24, mid34, no[4],
1487  attr, -1, fa[1], fa[2], fa[3], fa[4]);
1488 
1489  child[6] = NewTetrahedron(mid01, midf0, mid04, mid14,
1490  attr, -1, -1, -1, fa[1]);
1491 
1492  child[7] = NewTetrahedron(midf0, mid14, mid12, mid24,
1493  attr, -1, -1, fa[2], -1);
1494 
1495  child[8] = NewTetrahedron(midf0, mid23, mid34, mid24,
1496  attr, -1, -1, fa[3], -1);
1497 
1498  child[9] = NewTetrahedron(mid03, mid04, midf0, mid34,
1499  attr, -1, fa[4], -1, -1);
1500 
1501  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid03, midf0);
1502  }
1503  else if (el.Geom() == Geometry::SQUARE)
1504  {
1505  ref_type &= 0x3; // ignore Z bit
1506 
1507  if (ref_type == Refinement::X) // X split
1508  {
1509  int mid01 = nodes.GetId(no[0], no[1]);
1510  int mid23 = nodes.GetId(no[2], no[3]);
1511 
1512  child[0] = NewQuadrilateral(no[0], mid01, mid23, no[3],
1513  attr, fa[0], -1, fa[2], fa[3]);
1514 
1515  child[1] = NewQuadrilateral(mid01, no[1], no[2], mid23,
1516  attr, fa[0], fa[1], fa[2], -1);
1517  }
1518  else if (ref_type == Refinement::Y) // Y split
1519  {
1520  int mid12 = nodes.GetId(no[1], no[2]);
1521  int mid30 = nodes.GetId(no[3], no[0]);
1522 
1523  child[0] = NewQuadrilateral(no[0], no[1], mid12, mid30,
1524  attr, fa[0], fa[1], -1, fa[3]);
1525 
1526  child[1] = NewQuadrilateral(mid30, mid12, no[2], no[3],
1527  attr, -1, fa[1], fa[2], fa[3]);
1528  }
1529  else if (ref_type == Refinement::XY) // iso split
1530  {
1531  int mid01 = nodes.GetId(no[0], no[1]);
1532  int mid12 = nodes.GetId(no[1], no[2]);
1533  int mid23 = nodes.GetId(no[2], no[3]);
1534  int mid30 = nodes.GetId(no[3], no[0]);
1535 
1536  int midel = nodes.GetId(mid01, mid23);
1537 
1538  child[0] = NewQuadrilateral(no[0], mid01, midel, mid30,
1539  attr, fa[0], -1, -1, fa[3]);
1540 
1541  child[1] = NewQuadrilateral(mid01, no[1], mid12, midel,
1542  attr, fa[0], fa[1], -1, -1);
1543 
1544  child[2] = NewQuadrilateral(midel, mid12, no[2], mid23,
1545  attr, -1, fa[1], fa[2], -1);
1546 
1547  child[3] = NewQuadrilateral(mid30, midel, mid23, no[3],
1548  attr, -1, -1, fa[2], fa[3]);
1549  }
1550  else
1551  {
1552  MFEM_ABORT("Invalid refinement type.");
1553  }
1554 
1555  if (ref_type != Refinement::XY) { Iso = false; }
1556  }
1557  else if (el.Geom() == Geometry::TRIANGLE)
1558  {
1559  ref_type = Refinement::XY; // for consistence
1560 
1561  // isotropic split - the only ref_type available for triangles
1562  int mid01 = nodes.GetId(no[0], no[1]);
1563  int mid12 = nodes.GetId(no[1], no[2]);
1564  int mid20 = nodes.GetId(no[2], no[0]);
1565 
1566  child[0] = NewTriangle(no[0], mid01, mid20, attr, fa[0], -1, fa[2]);
1567  child[1] = NewTriangle(mid01, no[1], mid12, attr, fa[0], fa[1], -1);
1568  child[2] = NewTriangle(mid20, mid12, no[2], attr, -1, fa[1], fa[2]);
1569  child[3] = NewTriangle(mid12, mid20, mid01, attr, -1, -1, -1);
1570  }
1571  else if (el.Geom() == Geometry::SEGMENT)
1572  {
1573  ref_type = Refinement::X; // for consistence
1574 
1575  int mid = nodes.GetId(no[0], no[1]);
1576  child[0] = NewSegment(no[0], mid, attr, fa[0], -1);
1577  child[1] = NewSegment(mid, no[1], attr, -1, fa[1]);
1578  }
1579  else
1580  {
1581  MFEM_ABORT("Unsupported element geometry.");
1582  }
1583 
1584  // start using the nodes of the children, create edges & faces
1585  for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1586  {
1587  ReferenceElement(child[i]);
1588  }
1589 
1590  int buf[MaxElemFaces];
1591  Array<int> parentFaces(buf, MaxElemFaces);
1592  parentFaces.SetSize(0);
1593 
1594  // sign off of all nodes of the parent, clean up unused nodes, but keep faces
1595  UnreferenceElement(elem, parentFaces);
1596 
1597  // register the children in their faces
1598  for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1599  {
1600  RegisterFaces(child[i]);
1601  }
1602 
1603  // clean up parent faces, if unused
1604  DeleteUnusedFaces(parentFaces);
1605 
1606  // make the children inherit our rank; set the parent element
1607  for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1608  {
1609  Element &ch = elements[child[i]];
1610  ch.rank = el.rank;
1611  ch.parent = elem;
1612  }
1613 
1614  // finish the refinement
1615  el.ref_type = ref_type;
1616  std::memcpy(el.child, child, sizeof(el.child));
1617 }
1618 
1619 
1620 void NCMesh::Refine(const Array<Refinement>& refinements)
1621 {
1622  // push all refinements on the stack in reverse order
1623  ref_stack.Reserve(refinements.Size());
1624  for (int i = refinements.Size()-1; i >= 0; i--)
1625  {
1626  const Refinement& ref = refinements[i];
1627  ref_stack.Append(Refinement(leaf_elements[ref.index], ref.ref_type));
1628  }
1629 
1630  // keep refining as long as the stack contains something
1631  int nforced = 0;
1632  while (ref_stack.Size())
1633  {
1634  Refinement ref = ref_stack.Last();
1635  ref_stack.DeleteLast();
1636 
1637  int size = ref_stack.Size();
1638  RefineElement(ref.index, ref.ref_type);
1639  nforced += ref_stack.Size() - size;
1640  }
1641 
1642  /* TODO: the current algorithm of forced refinements is not optimal. As
1643  forced refinements spread through the mesh, some may not be necessary
1644  in the end, since the affected elements may still be scheduled for
1645  refinement that could stop the propagation. We should introduce the
1646  member Element::ref_pending that would show the intended refinement in
1647  the batch. A forced refinement would be combined with ref_pending to
1648  (possibly) stop the propagation earlier.
1649 
1650  Update: what about a FIFO instead of ref_stack? */
1651 
1652 #if defined(MFEM_DEBUG) && !defined(MFEM_USE_MPI)
1653  mfem::out << "Refined " << refinements.Size() << " + " << nforced
1654  << " elements" << std::endl;
1655 #else
1656  MFEM_CONTRACT_VAR(nforced);
1657 #endif
1658 
1659  ref_stack.DeleteAll();
1660  shadow.DeleteAll();
1661 
1662  Update();
1663 }
1664 
1665 
1666 //// Derefinement //////////////////////////////////////////////////////////////
1667 
1669 {
1670  if (!el.ref_type) { return el.node[index]; }
1671 
1672  // need to retrieve node from a child element (there is always a child
1673  // that inherited the parent's corner under the same index)
1674  int ch;
1675  switch (el.Geom())
1676  {
1677  case Geometry::CUBE:
1678  ch = el.child[hex_deref_table[el.ref_type - 1][index]];
1679  break;
1680 
1681  case Geometry::PRISM:
1682  ch = prism_deref_table[el.ref_type - 1][index];
1683  MFEM_ASSERT(ch != -1, "");
1684  ch = el.child[ch];
1685  break;
1686 
1687  case Geometry::PYRAMID:
1688  ch = pyramid_deref_table[el.ref_type - 1][index];
1689  MFEM_ASSERT(ch != -1, "");
1690  ch = el.child[ch];
1691  break;
1692 
1693  case Geometry::SQUARE:
1694  ch = el.child[quad_deref_table[el.ref_type - 1][index]];
1695  break;
1696 
1697  case Geometry::TETRAHEDRON:
1698  case Geometry::TRIANGLE:
1699  ch = el.child[index];
1700  break;
1701 
1702  default:
1703  ch = 0; // suppress compiler warning
1704  MFEM_ABORT("Unsupported element geometry.");
1705  }
1706  return RetrieveNode(elements[ch], index);
1707 }
1708 
1709 
1711 {
1712  Element &el = elements[elem];
1713  if (!el.ref_type) { return; }
1714 
1715  int child[MaxElemChildren];
1716  std::memcpy(child, el.child, sizeof(child));
1717 
1718  // first make sure that all children are leaves, derefine them if not
1719  for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1720  {
1721  if (elements[child[i]].ref_type)
1722  {
1723  DerefineElement(child[i]);
1724  }
1725  }
1726 
1727  int faces_attribute[MaxElemFaces];
1728  int ref_type_key = el.ref_type - 1;
1729 
1730  for (int i = 0; i < MaxElemNodes; i++) { el.node[i] = -1; }
1731 
1732  // retrieve original corner nodes and face attributes from the children
1733  if (el.Geom() == Geometry::CUBE)
1734  {
1735  // Sets corner nodes from childs
1736  constexpr int nb_cube_childs = 8;
1737  for (int i = 0; i < nb_cube_childs; i++)
1738  {
1739  const int child_local_index = hex_deref_table[ref_type_key][i];
1740  const int child_global_index = child[child_local_index];
1741  Element &ch = elements[child_global_index];
1742  el.node[i] = ch.node[i];
1743  }
1744  // Sets faces attributes from childs' faces
1745  constexpr int nb_cube_faces = 6;
1746  for (int i = 0; i < nb_cube_faces; i++)
1747  {
1748  const int child_local_index = hex_deref_table[ref_type_key]
1749  [i + nb_cube_childs];
1750  const int child_global_index = child[child_local_index];
1751  Element &ch = elements[child_global_index];
1752  const int* fv = GI[el.Geom()].faces[i];
1753  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1754  ch.node[fv[2]], ch.node[fv[3]])
1755  ->attribute;
1756  }
1757  }
1758  else if (el.Geom() == Geometry::PRISM)
1759  {
1760  MFEM_ASSERT(prism_deref_table[ref_type_key][0] != -1,
1761  "invalid prism refinement");
1762  constexpr int nb_prism_childs = 6;
1763  for (int i = 0; i < nb_prism_childs; i++)
1764  {
1765  const int child_local_index = prism_deref_table[ref_type_key][i];
1766  const int child_global_index = child[child_local_index];
1767  Element &ch = elements[child_global_index];
1768  el.node[i] = ch.node[i];
1769  }
1770  el.node[6] = el.node[7] = -1;
1771 
1772  constexpr int nb_prism_faces = 5;
1773  for (int i = 0; i < nb_prism_faces; i++)
1774  {
1775  const int child_local_index = prism_deref_table[ref_type_key]
1776  [i + nb_prism_childs];
1777  const int child_global_index = child[child_local_index];
1778  Element &ch = elements[child_global_index];
1779  const int* fv = GI[el.Geom()].faces[i];
1780  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1781  ch.node[fv[2]], ch.node[fv[3]])
1782  ->attribute;
1783  }
1784  }
1785  else if (el.Geom() == Geometry::PYRAMID)
1786  {
1787  MFEM_ASSERT(pyramid_deref_table[ref_type_key][0] != -1,
1788  "invalid pyramid refinement");
1789  constexpr int nb_pyramid_childs = 5;
1790  for (int i = 0; i < nb_pyramid_childs; i++)
1791  {
1792  const int child_local_index = pyramid_deref_table[ref_type_key][i];
1793  const int child_global_index = child[child_local_index];
1794  Element &ch = elements[child_global_index];
1795  el.node[i] = ch.node[i];
1796  }
1797  el.node[5] = el.node[6] = el.node[7] = -1;
1798 
1799 
1800  constexpr int nb_pyramid_faces = 5;
1801  for (int i = 0; i < nb_pyramid_faces; i++)
1802  {
1803  const int child_local_index = pyramid_deref_table[ref_type_key]
1804  [i + nb_pyramid_childs];
1805  const int child_global_index = child[child_local_index];
1806  Element &ch = elements[child_global_index];
1807  const int* fv = GI[el.Geom()].faces[i];
1808  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1809  ch.node[fv[2]], ch.node[fv[3]])
1810  ->attribute;
1811  }
1812  }
1813  else if (el.Geom() == Geometry::TETRAHEDRON)
1814  {
1815  for (int i = 0; i < 4; i++)
1816  {
1817  Element& ch1 = elements[child[i]];
1818  Element& ch2 = elements[child[(i+1) & 0x3]];
1819  el.node[i] = ch1.node[i];
1820  const int* fv = GI[el.Geom()].faces[i];
1821  faces_attribute[i] = faces.Find(ch2.node[fv[0]], ch2.node[fv[1]],
1822  ch2.node[fv[2]], ch2.node[fv[3]])
1823  ->attribute;
1824  }
1825  }
1826  else if (el.Geom() == Geometry::SQUARE)
1827  {
1828  constexpr int nb_square_childs = 4;
1829  for (int i = 0; i < nb_square_childs; i++)
1830  {
1831  const int child_local_index = quad_deref_table[ref_type_key][i];
1832  const int child_global_index = child[child_local_index];
1833  Element &ch = elements[child_global_index];
1834  el.node[i] = ch.node[i];
1835  }
1836  constexpr int nb_square_faces = 4;
1837  for (int i = 0; i < nb_square_faces; i++)
1838  {
1839  const int child_local_index = quad_deref_table[ref_type_key]
1840  [i + nb_square_childs];
1841  const int child_global_index = child[child_local_index];
1842  Element &ch = elements[child_global_index];
1843  const int* fv = GI[el.Geom()].faces[i];
1844  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1845  ch.node[fv[2]], ch.node[fv[3]])
1846  ->attribute;
1847  }
1848  }
1849  else if (el.Geom() == Geometry::TRIANGLE)
1850  {
1851  constexpr int nb_triangle_childs = 3;
1852  for (int i = 0; i < nb_triangle_childs; i++)
1853  {
1854  Element& ch = elements[child[i]];
1855  el.node[i] = ch.node[i];
1856  const int* fv = GI[el.Geom()].faces[i];
1857  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1858  ch.node[fv[2]], ch.node[fv[3]])
1859  ->attribute;
1860  }
1861  }
1862  else if (el.Geom() == Geometry::SEGMENT)
1863  {
1864  constexpr int nb_segment_childs = 2;
1865  for (int i = 0; i < nb_segment_childs; i++)
1866  {
1867  int ni = elements[child[i]].node[i];
1868  el.node[i] = ni;
1869  faces_attribute[i] = faces.Find(ni, ni, ni, ni)->attribute;
1870  }
1871  }
1872  else
1873  {
1874  MFEM_ABORT("Unsupported element geometry.");
1875  }
1876 
1877  // sign in to all nodes
1878  ReferenceElement(elem);
1879 
1880  int buf[MaxElemChildren*MaxElemFaces];
1881  Array<int> childFaces(buf, MaxElemChildren*MaxElemFaces);
1882  childFaces.SetSize(0);
1883 
1884  // delete children, determine rank
1885  el.rank = std::numeric_limits<int>::max();
1886  for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
1887  {
1888  el.rank = std::min(el.rank, elements[child[i]].rank);
1889  UnreferenceElement(child[i], childFaces);
1890  FreeElement(child[i]);
1891  }
1892 
1893  RegisterFaces(elem, faces_attribute);
1894 
1895  // delete unused faces
1896  childFaces.Sort();
1897  childFaces.Unique();
1898  DeleteUnusedFaces(childFaces);
1899 
1900  el.ref_type = 0;
1901 }
1902 
1903 
1905 {
1906  Element &el = elements[elem];
1907  if (!el.ref_type) { return; }
1908 
1909  int total = 0, ref = 0, ghost = 0;
1910  for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
1911  {
1912  total++;
1913  Element &ch = elements[el.child[i]];
1914  if (ch.ref_type) { ref++; break; }
1915  if (IsGhost(ch)) { ghost++; }
1916  }
1917 
1918  if (!ref && ghost < total)
1919  {
1920  // can be derefined, add to list
1921  int next_row = list.Size() ? (list.Last().from + 1) : 0;
1922  for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
1923  {
1924  Element &ch = elements[el.child[i]];
1925  list.Append(Connection(next_row, ch.index));
1926  }
1927  }
1928  else
1929  {
1930  for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
1931  {
1932  CollectDerefinements(el.child[i], list);
1933  }
1934  }
1935 }
1936 
1938 {
1939  Array<Connection> list;
1940  list.Reserve(leaf_elements.Size());
1941 
1942  for (int i = 0; i < root_state.Size(); i++)
1943  {
1944  CollectDerefinements(i, list);
1945  }
1946 
1947  int size = list.Size() ? (list.Last().from + 1) : 0;
1948  derefinements.MakeFromList(size, list);
1949  return derefinements;
1950 }
1951 
1952 void NCMesh::CheckDerefinementNCLevel(const Table &deref_table,
1953  Array<int> &level_ok, int max_nc_level)
1954 {
1955  level_ok.SetSize(deref_table.Size());
1956  for (int i = 0; i < deref_table.Size(); i++)
1957  {
1958  const int* fine = deref_table.GetRow(i), size = deref_table.RowSize(i);
1959  Element &parent = elements[elements[leaf_elements[fine[0]]].parent];
1960 
1961  int ok = 1;
1962  for (int j = 0; j < size; j++)
1963  {
1964  int splits[3];
1965  CountSplits(leaf_elements[fine[j]], splits);
1966 
1967  for (int k = 0; k < Dim; k++)
1968  {
1969  if ((parent.ref_type & (1 << k)) &&
1970  splits[k] >= max_nc_level)
1971  {
1972  ok = 0; break;
1973  }
1974  }
1975  if (!ok) { break; }
1976  }
1977  level_ok[i] = ok;
1978  }
1979 }
1980 
1981 void NCMesh::Derefine(const Array<int> &derefs)
1982 {
1983  MFEM_VERIFY(Dim < 3 || Iso,
1984  "derefinement of 3D anisotropic meshes not implemented yet.");
1985 
1987 
1988  Array<int> fine_coarse;
1989  leaf_elements.Copy(fine_coarse);
1990 
1991  // perform the derefinements
1992  for (int i = 0; i < derefs.Size(); i++)
1993  {
1994  int row = derefs[i];
1995  MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
1996  "invalid derefinement number.");
1997 
1998  const int* fine = derefinements.GetRow(row);
1999  int parent = elements[leaf_elements[fine[0]]].parent;
2000 
2001  // record the relation of the fine elements to their parent
2002  SetDerefMatrixCodes(parent, fine_coarse);
2003 
2004  DerefineElement(parent);
2005  }
2006 
2007  // update leaf_elements, Element::index etc.
2008  Update();
2009 
2010  // link old fine elements to the new coarse elements
2011  for (int i = 0; i < fine_coarse.Size(); i++)
2012  {
2013  transforms.embeddings[i].parent = elements[fine_coarse[i]].index;
2014  }
2015 }
2016 
2018 {
2019  int nfine = leaf_elements.Size();
2020 
2021  // this will tell GetDerefinementTransforms that transforms are not finished
2022  transforms.Clear();
2023 
2024  transforms.embeddings.SetSize(nfine);
2025  for (int i = 0; i < nfine; i++)
2026  {
2027  Embedding &emb = transforms.embeddings[i];
2028  emb.parent = -1;
2029  emb.matrix = 0;
2030  Element &el = elements[leaf_elements[i]];
2031  emb.geom = el.Geom();
2032  emb.ghost = IsGhost(el);
2033  }
2034 }
2035 
2036 void NCMesh::SetDerefMatrixCodes(int parent, Array<int> &fine_coarse)
2037 {
2038  // encode the ref_type and child number for GetDerefinementTransforms()
2039  Element &prn = elements[parent];
2040  for (int i = 0; i < MaxElemChildren && prn.child[i] >= 0; i++)
2041  {
2042  Element &ch = elements[prn.child[i]];
2043  if (ch.index >= 0)
2044  {
2045  int code = (prn.ref_type << 4) | i;
2046  transforms.embeddings[ch.index].matrix = code;
2047  fine_coarse[ch.index] = parent;
2048  }
2049  }
2050 }
2051 
2052 
2053 //// Mesh Interface ////////////////////////////////////////////////////////////
2054 
2055 void NCMesh::CollectLeafElements(int elem, int state, Array<int> &ghosts,
2056  int &counter)
2057 {
2058  Element &el = elements[elem];
2059  if (!el.ref_type)
2060  {
2061  if (el.rank >= 0) // skip elements beyond the ghost layer in parallel
2062  {
2063  if (!IsGhost(el))
2064  {
2065  leaf_elements.Append(elem);
2066  }
2067  else
2068  {
2069  // in parallel (or in serial loading a parallel file), collect
2070  // elements of neighboring ranks in a separate array
2071  ghosts.Append(elem);
2072  }
2073 
2074  // assign the SFC index (temporarily, will be replaced by Mesh index)
2075  el.index = counter++;
2076  }
2077  else
2078  {
2079  // elements beyond the ghost layer are invalid and don't appear in
2080  // 'leaf_elements' (also for performance reasons)
2081  el.index = -1;
2082  }
2083  }
2084  else // Refined element
2085  {
2086  // in non-leaf elements, the 'rank' and 'index' members have no meaning
2087  el.rank = -1;
2088  el.index = -1;
2089 
2090  // recurse to subtrees; try to order leaf elements along a space-filling
2091  // curve by changing the order the children are visited at each level
2092  if (el.Geom() == Geometry::SQUARE && el.ref_type == Refinement::XY)
2093  {
2094  for (int i = 0; i < 4; i++)
2095  {
2096  int ch = quad_hilbert_child_order[state][i];
2097  int st = quad_hilbert_child_state[state][i];
2098  CollectLeafElements(el.child[ch], st, ghosts, counter);
2099  }
2100  }
2101  else if (el.Geom() == Geometry::CUBE && el.ref_type == Refinement::XYZ)
2102  {
2103  for (int i = 0; i < 8; i++)
2104  {
2105  int ch = hex_hilbert_child_order[state][i];
2106  int st = hex_hilbert_child_state[state][i];
2107  CollectLeafElements(el.child[ch], st, ghosts, counter);
2108  }
2109  }
2110  else // no space filling curve tables yet for remaining cases
2111  {
2112  for (int i = 0; i < MaxElemChildren; i++)
2113  {
2114  if (el.child[i] >= 0)
2115  {
2116  CollectLeafElements(el.child[i], state, ghosts, counter);
2117  }
2118  }
2119  }
2120  }
2121 }
2122 
2124 {
2125  Array<int> ghosts;
2126 
2127  // collect leaf elements in leaf_elements and ghosts elements in ghosts from
2128  // all roots
2130  for (int i = 0, counter = 0; i < root_state.Size(); i++)
2131  {
2132  CollectLeafElements(i, root_state[i], ghosts, counter);
2133  }
2134 
2136  NGhostElements = ghosts.Size();
2137 
2138  // append ghost elements at the end of 'leaf_element' (if any)
2139  // and assign the final (Mesh) indices of leaves
2140  leaf_elements.Append(ghosts);
2142 
2143  for (int i = 0; i < leaf_elements.Size(); i++)
2144  {
2145  Element &el = elements[leaf_elements[i]];
2146  leaf_sfc_index[i] = el.index;
2147  el.index = i;
2148  }
2149 }
2150 
2152 {
2153 #ifndef MFEM_NCMESH_OLD_VERTEX_ORDERING
2154  // This method assigns indices to vertices (Node::vert_index) that will
2155  // be seen by the Mesh class and the rest of MFEM. We must be careful to:
2156  //
2157  // 1. Stay compatible with the conforming code, which expects top-level
2158  // (original) vertices to be indexed first, otherwise GridFunctions
2159  // defined on a conforming mesh would no longer be valid when the
2160  // mesh is converted to an NC mesh.
2161  //
2162  // 2. Make sure serial NCMesh is compatible with the parallel ParNCMesh,
2163  // so it is possible to read parallel partial solutions in serial code
2164  // (e.g., serial GLVis). This means handling ghost elements, if present.
2165  //
2166  // 3. Assign vertices in a globally consistent order for parallel meshes:
2167  // if two vertices i,j are shared by two ranks r1,r2, and i<j on r1,
2168  // then i<j on r2 as well. This is true for top-level vertices but also
2169  // for the remaining shared vertices thanks to the globally consistent
2170  // SFC ordering of the leaf elements. This property reduces communication
2171  // and simplifies ParNCMesh.
2172 
2173  // STEP 1: begin by splitting vertices into 4 classes:
2174  // - local top-level vertices (code -1)
2175  // - local non-top level vertices (code -2)
2176  // - ghost (non-local) vertices (code -3)
2177  // - vertices beyond the ghost layer (code -4)
2178 
2179  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2180  {
2181  node->vert_index = -4; // assume beyond ghost layer
2182  }
2183 
2184  for (int i = 0; i < leaf_elements.Size(); i++)
2185  {
2186  Element &el = elements[leaf_elements[i]];
2187  for (int j = 0; j < GI[el.Geom()].nv; j++)
2188  {
2189  Node &nd = nodes[el.node[j]];
2190  if (el.rank == MyRank)
2191  {
2192  if (nd.p1 == nd.p2) // local top-level vertex
2193  {
2194  if (nd.vert_index < -1) { nd.vert_index = -1; }
2195  }
2196  else // local non-top-level vertex
2197  {
2198  if (nd.vert_index < -2) { nd.vert_index = -2; }
2199  }
2200  }
2201  else // ghost vertex
2202  {
2203  if (nd.vert_index < -3) { nd.vert_index = -3; }
2204  }
2205  }
2206  }
2207 
2208  // STEP 2: assign indices of top-level local vertices, in original order
2209 
2210  NVertices = 0;
2211  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2212  {
2213  if (node->vert_index == -1)
2214  {
2215  node->vert_index = NVertices++;
2216  }
2217  }
2218 
2219  // STEP 3: go over all elements (local and ghost) in SFC order and assign
2220  // remaining local vertices in that order.
2221 
2222  Array<int> sfc_order(leaf_elements.Size());
2223  for (int i = 0; i < sfc_order.Size(); i++)
2224  {
2225  sfc_order[leaf_sfc_index[i]] = leaf_elements[i];
2226  }
2227 
2228  for (int i = 0; i < sfc_order.Size(); i++)
2229  {
2230  const Element &el = elements[sfc_order[i]];
2231  for (int j = 0; j < GI[el.Geom()].nv; j++)
2232  {
2233  Node &nd = nodes[el.node[j]];
2234  if (nd.vert_index == -2) { nd.vert_index = NVertices++; }
2235  }
2236  }
2237 
2238  // STEP 4: create the mapping from Mesh vertex index to NCMesh node index
2239 
2241  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2242  {
2243  if (node->HasVertex() && node->vert_index >= 0)
2244  {
2245  MFEM_ASSERT(node->vert_index < vertex_nodeId.Size(), "");
2246  vertex_nodeId[node->vert_index] = node.index();
2247  }
2248  }
2249 
2250  // STEP 5: assign remaining ghost vertices, ignore vertices beyond the
2251  // ghost layer
2252 
2253  NGhostVertices = 0;
2254  for (int i = 0; i < sfc_order.Size(); i++)
2255  {
2256  const Element &el = elements[sfc_order[i]];
2257  for (int j = 0; j < GI[el.Geom()].nv; j++)
2258  {
2259  Node &nd = nodes[el.node[j]];
2260  if (nd.vert_index == -3)
2261  {
2263  }
2264  }
2265  }
2266 
2267 #else // old ordering for debugging/testing only
2268  bool parallel = false;
2269 #ifdef MFEM_USE_MPI
2270  if (dynamic_cast<ParNCMesh*>(this)) { parallel = true; }
2271 #endif
2272 
2273  if (!parallel)
2274  {
2275  NVertices = 0;
2276  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2277  {
2278  if (node->HasVertex()) { node->vert_index = NVertices++; }
2279  }
2280 
2282 
2283  NVertices = 0;
2284  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2285  {
2286  if (node->HasVertex()) { vertex_nodeId[NVertices++] = node.index(); }
2287  }
2288  }
2289  else
2290  {
2291  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2292  {
2293  if (node->HasVertex()) { node->vert_index = -1; }
2294  }
2295 
2296  NVertices = 0;
2297  for (int i = 0; i < leaf_elements.Size(); i++)
2298  {
2299  Element &el = elements[leaf_elements[i]];
2300  if (el.rank == MyRank)
2301  {
2302  for (int j = 0; j < GI[el.Geom()].nv; j++)
2303  {
2304  int &vindex = nodes[el.node[j]].vert_index;
2305  if (vindex < 0) { vindex = NVertices++; }
2306  }
2307  }
2308  }
2309 
2311  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2312  {
2313  if (node->HasVertex() && node->vert_index >= 0)
2314  {
2315  vertex_nodeId[node->vert_index] = node.index();
2316  }
2317  }
2318 
2319  NGhostVertices = 0;
2320  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2321  {
2322  if (node->HasVertex() && node->vert_index < 0)
2323  {
2324  node->vert_index = NVertices + (NGhostVertices++);
2325  }
2326  }
2327  }
2328 #endif
2329 }
2330 
2331 void NCMesh::InitRootState(int root_count)
2332 {
2333  root_state.SetSize(root_count);
2334  root_state = 0;
2335 
2336  char* node_order;
2337  int nch;
2338 
2339  switch (elements[0].Geom()) // TODO: mixed meshes
2340  {
2341  case Geometry::SQUARE:
2342  nch = 4;
2343  node_order = (char*) quad_hilbert_child_order;
2344  break;
2345 
2346  case Geometry::CUBE:
2347  nch = 8;
2348  node_order = (char*) hex_hilbert_child_order;
2349  break;
2350 
2351  default:
2352  return; // do nothing, all states stay zero
2353  }
2354 
2355  int entry_node = -2;
2356 
2357  // process the root element sequence
2358  for (int i = 0; i < root_count; i++)
2359  {
2360  Element &el = elements[i];
2361 
2362  int v_in = FindNodeExt(el, entry_node, false);
2363  if (v_in < 0) { v_in = 0; }
2364 
2365  // determine which nodes are shared with the next element
2366  bool shared[MaxElemNodes];
2367  for (int ni = 0; ni < MaxElemNodes; ++ni) { shared[ni] = 0; }
2368  if (i+1 < root_count)
2369  {
2370  Element &next = elements[i+1];
2371  for (int j = 0; j < nch; j++)
2372  {
2373  int node = FindNodeExt(el, RetrieveNode(next, j), false);
2374  if (node >= 0) { shared[node] = true; }
2375  }
2376  }
2377 
2378  // select orientation that starts in v_in and exits in shared node
2379  int state = Dim*v_in;
2380  for (int j = 0; j < Dim; j++)
2381  {
2382  if (shared[(int) node_order[nch*(state + j) + nch-1]])
2383  {
2384  state += j;
2385  break;
2386  }
2387  }
2388 
2389  root_state[i] = state;
2390 
2391  entry_node = RetrieveNode(el, node_order[nch*state + nch-1]);
2392  }
2393 }
2394 
2396 {
2397  switch (geom)
2398  {
2399  case Geometry::CUBE: return new mfem::Hexahedron;
2400  case Geometry::PRISM: return new mfem::Wedge;
2401  case Geometry::PYRAMID: return new mfem::Pyramid;
2402  case Geometry::TETRAHEDRON: return new mfem::Tetrahedron;
2403  case Geometry::SQUARE: return new mfem::Quadrilateral;
2404  case Geometry::TRIANGLE: return new mfem::Triangle;
2405  }
2406  MFEM_ABORT("invalid geometry");
2407  return NULL;
2408 }
2409 
2410 const double* NCMesh::CalcVertexPos(int node) const
2411 {
2412  const Node &nd = nodes[node];
2413  if (nd.p1 == nd.p2) // top-level vertex
2414  {
2415  return &coordinates[3*nd.p1];
2416  }
2417 
2418  TmpVertex &tv = tmp_vertex[node];
2419  if (tv.valid) { return tv.pos; }
2420 
2421  MFEM_VERIFY(tv.visited == false, "cyclic vertex dependencies.");
2422  tv.visited = true;
2423 
2424  const double* pos1 = CalcVertexPos(nd.p1);
2425  const double* pos2 = CalcVertexPos(nd.p2);
2426 
2427  for (int i = 0; i < 3; i++)
2428  {
2429  tv.pos[i] = (pos1[i] + pos2[i]) * 0.5;
2430  }
2431  tv.valid = true;
2432  return tv.pos;
2433 }
2434 
2436 {
2437  mesh.vertices.SetSize(vertex_nodeId.Size());
2438  if (coordinates.Size())
2439  {
2440  // calculate vertex positions from stored top-level vertex coordinates
2441  tmp_vertex = new TmpVertex[nodes.NumIds()];
2442  for (int i = 0; i < mesh.vertices.Size(); i++)
2443  {
2444  mesh.vertices[i].SetCoords(spaceDim, CalcVertexPos(vertex_nodeId[i]));
2445  }
2446  delete [] tmp_vertex;
2447  }
2448  // NOTE: if the mesh is curved ('coordinates' is empty), mesh.vertices are
2449  // left uninitialized here; they will be initialized later by the Mesh from
2450  // Nodes -- here we just make sure mesh.vertices has the correct size.
2451 
2452  mesh.elements.SetSize(0);
2453 
2454  mesh.boundary.SetSize(0);
2455 
2456  // create an mfem::Element for each leaf Element
2457  for (int i = 0; i < NElements; i++)
2458  {
2459  const Element &nc_elem = elements[leaf_elements[i]];
2460 
2461  const int* node = nc_elem.node;
2462  GeomInfo& gi = GI[(int) nc_elem.geom];
2463 
2464  mfem::Element* elem = mesh.NewElement(nc_elem.geom);
2465  mesh.elements.Append(elem);
2466 
2467  elem->SetAttribute(nc_elem.attribute);
2468  for (int j = 0; j < gi.nv; j++)
2469  {
2470  elem->GetVertices()[j] = nodes[node[j]].vert_index;
2471  }
2472 
2473  // create boundary elements
2474  // TODO: use boundary_faces?
2475  for (int k = 0; k < gi.nf; k++)
2476  {
2477  const int* fv = gi.faces[k];
2478  const int nfv = gi.nfv[k];
2479  const Face* face = faces.Find(node[fv[0]], node[fv[1]],
2480  node[fv[2]], node[fv[3]]);
2481  if (face->Boundary())
2482  {
2483  if ((nc_elem.geom == Geometry::CUBE) ||
2484  ((nc_elem.geom == Geometry::PRISM ||
2485  nc_elem.geom == Geometry::PYRAMID) && nfv == 4))
2486  {
2487  auto* quad = (Quadrilateral*) mesh.NewElement(Geometry::SQUARE);
2488  quad->SetAttribute(face->attribute);
2489  for (int j = 0; j < 4; j++)
2490  {
2491  quad->GetVertices()[j] = nodes[node[fv[j]]].vert_index;
2492  }
2493  mesh.boundary.Append(quad);
2494  }
2495  else if (nc_elem.geom == Geometry::PRISM ||
2496  nc_elem.geom == Geometry::PYRAMID ||
2497  nc_elem.geom == Geometry::TETRAHEDRON)
2498  {
2499  MFEM_ASSERT(nfv == 3, "");
2500  auto* tri = (Triangle*) mesh.NewElement(Geometry::TRIANGLE);
2501  tri->SetAttribute(face->attribute);
2502  for (int j = 0; j < 3; j++)
2503  {
2504  tri->GetVertices()[j] = nodes[node[fv[j]]].vert_index;
2505  }
2506  mesh.boundary.Append(tri);
2507  }
2508  else if (nc_elem.geom == Geometry::SQUARE ||
2509  nc_elem.geom == Geometry::TRIANGLE)
2510  {
2511  auto* segment = (Segment*) mesh.NewElement(Geometry::SEGMENT);
2512  segment->SetAttribute(face->attribute);
2513  for (int j = 0; j < 2; j++)
2514  {
2515  segment->GetVertices()[j] = nodes[node[fv[2*j]]].vert_index;
2516  }
2517  mesh.boundary.Append(segment);
2518  }
2519  else
2520  {
2521  MFEM_ASSERT(nc_elem.geom == Geometry::SEGMENT, "");
2522  auto* point = (mfem::Point*) mesh.NewElement(Geometry::POINT);
2523  point->SetAttribute(face->attribute);
2524  point->GetVertices()[0] = nodes[node[fv[0]]].vert_index;
2525  mesh.boundary.Append(point);
2526  }
2527  }
2528  }
2529  }
2530 }
2531 
2533 {
2534  //// PART 1: pull indices of regular edges/faces from the Mesh
2535 
2536  NEdges = mesh->GetNEdges();
2537  NFaces = mesh->GetNumFaces();
2538  if (Dim < 2) { NFaces = 0; }
2539  // clear Node::edge_index and Face::index
2540  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2541  {
2542  if (node->HasEdge()) { node->edge_index = -1; }
2543  }
2544  for (auto face = faces.begin(); face != faces.end(); ++face)
2545  {
2546  face->index = -1;
2547  }
2548 
2549  // get edge enumeration from the Mesh
2550  Table *edge_vertex = mesh->GetEdgeVertexTable();
2551  for (int i = 0; i < edge_vertex->Size(); i++)
2552  {
2553  const int *ev = edge_vertex->GetRow(i);
2554  Node* node = nodes.Find(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
2555 
2556  MFEM_ASSERT(node && node->HasEdge(),
2557  "edge (" << ev[0] << "," << ev[1] << ") not found, "
2558  "node = " << node);
2559 
2560  node->edge_index = i;
2561  }
2562 
2563  // get face enumeration from the Mesh, initialize 'face_geom'
2565  for (int i = 0; i < NFaces; i++)
2566  {
2567  const int* fv = mesh->GetFace(i)->GetVertices();
2568  const int nfv = mesh->GetFace(i)->GetNVertices();
2569 
2570  Face* face;
2571  if (Dim == 3)
2572  {
2573  if (nfv == 4)
2574  {
2576  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2577  vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
2578  }
2579  else
2580  {
2581  MFEM_ASSERT(nfv == 3, "");
2583  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2584  vertex_nodeId[fv[2]]);
2585  }
2586  }
2587  else
2588  {
2589  MFEM_ASSERT(nfv == 2, "");
2591  int n0 = vertex_nodeId[fv[0]], n1 = vertex_nodeId[fv[1]];
2592  face = faces.Find(n0, n0, n1, n1); // look up degenerate face
2593 
2594 #ifdef MFEM_DEBUG
2595  // (non-ghost) edge and face numbers must match in 2D
2596  const int *ev = edge_vertex->GetRow(i);
2597  MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
2598  (ev[1] == fv[0] && ev[0] == fv[1]), "");
2599 #endif
2600  }
2601 
2602  MFEM_VERIFY(face, "face not found.");
2603  face->index = i;
2604  }
2605 
2606  //// PART 2: assign indices of ghost edges/faces, if any
2607 
2608  // count ghost edges and assign their indices
2609  NGhostEdges = 0;
2610  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2611  {
2612  if (node->HasEdge() && node->edge_index < 0)
2613  {
2614  node->edge_index = NEdges + (NGhostEdges++);
2615  }
2616  }
2617 
2618  // count ghost faces
2619  NGhostFaces = 0;
2620  for (auto face = faces.begin(); face != faces.end(); ++face)
2621  {
2622  if (face->index < 0) { NGhostFaces++; }
2623  }
2624 
2625  if (Dim == 2)
2626  {
2627  // in 2D we have fake faces because of DG
2628  MFEM_ASSERT(NFaces == NEdges, "");
2629  MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
2630  }
2631 
2632  // resize face_geom (default_geom is for slave faces beyond the ghost layer)
2633  Geometry::Type default_geom = Geometry::SQUARE;
2634  face_geom.SetSize(NFaces + NGhostFaces, default_geom);
2635 
2636  // update 'face_geom' for ghost faces, assign ghost face indices
2637  int nghosts = 0;
2638  for (int i = 0; i < NGhostElements; i++)
2639  {
2640  Element &el = elements[leaf_elements[NElements + i]]; // ghost element
2641  GeomInfo &gi = GI[el.Geom()];
2642 
2643  for (int j = 0; j < gi.nf; j++)
2644  {
2645  const int *fv = gi.faces[j];
2646  Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
2647  el.node[fv[2]], el.node[fv[3]]);
2648  MFEM_ASSERT(face, "face not found!");
2649 
2650  if (face->index < 0)
2651  {
2652  face->index = NFaces + (nghosts++);
2653 
2654  // store the face geometry
2655  static const Geometry::Type types[5] =
2656  {
2659  };
2660  face_geom[face->index] = types[gi.nfv[j]];
2661  }
2662  }
2663  }
2664 
2665  // assign valid indices also to faces beyond the ghost layer
2666  for (auto face = faces.begin(); face != faces.end(); ++face)
2667  {
2668  if (face->index < 0) { face->index = NFaces + (nghosts++); }
2669  }
2670  MFEM_ASSERT(nghosts == NGhostFaces, "");
2671 }
2672 
2673 
2674 //// Face/edge lists ///////////////////////////////////////////////////////////
2675 
2676 int NCMesh::QuadFaceSplitType(int v1, int v2, int v3, int v4,
2677  int mid[5]) const
2678 {
2679  MFEM_ASSERT(Dim >= 3, "");
2680 
2681  // find edge nodes
2682  int e1 = FindMidEdgeNode(v1, v2);
2683  int e2 = FindMidEdgeNode(v2, v3);
2684  int e3 = (e1 >= 0 && nodes[e1].HasVertex()) ? FindMidEdgeNode(v3, v4) : -1;
2685  int e4 = (e2 >= 0 && nodes[e2].HasVertex()) ? FindMidEdgeNode(v4, v1) : -1;
2686 
2687  // optional: return the mid-edge nodes if requested
2688  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
2689 
2690  // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
2691  int midf1 = -1, midf2 = -1;
2692  if (e1 >= 0 && e3 >= 0) { midf1 = FindMidEdgeNode(e1, e3); }
2693  if (e2 >= 0 && e4 >= 0) { midf2 = FindMidEdgeNode(e2, e4); }
2694 
2695  // get proper node if shadow node exists
2696  if (midf1 >= 0 && midf1 == midf2)
2697  {
2698  const Node &nd = nodes[midf1];
2699  if (nd.p1 != e1 && nd.p2 != e1) { midf1 = -1; }
2700  if (nd.p1 != e2 && nd.p2 != e2) { midf2 = -1; }
2701  }
2702 
2703  // only one way to access the mid-face node must always exist
2704  MFEM_ASSERT(!(midf1 >= 0 && midf2 >= 0), "incorrectly split face!");
2705 
2706  if (midf1 < 0 && midf2 < 0) // face not split
2707  {
2708  if (mid) { mid[4] = -1; }
2709  return 0;
2710  }
2711  else if (midf1 >= 0) // face split "vertically"
2712  {
2713  if (mid) { mid[4] = midf1; }
2714  return 1;
2715  }
2716  else // face split "horizontally"
2717  {
2718  if (mid) { mid[4] = midf2; }
2719  return 2;
2720  }
2721 }
2722 
2723 bool NCMesh::TriFaceSplit(int v1, int v2, int v3, int mid[3]) const
2724 {
2725  int e1 = nodes.FindId(v1, v2);
2726  if (e1 < 0 || !nodes[e1].HasVertex()) { return false; }
2727 
2728  int e2 = nodes.FindId(v2, v3);
2729  if (e2 < 0 || !nodes[e2].HasVertex()) { return false; }
2730 
2731  int e3 = nodes.FindId(v3, v1);
2732  if (e3 < 0 || !nodes[e3].HasVertex()) { return false; }
2733 
2734  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3; }
2735 
2736  // NOTE: face (v1, v2, v3) still needs to be checked
2737  return true;
2738 }
2739 
2740 int NCMesh::find_node(const Element &el, int node)
2741 {
2742  for (int i = 0; i < MaxElemNodes; i++)
2743  {
2744  if (el.node[i] == node) { return i; }
2745  }
2746  MFEM_ABORT("Node not found.");
2747  return -1;
2748 }
2749 
2750 int NCMesh::FindNodeExt(const Element &el, int node, bool abort)
2751 {
2752  for (int i = 0; i < GI[el.Geom()].nv; i++)
2753  {
2754  if (RetrieveNode(el, i) == node) { return i; }
2755  }
2756  if (abort) { MFEM_ABORT("Node not found."); }
2757  return -1;
2758 }
2759 
2760 int NCMesh::find_element_edge(const Element &el, int vn0, int vn1, bool abort)
2761 {
2762  MFEM_ASSERT(!el.ref_type, "");
2763 
2764  GeomInfo &gi = GI[el.Geom()];
2765  for (int i = 0; i < gi.ne; i++)
2766  {
2767  const int* ev = gi.edges[i];
2768  int n0 = el.node[ev[0]];
2769  int n1 = el.node[ev[1]];
2770  if ((n0 == vn0 && n1 == vn1) ||
2771  (n0 == vn1 && n1 == vn0)) { return i; }
2772  }
2773 
2774  if (abort) { MFEM_ABORT("Edge (" << vn0 << ", " << vn1 << ") not found"); }
2775  return -1;
2776 }
2777 
2778 int NCMesh::find_local_face(int geom, int a, int b, int c)
2779 {
2780  GeomInfo &gi = GI[geom];
2781  for (int i = 0; i < gi.nf; i++)
2782  {
2783  const int* fv = gi.faces[i];
2784  if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
2785  (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
2786  (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
2787  {
2788  return i;
2789  }
2790  }
2791  MFEM_ABORT("Face not found.");
2792  return -1;
2793 }
2794 
2795 
2796 /// Hash function for a PointMatrix, used in MatrixMap::map.
2797 struct PointMatrixHash
2798 {
2799  std::size_t operator()(const NCMesh::PointMatrix &pm) const
2800  {
2801  MFEM_ASSERT(sizeof(double) == sizeof(std::uint64_t), "");
2802 
2803  // This is a variation on "Hashing an array of floats" from here:
2804  // https://cs.stackexchange.com/questions/37952
2805  std::uint64_t hash = 0xf9ca9ba106acbba9; // random initial value
2806  for (int i = 0; i < pm.np; i++)
2807  {
2808  for (int j = 0; j < pm.points[i].dim; j++)
2809  {
2810  // mix the doubles by adding their binary representations
2811  // many times over (note: 31 is 11111 in binary)
2812  double coord = pm.points[i].coord[j];
2813  hash = 31*hash + *((std::uint64_t*) &coord);
2814  }
2815  }
2816  return hash; // return the lowest bits of the huge sum
2817  }
2818 };
2819 
2820 /** Helper container to keep track of point matrices encountered during
2821  * face/edge traversal and to assign unique indices to them.
2822  */
2823 struct MatrixMap
2824 {
2825  int GetIndex(const NCMesh::PointMatrix &pm)
2826  {
2827  int &index = map[pm];
2828  if (!index) { index = map.size(); }
2829  return index - 1;
2830  }
2831 
2832  void ExportMatrices(Array<DenseMatrix*> &point_matrices) const
2833  {
2834  point_matrices.SetSize(map.size());
2835  for (const auto &pair : map)
2836  {
2837  DenseMatrix* mat = new DenseMatrix();
2838  pair.first.GetMatrix(*mat);
2839  point_matrices[pair.second - 1] = mat;
2840  }
2841  }
2842 
2843  void DumpBucketSizes() const
2844  {
2845  for (unsigned i = 0; i < map.bucket_count(); i++)
2846  {
2847  mfem::out << map.bucket_size(i) << " ";
2848  }
2849  }
2850 
2851 private:
2852  std::unordered_map<NCMesh::PointMatrix, int, PointMatrixHash> map;
2853 };
2854 
2855 
2856 int NCMesh::ReorderFacePointMat(int v0, int v1, int v2, int v3,
2857  int elem, const PointMatrix &pm,
2858  PointMatrix &reordered) const
2859 {
2860  const Element &el = elements[elem];
2861  int master[4] =
2862  {
2863  find_node(el, v0), find_node(el, v1), find_node(el, v2),
2864  (v3 >= 0) ? find_node(el, v3) : -1
2865  };
2866  int nfv = (v3 >= 0) ? 4 : 3;
2867 
2868  int local = find_local_face(el.Geom(), master[0], master[1], master[2]);
2869  const int* fv = GI[el.Geom()].faces[local];
2870 
2871  reordered.np = pm.np;
2872  for (int i = 0, j; i < nfv; i++)
2873  {
2874  for (j = 0; j < nfv; j++)
2875  {
2876  if (fv[i] == master[j])
2877  {
2878  reordered.points[i] = pm.points[j];
2879  break;
2880  }
2881  }
2882  MFEM_ASSERT(j != nfv, "node not found.");
2883  }
2884  return local;
2885 }
2886 
2887 void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
2888  const PointMatrix& pm, int level,
2889  Face* eface[4], MatrixMap &matrix_map)
2890 {
2891  if (level > 0)
2892  {
2893  // check if we made it to a face that is not split further
2894  Face* fa = faces.Find(vn0, vn1, vn2, vn3);
2895  if (fa)
2896  {
2897  // we have a slave face, add it to the list
2898  int elem = fa->GetSingleElement();
2899  face_list.slaves.Append(
2900  Slave(fa->index, elem, -1, Geometry::SQUARE));
2901  Slave &sl = face_list.slaves.Last();
2902 
2903  // reorder the point matrix according to slave face orientation
2904  PointMatrix pm_r;
2905  sl.local = ReorderFacePointMat(vn0, vn1, vn2, vn3, elem, pm, pm_r);
2906  sl.matrix = matrix_map.GetIndex(pm_r);
2907 
2908  eface[0] = eface[2] = fa;
2909  eface[1] = eface[3] = fa;
2910 
2911  return;
2912  }
2913  }
2914 
2915  // we need to recurse deeper
2916  int mid[5];
2917  int split = QuadFaceSplitType(vn0, vn1, vn2, vn3, mid);
2918 
2919  Face *ef[2][4];
2920  if (split == 1) // "X" split face
2921  {
2922  Point pmid0(pm(0), pm(1)), pmid2(pm(2), pm(3));
2923 
2924  TraverseQuadFace(vn0, mid[0], mid[2], vn3,
2925  PointMatrix(pm(0), pmid0, pmid2, pm(3)),
2926  level+1, ef[0], matrix_map);
2927 
2928  TraverseQuadFace(mid[0], vn1, vn2, mid[2],
2929  PointMatrix(pmid0, pm(1), pm(2), pmid2),
2930  level+1, ef[1], matrix_map);
2931 
2932  eface[1] = ef[1][1];
2933  eface[3] = ef[0][3];
2934  eface[0] = eface[2] = NULL;
2935  }
2936  else if (split == 2) // "Y" split face
2937  {
2938  Point pmid1(pm(1), pm(2)), pmid3(pm(3), pm(0));
2939 
2940  TraverseQuadFace(vn0, vn1, mid[1], mid[3],
2941  PointMatrix(pm(0), pm(1), pmid1, pmid3),
2942  level+1, ef[0], matrix_map);
2943 
2944  TraverseQuadFace(mid[3], mid[1], vn2, vn3,
2945  PointMatrix(pmid3, pmid1, pm(2), pm(3)),
2946  level+1, ef[1], matrix_map);
2947 
2948  eface[0] = ef[0][0];
2949  eface[2] = ef[1][2];
2950  eface[1] = eface[3] = NULL;
2951  }
2952 
2953  // check for a prism edge constrained by the master face
2954  if (HavePrisms() && mid[4] >= 0)
2955  {
2956  Node& enode = nodes[mid[4]];
2957  if (enode.HasEdge())
2958  {
2959  // process the edge only if it's not shared by slave faces
2960  // within this master face (i.e. the edge is "hidden")
2961  const int fi[3][2] = {{0, 0}, {1, 3}, {2, 0}};
2962  if (!ef[0][fi[split][0]] && !ef[1][fi[split][1]])
2963  {
2964  MFEM_ASSERT(enode.edge_refc == 1, "");
2965 
2966  MeshId buf[4];
2967  Array<MeshId> eid(buf, 4);
2968 
2969  (split == 1) ? FindEdgeElements(mid[0], vn1, vn2, mid[2], eid)
2970  /* */ : FindEdgeElements(mid[3], vn0, vn1, mid[1], eid);
2971 
2972  MFEM_ASSERT(eid.Size() > 0, "edge prism not found");
2973  MFEM_ASSERT(eid.Size() < 2, "non-unique edge prism");
2974 
2975  // create a slave face record with a degenerate point matrix
2976  face_list.slaves.Append(
2977  Slave(-1 - enode.edge_index,
2978  eid[0].element, eid[0].local, Geometry::SQUARE));
2979  Slave &sl = face_list.slaves.Last();
2980 
2981  if (split == 1)
2982  {
2983  Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
2984  int v1 = nodes[mid[0]].vert_index;
2985  int v2 = nodes[mid[2]].vert_index;
2986  sl.matrix =
2987  matrix_map.GetIndex(
2988  (v1 < v2) ? PointMatrix(mid0, mid2, mid2, mid0) :
2989  /* */ PointMatrix(mid2, mid0, mid0, mid2));
2990  }
2991  else
2992  {
2993  Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
2994  int v1 = nodes[mid[1]].vert_index;
2995  int v2 = nodes[mid[3]].vert_index;
2996  sl.matrix =
2997  matrix_map.GetIndex(
2998  (v1 < v2) ? PointMatrix(mid1, mid3, mid3, mid1) :
2999  /* */ PointMatrix(mid3, mid1, mid1, mid3));
3000  }
3001  }
3002  }
3003  }
3004 }
3005 
3006 void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
3007  MatrixMap &matrix_map)
3008 {
3009  int mid = nodes.FindId(vn0, vn1);
3010  if (mid < 0) { return; }
3011 
3012  const Node &nd = nodes[mid];
3013  if (nd.HasEdge())
3014  {
3015  // check if the edge is already a master in 'edge_list'
3016  int type;
3017  const MeshId &eid = edge_list.LookUp(nd.edge_index, &type);
3018  if (type == 1)
3019  {
3020  // in this case we need to add an edge-face constraint, because the
3021  // master edge is really a (face-)slave itself
3022 
3023  face_list.slaves.Append(
3024  Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
3025 
3026  int v0index = nodes[vn0].vert_index;
3027  int v1index = nodes[vn1].vert_index;
3028 
3029  face_list.slaves.Last().matrix =
3030  matrix_map.GetIndex((v0index < v1index) ? PointMatrix(p0, p1, p0)
3031  /* */ : PointMatrix(p1, p0, p1));
3032 
3033  return; // no need to continue deeper
3034  }
3035  }
3036 
3037  // recurse deeper
3038  Point pmid(p0, p1);
3039  TraverseTetEdge(vn0, mid, p0, pmid, matrix_map);
3040  TraverseTetEdge(mid, vn1, pmid, p1, matrix_map);
3041 }
3042 
3043 bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
3044  const PointMatrix& pm, int level,
3045  MatrixMap &matrix_map)
3046 {
3047  if (level > 0)
3048  {
3049  // check if we made it to a face that is not split further
3050  Face* fa = faces.Find(vn0, vn1, vn2);
3051  if (fa)
3052  {
3053  // we have a slave face, add it to the list
3054  int elem = fa->GetSingleElement();
3055  face_list.slaves.Append(
3056  Slave(fa->index, elem, -1, Geometry::TRIANGLE));
3057  Slave &sl = face_list.slaves.Last();
3058 
3059  // reorder the point matrix according to slave face orientation
3060  PointMatrix pm_r;
3061  sl.local = ReorderFacePointMat(vn0, vn1, vn2, -1, elem, pm, pm_r);
3062  sl.matrix = matrix_map.GetIndex(pm_r);
3063 
3064  return true;
3065  }
3066  }
3067 
3068  int mid[3];
3069  if (TriFaceSplit(vn0, vn1, vn2, mid))
3070  {
3071  Point pmid0(pm(0), pm(1)), pmid1(pm(1), pm(2)), pmid2(pm(2), pm(0));
3072  bool b[4];
3073 
3074  b[0] = TraverseTriFace(vn0, mid[0], mid[2],
3075  PointMatrix(pm(0), pmid0, pmid2),
3076  level+1, matrix_map);
3077 
3078  b[1] = TraverseTriFace(mid[0], vn1, mid[1],
3079  PointMatrix(pmid0, pm(1), pmid1),
3080  level+1, matrix_map);
3081 
3082  b[2] = TraverseTriFace(mid[2], mid[1], vn2,
3083  PointMatrix(pmid2, pmid1, pm(2)),
3084  level+1, matrix_map);
3085 
3086  b[3] = TraverseTriFace(mid[1], mid[2], mid[0],
3087  PointMatrix(pmid1, pmid2, pmid0),
3088  level+1, matrix_map);
3089 
3090  // traverse possible tet edges constrained by the master face
3091  if (HaveTets() && !b[3])
3092  {
3093  if (!b[1]) { TraverseTetEdge(mid[0],mid[1], pmid0,pmid1, matrix_map); }
3094  if (!b[2]) { TraverseTetEdge(mid[1],mid[2], pmid1,pmid2, matrix_map); }
3095  if (!b[0]) { TraverseTetEdge(mid[2],mid[0], pmid2,pmid0, matrix_map); }
3096  }
3097  }
3098 
3099  return false;
3100 }
3101 
3103 {
3104  face_list.Clear();
3105  if (Dim < 3) { return; }
3106 
3107  if (HaveTets()) { GetEdgeList(); } // needed by TraverseTetEdge()
3108 
3110 
3111  Array<char> processed_faces(faces.NumIds());
3112  processed_faces = 0;
3113 
3114  MatrixMap matrix_maps[Geometry::NumGeom];
3115 
3116  // visit faces of leaf elements
3117  for (int i = 0; i < leaf_elements.Size(); i++)
3118  {
3119  int elem = leaf_elements[i];
3120  Element &el = elements[elem];
3121  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3122 
3123  GeomInfo& gi = GI[el.Geom()];
3124  for (int j = 0; j < gi.nf; j++)
3125  {
3126  // get nodes for this face
3127  int node[4];
3128  for (int k = 0; k < 4; k++)
3129  {
3130  node[k] = el.node[gi.faces[j][k]];
3131  }
3132 
3133  int face = faces.FindId(node[0], node[1], node[2], node[3]);
3134  MFEM_ASSERT(face >= 0, "face not found!");
3135 
3136  // tell ParNCMesh about the face
3137  ElementSharesFace(elem, j, face);
3138 
3139  // have we already processed this face? skip if yes
3140  if (processed_faces[face]) { continue; }
3141  processed_faces[face] = 1;
3142 
3143  int fgeom = (node[3] >= 0) ? Geometry::SQUARE : Geometry::TRIANGLE;
3144 
3145  Face &fa = faces[face];
3146  if (fa.elem[0] >= 0 && fa.elem[1] >= 0)
3147  {
3148  // this is a conforming face, add it to the list
3149  face_list.conforming.Append(MeshId(fa.index, elem, j, fgeom));
3150  }
3151  else
3152  {
3153  // this is either a master face or a slave face, but we can't
3154  // tell until we traverse the face refinement 'tree'...
3155  int sb = face_list.slaves.Size();
3156  if (fgeom == Geometry::SQUARE)
3157  {
3158  Face* dummy[4];
3159  TraverseQuadFace(node[0], node[1], node[2], node[3],
3160  pm_quad_identity, 0, dummy, matrix_maps[fgeom]);
3161  }
3162  else
3163  {
3164  TraverseTriFace(node[0], node[1], node[2],
3165  pm_tri_identity, 0, matrix_maps[fgeom]);
3166  }
3167 
3168  int se = face_list.slaves.Size();
3169  if (sb < se)
3170  {
3171  // found slaves, so this is a master face; add it to the list
3172  face_list.masters.Append(
3173  Master(fa.index, elem, j, fgeom, sb, se));
3174 
3175  // also, set the master index for the slaves
3176  for (int ii = sb; ii < se; ii++)
3177  {
3178  face_list.slaves[ii].master = fa.index;
3179  }
3180  }
3181  }
3182 
3183  if (fa.Boundary()) { boundary_faces.Append(face); }
3184  }
3185  }
3186 
3187  // export unique point matrices
3188  for (int i = 0; i < Geometry::NumGeom; i++)
3189  {
3190  matrix_maps[i].ExportMatrices(face_list.point_matrices[i]);
3191  }
3192 }
3193 
3194 void NCMesh::TraverseEdge(int vn0, int vn1, double t0, double t1, int flags,
3195  int level, MatrixMap &matrix_map)
3196 {
3197  int mid = nodes.FindId(vn0, vn1);
3198  if (mid < 0) { return; }
3199 
3200  Node &nd = nodes[mid];
3201  if (nd.HasEdge() && level > 0)
3202  {
3203  // we have a slave edge, add it to the list
3204  edge_list.slaves.Append(Slave(nd.edge_index, -1, -1, Geometry::SEGMENT));
3205 
3206  Slave &sl = edge_list.slaves.Last();
3207  sl.matrix = matrix_map.GetIndex(PointMatrix(Point(t0), Point(t1)));
3208 
3209  // handle slave edge orientation
3210  sl.edge_flags = flags;
3211  int v0index = nodes[vn0].vert_index;
3212  int v1index = nodes[vn1].vert_index;
3213  if (v0index > v1index) { sl.edge_flags |= 2; }
3214  }
3215 
3216  // recurse deeper
3217  double tmid = (t0 + t1) / 2;
3218  TraverseEdge(vn0, mid, t0, tmid, flags, level+1, matrix_map);
3219  TraverseEdge(mid, vn1, tmid, t1, flags, level+1, matrix_map);
3220 }
3221 
3223 {
3224  edge_list.Clear();
3225  if (Dim < 3) { boundary_faces.SetSize(0); }
3226 
3227  Array<char> processed_edges(nodes.NumIds());
3228  processed_edges = 0;
3229 
3230  Array<int> edge_element(nodes.NumIds());
3231  Array<signed char> edge_local(nodes.NumIds());
3232  edge_local = -1;
3233 
3234  MatrixMap matrix_map;
3235 
3236  // visit edges of leaf elements
3237  for (int i = 0; i < leaf_elements.Size(); i++)
3238  {
3239  int elem = leaf_elements[i];
3240  Element &el = elements[elem];
3241  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3242 
3243  GeomInfo& gi = GI[el.Geom()];
3244  for (int j = 0; j < gi.ne; j++)
3245  {
3246  // get nodes for this edge
3247  const int* ev = gi.edges[j];
3248  int node[2] = { el.node[ev[0]], el.node[ev[1]] };
3249 
3250  int enode = nodes.FindId(node[0], node[1]);
3251  MFEM_ASSERT(enode >= 0, "edge node not found!");
3252 
3253  Node &nd = nodes[enode];
3254  MFEM_ASSERT(nd.HasEdge(), "edge not found!");
3255 
3256  // tell ParNCMesh about the edge
3257  ElementSharesEdge(elem, j, enode);
3258 
3259  // (2D only, store boundary faces)
3260  if (Dim <= 2)
3261  {
3262  int face = faces.FindId(node[0], node[0], node[1], node[1]);
3263  MFEM_ASSERT(face >= 0, "face not found!");
3264  if (faces[face].Boundary()) { boundary_faces.Append(face); }
3265  }
3266 
3267  // store element/local for later
3268  edge_element[nd.edge_index] = elem;
3269  edge_local[nd.edge_index] = j;
3270 
3271  // skip slave edges here, they will be reached from their masters
3272  if (GetEdgeMaster(enode) >= 0) { continue; }
3273 
3274  // have we already processed this edge? skip if yes
3275  if (processed_edges[enode]) { continue; }
3276  processed_edges[enode] = 1;
3277 
3278  // prepare edge interval for slave traversal, handle orientation
3279  double t0 = 0.0, t1 = 1.0;
3280  int v0index = nodes[node[0]].vert_index;
3281  int v1index = nodes[node[1]].vert_index;
3282  int flags = (v0index > v1index) ? 1 : 0;
3283 
3284  // try traversing the edge to find slave edges
3285  int sb = edge_list.slaves.Size();
3286  TraverseEdge(node[0], node[1], t0, t1, flags, 0, matrix_map);
3287 
3288  int se = edge_list.slaves.Size();
3289  if (sb < se)
3290  {
3291  // found slaves, this is a master face; add it to the list
3292  edge_list.masters.Append(
3293  Master(nd.edge_index, elem, j, Geometry::SEGMENT, sb, se));
3294 
3295  // also, set the master index for the slaves
3296  for (int ii = sb; ii < se; ii++)
3297  {
3298  edge_list.slaves[ii].master = nd.edge_index;
3299  }
3300  }
3301  else
3302  {
3303  // no slaves, this is a conforming edge
3304  edge_list.conforming.Append(MeshId(nd.edge_index, elem, j));
3305  }
3306  }
3307  }
3308 
3309  // fix up slave edge element/local
3310  for (int i = 0; i < edge_list.slaves.Size(); i++)
3311  {
3312  Slave &sl = edge_list.slaves[i];
3313  int local = edge_local[sl.index];
3314  if (local >= 0)
3315  {
3316  sl.local = local;
3317  sl.element = edge_element[sl.index];
3318  }
3319  }
3320 
3321  // export unique point matrices
3322  matrix_map.ExportMatrices(edge_list.point_matrices[Geometry::SEGMENT]);
3323 }
3324 
3326 {
3327  int total = NVertices + NGhostVertices;
3328 
3329  vertex_list.Clear();
3330  vertex_list.conforming.Reserve(total);
3331 
3332  Array<char> processed_vertices(total);
3333  processed_vertices = 0;
3334 
3335  // analogously to above, visit vertices of leaf elements
3336  for (int i = 0; i < leaf_elements.Size(); i++)
3337  {
3338  int elem = leaf_elements[i];
3339  Element &el = elements[elem];
3340 
3341  for (int j = 0; j < GI[el.Geom()].nv; j++)
3342  {
3343  int node = el.node[j];
3344  Node &nd = nodes[node];
3345 
3346  int index = nd.vert_index;
3347  if (index >= 0)
3348  {
3349  ElementSharesVertex(elem, j, node);
3350 
3351  if (processed_vertices[index]) { continue; }
3352  processed_vertices[index] = 1;
3353 
3354  vertex_list.conforming.Append(MeshId(index, elem, j));
3355  }
3356  }
3357  }
3358 }
3359 
3361  DenseMatrix &oriented_matrix) const
3362 {
3363  oriented_matrix = *(point_matrices[slave.Geom()][slave.matrix]);
3364 
3365  if (slave.edge_flags)
3366  {
3367  MFEM_ASSERT(oriented_matrix.Height() == 1 &&
3368  oriented_matrix.Width() == 2, "not an edge point matrix");
3369 
3370  if (slave.edge_flags & 1) // master inverted
3371  {
3372  oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
3373  oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
3374  }
3375  if (slave.edge_flags & 2) // slave inverted
3376  {
3377  std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
3378  }
3379  }
3380 }
3381 
3383 {
3384  conforming.DeleteAll();
3385  masters.DeleteAll();
3386  slaves.DeleteAll();
3387 
3388  for (int i = 0; i < Geometry::NumGeom; i++)
3389  {
3390  for (int j = 0; j < point_matrices[i].Size(); j++)
3391  {
3392  delete point_matrices[i][j];
3393  }
3394  point_matrices[i].DeleteAll();
3395  }
3396 
3397  inv_index.DeleteAll();
3398 }
3399 
3401 {
3402  return conforming.Size() + masters.Size() + slaves.Size();
3403 }
3404 
3405 const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
3406 {
3407  if (!inv_index.Size())
3408  {
3409  int max_index = -1;
3410  for (int i = 0; i < conforming.Size(); i++)
3411  {
3412  max_index = std::max(conforming[i].index, max_index);
3413  }
3414  for (int i = 0; i < masters.Size(); i++)
3415  {
3416  max_index = std::max(masters[i].index, max_index);
3417  }
3418  for (int i = 0; i < slaves.Size(); i++)
3419  {
3420  if (slaves[i].index < 0) { continue; }
3421  max_index = std::max(slaves[i].index, max_index);
3422  }
3423 
3424  inv_index.SetSize(max_index + 1);
3425  inv_index = -1;
3426 
3427  for (int i = 0; i < conforming.Size(); i++)
3428  {
3429  inv_index[conforming[i].index] = (i << 2);
3430  }
3431  for (int i = 0; i < masters.Size(); i++)
3432  {
3433  inv_index[masters[i].index] = (i << 2) + 1;
3434  }
3435  for (int i = 0; i < slaves.Size(); i++)
3436  {
3437  if (slaves[i].index < 0) { continue; }
3438  inv_index[slaves[i].index] = (i << 2) + 2;
3439  }
3440  }
3441 
3442  MFEM_ASSERT(index >= 0 && index < inv_index.Size(), "");
3443  int key = inv_index[index];
3444 
3445  if (!type)
3446  {
3447  MFEM_VERIFY(key >= 0, "entity not found.");
3448  }
3449  else // return entity type if requested, don't abort when not found
3450  {
3451  *type = (key >= 0) ? (key & 0x3) : -1;
3452 
3453  static MeshId invalid;
3454  if (*type < 0) { return invalid; } // not found
3455  }
3456 
3457  // return found entity MeshId
3458  switch (key & 0x3)
3459  {
3460  case 0: return conforming[key >> 2];
3461  case 1: return masters[key >> 2];
3462  case 2: return slaves[key >> 2];
3463  default: MFEM_ABORT("internal error"); return conforming[0];
3464  }
3465 }
3466 
3467 
3468 //// Neighbors /////////////////////////////////////////////////////////////////
3469 
3470 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
3471 {
3472  int mid = nodes.FindId(v0, v1);
3473  if (mid >= 0 && nodes[mid].HasVertex())
3474  {
3475  indices.Append(mid);
3476 
3477  CollectEdgeVertices(v0, mid, indices);
3478  CollectEdgeVertices(mid, v1, indices);
3479  }
3480 }
3481 
3482 void NCMesh::CollectTriFaceVertices(int v0, int v1, int v2, Array<int> &indices)
3483 {
3484  int mid[3];
3485  if (TriFaceSplit(v0, v1, v2, mid))
3486  {
3487  for (int i = 0; i < 3; i++)
3488  {
3489  indices.Append(mid[i]);
3490  }
3491 
3492  CollectTriFaceVertices(v0, mid[0], mid[2], indices);
3493  CollectTriFaceVertices(mid[0], v1, mid[1], indices);
3494  CollectTriFaceVertices(mid[2], mid[1], v2, indices);
3495  CollectTriFaceVertices(mid[0], mid[1], mid[2], indices);
3496 
3497  if (HaveTets()) // possible edge-face contact
3498  {
3499  CollectEdgeVertices(mid[0], mid[1], indices);
3500  CollectEdgeVertices(mid[1], mid[2], indices);
3501  CollectEdgeVertices(mid[2], mid[0], indices);
3502  }
3503  }
3504 }
3505 
3506 void NCMesh::CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
3507  Array<int> &indices)
3508 {
3509  int mid[5];
3510  switch (QuadFaceSplitType(v0, v1, v2, v3, mid))
3511  {
3512  case 1:
3513  indices.Append(mid[0]);
3514  indices.Append(mid[2]);
3515 
3516  CollectQuadFaceVertices(v0, mid[0], mid[2], v3, indices);
3517  CollectQuadFaceVertices(mid[0], v1, v2, mid[2], indices);
3518 
3519  if (HavePrisms()) // possible edge-face contact
3520  {
3521  CollectEdgeVertices(mid[0], mid[2], indices);
3522  }
3523  break;
3524 
3525  case 2:
3526  indices.Append(mid[1]);
3527  indices.Append(mid[3]);
3528 
3529  CollectQuadFaceVertices(v0, v1, mid[1], mid[3], indices);
3530  CollectQuadFaceVertices(mid[3], mid[1], v2, v3, indices);
3531 
3532  if (HavePrisms()) // possible edge-face contact
3533  {
3534  CollectEdgeVertices(mid[1], mid[3], indices);
3535  }
3536  break;
3537  }
3538 }
3539 
3541 {
3542  int nrows = leaf_elements.Size();
3543  int* I = Memory<int>(nrows + 1);
3544  int** JJ = new int*[nrows];
3545 
3546  Array<int> indices;
3547  indices.Reserve(128);
3548 
3549  // collect vertices coinciding with each element, including hanging vertices
3550  for (int i = 0; i < leaf_elements.Size(); i++)
3551  {
3552  int elem = leaf_elements[i];
3553  Element &el = elements[elem];
3554  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3555 
3556  GeomInfo& gi = GI[el.Geom()];
3557  int* node = el.node;
3558 
3559  indices.SetSize(0);
3560  for (int j = 0; j < gi.ne; j++)
3561  {
3562  const int* ev = gi.edges[j];
3563  CollectEdgeVertices(node[ev[0]], node[ev[1]], indices);
3564  }
3565 
3566  if (Dim >= 3)
3567  {
3568  for (int j = 0; j < gi.nf; j++)
3569  {
3570  const int* fv = gi.faces[j];
3571  if (gi.nfv[j] == 4)
3572  {
3573  CollectQuadFaceVertices(node[fv[0]], node[fv[1]],
3574  node[fv[2]], node[fv[3]], indices);
3575  }
3576  else
3577  {
3578  CollectTriFaceVertices(node[fv[0]], node[fv[1]], node[fv[2]],
3579  indices);
3580  }
3581  }
3582  }
3583 
3584  // temporarily store one row of the table
3585  indices.Sort();
3586  indices.Unique();
3587  int size = indices.Size();
3588  I[i] = size;
3589  JJ[i] = new int[size];
3590  std::memcpy(JJ[i], indices.GetData(), size * sizeof(int));
3591  }
3592 
3593  // finalize the I array of the table
3594  int nnz = 0;
3595  for (int i = 0; i < nrows; i++)
3596  {
3597  int cnt = I[i];
3598  I[i] = nnz;
3599  nnz += cnt;
3600  }
3601  I[nrows] = nnz;
3602 
3603  // copy the temporarily stored rows into one J array
3604  int *J = Memory<int>(nnz);
3605  nnz = 0;
3606  for (int i = 0; i < nrows; i++)
3607  {
3608  int cnt = I[i+1] - I[i];
3609  std::memcpy(J+nnz, JJ[i], cnt * sizeof(int));
3610  delete [] JJ[i];
3611  nnz += cnt;
3612  }
3613 
3614  element_vertex.SetIJ(I, J, nrows);
3615 
3616  delete [] JJ;
3617 }
3618 
3619 
3621  Array<int> *neighbors,
3622  Array<char> *neighbor_set)
3623 {
3624  // If A is the element-to-vertex table (see 'element_vertex') listing all
3625  // vertices touching each element, including hanging vertices, then A*A^T is
3626  // the element-to-neighbor table. Multiplying the element set with A*A^T
3627  // gives the neighbor set. To save memory, this function only computes the
3628  // action of A*A^T, the product itself is not stored anywhere.
3629 
3630  // Optimization: the 'element_vertex' table does not store the obvious
3631  // corner nodes in it. The table is therefore empty for conforming meshes.
3632 
3634 
3635  int nleaves = leaf_elements.Size();
3636  MFEM_VERIFY(elem_set.Size() == nleaves, "");
3637  MFEM_ASSERT(element_vertex.Size() == nleaves, "");
3638 
3639  // step 1: vertices = A^T * elem_set, i.e, find all vertices touching the
3640  // element set
3641 
3642  Array<char> vmark(nodes.NumIds());
3643  vmark = 0;
3644 
3645  for (int i = 0; i < nleaves; i++)
3646  {
3647  if (elem_set[i])
3648  {
3649  int *v = element_vertex.GetRow(i);
3650  int nv = element_vertex.RowSize(i);
3651  for (int j = 0; j < nv; j++)
3652  {
3653  vmark[v[j]] = 1;
3654  }
3655 
3656  Element &el = elements[leaf_elements[i]];
3657  nv = GI[el.Geom()].nv;
3658  for (int j = 0; j < nv; j++)
3659  {
3660  vmark[el.node[j]] = 1;
3661  }
3662  }
3663  }
3664 
3665  // step 2: neighbors = A * vertices, i.e., find all elements coinciding with
3666  // vertices from step 1; NOTE: in the result we don't include elements from
3667  // the original set
3668 
3669  if (neighbor_set)
3670  {
3671  neighbor_set->SetSize(nleaves);
3672  *neighbor_set = 0;
3673  }
3674 
3675  for (int i = 0; i < nleaves; i++)
3676  {
3677  if (!elem_set[i])
3678  {
3679  bool hit = false;
3680 
3681  int *v = element_vertex.GetRow(i);
3682  int nv = element_vertex.RowSize(i);
3683  for (int j = 0; j < nv; j++)
3684  {
3685  if (vmark[v[j]]) { hit = true; break; }
3686  }
3687 
3688  if (!hit)
3689  {
3690  Element &el = elements[leaf_elements[i]];
3691  nv = GI[el.Geom()].nv;
3692  for (int j = 0; j < nv; j++)
3693  {
3694  if (vmark[el.node[j]]) { hit = true; break; }
3695  }
3696  }
3697 
3698  if (hit)
3699  {
3700  if (neighbors) { neighbors->Append(leaf_elements[i]); }
3701  if (neighbor_set) { (*neighbor_set)[i] = 1; }
3702  }
3703  }
3704  }
3705 }
3706 
3707 static bool sorted_lists_intersect(const int* a, const int* b, int na, int nb)
3708 {
3709  // pointers to "end" sentinel, not last entry. Not for dereferencing.
3710  const int * const a_end = a + na;
3711  const int * const b_end = b + nb;
3712  while (a != a_end && b != b_end)
3713  {
3714  if (*a < *b)
3715  {
3716  ++a;
3717  }
3718  else if (*b < *a)
3719  {
3720  ++b;
3721  }
3722  else
3723  {
3724  return true; // neither *a < *b nor *b < *a thus a == b
3725  }
3726  }
3727  return false; // no common element found
3728 }
3729 
3730 
3731 void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
3732  const Array<int> *search_set)
3733 {
3734  // TODO future: this function is inefficient. For a single element, an
3735  // octree neighbor search algorithm would be better. However, the octree
3736  // neighbor algorithm is hard to get right in the multi-octree case due to
3737  // the different orientations of the octrees (i.e., the root elements).
3738 
3740 
3741  // sorted list of all vertex nodes touching 'elem'
3742  Array<int> vert;
3743  vert.Reserve(128);
3744 
3745  // support for non-leaf 'elem', collect vertices of all children
3746  Array<int> stack;
3747  stack.Reserve(64);
3748  stack.Append(elem);
3749 
3750  while (stack.Size())
3751  {
3752  Element &el = elements[stack.Last()];
3753  stack.DeleteLast();
3754 
3755  if (!el.ref_type)
3756  {
3757  int *v = element_vertex.GetRow(el.index);
3758  int nv = element_vertex.RowSize(el.index);
3759  for (int i = 0; i < nv; i++)
3760  {
3761  vert.Append(v[i]);
3762  }
3763 
3764  nv = GI[el.Geom()].nv;
3765  for (int i = 0; i < nv; i++)
3766  {
3767  vert.Append(el.node[i]);
3768  }
3769  }
3770  else
3771  {
3772  for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
3773  {
3774  stack.Append(el.child[i]);
3775  }
3776  }
3777  }
3778 
3779  vert.Sort();
3780  vert.Unique();
3781 
3782  int *v1 = vert.GetData();
3783  int nv1 = vert.Size();
3784 
3785  if (!search_set) { search_set = &leaf_elements; }
3786 
3787  // test *all* potential neighbors from the search set
3788  for (int i = 0; i < search_set->Size(); i++)
3789  {
3790  int testme = (*search_set)[i];
3791  if (testme != elem)
3792  {
3793  Element &el = elements[testme];
3794  int *v2 = element_vertex.GetRow(el.index);
3795  int nv2 = element_vertex.RowSize(el.index);
3796 
3797  bool hit = sorted_lists_intersect(v1, v2, nv1, nv2);
3798 
3799  if (!hit)
3800  {
3801  int nv = GI[el.Geom()].nv;
3802  for (int j = 0; j < nv; j++)
3803  {
3804  hit = sorted_lists_intersect(&el.node[j], v1, 1, nv1);
3805  if (hit) { break; }
3806  }
3807  }
3808 
3809  if (hit) { neighbors.Append(testme); }
3810  }
3811  }
3812 }
3813 
3815  Array<int> &expanded,
3816  const Array<int> *search_set)
3817 {
3819 
3820  Array<char> vmark(nodes.NumIds());
3821  vmark = 0;
3822 
3823  for (int i = 0; i < elems.Size(); i++)
3824  {
3825  Element &el = elements[elems[i]];
3826 
3827  int *v = element_vertex.GetRow(el.index);
3828  int nv = element_vertex.RowSize(el.index);
3829  for (int j = 0; j < nv; j++)
3830  {
3831  vmark[v[j]] = 1;
3832  }
3833 
3834  nv = GI[el.Geom()].nv;
3835  for (int j = 0; j < nv; j++)
3836  {
3837  vmark[el.node[j]] = 1;
3838  }
3839  }
3840 
3841  if (!search_set)
3842  {
3843  search_set = &leaf_elements;
3844  }
3845 
3846  expanded.SetSize(0);
3847  for (int i = 0; i < search_set->Size(); i++)
3848  {
3849  int testme = (*search_set)[i];
3850  Element &el = elements[testme];
3851  bool hit = false;
3852 
3853  int *v = element_vertex.GetRow(el.index);
3854  int nv = element_vertex.RowSize(el.index);
3855  for (int j = 0; j < nv; j++)
3856  {
3857  if (vmark[v[j]]) { hit = true; break; }
3858  }
3859 
3860  if (!hit)
3861  {
3862  nv = GI[el.Geom()].nv;
3863  for (int j = 0; j < nv; j++)
3864  {
3865  if (vmark[el.node[j]]) { hit = true; break; }
3866  }
3867  }
3868 
3869  if (hit) { expanded.Append(testme); }
3870  }
3871 }
3872 
3873 int NCMesh::GetVertexRootCoord(int elem, RefCoord coord[3]) const
3874 {
3875  while (1)
3876  {
3877  const Element &el = elements[elem];
3878  if (el.parent < 0) { return elem; }
3879 
3880  const Element &pa = elements[el.parent];
3881  MFEM_ASSERT(pa.ref_type, "internal error");
3882 
3883  int ch = 0;
3884  while (ch < MaxElemChildren && pa.child[ch] != elem) { ch++; }
3885  MFEM_ASSERT(ch < MaxElemChildren, "internal error");
3886 
3887  MFEM_ASSERT(geom_parent[el.Geom()], "unsupported geometry");
3888  const RefTrf &tr = geom_parent[el.Geom()][(int) pa.ref_type][ch];
3889  tr.Apply(coord, coord);
3890 
3891  elem = el.parent;
3892  }
3893 }
3894 
3895 static bool RefPointInside(Geometry::Type geom, const RefCoord pt[3])
3896 {
3897  switch (geom)
3898  {
3899  case Geometry::SQUARE:
3900  return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
3901  (pt[1] >= 0) && (pt[1] <= T_ONE);
3902 
3903  case Geometry::CUBE:
3904  return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
3905  (pt[1] >= 0) && (pt[1] <= T_ONE) &&
3906  (pt[2] >= 0) && (pt[2] <= T_ONE);
3907 
3908  case Geometry::TRIANGLE:
3909  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE);
3910 
3911  case Geometry::PRISM:
3912  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE) &&
3913  (pt[2] >= 0) && (pt[2] <= T_ONE);
3914 
3915  case Geometry::PYRAMID:
3916  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[2] >= 0.0) &&
3917  (pt[0] + pt[2] <= T_ONE) && (pt[1] + pt[2] <= T_ONE) &&
3918  (pt[2] <= T_ONE);
3919 
3920  default:
3921  MFEM_ABORT("unsupported geometry");
3922  return false;
3923  }
3924 }
3925 
3926 void NCMesh::CollectIncidentElements(int elem, const RefCoord coord[3],
3927  Array<int> &list) const
3928 {
3929  const Element &el = elements[elem];
3930  if (!el.ref_type)
3931  {
3932  list.Append(elem);
3933  return;
3934  }
3935 
3936  RefCoord tcoord[3];
3937  for (int ch = 0; ch < MaxElemChildren && el.child[ch] >= 0; ch++)
3938  {
3939  const RefTrf &tr = geom_child[el.Geom()][(int) el.ref_type][ch];
3940  tr.Apply(coord, tcoord);
3941 
3942  if (RefPointInside(el.Geom(), tcoord))
3943  {
3944  CollectIncidentElements(el.child[ch], tcoord, list);
3945  }
3946  }
3947 }
3948 
3949 void NCMesh::FindVertexCousins(int elem, int local, Array<int> &cousins) const
3950 {
3951  const Element &el = elements[elem];
3952 
3953  RefCoord coord[3];
3954  MFEM_ASSERT(geom_corners[el.Geom()], "unsupported geometry");
3955  std::memcpy(coord, geom_corners[el.Geom()][local], sizeof(coord));
3956 
3957  int root = GetVertexRootCoord(elem, coord);
3958 
3959  cousins.SetSize(0);
3960  CollectIncidentElements(root, coord, cousins);
3961 }
3962 
3963 
3964 //// Coarse/fine transformations ///////////////////////////////////////////////
3965 
3967 {
3968  MFEM_ASSERT(np == pm.np, "");
3969  for (int i = 0; i < np; i++)
3970  {
3971  MFEM_ASSERT(points[i].dim == pm.points[i].dim, "");
3972  for (int j = 0; j < points[i].dim; j++)
3973  {
3974  if (points[i].coord[j] != pm.points[i].coord[j]) { return false; }
3975  }
3976  }
3977  return true;
3978 }
3979 
3981 {
3982  point_matrix.SetSize(points[0].dim, np);
3983  for (int i = 0; i < np; i++)
3984  {
3985  for (int j = 0; j < points[0].dim; j++)
3986  {
3987  point_matrix(j, i) = points[i].coord[j];
3988  }
3989  }
3990 }
3991 
3993  Point(0), Point(1)
3994 );
3996  Point(0, 0), Point(1, 0), Point(0, 1)
3997 );
3999  Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)
4000 );
4002  Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)
4003 );
4005  Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0),
4006  Point(0, 0, 1), Point(1, 0, 1), Point(0, 1, 1)
4007 );
4009  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0),
4010  Point(0, 1, 0), Point(0, 0, 1)
4011 );
4013  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
4014  Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
4015 );
4016 
4018 {
4019  switch (geom)
4020  {
4021  case Geometry::SEGMENT: return pm_seg_identity;
4022  case Geometry::TRIANGLE: return pm_tri_identity;
4023  case Geometry::SQUARE: return pm_quad_identity;
4025  case Geometry::PRISM: return pm_prism_identity;
4027  case Geometry::CUBE: return pm_hex_identity;
4028  default:
4029  MFEM_ABORT("unsupported geometry " << geom);
4030  return pm_tri_identity;
4031  }
4032 }
4033 
4034 void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
4035  DenseMatrix& matrix)
4036 {
4037  PointMatrix pm = GetGeomIdentity(geom);
4038 
4039  while (*ref_path)
4040  {
4041  int ref_type = *ref_path++;
4042  int child = *ref_path++;
4043 
4044  // TODO: do this with the new child transform tables
4045 
4046  if (geom == Geometry::CUBE)
4047  {
4048  if (ref_type == 1) // split along X axis
4049  {
4050  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4051  Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
4052 
4053  if (child == 0)
4054  {
4055  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
4056  pm(4), mid45, mid67, pm(7));
4057  }
4058  else if (child == 1)
4059  {
4060  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
4061  mid45, pm(5), pm(6), mid67);
4062  }
4063  }
4064  else if (ref_type == 2) // split along Y axis
4065  {
4066  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4067  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4068 
4069  if (child == 0)
4070  {
4071  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4072  pm(4), pm(5), mid56, mid74);
4073  }
4074  else if (child == 1)
4075  {
4076  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4077  mid74, mid56, pm(6), pm(7));
4078  }
4079  }
4080  else if (ref_type == 4) // split along Z axis
4081  {
4082  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4083  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4084 
4085  if (child == 0)
4086  {
4087  pm = PointMatrix(pm(0), pm(1), pm(2), pm(3),
4088  mid04, mid15, mid26, mid37);
4089  }
4090  else if (child == 1)
4091  {
4092  pm = PointMatrix(mid04, mid15, mid26, mid37,
4093  pm(4), pm(5), pm(6), pm(7));
4094  }
4095  }
4096  else if (ref_type == 3) // XY split
4097  {
4098  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4099  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4100  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4101  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4102 
4103  Point midf0(mid23, mid12, mid01, mid30);
4104  Point midf5(mid45, mid56, mid67, mid74);
4105 
4106  if (child == 0)
4107  {
4108  pm = PointMatrix(pm(0), mid01, midf0, mid30,
4109  pm(4), mid45, midf5, mid74);
4110  }
4111  else if (child == 1)
4112  {
4113  pm = PointMatrix(mid01, pm(1), mid12, midf0,
4114  mid45, pm(5), mid56, midf5);
4115  }
4116  else if (child == 2)
4117  {
4118  pm = PointMatrix(midf0, mid12, pm(2), mid23,
4119  midf5, mid56, pm(6), mid67);
4120  }
4121  else if (child == 3)
4122  {
4123  pm = PointMatrix(mid30, midf0, mid23, pm(3),
4124  mid74, midf5, mid67, pm(7));
4125  }
4126  }
4127  else if (ref_type == 5) // XZ split
4128  {
4129  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4130  Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
4131  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4132  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4133 
4134  Point midf1(mid01, mid15, mid45, mid04);
4135  Point midf3(mid23, mid37, mid67, mid26);
4136 
4137  if (child == 0)
4138  {
4139  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
4140  mid04, midf1, midf3, mid37);
4141  }
4142  else if (child == 1)
4143  {
4144  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
4145  midf1, mid15, mid26, midf3);
4146  }
4147  else if (child == 2)
4148  {
4149  pm = PointMatrix(midf1, mid15, mid26, midf3,
4150  mid45, pm(5), pm(6), mid67);
4151  }
4152  else if (child == 3)
4153  {
4154  pm = PointMatrix(mid04, midf1, midf3, mid37,
4155  pm(4), mid45, mid67, pm(7));
4156  }
4157  }
4158  else if (ref_type == 6) // YZ split
4159  {
4160  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4161  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4162  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4163  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4164 
4165  Point midf2(mid12, mid26, mid56, mid15);
4166  Point midf4(mid30, mid04, mid74, mid37);
4167 
4168  if (child == 0)
4169  {
4170  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4171  mid04, mid15, midf2, midf4);
4172  }
4173  else if (child == 1)
4174  {
4175  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4176  midf4, midf2, mid26, mid37);
4177  }
4178  else if (child == 2)
4179  {
4180  pm = PointMatrix(mid04, mid15, midf2, midf4,
4181  pm(4), pm(5), mid56, mid74);
4182  }
4183  else if (child == 3)
4184  {
4185  pm = PointMatrix(midf4, midf2, mid26, mid37,
4186  mid74, mid56, pm(6), pm(7));
4187  }
4188  }
4189  else if (ref_type == 7) // full isotropic refinement
4190  {
4191  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4192  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4193  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4194  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4195  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4196  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4197 
4198  Point midf0(mid23, mid12, mid01, mid30);
4199  Point midf1(mid01, mid15, mid45, mid04);
4200  Point midf2(mid12, mid26, mid56, mid15);
4201  Point midf3(mid23, mid37, mid67, mid26);
4202  Point midf4(mid30, mid04, mid74, mid37);
4203  Point midf5(mid45, mid56, mid67, mid74);
4204 
4205  Point midel(midf1, midf3);
4206 
4207  if (child == 0)
4208  {
4209  pm = PointMatrix(pm(0), mid01, midf0, mid30,
4210  mid04, midf1, midel, midf4);
4211  }
4212  else if (child == 1)
4213  {
4214  pm = PointMatrix(mid01, pm(1), mid12, midf0,
4215  midf1, mid15, midf2, midel);
4216  }
4217  else if (child == 2)
4218  {
4219  pm = PointMatrix(midf0, mid12, pm(2), mid23,
4220  midel, midf2, mid26, midf3);
4221  }
4222  else if (child == 3)
4223  {
4224  pm = PointMatrix(mid30, midf0, mid23, pm(3),
4225  midf4, midel, midf3, mid37);
4226  }
4227  else if (child == 4)
4228  {
4229  pm = PointMatrix(mid04, midf1, midel, midf4,
4230  pm(4), mid45, midf5, mid74);
4231  }
4232  else if (child == 5)
4233  {
4234  pm = PointMatrix(midf1, mid15, midf2, midel,
4235  mid45, pm(5), mid56, midf5);
4236  }
4237  else if (child == 6)
4238  {
4239  pm = PointMatrix(midel, midf2, mid26, midf3,
4240  midf5, mid56, pm(6), mid67);
4241  }
4242  else if (child == 7)
4243  {
4244  pm = PointMatrix(midf4, midel, midf3, mid37,
4245  mid74, midf5, mid67, pm(7));
4246  }
4247  }
4248  }
4249  else if (geom == Geometry::PRISM)
4250  {
4251  if (ref_type < 4) // XY split
4252  {
4253  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4254  Point mid20(pm(2), pm(0)), mid34(pm(3), pm(4));
4255  Point mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4256 
4257  if (child == 0)
4258  {
4259  pm = PointMatrix(pm(0), mid01, mid20, pm(3), mid34, mid53);
4260  }
4261  else if (child == 1)
4262  {
4263  pm = PointMatrix(mid01, pm(1), mid12, mid34, pm(4), mid45);
4264  }
4265  else if (child == 2)
4266  {
4267  pm = PointMatrix(mid20, mid12, pm(2), mid53, mid45, pm(5));
4268  }
4269  else if (child == 3)
4270  {
4271  pm = PointMatrix(mid12, mid20, mid01, mid45, mid53, mid34);
4272  }
4273  }
4274  else if (ref_type == 4) // Z split
4275  {
4276  Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4277 
4278  if (child == 0)
4279  {
4280  pm = PointMatrix(pm(0), pm(1), pm(2), mid03, mid14, mid25);
4281  }
4282  else if (child == 1)
4283  {
4284  pm = PointMatrix(mid03, mid14, mid25, pm(3), pm(4), pm(5));
4285  }
4286  }
4287  else // ref_type > 4, iso split
4288  {
4289  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4290  Point mid34(pm(3), pm(4)), mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4291  Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4292 
4293  Point midf2(mid01, mid14, mid34, mid03);
4294  Point midf3(mid12, mid25, mid45, mid14);
4295  Point midf4(mid20, mid03, mid53, mid25);
4296 
4297  if (child == 0)
4298  {
4299  pm = PointMatrix(pm(0), mid01, mid20, mid03, midf2, midf4);
4300  }
4301  else if (child == 1)
4302  {
4303  pm = PointMatrix(mid01, pm(1), mid12, midf2, mid14, midf3);
4304  }
4305  else if (child == 2)
4306  {
4307  pm = PointMatrix(mid20, mid12, pm(2), midf4, midf3, mid25);
4308  }
4309  else if (child == 3)
4310  {
4311  pm = PointMatrix(mid12, mid20, mid01, midf3, midf4, midf2);
4312  }
4313  else if (child == 4)
4314  {
4315  pm = PointMatrix(mid03, midf2, midf4, pm(3), mid34, mid53);
4316  }
4317  else if (child == 5)
4318  {
4319  pm = PointMatrix(midf2, mid14, midf3, mid34, pm(4), mid45);
4320  }
4321  else if (child == 6)
4322  {
4323  pm = PointMatrix(midf4, midf3, mid25, mid53, mid45, pm(5));
4324  }
4325  else if (child == 7)
4326  {
4327  pm = PointMatrix(midf3, midf4, midf2, mid45, mid53, mid34);
4328  }
4329  }
4330  }
4331  else if (geom == Geometry::PYRAMID)
4332  {
4333  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4334  Point mid03(pm(0), pm(3)), mid12(pm(1), pm(2));
4335  Point mid04(pm(0), pm(4)), mid14(pm(1), pm(4));
4336  Point mid24(pm(2), pm(4)), mid34(pm(3), pm(4));
4337  Point midf0(mid23, mid12, mid01, mid03);
4338 
4339  if (child == 0) // Pyramid
4340  {
4341  pm = PointMatrix(pm(0), mid01, midf0, mid03, mid04);
4342  }
4343  else if (child == 1) // Pyramid
4344  {
4345  pm = PointMatrix(mid01, pm(1), mid12, midf0, mid14);
4346  }
4347  else if (child == 2) // Pyramid
4348  {
4349  pm = PointMatrix(midf0, mid12, pm(2), mid23, mid24);
4350  }
4351  else if (child == 3) // Pyramid
4352  {
4353  pm = PointMatrix(mid03, midf0, mid23, pm(3), mid34);
4354  }
4355  else if (child == 4) // Pyramid
4356  {
4357  pm = PointMatrix(mid24, mid14, mid04, mid34, midf0);
4358  }
4359  else if (child == 5) // Pyramid
4360  {
4361  pm = PointMatrix(mid04, mid14, mid24, mid34, pm(4));
4362  }
4363  else if (child == 6) // Tet
4364  {
4365  pm = PointMatrix(mid01, midf0, mid04, mid14);
4366  }
4367  else if (child == 7) // Tet
4368  {
4369  pm = PointMatrix(midf0, mid14, mid12, mid24);
4370  }
4371  else if (child == 8) // Tet
4372  {
4373  pm = PointMatrix(midf0, mid23, mid34, mid24);
4374  }
4375  else if (child == 9) // Tet
4376  {
4377  pm = PointMatrix(mid03, mid04, midf0, mid34);
4378  }
4379  }
4380  else if (geom == Geometry::TETRAHEDRON)
4381  {
4382  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid02(pm(2), pm(0));
4383  Point mid03(pm(0), pm(3)), mid13(pm(1), pm(3)), mid23(pm(2), pm(3));
4384 
4385  if (child == 0)
4386  {
4387  pm = PointMatrix(pm(0), mid01, mid02, mid03);
4388  }
4389  else if (child == 1)
4390  {
4391  pm = PointMatrix(mid01, pm(1), mid12, mid13);
4392  }
4393  else if (child == 2)
4394  {
4395  pm = PointMatrix(mid02, mid12, pm(2), mid23);
4396  }
4397  else if (child == 3)
4398  {
4399  pm = PointMatrix(mid03, mid13, mid23, pm(3));
4400  }
4401  else if (child == 4)
4402  {
4403  pm = PointMatrix(mid01, mid23, mid02, mid03);
4404  }
4405  else if (child == 5)
4406  {
4407  pm = PointMatrix(mid01, mid23, mid03, mid13);
4408  }
4409  else if (child == 6)
4410  {
4411  pm = PointMatrix(mid01, mid23, mid13, mid12);
4412  }
4413  else if (child == 7)
4414  {
4415  pm = PointMatrix(mid01, mid23, mid12, mid02);
4416  }
4417  }
4418  else if (geom == Geometry::SQUARE)
4419  {
4420  if (ref_type == 1) // X split
4421  {
4422  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4423 
4424  if (child == 0)
4425  {
4426  pm = PointMatrix(pm(0), mid01, mid23, pm(3));
4427  }
4428  else if (child == 1)
4429  {
4430  pm = PointMatrix(mid01, pm(1), pm(2), mid23);
4431  }
4432  }
4433  else if (ref_type == 2) // Y split
4434  {
4435  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4436 
4437  if (child == 0)
4438  {
4439  pm = PointMatrix(pm(0), pm(1), mid12, mid30);
4440  }
4441  else if (child == 1)
4442  {
4443  pm = PointMatrix(mid30, mid12, pm(2), pm(3));
4444  }
4445  }
4446  else if (ref_type == 3) // iso split
4447  {
4448  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4449  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4450  Point midel(mid01, mid23);
4451 
4452  if (child == 0)
4453  {
4454  pm = PointMatrix(pm(0), mid01, midel, mid30);
4455  }
4456  else if (child == 1)
4457  {
4458  pm = PointMatrix(mid01, pm(1), mid12, midel);
4459  }
4460  else if (child == 2)
4461  {
4462  pm = PointMatrix(midel, mid12, pm(2), mid23);
4463  }
4464  else if (child == 3)
4465  {
4466  pm = PointMatrix(mid30, midel, mid23, pm(3));
4467  }
4468  }
4469  }
4470  else if (geom == Geometry::TRIANGLE)
4471  {
4472  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4473 
4474  if (child == 0)
4475  {
4476  pm = PointMatrix(pm(0), mid01, mid20);
4477  }
4478  else if (child == 1)
4479  {
4480  pm = PointMatrix(mid01, pm(1), mid12);
4481  }
4482  else if (child == 2)
4483  {
4484  pm = PointMatrix(mid20, mid12, pm(2));
4485  }
4486  else if (child == 3)
4487  {
4488  pm = PointMatrix(mid12, mid20, mid01);
4489  }
4490  }
4491  else if (geom == Geometry::SEGMENT)
4492  {
4493  Point mid01(pm(0), pm(1));
4494 
4495  if (child == 0)
4496  {
4497  pm = PointMatrix(pm(0), mid01);
4498  }
4499  else if (child == 1)
4500  {
4501  pm = PointMatrix(mid01, pm(1));
4502  }
4503  }
4504  }
4505 
4506  // write the points to the matrix
4507  for (int i = 0; i < pm.np; i++)
4508  {
4509  for (int j = 0; j < pm(i).dim; j++)
4510  {
4511  matrix(j, i) = pm(i).coord[j];
4512  }
4513  }
4514 }
4515 
4517 {
4520 
4521  for (int i = 0; i < leaf_elements.Size(); i++)
4522  {
4523  int elem = leaf_elements[i];
4524  if (!IsGhost(elements[elem])) { coarse_elements.Append(elem); }
4525  }
4526 
4527  transforms.embeddings.DeleteAll();
4528 }
4529 
4530 void NCMesh::TraverseRefinements(int elem, int coarse_index,
4531  std::string &ref_path, RefPathMap &map)
4532 {
4533  Element &el = elements[elem];
4534  if (!el.ref_type)
4535  {
4536  int &matrix = map[ref_path];
4537  if (!matrix) { matrix = map.size(); }
4538 
4539  Embedding &emb = transforms.embeddings[el.index];
4540  emb.parent = coarse_index;
4541  emb.matrix = matrix - 1;
4542  emb.geom = el.Geom();
4543  emb.ghost = IsGhost(el);
4544  }
4545  else
4546  {
4547  MFEM_ASSERT(el.tet_type == 0, "not implemented");
4548 
4549  ref_path.push_back(el.ref_type);
4550  ref_path.push_back(0);
4551 
4552  for (int i = 0; i < MaxElemChildren; i++)
4553  {
4554  if (el.child[i] >= 0)
4555  {
4556  ref_path[ref_path.length()-1] = i;
4557  TraverseRefinements(el.child[i], coarse_index, ref_path, map);
4558  }
4559  }
4560  ref_path.resize(ref_path.length()-2);
4561  }
4562 }
4563 
4565 {
4566  MFEM_VERIFY(coarse_elements.Size() || !leaf_elements.Size(),
4567  "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
4568  " and Refine().");
4569 
4570  if (!transforms.embeddings.Size())
4571  {
4572  transforms.Clear();
4573  transforms.embeddings.SetSize(NElements);
4574 
4575  std::string ref_path;
4576  ref_path.reserve(100);
4577 
4578  RefPathMap path_map[Geometry::NumGeom];
4579  for (int g = 0; g < Geometry::NumGeom; g++)
4580  {
4581  path_map[g][ref_path] = 1; // empty path == identity
4582  }
4583 
4584  int used_geoms = 0;
4585  for (int i = 0; i < coarse_elements.Size(); i++)
4586  {
4587  int geom = elements[coarse_elements[i]].geom;
4588  TraverseRefinements(coarse_elements[i], i, ref_path, path_map[geom]);
4589  used_geoms |= (1 << geom);
4590  }
4591 
4592  for (int g = 0; g < Geometry::NumGeom; g++)
4593  {
4594  if (used_geoms & (1 << g))
4595  {
4596  Geometry::Type geom = Geometry::Type(g);
4597  const PointMatrix &identity = GetGeomIdentity(geom);
4598 
4600  .SetSize(Dim, identity.np, path_map[g].size());
4601 
4602  // calculate the point matrices
4603  RefPathMap::iterator it;
4604  for (it = path_map[g].begin(); it != path_map[g].end(); ++it)
4605  {
4606  GetPointMatrix(geom, it->first.c_str(),
4607  transforms.point_matrices[g](it->second-1));
4608  }
4609  }
4610  }
4611  }
4612  return transforms;
4613 }
4614 
4616 {
4617  MFEM_VERIFY(transforms.embeddings.Size() || !leaf_elements.Size(),
4618  "GetDerefinementTransforms() must be preceded by Derefine().");
4619 
4620  if (!transforms.IsInitialized())
4621  {
4622  std::map<int, int> mat_no[Geometry::NumGeom];
4623  for (int g = 0; g < Geometry::NumGeom; g++)
4624  {
4625  mat_no[g][0] = 1; // 0 == identity
4626  }
4627 
4628  // assign numbers to the different matrices used
4629  for (int i = 0; i < transforms.embeddings.Size(); i++)
4630  {
4631  Embedding &emb = transforms.embeddings[i];
4632  int code = emb.matrix; // see SetDerefMatrixCodes()
4633  if (code)
4634  {
4635  int &matrix = mat_no[emb.geom][code];
4636  if (!matrix) { matrix = mat_no[emb.geom].size(); }
4637 
4638  emb.matrix = matrix - 1;
4639  }
4640  }
4641 
4642  for (int g = 0; g < Geometry::NumGeom; g++)
4643  {
4644  if (Geoms & (1 << g))
4645  {
4646  Geometry::Type geom = Geometry::Type(g);
4647  const PointMatrix &identity = GetGeomIdentity(geom);
4648 
4650  .SetSize(Dim, identity.np, mat_no[geom].size());
4651 
4652  // calculate point matrices
4653  for (auto it = mat_no[geom].begin(); it != mat_no[geom].end(); ++it)
4654  {
4655  char path[3] = { 0, 0, 0 };
4656 
4657  int code = it->first;
4658  if (code)
4659  {
4660  path[0] = code >> 4; // ref_type (see SetDerefMatrixCodes())
4661  path[1] = code & 0xf; // child
4662  }
4663 
4664  GetPointMatrix(geom, path,
4665  transforms.point_matrices[geom](it->second-1));
4666  }
4667  }
4668  }
4669  }
4670  return transforms;
4671 }
4672 
4674  bool want_ghosts) const
4675 {
4676  Array<Connection> conn;
4677  conn.Reserve(embeddings.Size());
4678 
4679  int max_parent = -1;
4680  for (int i = 0; i < embeddings.Size(); i++)
4681  {
4682  const Embedding &emb = embeddings[i];
4683  if ((emb.parent >= 0) &&
4684  (!emb.ghost || want_ghosts))
4685  {
4686  conn.Append(Connection(emb.parent, i));
4687  max_parent = std::max(emb.parent, max_parent);
4688  }
4689  }
4690 
4691  conn.Sort(); // NOTE: unique is not necessary
4692  coarse_to_fine.MakeFromList(max_parent+1, conn);
4693 }
4694 
4696 {
4698  transforms.Clear();
4699 }
4700 
4702 {
4703  for (int i = 0; i < Geometry::NumGeom; i++)
4704  {
4705  point_matrices[i].SetSize(0, 0, 0);
4706  }
4707  embeddings.DeleteAll();
4708 }
4709 
4711 {
4712  // return true if point matrices are present for any geometry
4713  for (int i = 0; i < Geometry::NumGeom; i++)
4714  {
4715  if (point_matrices[i].SizeK()) { return true; }
4716  }
4717  return false;
4718 }
4719 
4721 {
4722  for (int g = 0; g < Geometry::NumGeom; ++g)
4723  {
4724  a.point_matrices[g].Swap(b.point_matrices[g]);
4725  }
4726  Swap(a.embeddings, b.embeddings);
4727 }
4728 
4729 
4730 //// SFC Ordering //////////////////////////////////////////////////////////////
4731 
4732 static int sgn(int x)
4733 {
4734  return (x < 0) ? -1 : (x > 0) ? 1 : 0;
4735 }
4736 
4737 static void HilbertSfc2D(int x, int y, int ax, int ay, int bx, int by,
4738  Array<int> &coords)
4739 {
4740  int w = std::abs(ax + ay);
4741  int h = std::abs(bx + by);
4742 
4743  int dax = sgn(ax), day = sgn(ay); // unit major direction ("right")
4744  int dbx = sgn(bx), dby = sgn(by); // unit orthogonal direction ("up")
4745 
4746  if (h == 1) // trivial row fill
4747  {
4748  for (int i = 0; i < w; i++, x += dax, y += day)
4749  {
4750  coords.Append(x);
4751  coords.Append(y);
4752  }
4753  return;
4754  }
4755  if (w == 1) // trivial column fill
4756  {
4757  for (int i = 0; i < h; i++, x += dbx, y += dby)
4758  {
4759  coords.Append(x);
4760  coords.Append(y);
4761  }
4762  return;
4763  }
4764 
4765  int ax2 = ax/2, ay2 = ay/2;
4766  int bx2 = bx/2, by2 = by/2;
4767 
4768  int w2 = std::abs(ax2 + ay2);
4769  int h2 = std::abs(bx2 + by2);
4770 
4771  if (2*w > 3*h) // long case: split in two parts only
4772  {
4773  if ((w2 & 0x1) && (w > 2))
4774  {
4775  ax2 += dax, ay2 += day; // prefer even steps
4776  }
4777 
4778  HilbertSfc2D(x, y, ax2, ay2, bx, by, coords);
4779  HilbertSfc2D(x+ax2, y+ay2, ax-ax2, ay-ay2, bx, by, coords);
4780  }
4781  else // standard case: one step up, one long horizontal step, one step down
4782  {
4783  if ((h2 & 0x1) && (h > 2))
4784  {
4785  bx2 += dbx, by2 += dby; // prefer even steps
4786  }
4787 
4788  HilbertSfc2D(x, y, bx2, by2, ax2, ay2, coords);
4789  HilbertSfc2D(x+bx2, y+by2, ax, ay, bx-bx2, by-by2, coords);
4790  HilbertSfc2D(x+(ax-dax)+(bx2-dbx), y+(ay-day)+(by2-dby),
4791  -bx2, -by2, -(ax-ax2), -(ay-ay2), coords);
4792  }
4793 }
4794 
4795 static void HilbertSfc3D(int x, int y, int z,
4796  int ax, int ay, int az,
4797  int bx, int by, int bz,
4798  int cx, int cy, int cz,
4799  Array<int> &coords)
4800 {
4801  int w = std::abs(ax + ay + az);
4802  int h = std::abs(bx + by + bz);
4803  int d = std::abs(cx + cy + cz);
4804 
4805  int dax = sgn(ax), day = sgn(ay), daz = sgn(az); // unit major dir ("right")
4806  int dbx = sgn(bx), dby = sgn(by), dbz = sgn(bz); // unit ortho dir ("forward")
4807  int dcx = sgn(cx), dcy = sgn(cy), dcz = sgn(cz); // unit ortho dir ("up")
4808 
4809  // trivial row/column fills
4810  if (h == 1 && d == 1)
4811  {
4812  for (int i = 0; i < w; i++, x += dax, y += day, z += daz)
4813  {
4814  coords.Append(x);
4815  coords.Append(y);
4816  coords.Append(z);
4817  }
4818  return;
4819  }
4820  if (w == 1 && d == 1)
4821  {
4822  for (int i = 0; i < h; i++, x += dbx, y += dby, z += dbz)
4823  {
4824  coords.Append(x);
4825  coords.Append(y);
4826  coords.Append(z);
4827  }
4828  return;
4829  }
4830  if (w == 1 && h == 1)
4831  {
4832  for (int i = 0; i < d; i++, x += dcx, y += dcy, z += dcz)
4833  {
4834  coords.Append(x);
4835  coords.Append(y);
4836  coords.Append(z);
4837  }
4838  return;
4839  }
4840 
4841  int ax2 = ax/2, ay2 = ay/2, az2 = az/2;
4842  int bx2 = bx/2, by2 = by/2, bz2 = bz/2;
4843  int cx2 = cx/2, cy2 = cy/2, cz2 = cz/2;
4844 
4845  int w2 = std::abs(ax2 + ay2 + az2);
4846  int h2 = std::abs(bx2 + by2 + bz2);
4847  int d2 = std::abs(cx2 + cy2 + cz2);
4848 
4849  // prefer even steps
4850  if ((w2 & 0x1) && (w > 2))
4851  {
4852  ax2 += dax, ay2 += day, az2 += daz;
4853  }
4854  if ((h2 & 0x1) && (h > 2))
4855  {
4856  bx2 += dbx, by2 += dby, bz2 += dbz;
4857  }
4858  if ((d2 & 0x1) && (d > 2))
4859  {
4860  cx2 += dcx, cy2 += dcy, cz2 += dcz;
4861  }
4862 
4863  // wide case, split in w only
4864  if ((2*w > 3*h) && (2*w > 3*d))
4865  {
4866  HilbertSfc3D(x, y, z,
4867  ax2, ay2, az2,
4868  bx, by, bz,
4869  cx, cy, cz, coords);
4870 
4871  HilbertSfc3D(x+ax2, y+ay2, z+az2,
4872  ax-ax2, ay-ay2, az-az2,
4873  bx, by, bz,
4874  cx, cy, cz, coords);
4875  }
4876  // do not split in d
4877  else if (3*h > 4*d)
4878  {
4879  HilbertSfc3D(x, y, z,
4880  bx2, by2, bz2,
4881  cx, cy, cz,
4882  ax2, ay2, az2, coords);
4883 
4884  HilbertSfc3D(x+bx2, y+by2, z+bz2,
4885  ax, ay, az,
4886  bx-bx2, by-by2, bz-bz2,
4887  cx, cy, cz, coords);
4888 
4889  HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
4890  y+(ay-day)+(by2-dby),
4891  z+(az-daz)+(bz2-dbz),
4892  -bx2, -by2, -bz2,
4893  cx, cy, cz,
4894  -(ax-ax2), -(ay-ay2), -(az-az2), coords);
4895  }
4896  // do not split in h
4897  else if (3*d > 4*h)
4898  {
4899  HilbertSfc3D(x, y, z,
4900  cx2, cy2, cz2,
4901  ax2, ay2, az2,
4902  bx, by, bz, coords);
4903 
4904  HilbertSfc3D(x+cx2, y+cy2, z+cz2,
4905  ax, ay, az,
4906  bx, by, bz,
4907  cx-cx2, cy-cy2, cz-cz2, coords);
4908 
4909  HilbertSfc3D(x+(ax-dax)+(cx2-dcx),
4910  y+(ay-day)+(cy2-dcy),
4911  z+(az-daz)+(cz2-dcz),
4912  -cx2, -cy2, -cz2,
4913  -(ax-ax2), -(ay-ay2), -(az-az2),
4914  bx, by, bz, coords);
4915  }
4916  // regular case, split in all w/h/d
4917  else
4918  {
4919  HilbertSfc3D(x, y, z,
4920  bx2, by2, bz2,
4921  cx2, cy2, cz2,
4922  ax2, ay2, az2, coords);
4923 
4924  HilbertSfc3D(x+bx2, y+by2, z+bz2,
4925  cx, cy, cz,
4926  ax2, ay2, az2,
4927  bx-bx2, by-by2, bz-bz2, coords);
4928 
4929  HilbertSfc3D(x+(bx2-dbx)+(cx-dcx),
4930  y+(by2-dby)+(cy-dcy),
4931  z+(bz2-dbz)+(cz-dcz),
4932  ax, ay, az,
4933  -bx2, -by2, -bz2,
4934  -(cx-cx2), -(cy-cy2), -(cz-cz2), coords);
4935 
4936  HilbertSfc3D(x+(ax-dax)+bx2+(cx-dcx),
4937  y+(ay-day)+by2+(cy-dcy),
4938  z+(az-daz)+bz2+(cz-dcz),
4939  -cx, -cy, -cz,
4940  -(ax-ax2), -(ay-ay2), -(az-az2),
4941  bx-bx2, by-by2, bz-bz2, coords);
4942 
4943  HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
4944  y+(ay-day)+(by2-dby),
4945  z+(az-daz)+(bz2-dbz),
4946  -bx2, -by2, -bz2,
4947  cx2, cy2, cz2,
4948  -(ax-ax2), -(ay-ay2), -(az-az2), coords);
4949  }
4950 }
4951 
4952 void NCMesh::GridSfcOrdering2D(int width, int height, Array<int> &coords)
4953 {
4954  coords.SetSize(0);
4955  coords.Reserve(2*width*height);
4956 
4957  if (width >= height)
4958  {
4959  HilbertSfc2D(0, 0, width, 0, 0, height, coords);
4960  }
4961  else
4962  {
4963  HilbertSfc2D(0, 0, 0, height, width, 0, coords);
4964  }
4965 }
4966 
4967 void NCMesh::GridSfcOrdering3D(int width, int height, int depth,
4968  Array<int> &coords)
4969 {
4970  coords.SetSize(0);
4971  coords.Reserve(3*width*height*depth);
4972 
4973  if (width >= height && width >= depth)
4974  {
4975  HilbertSfc3D(0, 0, 0,
4976  width, 0, 0,
4977  0, height, 0,
4978  0, 0, depth, coords);
4979  }
4980  else if (height >= width && height >= depth)
4981  {
4982  HilbertSfc3D(0, 0, 0,
4983  0, height, 0,
4984  width, 0, 0,
4985  0, 0, depth, coords);
4986  }
4987  else // depth >= width && depth >= height
4988  {
4989  HilbertSfc3D(0, 0, 0,
4990  0, 0, depth,
4991  width, 0, 0,
4992  0, height, 0, coords);
4993  }
4994 }
4995 
4996 
4997 //// Utility ///////////////////////////////////////////////////////////////////
4998 
4999 void NCMesh::GetEdgeVertices(const MeshId &edge_id, int vert_index[2],
5000  bool oriented) const
5001 {
5002  const Element &el = elements[edge_id.element];
5003  const GeomInfo& gi = GI[el.Geom()];
5004  const int* ev = gi.edges[(int) edge_id.local];
5005 
5006  int n0 = el.node[ev[0]], n1 = el.node[ev[1]];
5007  if (n0 > n1) { std::swap(n0, n1); }
5008 
5009  vert_index[0] = nodes[n0].vert_index;
5010  vert_index[1] = nodes[n1].vert_index;
5011 
5012  if (oriented && vert_index[0] > vert_index[1])
5013  {
5014  std::swap(vert_index[0], vert_index[1]);
5015  }
5016 }
5017 
5019 {
5020  const Element &el = elements[edge_id.element];
5021  const GeomInfo& gi = GI[el.Geom()];
5022  const int* ev = gi.edges[(int) edge_id.local];
5023 
5024  int v0 = nodes[el.node[ev[0]]].vert_index;
5025  int v1 = nodes[el.node[ev[1]]].vert_index;
5026 
5027  return ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1])) ? -1 : 1;
5028 }
5029 
5031  int vert_index[4], int edge_index[4],
5032  int edge_orientation[4]) const
5033 {
5034  MFEM_ASSERT(Dim >= 3, "");
5035 
5036  const Element &el = elements[face_id.element];
5037  const GeomInfo& gi = GI[el.Geom()];
5038 
5039  const int *fv = gi.faces[(int) face_id.local];
5040  const int nfv = gi.nfv[(int) face_id.local];
5041 
5042  vert_index[3] = edge_index[3] = -1;
5043  edge_orientation[3] = 0;
5044 
5045  for (int i = 0; i < nfv; i++)
5046  {
5047  vert_index[i] = nodes[el.node[fv[i]]].vert_index;
5048  }
5049 
5050  for (int i = 0; i < nfv; i++)
5051  {
5052  int j = i+1;
5053  if (j >= nfv) { j = 0; }
5054 
5055  int n1 = el.node[fv[i]];
5056  int n2 = el.node[fv[j]];
5057 
5058  const Node* en = nodes.Find(n1, n2);
5059  MFEM_ASSERT(en != NULL, "edge not found.");
5060 
5061  edge_index[i] = en->edge_index;
5062  edge_orientation[i] = (vert_index[i] < vert_index[j]) ? 1 : -1;
5063  }
5064 
5065  return nfv;
5066 }
5067 
5068 int NCMesh::GetEdgeMaster(int node) const
5069 {
5070  MFEM_ASSERT(node >= 0, "edge node not found.");
5071  const Node &nd = nodes[node];
5072 
5073  int p1 = nd.p1, p2 = nd.p2;
5074  MFEM_ASSERT(p1 != p2, "invalid edge node.");
5075 
5076  const Node &n1 = nodes[p1], &n2 = nodes[p2];
5077 
5078  int n1p1 = n1.p1, n1p2 = n1.p2;
5079  int n2p1 = n2.p1, n2p2 = n2.p2;
5080 
5081  if ((n2p1 != n2p2) && (p1 == n2p1 || p1 == n2p2))
5082  {
5083  // n1 is parent of n2:
5084  // (n1)--(nd)--(n2)------(*)
5085  if (n2.HasEdge()) { return p2; }
5086  else { return GetEdgeMaster(p2); }
5087  }
5088 
5089  if ((n1p1 != n1p2) && (p2 == n1p1 || p2 == n1p2))
5090  {
5091  // n2 is parent of n1:
5092  // (n2)--(nd)--(n1)------(*)
5093  if (n1.HasEdge()) { return p1; }
5094  else { return GetEdgeMaster(p1); }
5095  }
5096 
5097  return -1;
5098 }
5099 
5100 int NCMesh::GetEdgeMaster(int v1, int v2) const
5101 {
5102  int node = nodes.FindId(vertex_nodeId[v1], vertex_nodeId[v2]);
5103  MFEM_ASSERT(node >= 0 && nodes[node].HasEdge(), "(v1, v2) is not an edge.");
5104 
5105  int master = GetEdgeMaster(node);
5106  return (master >= 0) ? nodes[master].edge_index : -1;
5107 }
5108 
5109 int NCMesh::GetElementDepth(int i) const
5110 {
5111  int elem = leaf_elements[i];
5112  int depth = 0, parent;
5113  while ((parent = elements[elem].parent) != -1)
5114  {
5115  elem = parent;
5116  depth++;
5117  }
5118  return depth;
5119 }
5120 
5122 {
5123  int elem = leaf_elements[i];
5124  int parent, reduction = 1;
5125  while ((parent = elements[elem].parent) != -1)
5126  {
5127  if (elements[parent].ref_type & 1) { reduction *= 2; }
5128  if (elements[parent].ref_type & 2) { reduction *= 2; }
5129  if (elements[parent].ref_type & 4) { reduction *= 2; }
5130  elem = parent;
5131  }
5132  return reduction;
5133 }
5134 
5136  Array<int> &face_indices,
5137  Array<int> &face_attribs) const
5138 {
5139  const Element &el = elements[leaf_elements[leaf_elem]];
5140  const GeomInfo& gi = GI[el.Geom()];
5141 
5142  face_indices.SetSize(gi.nf);
5143  face_attribs.SetSize(gi.nf);
5144 
5145  for (int i = 0; i < gi.nf; i++)
5146  {
5147  const int* fv = gi.faces[i];
5148  const Face *face = faces.Find(el.node[fv[0]], el.node[fv[1]],
5149  el.node[fv[2]], el.node[fv[3]]);
5150  MFEM_ASSERT(face, "face not found");
5151  face_indices[i] = face->index;
5152  face_attribs[i] = face->attribute;
5153  }
5154 }
5155 
5156 void NCMesh::FindFaceNodes(int face, int node[4])
5157 {
5158  // Obtain face nodes from one of its elements (note that face->p1, p2, p3
5159  // cannot be used directly since they are not in order and p4 is missing).
5160 
5161  Face &fa = faces[face];
5162 
5163  int elem = fa.elem[0];
5164  if (elem < 0) { elem = fa.elem[1]; }
5165  MFEM_ASSERT(elem >= 0, "Face has no elements?");
5166 
5167  Element &el = elements[elem];
5168  int f = find_local_face(el.Geom(),
5169  find_node(el, fa.p1),
5170  find_node(el, fa.p2),
5171  find_node(el, fa.p3));
5172 
5173  const int* fv = GI[el.Geom()].faces[f];
5174  for (int i = 0; i < 4; i++)
5175  {
5176  node[i] = el.node[fv[i]];
5177  }
5178 }
5179 
5180 void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
5181  Array<int> &bdr_vertices, Array<int> &bdr_edges)
5182 {
5183  bdr_vertices.SetSize(0);
5184  bdr_edges.SetSize(0);
5185 
5186  if (Dim == 3)
5187  {
5188  GetFaceList(); // make sure 'boundary_faces' is up to date
5189 
5190  for (int i = 0; i < boundary_faces.Size(); i++)
5191  {
5192  int face = boundary_faces[i];
5193  if (bdr_attr_is_ess[faces[face].attribute - 1])
5194  {
5195  int node[4];
5196  FindFaceNodes(face, node);
5197  int nfv = (node[3] < 0) ? 3 : 4;
5198 
5199  for (int j = 0; j < nfv; j++)
5200  {
5201  bdr_vertices.Append(nodes[node[j]].vert_index);
5202 
5203  int enode = nodes.FindId(node[j], node[(j+1) % nfv]);
5204  MFEM_ASSERT(enode >= 0 && nodes[enode].HasEdge(), "Edge not found.");
5205  bdr_edges.Append(nodes[enode].edge_index);
5206 
5207  while ((enode = GetEdgeMaster(enode)) >= 0)
5208  {
5209  // append master edges that may not be accessible from any
5210  // boundary element, this happens in 3D in re-entrant corners
5211  bdr_edges.Append(nodes[enode].edge_index);
5212  }
5213  }
5214  }
5215  }
5216  }
5217  else if (Dim == 2)
5218  {
5219  GetEdgeList(); // make sure 'boundary_faces' is up to date
5220 
5221  for (int i = 0; i < boundary_faces.Size(); i++)
5222  {
5223  int face = boundary_faces[i];
5224  Face &fc = faces[face];
5225  if (bdr_attr_is_ess[fc.attribute - 1])
5226  {
5227  bdr_vertices.Append(nodes[fc.p1].vert_index);
5228  bdr_vertices.Append(nodes[fc.p3].vert_index);
5229  }
5230  }
5231  }
5232 
5233  bdr_vertices.Sort();
5234  bdr_vertices.Unique();
5235 
5236  bdr_edges.Sort();
5237  bdr_edges.Unique();
5238 }
5239 
5240 static int max4(int a, int b, int c, int d)
5241 {
5242  return std::max(std::max(a, b), std::max(c, d));
5243 }
5244 static int max6(int a, int b, int c, int d, int e, int f)
5245 {
5246  return std::max(max4(a, b, c, d), std::max(e, f));
5247 }
5248 static int max8(int a, int b, int c, int d, int e, int f, int g, int h)
5249 {
5250  return std::max(max4(a, b, c, d), max4(e, f, g, h));
5251 }
5252 
5253 int NCMesh::EdgeSplitLevel(int vn1, int vn2) const
5254 {
5255  int mid = nodes.FindId(vn1, vn2);
5256  if (mid < 0 || !nodes[mid].HasVertex()) { return 0; }
5257  return 1 + std::max(EdgeSplitLevel(vn1, mid), EdgeSplitLevel(mid, vn2));
5258 }
5259 
5260 int NCMesh::TriFaceSplitLevel(int vn1, int vn2, int vn3) const
5261 {
5262  int mid[3];
5263  if (TriFaceSplit(vn1, vn2, vn3, mid) &&
5264  faces.FindId(vn1, vn2, vn3) < 0)
5265  {
5266  return 1 + max4(TriFaceSplitLevel(vn1, mid[0], mid[2]),
5267  TriFaceSplitLevel(mid[0], vn2, mid[1]),
5268  TriFaceSplitLevel(mid[2], mid[1], vn3),
5269  TriFaceSplitLevel(mid[0], mid[1], mid[2]));
5270  }
5271  else // not split
5272  {
5273  return 0;
5274  }
5275 }
5276 
5277 void NCMesh::QuadFaceSplitLevel(int vn1, int vn2, int vn3, int vn4,
5278  int& h_level, int& v_level) const
5279 {
5280  int hl1, hl2, vl1, vl2;
5281  int mid[5];
5282 
5283  switch (QuadFaceSplitType(vn1, vn2, vn3, vn4, mid))
5284  {
5285  case 0: // not split
5286  h_level = v_level = 0;
5287  break;
5288 
5289  case 1: // vertical
5290  QuadFaceSplitLevel(vn1, mid[0], mid[2], vn4, hl1, vl1);
5291  QuadFaceSplitLevel(mid[0], vn2, vn3, mid[2], hl2, vl2);
5292  h_level = std::max(hl1, hl2);
5293  v_level = std::max(vl1, vl2) + 1;
5294  break;
5295 
5296  default: // horizontal
5297  QuadFaceSplitLevel(vn1, vn2, mid[1], mid[3], hl1, vl1);
5298  QuadFaceSplitLevel(mid[3], mid[1], vn3, vn4, hl2, vl2);
5299  h_level = std::max(hl1, hl2) + 1;
5300  v_level = std::max(vl1, vl2);
5301  }
5302 }
5303 
5304 void NCMesh::CountSplits(int elem, int splits[3]) const
5305 {
5306  const Element &el = elements[elem];
5307  const int* node = el.node;
5308  GeomInfo& gi = GI[el.Geom()];
5309 
5310  int elevel[MaxElemEdges];
5311  for (int i = 0; i < gi.ne; i++)
5312  {
5313  const int* ev = gi.edges[i];
5314  elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
5315  }
5316 
5317  int flevel[MaxElemFaces][2];
5318  if (Dim >= 3)
5319  {
5320  for (int i = 0; i < gi.nf; i++)
5321  {
5322  const int* fv = gi.faces[i];
5323  if (gi.nfv[i] == 4)
5324  {
5325  QuadFaceSplitLevel(node[fv[0]], node[fv[1]],
5326  node[fv[2]], node[fv[3]],
5327  flevel[i][1], flevel[i][0]);
5328  }
5329  else
5330  {
5331  flevel[i][1] = 0;
5332  flevel[i][0] =
5333  TriFaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]]);
5334  }
5335  }
5336  }
5337 
5338  if (el.Geom() == Geometry::CUBE)
5339  {
5340  splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
5341  elevel[0], elevel[2], elevel[4], elevel[6]);
5342 
5343  splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
5344  elevel[1], elevel[3], elevel[5], elevel[7]);
5345 
5346  splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
5347  elevel[8], elevel[9], elevel[10], elevel[11]);
5348  }
5349  else if (el.Geom() == Geometry::PRISM)
5350  {
5351  splits[0] = splits[1] =
5352  std::max(
5353  max6(flevel[0][0], flevel[1][0], 0,
5354  flevel[2][0], flevel[3][0], flevel[4][0]),
5355  max6(elevel[0], elevel[1], elevel[2],
5356  elevel[3], elevel[4], elevel[5]));
5357 
5358  splits[2] = max6(flevel[2][1], flevel[3][1], flevel[4][1],
5359  elevel[6], elevel[7], elevel[8]);
5360  }
5361  else if (el.Geom() == Geometry::PYRAMID)
5362  {
5363  splits[0] = std::max(
5364  max6(flevel[0][0], flevel[1][0], 0,
5365  flevel[2][0], flevel[3][0], flevel[4][0]),
5366  max8(elevel[0], elevel[1], elevel[2],
5367  elevel[3], elevel[4], elevel[5],
5368  elevel[6], elevel[7]));
5369 
5370  splits[1] = splits[0];
5371  splits[2] = splits[0];
5372  }
5373  else if (el.Geom() == Geometry::TETRAHEDRON)
5374  {
5375  splits[0] = std::max(
5376  max4(flevel[0][0], flevel[1][0], flevel[2][0], flevel[3][0]),
5377  max6(elevel[0], elevel[1], elevel[2],
5378  elevel[3], elevel[4], elevel[5]));
5379 
5380  splits[1] = splits[0];
5381  splits[2] = splits[0];
5382  }
5383  else if (el.Geom() == Geometry::SQUARE)
5384  {
5385  splits[0] = std::max(elevel[0], elevel[2]);
5386  splits[1] = std::max(elevel[1], elevel[3]);
5387  }
5388  else if (el.Geom() == Geometry::TRIANGLE)
5389  {
5390  splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
5391  splits[1] = splits[0];
5392  }
5393  else
5394  {
5395  MFEM_ABORT("Unsupported element geometry.");
5396  }
5397 }
5398 
5399 void NCMesh::GetLimitRefinements(Array<Refinement> &refinements, int max_level)
5400 {
5401  for (int i = 0; i < leaf_elements.Size(); i++)
5402  {
5403  if (IsGhost(elements[leaf_elements[i]])) { break; } // TODO: NElements
5404 
5405  int splits[3];
5406  CountSplits(leaf_elements[i], splits);
5407 
5408  char ref_type = 0;
5409  for (int k = 0; k < Dim; k++)
5410  {
5411  if (splits[k] > max_level)
5412  {
5413  ref_type |= (1 << k);
5414  }
5415  }
5416 
5417  if (ref_type)
5418  {
5419  if (Iso)
5420  {
5421  // iso meshes should only be modified by iso refinements
5422  ref_type = 7;
5423  }
5424  refinements.Append(Refinement(i, ref_type));
5425  }
5426  }
5427 }
5428 
5429 void NCMesh::LimitNCLevel(int max_nc_level)
5430 {
5431  MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
5432 
5433  while (1)
5434  {
5435  Array<Refinement> refinements;
5436  GetLimitRefinements(refinements, max_nc_level);
5437 
5438  if (!refinements.Size()) { break; }
5439 
5440  Refine(refinements);
5441  }
5442 }
5443 
5444 
5445 //// I/O ////////////////////////////////////////////////////////////////////////
5446 
5447 int NCMesh::PrintVertexParents(std::ostream *os) const
5448 {
5449  if (!os)
5450  {
5451  // count vertex nodes with parents
5452  int nv = 0;
5453  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5454  {
5455  if (node->HasVertex() && node->p1 != node->p2) { nv++; }
5456  }
5457  return nv;
5458  }
5459  else
5460  {
5461  // print the relations
5462  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5463  {
5464  if (node->HasVertex() && node->p1 != node->p2)
5465  {
5466  MFEM_ASSERT(nodes[node->p1].HasVertex(), "");
5467  MFEM_ASSERT(nodes[node->p2].HasVertex(), "");
5468 
5469  (*os) << node.index() << " " << node->p1 << " " << node->p2 << "\n";
5470  }
5471  }
5472  return 0;
5473  }
5474 }
5475 
5476 void NCMesh::LoadVertexParents(std::istream &input)
5477 {
5478  int nv;
5479  input >> nv;
5480  while (nv--)
5481  {
5482  int id, p1, p2;
5483  input >> id >> p1 >> p2;
5484  MFEM_VERIFY(input, "problem reading vertex parents.");
5485 
5486  MFEM_VERIFY(nodes.IdExists(id), "vertex " << id << " not found.");
5487  MFEM_VERIFY(nodes.IdExists(p1), "parent " << p1 << " not found.");
5488  MFEM_VERIFY(nodes.IdExists(p2), "parent " << p2 << " not found.");
5489 
5490  int check = nodes.FindId(p1, p2);
5491  MFEM_VERIFY(check < 0, "parents (" << p1 << ", " << p2 << ") already "
5492  "assigned to node " << check);
5493 
5494  // assign new parents for the node
5495  nodes.Reparent(id, p1, p2);
5496  }
5497 }
5498 
5499 int NCMesh::PrintBoundary(std::ostream *os) const
5500 {
5501  static const int nfv2geom[5] =
5502  {
5505  };
5506  int deg = (Dim == 2) ? 2 : 1; // for degenerate faces in 2D
5507 
5508  int count = 0;
5509  for (int i = 0; i < elements.Size(); i++)
5510  {
5511  const Element &el = elements[i];
5512  if (!el.IsLeaf()) { continue; }
5513 
5514  GeomInfo& gi = GI[el.Geom()];
5515  for (int k = 0; k < gi.nf; k++)
5516  {
5517  const int* fv = gi.faces[k];
5518  const int nfv = gi.nfv[k];
5519  const Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
5520  el.node[fv[2]], el.node[fv[3]]);
5521  MFEM_ASSERT(face != NULL, "face not found");
5522  if (face->Boundary())
5523  {
5524  if (!os) { count++; continue; }
5525 
5526  (*os) << face->attribute << " " << nfv2geom[nfv];
5527  for (int j = 0; j < nfv; j++)
5528  {
5529  (*os) << " " << el.node[fv[j*deg]];
5530  }
5531  (*os) << "\n";
5532  }
5533  }
5534  }
5535  return count;
5536 }
5537 
5538 void NCMesh::LoadBoundary(std::istream &input)
5539 {
5540  int nb, attr, geom;
5541  input >> nb;
5542  for (int i = 0; i < nb; i++)
5543  {
5544  input >> attr >> geom;
5545 
5546  int v1, v2, v3, v4;
5547  if (geom == Geometry::SQUARE)
5548  {
5549  input >> v1 >> v2 >> v3 >> v4;
5550  Face* face = faces.Get(v1, v2, v3, v4);
5551  face->attribute = attr;
5552  }
5553  else if (geom == Geometry::TRIANGLE)
5554  {
5555  input >> v1 >> v2 >> v3;
5556  Face* face = faces.Get(v1, v2, v3);
5557  face->attribute = attr;
5558  }
5559  else if (geom == Geometry::SEGMENT)
5560  {
5561  input >> v1 >> v2;
5562  Face* face = faces.Get(v1, v1, v2, v2);
5563  face->attribute = attr;
5564  }
5565  else if (geom == Geometry::POINT)
5566  {
5567  input >> v1;
5568  Face* face = faces.Get(v1, v1, v1, v1);
5569  face->attribute = attr;
5570  }
5571  else
5572  {
5573  MFEM_ABORT("unsupported boundary element geometry: " << geom);
5574  }
5575  }
5576 }
5577 
5578 void NCMesh::PrintCoordinates(std::ostream &os) const
5579 {
5580  int nv = coordinates.Size()/3;
5581  os << nv << "\n";
5582  if (!nv) { return; }
5583 
5584  os << spaceDim << "\n";
5585  for (int i = 0; i < nv; i++)
5586  {
5587  os << coordinates[3*i];
5588  for (int j = 1; j < spaceDim; j++)
5589  {
5590  os << " " << coordinates[3*i + j];
5591  }
5592  os << "\n";
5593  }
5594 }
5595 
5596 void NCMesh::LoadCoordinates(std::istream &input)
5597 {
5598  int nv;
5599  input >> nv;
5600  if (!nv) { return; }
5601 
5602  input >> spaceDim;
5603 
5604  coordinates.SetSize(nv * 3);
5605  coordinates = 0.0;
5606 
5607  for (int i = 0; i < nv; i++)
5608  {
5609  for (int j = 0; j < spaceDim; j++)
5610  {
5611  input >> coordinates[3*i + j];
5612  MFEM_VERIFY(input.good(), "unexpected EOF");
5613  }
5614  }
5615 }
5616 
5618 {
5619  for (int i = 0; i < root_state.Size(); i++)
5620  {
5621  if (root_state[i]) { return false; }
5622  }
5623  return true;
5624 }
5625 
5626 void NCMesh::Print(std::ostream &os) const
5627 {
5628  os << "MFEM NC mesh v1.0\n\n"
5629  "# NCMesh supported geometry types:\n"
5630  "# SEGMENT = 1\n"
5631  "# TRIANGLE = 2\n"
5632  "# SQUARE = 3\n"
5633  "# TETRAHEDRON = 4\n"
5634  "# CUBE = 5\n"
5635  "# PRISM = 6\n"
5636  "# PYRAMID = 7\n";
5637 
5638  os << "\ndimension\n" << Dim << "\n";
5639 
5640 #ifndef MFEM_USE_MPI
5641  if (MyRank != 0) // don't print this section in serial: default rank is 0
5642 #endif
5643  {
5644  os << "\nrank\n" << MyRank << "\n";
5645  }
5646 
5647  os << "\n# rank attr geom ref_type nodes/children";
5648  os << "\nelements\n" << elements.Size() << "\n";
5649 
5650  for (int i = 0; i < elements.Size(); i++)
5651  {
5652  const Element &el = elements[i];
5653  os << el.rank << " " << el.attribute << " ";
5654  if (el.parent == -2) { os << "-1\n"; continue; } // unused element
5655 
5656  os << int(el.geom) << " " << int(el.ref_type);
5657  for (int j = 0; j < MaxElemNodes && el.node[j] >= 0; j++)
5658  {
5659  os << " " << el.node[j];
5660  }
5661  os << "\n";
5662  }
5663 
5664  int nb = PrintBoundary(NULL);
5665  if (nb)
5666  {
5667  os << "\n# attr geom nodes";
5668  os << "\nboundary\n" << nb << "\n";
5669 
5670  PrintBoundary(&os);
5671  }
5672 
5673  int nvp = PrintVertexParents(NULL);
5674  if (nvp)
5675  {
5676  os << "\n# vert_id p1 p2";
5677  os << "\nvertex_parents\n" << nvp << "\n";
5678 
5679  PrintVertexParents(&os);
5680  }
5681 
5682  if (!ZeroRootStates()) // root_state section is optional
5683  {
5684  os << "\n# root element orientation";
5685  os << "\nroot_state\n" << root_state.Size() << "\n";
5686 
5687  for (int i = 0; i < root_state.Size(); i++)
5688  {
5689  os << root_state[i] << "\n";
5690  }
5691  }
5692 
5693  if (coordinates.Size())
5694  {
5695  os << "\n# top-level node coordinates";
5696  os << "\ncoordinates\n";
5697 
5698  PrintCoordinates(os);
5699  }
5700  else
5701  {
5702  // 'nodes' will be printed one level up by Mesh::Printer()
5703  }
5704 }
5705 
5707 {
5708  // initialize Element::parent
5709  for (int i = 0; i < elements.Size(); i++)
5710  {
5711  Element &el = elements[i];
5712  if (el.ref_type)
5713  {
5714  for (int j = 0; j < MaxElemChildren && el.child[j] >= 0; j++)
5715  {
5716  int child = el.child[j];
5717  MFEM_VERIFY(child < elements.Size(), "invalid mesh file: "
5718  "element " << i << " references invalid child " << child);
5719  elements[child].parent = i;
5720  }
5721  }
5722  }
5723 
5724  // count the root elements
5725  int nroots = 0;
5726  while (nroots < elements.Size() &&
5727  elements[nroots].parent == -1)
5728  {
5729  nroots++;
5730  }
5731  MFEM_VERIFY(nroots, "invalid mesh file: no root elements found.");
5732 
5733  // check that only the first 'nroot' elements are roots (have no parent)
5734  for (int i = nroots; i < elements.Size(); i++)
5735  {
5736  MFEM_VERIFY(elements[i].parent != -1,
5737  "invalid mesh file: only the first M elements can be roots. "
5738  "Found element " << i << " with no parent.");
5739  }
5740 
5741  // default root state is zero
5742  root_state.SetSize(nroots);
5743  root_state = 0;
5744 }
5745 
5747 {
5748  int ntop = 0;
5749  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5750  {
5751  if (node->p1 == node->p2) { ntop = node.index() + 1; }
5752  }
5753  return ntop;
5754 }
5755 
5756 NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
5757  : spaceDim(0), MyRank(0), Iso(true), Legacy(false)
5758 {
5759  is_nc = 1;
5760  if (version == 1) // old MFEM mesh v1.1 format
5761  {
5762  LoadLegacyFormat(input, curved, is_nc);
5763  Legacy = true;
5764  return;
5765  }
5766 
5767  MFEM_ASSERT(version == 10, "");
5768  std::string ident;
5769  int count;
5770 
5771  // load dimension
5772  skip_comment_lines(input, '#');
5773  input >> ident;
5774  MFEM_VERIFY(ident == "dimension", "Invalid mesh file: " << ident);
5775  input >> Dim;
5776 
5777  // load rank, if present
5778  skip_comment_lines(input, '#');
5779  input >> ident;
5780  if (ident == "rank")
5781  {
5782  input >> MyRank;
5783  MFEM_VERIFY(MyRank >= 0, "Invalid rank");
5784 
5785  skip_comment_lines(input, '#');
5786  input >> ident;
5787  }
5788 
5789  // load file SFC version, if present (for future changes to SFC ordering)
5790  if (ident == "sfc_version")
5791  {
5792  int sfc_version; // TODO future: store as class member
5793  input >> sfc_version;
5794  MFEM_VERIFY(sfc_version == 0,
5795  "Unsupported mesh file SFC version (" << sfc_version << "). "
5796  "Please update MFEM.");
5797 
5798  skip_comment_lines(input, '#');
5799  input >> ident;
5800  }
5801 
5802  // load elements
5803  MFEM_VERIFY(ident == "elements", "Invalid mesh file: " << ident);
5804  input >> count;
5805  for (int i = 0; i < count; i++)
5806  {
5807  int rank, attr, geom, ref_type;
5808  input >> rank >> attr >> geom;
5809 
5810  Geometry::Type type = Geometry::Type(geom);
5811  elements.Append(Element(type, attr));
5812 
5813  MFEM_ASSERT(elements.Size() == i+1, "");
5814  Element &el = elements[i];
5815  el.rank = rank;
5816 
5817  if (geom >= 0)
5818  {
5819  CheckSupportedGeom(type);
5820  GI[geom].InitGeom(type);
5821 
5822  input >> ref_type;
5823  MFEM_VERIFY(ref_type >= 0 && ref_type < 8, "");
5824  el.ref_type = ref_type;
5825 
5826  if (ref_type) // refined element
5827  {
5828  for (int j = 0; j < ref_type_num_children[ref_type]; j++)
5829  {
5830  input >> el.child[j];
5831  }
5832  if (Dim == 3 && ref_type != 7) { Iso = false; }
5833  }
5834  else // leaf element
5835  {
5836  for (int j = 0; j < GI[geom].nv; j++)
5837  {
5838  int id;
5839  input >> id;