MFEM  v4.6.0
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)
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)
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)
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)
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  for (int i = 0; i < mesh.NumOfElements; i++)
2453  {
2454  mesh.FreeElement(mesh.elements[i]);
2455  }
2456  mesh.elements.SetSize(0);
2457 
2458  for (int i = 0; i < mesh.NumOfBdrElements; i++)
2459  {
2460  mesh.FreeElement(mesh.boundary[i]);
2461  }
2462  mesh.boundary.SetSize(0);
2463 
2464  // create an mfem::Element for each leaf Element
2465  for (int i = 0; i < NElements; i++)
2466  {
2467  const Element &nc_elem = elements[leaf_elements[i]];
2468 
2469  const int* node = nc_elem.node;
2470  GeomInfo& gi = GI[(int) nc_elem.geom];
2471 
2472  mfem::Element* elem = mesh.NewElement(nc_elem.geom);
2473  mesh.elements.Append(elem);
2474 
2475  elem->SetAttribute(nc_elem.attribute);
2476  for (int j = 0; j < gi.nv; j++)
2477  {
2478  elem->GetVertices()[j] = nodes[node[j]].vert_index;
2479  }
2480 
2481  // create boundary elements
2482  // TODO: use boundary_faces?
2483  for (int k = 0; k < gi.nf; k++)
2484  {
2485  const int* fv = gi.faces[k];
2486  const int nfv = gi.nfv[k];
2487  const Face* face = faces.Find(node[fv[0]], node[fv[1]],
2488  node[fv[2]], node[fv[3]]);
2489  if (face->Boundary())
2490  {
2491  if ((nc_elem.geom == Geometry::CUBE) ||
2492  ((nc_elem.geom == Geometry::PRISM ||
2493  nc_elem.geom == Geometry::PYRAMID) && nfv == 4))
2494  {
2495  auto* quad = (Quadrilateral*) mesh.NewElement(Geometry::SQUARE);
2496  quad->SetAttribute(face->attribute);
2497  for (int j = 0; j < 4; j++)
2498  {
2499  quad->GetVertices()[j] = nodes[node[fv[j]]].vert_index;
2500  }
2501  mesh.boundary.Append(quad);
2502  }
2503  else if (nc_elem.geom == Geometry::PRISM ||
2504  nc_elem.geom == Geometry::PYRAMID ||
2505  nc_elem.geom == Geometry::TETRAHEDRON)
2506  {
2507  MFEM_ASSERT(nfv == 3, "");
2508  auto* tri = (Triangle*) mesh.NewElement(Geometry::TRIANGLE);
2509  tri->SetAttribute(face->attribute);
2510  for (int j = 0; j < 3; j++)
2511  {
2512  tri->GetVertices()[j] = nodes[node[fv[j]]].vert_index;
2513  }
2514  mesh.boundary.Append(tri);
2515  }
2516  else if (nc_elem.geom == Geometry::SQUARE ||
2517  nc_elem.geom == Geometry::TRIANGLE)
2518  {
2519  auto* segment = (Segment*) mesh.NewElement(Geometry::SEGMENT);
2520  segment->SetAttribute(face->attribute);
2521  for (int j = 0; j < 2; j++)
2522  {
2523  segment->GetVertices()[j] = nodes[node[fv[2*j]]].vert_index;
2524  }
2525  mesh.boundary.Append(segment);
2526  }
2527  else
2528  {
2529  MFEM_ASSERT(nc_elem.geom == Geometry::SEGMENT, "");
2530  auto* point = (mfem::Point*) mesh.NewElement(Geometry::POINT);
2531  point->SetAttribute(face->attribute);
2532  point->GetVertices()[0] = nodes[node[fv[0]]].vert_index;
2533  mesh.boundary.Append(point);
2534  }
2535  }
2536  }
2537  }
2538 }
2539 
2541 {
2542  //// PART 1: pull indices of regular edges/faces from the Mesh
2543 
2544  NEdges = mesh->GetNEdges();
2545  NFaces = mesh->GetNumFaces();
2546  if (Dim < 2) { NFaces = 0; }
2547  // clear Node::edge_index and Face::index
2548  for (auto &node : nodes)
2549  {
2550  if (node.HasEdge()) { node.edge_index = -1; }
2551  }
2552  for (auto &face : faces)
2553  {
2554  face.index = -1;
2555  }
2556 
2557  // get edge enumeration from the Mesh
2558  Table *edge_vertex = mesh->GetEdgeVertexTable();
2559  for (int i = 0; i < edge_vertex->Size(); i++)
2560  {
2561  const int *ev = edge_vertex->GetRow(i);
2562  Node* node = nodes.Find(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
2563 
2564  MFEM_ASSERT(node && node->HasEdge(),
2565  "edge (" << ev[0] << "," << ev[1] << ") not found, "
2566  "node = " << node);
2567 
2568  node->edge_index = i;
2569  }
2570 
2571  // get face enumeration from the Mesh, initialize 'face_geom'
2573  for (int i = 0; i < NFaces; i++)
2574  {
2575  const int* fv = mesh->GetFace(i)->GetVertices();
2576  const int nfv = mesh->GetFace(i)->GetNVertices();
2577 
2578  Face* face;
2579  if (Dim == 3)
2580  {
2581  if (nfv == 4)
2582  {
2584  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2585  vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
2586  }
2587  else
2588  {
2589  MFEM_ASSERT(nfv == 3, "");
2591  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2592  vertex_nodeId[fv[2]]);
2593  }
2594  }
2595  else
2596  {
2597  MFEM_ASSERT(nfv == 2, "");
2599  int n0 = vertex_nodeId[fv[0]], n1 = vertex_nodeId[fv[1]];
2600  face = faces.Find(n0, n0, n1, n1); // look up degenerate face
2601 
2602 #ifdef MFEM_DEBUG
2603  // (non-ghost) edge and face numbers must match in 2D
2604  const int *ev = edge_vertex->GetRow(i);
2605  MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
2606  (ev[1] == fv[0] && ev[0] == fv[1]), "");
2607 #endif
2608  }
2609 
2610  MFEM_VERIFY(face, "face not found.");
2611  face->index = i;
2612  }
2613 
2614  //// PART 2: assign indices of ghost edges/faces, if any
2615 
2616  // count ghost edges and assign their indices
2617  NGhostEdges = 0;
2618  for (auto &node : nodes)
2619  {
2620  if (node.HasEdge() && node.edge_index < 0)
2621  {
2622  node.edge_index = NEdges + (NGhostEdges++);
2623  }
2624  }
2625 
2626  // count ghost faces
2627  NGhostFaces = 0;
2628  for (auto &face : faces)
2629  {
2630  if (face.index < 0) { NGhostFaces++; }
2631  }
2632 
2633  if (Dim == 2)
2634  {
2635  // in 2D we have fake faces because of DG
2636  MFEM_ASSERT(NFaces == NEdges, "");
2637  MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
2638  }
2639 
2640  // resize face_geom (default_geom is for slave faces beyond the ghost layer)
2641  Geometry::Type default_geom = Geometry::SQUARE;
2642  face_geom.SetSize(NFaces + NGhostFaces, default_geom);
2643 
2644  // update 'face_geom' for ghost faces, assign ghost face indices
2645  int nghosts = 0;
2646  for (int i = 0; i < NGhostElements; i++)
2647  {
2648  Element &el = elements[leaf_elements[NElements + i]]; // ghost element
2649  GeomInfo &gi = GI[el.Geom()];
2650 
2651  for (int j = 0; j < gi.nf; j++)
2652  {
2653  const int *fv = gi.faces[j];
2654  Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
2655  el.node[fv[2]], el.node[fv[3]]);
2656  MFEM_ASSERT(face, "face not found!");
2657 
2658  if (face->index < 0)
2659  {
2660  face->index = NFaces + (nghosts++);
2661 
2662  // store the face geometry
2663  static const Geometry::Type types[5] =
2664  {
2667  };
2668  face_geom[face->index] = types[gi.nfv[j]];
2669  }
2670  }
2671  }
2672 
2673  // assign valid indices also to faces beyond the ghost layer
2674  for (auto &face : faces)
2675  {
2676  if (face.index < 0) { face.index = NFaces + (nghosts++); }
2677  }
2678  MFEM_ASSERT(nghosts == NGhostFaces, "");
2679 }
2680 
2681 
2682 //// Face/edge lists ///////////////////////////////////////////////////////////
2683 
2684 int NCMesh::QuadFaceSplitType(int v1, int v2, int v3, int v4,
2685  int mid[5]) const
2686 {
2687  MFEM_ASSERT(Dim >= 3, "");
2688 
2689  // find edge nodes
2690  int e1 = FindMidEdgeNode(v1, v2);
2691  int e2 = FindMidEdgeNode(v2, v3);
2692  int e3 = (e1 >= 0 && nodes[e1].HasVertex()) ? FindMidEdgeNode(v3, v4) : -1;
2693  int e4 = (e2 >= 0 && nodes[e2].HasVertex()) ? FindMidEdgeNode(v4, v1) : -1;
2694 
2695  // optional: return the mid-edge nodes if requested
2696  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
2697 
2698  // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
2699  int midf1 = -1, midf2 = -1;
2700  if (e1 >= 0 && e3 >= 0) { midf1 = FindMidEdgeNode(e1, e3); }
2701  if (e2 >= 0 && e4 >= 0) { midf2 = FindMidEdgeNode(e2, e4); }
2702 
2703  // get proper node if shadow node exists
2704  if (midf1 >= 0 && midf1 == midf2)
2705  {
2706  const Node &nd = nodes[midf1];
2707  if (nd.p1 != e1 && nd.p2 != e1) { midf1 = -1; }
2708  if (nd.p1 != e2 && nd.p2 != e2) { midf2 = -1; }
2709  }
2710 
2711  // only one way to access the mid-face node must always exist
2712  MFEM_ASSERT(!(midf1 >= 0 && midf2 >= 0), "incorrectly split face!");
2713 
2714  if (midf1 < 0 && midf2 < 0) // face not split
2715  {
2716  if (mid) { mid[4] = -1; }
2717  return 0;
2718  }
2719  else if (midf1 >= 0) // face split "vertically"
2720  {
2721  if (mid) { mid[4] = midf1; }
2722  return 1;
2723  }
2724  else // face split "horizontally"
2725  {
2726  if (mid) { mid[4] = midf2; }
2727  return 2;
2728  }
2729 }
2730 
2731 bool NCMesh::TriFaceSplit(int v1, int v2, int v3, int mid[3]) const
2732 {
2733  int e1 = nodes.FindId(v1, v2);
2734  if (e1 < 0 || !nodes[e1].HasVertex()) { return false; }
2735 
2736  int e2 = nodes.FindId(v2, v3);
2737  if (e2 < 0 || !nodes[e2].HasVertex()) { return false; }
2738 
2739  int e3 = nodes.FindId(v3, v1);
2740  if (e3 < 0 || !nodes[e3].HasVertex()) { return false; }
2741 
2742  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3; }
2743 
2744  // NOTE: face (v1, v2, v3) still needs to be checked
2745  return true;
2746 }
2747 
2748 int NCMesh::find_node(const Element &el, int node)
2749 {
2750  for (int i = 0; i < MaxElemNodes; i++)
2751  {
2752  if (el.node[i] == node) { return i; }
2753  }
2754  MFEM_ABORT("Node not found.");
2755  return -1;
2756 }
2757 
2758 int NCMesh::FindNodeExt(const Element &el, int node, bool abort)
2759 {
2760  for (int i = 0; i < GI[el.Geom()].nv; i++)
2761  {
2762  if (RetrieveNode(el, i) == node) { return i; }
2763  }
2764  if (abort) { MFEM_ABORT("Node not found."); }
2765  return -1;
2766 }
2767 
2768 int NCMesh::find_element_edge(const Element &el, int vn0, int vn1, bool abort)
2769 {
2770  MFEM_ASSERT(!el.ref_type, "");
2771 
2772  GeomInfo &gi = GI[el.Geom()];
2773  for (int i = 0; i < gi.ne; i++)
2774  {
2775  const int* ev = gi.edges[i];
2776  int n0 = el.node[ev[0]];
2777  int n1 = el.node[ev[1]];
2778  if ((n0 == vn0 && n1 == vn1) ||
2779  (n0 == vn1 && n1 == vn0)) { return i; }
2780  }
2781 
2782  if (abort) { MFEM_ABORT("Edge (" << vn0 << ", " << vn1 << ") not found"); }
2783  return -1;
2784 }
2785 
2786 int NCMesh::find_local_face(int geom, int a, int b, int c)
2787 {
2788  GeomInfo &gi = GI[geom];
2789  for (int i = 0; i < gi.nf; i++)
2790  {
2791  const int* fv = gi.faces[i];
2792  if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
2793  (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
2794  (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
2795  {
2796  return i;
2797  }
2798  }
2799  MFEM_ABORT("Face not found.");
2800  return -1;
2801 }
2802 
2803 
2804 /// Hash function for a PointMatrix, used in MatrixMap::map.
2805 struct PointMatrixHash
2806 {
2807  std::size_t operator()(const NCMesh::PointMatrix &pm) const
2808  {
2809  MFEM_ASSERT(sizeof(double) == sizeof(std::uint64_t), "");
2810 
2811  // This is a variation on "Hashing an array of floats" from here:
2812  // https://cs.stackexchange.com/questions/37952
2813  std::uint64_t hash = 0xf9ca9ba106acbba9; // random initial value
2814  for (int i = 0; i < pm.np; i++)
2815  {
2816  for (int j = 0; j < pm.points[i].dim; j++)
2817  {
2818  // mix the doubles by adding their binary representations
2819  // many times over (note: 31 is 11111 in binary)
2820  double coord = pm.points[i].coord[j];
2821  hash = 31*hash + *((std::uint64_t*) &coord);
2822  }
2823  }
2824  return hash; // return the lowest bits of the huge sum
2825  }
2826 };
2827 
2828 /** Helper container to keep track of point matrices encountered during
2829  * face/edge traversal and to assign unique indices to them.
2830  */
2831 struct MatrixMap
2832 {
2833  int GetIndex(const NCMesh::PointMatrix &pm)
2834  {
2835  int &index = map[pm];
2836  if (!index) { index = map.size(); }
2837  return index - 1;
2838  }
2839 
2840  void ExportMatrices(Array<DenseMatrix*> &point_matrices) const
2841  {
2842  point_matrices.SetSize(map.size());
2843  for (const auto &pair : map)
2844  {
2845  DenseMatrix* mat = new DenseMatrix();
2846  pair.first.GetMatrix(*mat);
2847  point_matrices[pair.second - 1] = mat;
2848  }
2849  }
2850 
2851  void DumpBucketSizes() const
2852  {
2853  for (unsigned i = 0; i < map.bucket_count(); i++)
2854  {
2855  mfem::out << map.bucket_size(i) << " ";
2856  }
2857  }
2858 
2859 private:
2860  std::unordered_map<NCMesh::PointMatrix, int, PointMatrixHash> map;
2861 };
2862 
2863 
2864 int NCMesh::ReorderFacePointMat(int v0, int v1, int v2, int v3,
2865  int elem, const PointMatrix &pm,
2866  PointMatrix &reordered) const
2867 {
2868  const Element &el = elements[elem];
2869  int master[4] =
2870  {
2871  find_node(el, v0), find_node(el, v1), find_node(el, v2),
2872  (v3 >= 0) ? find_node(el, v3) : -1
2873  };
2874  int nfv = (v3 >= 0) ? 4 : 3;
2875 
2876  int local = find_local_face(el.Geom(), master[0], master[1], master[2]);
2877  const int* fv = GI[el.Geom()].faces[local];
2878 
2879  reordered.np = pm.np;
2880  for (int i = 0, j; i < nfv; i++)
2881  {
2882  for (j = 0; j < nfv; j++)
2883  {
2884  if (fv[i] == master[j])
2885  {
2886  reordered.points[i] = pm.points[j];
2887  break;
2888  }
2889  }
2890  MFEM_ASSERT(j != nfv, "node not found.");
2891  }
2892  return local;
2893 }
2894 
2895 void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
2896  const PointMatrix& pm, int level,
2897  Face* eface[4], MatrixMap &matrix_map)
2898 {
2899  if (level > 0)
2900  {
2901  // check if we made it to a face that is not split further
2902  Face* fa = faces.Find(vn0, vn1, vn2, vn3);
2903  if (fa)
2904  {
2905  // we have a slave face, add it to the list
2906  int elem = fa->GetSingleElement();
2907  face_list.slaves.Append(
2908  Slave(fa->index, elem, -1, Geometry::SQUARE));
2909  Slave &sl = face_list.slaves.Last();
2910 
2911  // reorder the point matrix according to slave face orientation
2912  PointMatrix pm_r;
2913  sl.local = ReorderFacePointMat(vn0, vn1, vn2, vn3, elem, pm, pm_r);
2914  sl.matrix = matrix_map.GetIndex(pm_r);
2915 
2916  eface[0] = eface[2] = fa;
2917  eface[1] = eface[3] = fa;
2918 
2919  return;
2920  }
2921  }
2922 
2923  // we need to recurse deeper
2924  int mid[5];
2925  int split = QuadFaceSplitType(vn0, vn1, vn2, vn3, mid);
2926 
2927  Face *ef[2][4];
2928  if (split == 1) // "X" split face
2929  {
2930  Point pmid0(pm(0), pm(1)), pmid2(pm(2), pm(3));
2931 
2932  TraverseQuadFace(vn0, mid[0], mid[2], vn3,
2933  PointMatrix(pm(0), pmid0, pmid2, pm(3)),
2934  level+1, ef[0], matrix_map);
2935 
2936  TraverseQuadFace(mid[0], vn1, vn2, mid[2],
2937  PointMatrix(pmid0, pm(1), pm(2), pmid2),
2938  level+1, ef[1], matrix_map);
2939 
2940  eface[1] = ef[1][1];
2941  eface[3] = ef[0][3];
2942  eface[0] = eface[2] = NULL;
2943  }
2944  else if (split == 2) // "Y" split face
2945  {
2946  Point pmid1(pm(1), pm(2)), pmid3(pm(3), pm(0));
2947 
2948  TraverseQuadFace(vn0, vn1, mid[1], mid[3],
2949  PointMatrix(pm(0), pm(1), pmid1, pmid3),
2950  level+1, ef[0], matrix_map);
2951 
2952  TraverseQuadFace(mid[3], mid[1], vn2, vn3,
2953  PointMatrix(pmid3, pmid1, pm(2), pm(3)),
2954  level+1, ef[1], matrix_map);
2955 
2956  eface[0] = ef[0][0];
2957  eface[2] = ef[1][2];
2958  eface[1] = eface[3] = NULL;
2959  }
2960 
2961  // check for a prism edge constrained by the master face
2962  if (HavePrisms() && mid[4] >= 0)
2963  {
2964  Node& enode = nodes[mid[4]];
2965  if (enode.HasEdge())
2966  {
2967  // process the edge only if it's not shared by slave faces
2968  // within this master face (i.e. the edge is "hidden")
2969  const int fi[3][2] = {{0, 0}, {1, 3}, {2, 0}};
2970  if (!ef[0][fi[split][0]] && !ef[1][fi[split][1]])
2971  {
2972  MFEM_ASSERT(enode.edge_refc == 1, "");
2973 
2974  MeshId buf[4];
2975  Array<MeshId> eid(buf, 4);
2976 
2977  (split == 1) ? FindEdgeElements(mid[0], vn1, vn2, mid[2], eid)
2978  /* */ : FindEdgeElements(mid[3], vn0, vn1, mid[1], eid);
2979 
2980  MFEM_ASSERT(eid.Size() > 0, "edge prism not found");
2981  MFEM_ASSERT(eid.Size() < 2, "non-unique edge prism");
2982 
2983  // create a slave face record with a degenerate point matrix
2984  face_list.slaves.Append(
2985  Slave(-1 - enode.edge_index,
2986  eid[0].element, eid[0].local, Geometry::SQUARE));
2987  Slave &sl = face_list.slaves.Last();
2988 
2989  if (split == 1)
2990  {
2991  Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
2992  int v1 = nodes[mid[0]].vert_index;
2993  int v2 = nodes[mid[2]].vert_index;
2994  sl.matrix =
2995  matrix_map.GetIndex(
2996  (v1 < v2) ? PointMatrix(mid0, mid2, mid2, mid0) :
2997  /* */ PointMatrix(mid2, mid0, mid0, mid2));
2998  }
2999  else
3000  {
3001  Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
3002  int v1 = nodes[mid[1]].vert_index;
3003  int v2 = nodes[mid[3]].vert_index;
3004  sl.matrix =
3005  matrix_map.GetIndex(
3006  (v1 < v2) ? PointMatrix(mid1, mid3, mid3, mid1) :
3007  /* */ PointMatrix(mid3, mid1, mid1, mid3));
3008  }
3009  }
3010  }
3011  }
3012 }
3013 
3014 void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
3015  MatrixMap &matrix_map)
3016 {
3017  int mid = nodes.FindId(vn0, vn1);
3018  if (mid < 0) { return; }
3019 
3020  const Node &nd = nodes[mid];
3021  if (nd.HasEdge())
3022  {
3023  // check if the edge is already a master in 'edge_list'
3024  int type;
3025  const MeshId &eid = edge_list.LookUp(nd.edge_index, &type);
3026  if (type == 1)
3027  {
3028  // in this case we need to add an edge-face constraint, because the
3029  // master edge is really a (face-)slave itself
3030 
3031  face_list.slaves.Append(
3032  Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
3033 
3034  int v0index = nodes[vn0].vert_index;
3035  int v1index = nodes[vn1].vert_index;
3036 
3037  face_list.slaves.Last().matrix =
3038  matrix_map.GetIndex((v0index < v1index) ? PointMatrix(p0, p1, p0)
3039  /* */ : PointMatrix(p1, p0, p1));
3040 
3041  return; // no need to continue deeper
3042  }
3043  }
3044 
3045  // recurse deeper
3046  Point pmid(p0, p1);
3047  TraverseTetEdge(vn0, mid, p0, pmid, matrix_map);
3048  TraverseTetEdge(mid, vn1, pmid, p1, matrix_map);
3049 }
3050 
3051 bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
3052  const PointMatrix& pm, int level,
3053  MatrixMap &matrix_map)
3054 {
3055  if (level > 0)
3056  {
3057  // check if we made it to a face that is not split further
3058  Face* fa = faces.Find(vn0, vn1, vn2);
3059  if (fa)
3060  {
3061  // we have a slave face, add it to the list
3062  int elem = fa->GetSingleElement();
3063  face_list.slaves.Append(
3064  Slave(fa->index, elem, -1, Geometry::TRIANGLE));
3065  Slave &sl = face_list.slaves.Last();
3066 
3067  // reorder the point matrix according to slave face orientation
3068  PointMatrix pm_r;
3069  sl.local = ReorderFacePointMat(vn0, vn1, vn2, -1, elem, pm, pm_r);
3070  sl.matrix = matrix_map.GetIndex(pm_r);
3071 
3072  return true;
3073  }
3074  }
3075 
3076  int mid[3];
3077  if (TriFaceSplit(vn0, vn1, vn2, mid))
3078  {
3079  Point pmid0(pm(0), pm(1)), pmid1(pm(1), pm(2)), pmid2(pm(2), pm(0));
3080  bool b[4];
3081 
3082  b[0] = TraverseTriFace(vn0, mid[0], mid[2],
3083  PointMatrix(pm(0), pmid0, pmid2),
3084  level+1, matrix_map);
3085 
3086  b[1] = TraverseTriFace(mid[0], vn1, mid[1],
3087  PointMatrix(pmid0, pm(1), pmid1),
3088  level+1, matrix_map);
3089 
3090  b[2] = TraverseTriFace(mid[2], mid[1], vn2,
3091  PointMatrix(pmid2, pmid1, pm(2)),
3092  level+1, matrix_map);
3093 
3094  b[3] = TraverseTriFace(mid[1], mid[2], mid[0],
3095  PointMatrix(pmid1, pmid2, pmid0),
3096  level+1, matrix_map);
3097 
3098  // traverse possible tet edges constrained by the master face
3099  if (HaveTets() && !b[3])
3100  {
3101  if (!b[1]) { TraverseTetEdge(mid[0],mid[1], pmid0,pmid1, matrix_map); }
3102  if (!b[2]) { TraverseTetEdge(mid[1],mid[2], pmid1,pmid2, matrix_map); }
3103  if (!b[0]) { TraverseTetEdge(mid[2],mid[0], pmid2,pmid0, matrix_map); }
3104  }
3105  }
3106 
3107  return false;
3108 }
3109 
3111 {
3112  face_list.Clear();
3113  if (Dim < 3) { return; }
3114 
3115  if (HaveTets()) { GetEdgeList(); } // needed by TraverseTetEdge()
3116 
3118 
3119  Array<char> processed_faces(faces.NumIds());
3120  processed_faces = 0;
3121 
3122  MatrixMap matrix_maps[Geometry::NumGeom];
3123 
3124  // visit faces of leaf elements
3125  for (int i = 0; i < leaf_elements.Size(); i++)
3126  {
3127  int elem = leaf_elements[i];
3128  Element &el = elements[elem];
3129  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3130 
3131  GeomInfo& gi = GI[el.Geom()];
3132  for (int j = 0; j < gi.nf; j++)
3133  {
3134  // get nodes for this face
3135  int node[4];
3136  for (int k = 0; k < 4; k++)
3137  {
3138  node[k] = el.node[gi.faces[j][k]];
3139  }
3140 
3141  int face = faces.FindId(node[0], node[1], node[2], node[3]);
3142  MFEM_ASSERT(face >= 0, "face not found!");
3143 
3144  // tell ParNCMesh about the face
3145  ElementSharesFace(elem, j, face);
3146 
3147  // have we already processed this face? skip if yes
3148  if (processed_faces[face]) { continue; }
3149  processed_faces[face] = 1;
3150 
3151  int fgeom = (node[3] >= 0) ? Geometry::SQUARE : Geometry::TRIANGLE;
3152 
3153  Face &fa = faces[face];
3154  if (fa.elem[0] >= 0 && fa.elem[1] >= 0)
3155  {
3156  // this is a conforming face, add it to the list
3157  face_list.conforming.Append(MeshId(fa.index, elem, j, fgeom));
3158  }
3159  else
3160  {
3161  // this is either a master face or a slave face, but we can't
3162  // tell until we traverse the face refinement 'tree'...
3163  int sb = face_list.slaves.Size();
3164  if (fgeom == Geometry::SQUARE)
3165  {
3166  Face* dummy[4];
3167  TraverseQuadFace(node[0], node[1], node[2], node[3],
3168  pm_quad_identity, 0, dummy, matrix_maps[fgeom]);
3169  }
3170  else
3171  {
3172  TraverseTriFace(node[0], node[1], node[2],
3173  pm_tri_identity, 0, matrix_maps[fgeom]);
3174  }
3175 
3176  int se = face_list.slaves.Size();
3177  if (sb < se)
3178  {
3179  // found slaves, so this is a master face; add it to the list
3180  face_list.masters.Append(
3181  Master(fa.index, elem, j, fgeom, sb, se));
3182 
3183  // also, set the master index for the slaves
3184  for (int ii = sb; ii < se; ii++)
3185  {
3186  face_list.slaves[ii].master = fa.index;
3187  }
3188  }
3189  }
3190 
3191  if (fa.Boundary()) { boundary_faces.Append(face); }
3192  }
3193  }
3194 
3195  // export unique point matrices
3196  for (int i = 0; i < Geometry::NumGeom; i++)
3197  {
3198  matrix_maps[i].ExportMatrices(face_list.point_matrices[i]);
3199  }
3200 }
3201 
3202 void NCMesh::TraverseEdge(int vn0, int vn1, double t0, double t1, int flags,
3203  int level, MatrixMap &matrix_map)
3204 {
3205  int mid = nodes.FindId(vn0, vn1);
3206  if (mid < 0) { return; }
3207 
3208  Node &nd = nodes[mid];
3209  if (nd.HasEdge() && level > 0)
3210  {
3211  // we have a slave edge, add it to the list
3212  edge_list.slaves.Append(Slave(nd.edge_index, -1, -1, Geometry::SEGMENT));
3213 
3214  Slave &sl = edge_list.slaves.Last();
3215  sl.matrix = matrix_map.GetIndex(PointMatrix(Point(t0), Point(t1)));
3216 
3217  // handle slave edge orientation
3218  sl.edge_flags = flags;
3219  int v0index = nodes[vn0].vert_index;
3220  int v1index = nodes[vn1].vert_index;
3221  if (v0index > v1index) { sl.edge_flags |= 2; }
3222  }
3223 
3224  // recurse deeper
3225  double tmid = (t0 + t1) / 2;
3226  TraverseEdge(vn0, mid, t0, tmid, flags, level+1, matrix_map);
3227  TraverseEdge(mid, vn1, tmid, t1, flags, level+1, matrix_map);
3228 }
3229 
3231 {
3232  edge_list.Clear();
3233  if (Dim < 3) { boundary_faces.SetSize(0); }
3234 
3235  Array<char> processed_edges(nodes.NumIds());
3236  processed_edges = 0;
3237 
3238  Array<int> edge_element(nodes.NumIds());
3239  Array<signed char> edge_local(nodes.NumIds());
3240  edge_local = -1;
3241 
3242  MatrixMap matrix_map;
3243 
3244  // visit edges of leaf elements
3245  for (int i = 0; i < leaf_elements.Size(); i++)
3246  {
3247  int elem = leaf_elements[i];
3248  Element &el = elements[elem];
3249  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3250 
3251  GeomInfo& gi = GI[el.Geom()];
3252  for (int j = 0; j < gi.ne; j++)
3253  {
3254  // get nodes for this edge
3255  const int* ev = gi.edges[j];
3256  int node[2] = { el.node[ev[0]], el.node[ev[1]] };
3257 
3258  int enode = nodes.FindId(node[0], node[1]);
3259  MFEM_ASSERT(enode >= 0, "edge node not found!");
3260 
3261  Node &nd = nodes[enode];
3262  MFEM_ASSERT(nd.HasEdge(), "edge not found!");
3263 
3264  // tell ParNCMesh about the edge
3265  ElementSharesEdge(elem, j, enode);
3266 
3267  // (2D only, store boundary faces)
3268  if (Dim <= 2)
3269  {
3270  int face = faces.FindId(node[0], node[0], node[1], node[1]);
3271  MFEM_ASSERT(face >= 0, "face not found!");
3272  if (faces[face].Boundary()) { boundary_faces.Append(face); }
3273  }
3274 
3275  // store element/local for later
3276  edge_element[nd.edge_index] = elem;
3277  edge_local[nd.edge_index] = j;
3278 
3279  // skip slave edges here, they will be reached from their masters
3280  if (GetEdgeMaster(enode) >= 0) { continue; }
3281 
3282  // have we already processed this edge? skip if yes
3283  if (processed_edges[enode]) { continue; }
3284  processed_edges[enode] = 1;
3285 
3286  // prepare edge interval for slave traversal, handle orientation
3287  double t0 = 0.0, t1 = 1.0;
3288  int v0index = nodes[node[0]].vert_index;
3289  int v1index = nodes[node[1]].vert_index;
3290  int flags = (v0index > v1index) ? 1 : 0;
3291 
3292  // try traversing the edge to find slave edges
3293  int sb = edge_list.slaves.Size();
3294  TraverseEdge(node[0], node[1], t0, t1, flags, 0, matrix_map);
3295 
3296  int se = edge_list.slaves.Size();
3297  if (sb < se)
3298  {
3299  // found slaves, this is a master face; add it to the list
3300  edge_list.masters.Append(
3301  Master(nd.edge_index, elem, j, Geometry::SEGMENT, sb, se));
3302 
3303  // also, set the master index for the slaves
3304  for (int ii = sb; ii < se; ii++)
3305  {
3306  edge_list.slaves[ii].master = nd.edge_index;
3307  }
3308  }
3309  else
3310  {
3311  // no slaves, this is a conforming edge
3312  edge_list.conforming.Append(MeshId(nd.edge_index, elem, j));
3313  }
3314  }
3315  }
3316 
3317  // fix up slave edge element/local
3318  for (int i = 0; i < edge_list.slaves.Size(); i++)
3319  {
3320  Slave &sl = edge_list.slaves[i];
3321  int local = edge_local[sl.index];
3322  if (local >= 0)
3323  {
3324  sl.local = local;
3325  sl.element = edge_element[sl.index];
3326  }
3327  }
3328 
3329  // export unique point matrices
3330  matrix_map.ExportMatrices(edge_list.point_matrices[Geometry::SEGMENT]);
3331 }
3332 
3334 {
3335  int total = NVertices + NGhostVertices;
3336 
3337  vertex_list.Clear();
3338  vertex_list.conforming.Reserve(total);
3339 
3340  Array<char> processed_vertices(total);
3341  processed_vertices = 0;
3342 
3343  // analogously to above, visit vertices of leaf elements
3344  for (int i = 0; i < leaf_elements.Size(); i++)
3345  {
3346  int elem = leaf_elements[i];
3347  Element &el = elements[elem];
3348 
3349  for (int j = 0; j < GI[el.Geom()].nv; j++)
3350  {
3351  int node = el.node[j];
3352  Node &nd = nodes[node];
3353 
3354  int index = nd.vert_index;
3355  if (index >= 0)
3356  {
3357  ElementSharesVertex(elem, j, node);
3358 
3359  if (processed_vertices[index]) { continue; }
3360  processed_vertices[index] = 1;
3361 
3362  vertex_list.conforming.Append(MeshId(index, elem, j));
3363  }
3364  }
3365  }
3366 }
3367 
3369  DenseMatrix &oriented_matrix) const
3370 {
3371  oriented_matrix = *(point_matrices[slave.Geom()][slave.matrix]);
3372 
3373  if (slave.edge_flags)
3374  {
3375  MFEM_ASSERT(oriented_matrix.Height() == 1 &&
3376  oriented_matrix.Width() == 2, "not an edge point matrix");
3377 
3378  if (slave.edge_flags & 1) // master inverted
3379  {
3380  oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
3381  oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
3382  }
3383  if (slave.edge_flags & 2) // slave inverted
3384  {
3385  std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
3386  }
3387  }
3388 }
3389 
3391 {
3392  conforming.DeleteAll();
3393  masters.DeleteAll();
3394  slaves.DeleteAll();
3395 
3396  for (int i = 0; i < Geometry::NumGeom; i++)
3397  {
3398  for (int j = 0; j < point_matrices[i].Size(); j++)
3399  {
3400  delete point_matrices[i][j];
3401  }
3402  point_matrices[i].DeleteAll();
3403  }
3404 
3405  inv_index.DeleteAll();
3406 }
3407 
3409 {
3410  return conforming.Size() + masters.Size() + slaves.Size();
3411 }
3412 
3413 const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
3414 {
3415  if (!inv_index.Size())
3416  {
3417  int max_index = -1;
3418  for (int i = 0; i < conforming.Size(); i++)
3419  {
3420  max_index = std::max(conforming[i].index, max_index);
3421  }
3422  for (int i = 0; i < masters.Size(); i++)
3423  {
3424  max_index = std::max(masters[i].index, max_index);
3425  }
3426  for (int i = 0; i < slaves.Size(); i++)
3427  {
3428  if (slaves[i].index < 0) { continue; }
3429  max_index = std::max(slaves[i].index, max_index);
3430  }
3431 
3432  inv_index.SetSize(max_index + 1);
3433  inv_index = -1;
3434 
3435  for (int i = 0; i < conforming.Size(); i++)
3436  {
3437  inv_index[conforming[i].index] = (i << 2);
3438  }
3439  for (int i = 0; i < masters.Size(); i++)
3440  {
3441  inv_index[masters[i].index] = (i << 2) + 1;
3442  }
3443  for (int i = 0; i < slaves.Size(); i++)
3444  {
3445  if (slaves[i].index < 0) { continue; }
3446  inv_index[slaves[i].index] = (i << 2) + 2;
3447  }
3448  }
3449 
3450  MFEM_ASSERT(index >= 0 && index < inv_index.Size(), "");
3451  int key = inv_index[index];
3452 
3453  if (!type)
3454  {
3455  MFEM_VERIFY(key >= 0, "index " << index << " not found.");
3456  }
3457  else // return entity type if requested, don't abort when not found
3458  {
3459  *type = (key >= 0) ? (key & 0x3) : -1;
3460 
3461  static MeshId invalid;
3462  if (*type < 0) { return invalid; } // not found
3463  }
3464 
3465  // return found entity MeshId
3466  switch (key & 0x3)
3467  {
3468  case 0: return conforming[key >> 2];
3469  case 1: return masters[key >> 2];
3470  case 2: return slaves[key >> 2];
3471  default: MFEM_ABORT("internal error"); return conforming[0];
3472  }
3473 }
3474 
3475 
3476 //// Neighbors /////////////////////////////////////////////////////////////////
3477 
3478 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
3479 {
3480  int mid = nodes.FindId(v0, v1);
3481  if (mid >= 0 && nodes[mid].HasVertex())
3482  {
3483  indices.Append(mid);
3484 
3485  CollectEdgeVertices(v0, mid, indices);
3486  CollectEdgeVertices(mid, v1, indices);
3487  }
3488 }
3489 
3490 void NCMesh::CollectTriFaceVertices(int v0, int v1, int v2, Array<int> &indices)
3491 {
3492  int mid[3];
3493  if (TriFaceSplit(v0, v1, v2, mid))
3494  {
3495  for (int i = 0; i < 3; i++)
3496  {
3497  indices.Append(mid[i]);
3498  }
3499 
3500  CollectTriFaceVertices(v0, mid[0], mid[2], indices);
3501  CollectTriFaceVertices(mid[0], v1, mid[1], indices);
3502  CollectTriFaceVertices(mid[2], mid[1], v2, indices);
3503  CollectTriFaceVertices(mid[0], mid[1], mid[2], indices);
3504 
3505  if (HaveTets()) // possible edge-face contact
3506  {
3507  CollectEdgeVertices(mid[0], mid[1], indices);
3508  CollectEdgeVertices(mid[1], mid[2], indices);
3509  CollectEdgeVertices(mid[2], mid[0], indices);
3510  }
3511  }
3512 }
3513 
3514 void NCMesh::CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
3515  Array<int> &indices)
3516 {
3517  int mid[5];
3518  switch (QuadFaceSplitType(v0, v1, v2, v3, mid))
3519  {
3520  case 1:
3521  indices.Append(mid[0]);
3522  indices.Append(mid[2]);
3523 
3524  CollectQuadFaceVertices(v0, mid[0], mid[2], v3, indices);
3525  CollectQuadFaceVertices(mid[0], v1, v2, mid[2], indices);
3526 
3527  if (HavePrisms()) // possible edge-face contact
3528  {
3529  CollectEdgeVertices(mid[0], mid[2], indices);
3530  }
3531  break;
3532 
3533  case 2:
3534  indices.Append(mid[1]);
3535  indices.Append(mid[3]);
3536 
3537  CollectQuadFaceVertices(v0, v1, mid[1], mid[3], indices);
3538  CollectQuadFaceVertices(mid[3], mid[1], v2, v3, indices);
3539 
3540  if (HavePrisms()) // possible edge-face contact
3541  {
3542  CollectEdgeVertices(mid[1], mid[3], indices);
3543  }
3544  break;
3545  }
3546 }
3547 
3549 {
3550  int nrows = leaf_elements.Size();
3551  int* I = Memory<int>(nrows + 1);
3552  int** JJ = new int*[nrows];
3553 
3554  Array<int> indices;
3555  indices.Reserve(128);
3556 
3557  // collect vertices coinciding with each element, including hanging vertices
3558  for (int i = 0; i < leaf_elements.Size(); i++)
3559  {
3560  int elem = leaf_elements[i];
3561  Element &el = elements[elem];
3562  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3563 
3564  GeomInfo& gi = GI[el.Geom()];
3565  int* node = el.node;
3566 
3567  indices.SetSize(0);
3568  for (int j = 0; j < gi.ne; j++)
3569  {
3570  const int* ev = gi.edges[j];
3571  CollectEdgeVertices(node[ev[0]], node[ev[1]], indices);
3572  }
3573 
3574  if (Dim >= 3)
3575  {
3576  for (int j = 0; j < gi.nf; j++)
3577  {
3578  const int* fv = gi.faces[j];
3579  if (gi.nfv[j] == 4)
3580  {
3581  CollectQuadFaceVertices(node[fv[0]], node[fv[1]],
3582  node[fv[2]], node[fv[3]], indices);
3583  }
3584  else
3585  {
3586  CollectTriFaceVertices(node[fv[0]], node[fv[1]], node[fv[2]],
3587  indices);
3588  }
3589  }
3590  }
3591 
3592  // temporarily store one row of the table
3593  indices.Sort();
3594  indices.Unique();
3595  int size = indices.Size();
3596  I[i] = size;
3597  JJ[i] = new int[size];
3598  std::memcpy(JJ[i], indices.GetData(), size * sizeof(int));
3599  }
3600 
3601  // finalize the I array of the table
3602  int nnz = 0;
3603  for (int i = 0; i < nrows; i++)
3604  {
3605  int cnt = I[i];
3606  I[i] = nnz;
3607  nnz += cnt;
3608  }
3609  I[nrows] = nnz;
3610 
3611  // copy the temporarily stored rows into one J array
3612  int *J = Memory<int>(nnz);
3613  nnz = 0;
3614  for (int i = 0; i < nrows; i++)
3615  {
3616  int cnt = I[i+1] - I[i];
3617  std::memcpy(J+nnz, JJ[i], cnt * sizeof(int));
3618  delete [] JJ[i];
3619  nnz += cnt;
3620  }
3621 
3622  element_vertex.SetIJ(I, J, nrows);
3623 
3624  delete [] JJ;
3625 }
3626 
3627 
3629  Array<int> *neighbors,
3630  Array<char> *neighbor_set)
3631 {
3632  // If A is the element-to-vertex table (see 'element_vertex') listing all
3633  // vertices touching each element, including hanging vertices, then A*A^T is
3634  // the element-to-neighbor table. Multiplying the element set with A*A^T
3635  // gives the neighbor set. To save memory, this function only computes the
3636  // action of A*A^T, the product itself is not stored anywhere.
3637 
3638  // Optimization: the 'element_vertex' table does not store the obvious
3639  // corner nodes in it. The table is therefore empty for conforming meshes.
3640 
3642 
3643  int nleaves = leaf_elements.Size();
3644  MFEM_VERIFY(elem_set.Size() == nleaves, "");
3645  MFEM_ASSERT(element_vertex.Size() == nleaves, "");
3646 
3647  // step 1: vertices = A^T * elem_set, i.e, find all vertices touching the
3648  // element set
3649 
3650  Array<char> vmark(nodes.NumIds());
3651  vmark = 0;
3652 
3653  for (int i = 0; i < nleaves; i++)
3654  {
3655  if (elem_set[i])
3656  {
3657  int *v = element_vertex.GetRow(i);
3658  int nv = element_vertex.RowSize(i);
3659  for (int j = 0; j < nv; j++)
3660  {
3661  vmark[v[j]] = 1;
3662  }
3663 
3664  Element &el = elements[leaf_elements[i]];
3665  nv = GI[el.Geom()].nv;
3666  for (int j = 0; j < nv; j++)
3667  {
3668  vmark[el.node[j]] = 1;
3669  }
3670  }
3671  }
3672 
3673  // step 2: neighbors = A * vertices, i.e., find all elements coinciding with
3674  // vertices from step 1; NOTE: in the result we don't include elements from
3675  // the original set
3676 
3677  if (neighbor_set)
3678  {
3679  neighbor_set->SetSize(nleaves);
3680  *neighbor_set = 0;
3681  }
3682 
3683  for (int i = 0; i < nleaves; i++)
3684  {
3685  if (!elem_set[i])
3686  {
3687  bool hit = false;
3688 
3689  int *v = element_vertex.GetRow(i);
3690  int nv = element_vertex.RowSize(i);
3691  for (int j = 0; j < nv; j++)
3692  {
3693  if (vmark[v[j]]) { hit = true; break; }
3694  }
3695 
3696  if (!hit)
3697  {
3698  Element &el = elements[leaf_elements[i]];
3699  nv = GI[el.Geom()].nv;
3700  for (int j = 0; j < nv; j++)
3701  {
3702  if (vmark[el.node[j]]) { hit = true; break; }
3703  }
3704  }
3705 
3706  if (hit)
3707  {
3708  if (neighbors) { neighbors->Append(leaf_elements[i]); }
3709  if (neighbor_set) { (*neighbor_set)[i] = 1; }
3710  }
3711  }
3712  }
3713 }
3714 
3715 static bool sorted_lists_intersect(const int* a, const int* b, int na, int nb)
3716 {
3717  // pointers to "end" sentinel, not last entry. Not for dereferencing.
3718  const int * const a_end = a + na;
3719  const int * const b_end = b + nb;
3720  while (a != a_end && b != b_end)
3721  {
3722  if (*a < *b)
3723  {
3724  ++a;
3725  }
3726  else if (*b < *a)
3727  {
3728  ++b;
3729  }
3730  else
3731  {
3732  return true; // neither *a < *b nor *b < *a thus a == b
3733  }
3734  }
3735  return false; // no common element found
3736 }
3737 
3738 
3739 void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
3740  const Array<int> *search_set)
3741 {
3742  // TODO future: this function is inefficient. For a single element, an
3743  // octree neighbor search algorithm would be better. However, the octree
3744  // neighbor algorithm is hard to get right in the multi-octree case due to
3745  // the different orientations of the octrees (i.e., the root elements).
3746 
3748 
3749  // sorted list of all vertex nodes touching 'elem'
3750  Array<int> vert;
3751  vert.Reserve(128);
3752 
3753  // support for non-leaf 'elem', collect vertices of all children
3754  Array<int> stack;
3755  stack.Reserve(64);
3756  stack.Append(elem);
3757 
3758  while (stack.Size())
3759  {
3760  Element &el = elements[stack.Last()];
3761  stack.DeleteLast();
3762 
3763  if (!el.ref_type)
3764  {
3765  int *v = element_vertex.GetRow(el.index);
3766  int nv = element_vertex.RowSize(el.index);
3767  for (int i = 0; i < nv; i++)
3768  {
3769  vert.Append(v[i]);
3770  }
3771 
3772  nv = GI[el.Geom()].nv;
3773  for (int i = 0; i < nv; i++)
3774  {
3775  vert.Append(el.node[i]);
3776  }
3777  }
3778  else
3779  {
3780  for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
3781  {
3782  stack.Append(el.child[i]);
3783  }
3784  }
3785  }
3786 
3787  vert.Sort();
3788  vert.Unique();
3789 
3790  int *v1 = vert.GetData();
3791  int nv1 = vert.Size();
3792 
3793  if (!search_set) { search_set = &leaf_elements; }
3794 
3795  // test *all* potential neighbors from the search set
3796  for (int i = 0; i < search_set->Size(); i++)
3797  {
3798  int testme = (*search_set)[i];
3799  if (testme != elem)
3800  {
3801  Element &el = elements[testme];
3802  int *v2 = element_vertex.GetRow(el.index);
3803  int nv2 = element_vertex.RowSize(el.index);
3804 
3805  bool hit = sorted_lists_intersect(v1, v2, nv1, nv2);
3806 
3807  if (!hit)
3808  {
3809  int nv = GI[el.Geom()].nv;
3810  for (int j = 0; j < nv; j++)
3811  {
3812  hit = sorted_lists_intersect(&el.node[j], v1, 1, nv1);
3813  if (hit) { break; }
3814  }
3815  }
3816 
3817  if (hit) { neighbors.Append(testme); }
3818  }
3819  }
3820 }
3821 
3823  Array<int> &expanded,
3824  const Array<int> *search_set)
3825 {
3827 
3828  Array<char> vmark(nodes.NumIds());
3829  vmark = 0;
3830 
3831  for (int i = 0; i < elems.Size(); i++)
3832  {
3833  Element &el = elements[elems[i]];
3834 
3835  int *v = element_vertex.GetRow(el.index);
3836  int nv = element_vertex.RowSize(el.index);
3837  for (int j = 0; j < nv; j++)
3838  {
3839  vmark[v[j]] = 1;
3840  }
3841 
3842  nv = GI[el.Geom()].nv;
3843  for (int j = 0; j < nv; j++)
3844  {
3845  vmark[el.node[j]] = 1;
3846  }
3847  }
3848 
3849  if (!search_set)
3850  {
3851  search_set = &leaf_elements;
3852  }
3853 
3854  expanded.SetSize(0);
3855  for (int i = 0; i < search_set->Size(); i++)
3856  {
3857  int testme = (*search_set)[i];
3858  Element &el = elements[testme];
3859  bool hit = false;
3860 
3861  int *v = element_vertex.GetRow(el.index);
3862  int nv = element_vertex.RowSize(el.index);
3863  for (int j = 0; j < nv; j++)
3864  {
3865  if (vmark[v[j]]) { hit = true; break; }
3866  }
3867 
3868  if (!hit)
3869  {
3870  nv = GI[el.Geom()].nv;
3871  for (int j = 0; j < nv; j++)
3872  {
3873  if (vmark[el.node[j]]) { hit = true; break; }
3874  }
3875  }
3876 
3877  if (hit) { expanded.Append(testme); }
3878  }
3879 }
3880 
3881 int NCMesh::GetVertexRootCoord(int elem, RefCoord coord[3]) const
3882 {
3883  while (1)
3884  {
3885  const Element &el = elements[elem];
3886  if (el.parent < 0) { return elem; }
3887 
3888  const Element &pa = elements[el.parent];
3889  MFEM_ASSERT(pa.ref_type, "internal error");
3890 
3891  int ch = 0;
3892  while (ch < MaxElemChildren && pa.child[ch] != elem) { ch++; }
3893  MFEM_ASSERT(ch < MaxElemChildren, "internal error");
3894 
3895  MFEM_ASSERT(geom_parent[el.Geom()], "unsupported geometry");
3896  const RefTrf &tr = geom_parent[el.Geom()][(int) pa.ref_type][ch];
3897  tr.Apply(coord, coord);
3898 
3899  elem = el.parent;
3900  }
3901 }
3902 
3903 static bool RefPointInside(Geometry::Type geom, const RefCoord pt[3])
3904 {
3905  switch (geom)
3906  {
3907  case Geometry::SQUARE:
3908  return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
3909  (pt[1] >= 0) && (pt[1] <= T_ONE);
3910 
3911  case Geometry::CUBE:
3912  return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
3913  (pt[1] >= 0) && (pt[1] <= T_ONE) &&
3914  (pt[2] >= 0) && (pt[2] <= T_ONE);
3915 
3916  case Geometry::TRIANGLE:
3917  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE);
3918 
3919  case Geometry::PRISM:
3920  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE) &&
3921  (pt[2] >= 0) && (pt[2] <= T_ONE);
3922 
3923  case Geometry::PYRAMID:
3924  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[2] >= 0.0) &&
3925  (pt[0] + pt[2] <= T_ONE) && (pt[1] + pt[2] <= T_ONE) &&
3926  (pt[2] <= T_ONE);
3927 
3928  default:
3929  MFEM_ABORT("unsupported geometry");
3930  return false;
3931  }
3932 }
3933 
3934 void NCMesh::CollectIncidentElements(int elem, const RefCoord coord[3],
3935  Array<int> &list) const
3936 {
3937  const Element &el = elements[elem];
3938  if (!el.ref_type)
3939  {
3940  list.Append(elem);
3941  return;
3942  }
3943 
3944  RefCoord tcoord[3];
3945  for (int ch = 0; ch < MaxElemChildren && el.child[ch] >= 0; ch++)
3946  {
3947  const RefTrf &tr = geom_child[el.Geom()][(int) el.ref_type][ch];
3948  tr.Apply(coord, tcoord);
3949 
3950  if (RefPointInside(el.Geom(), tcoord))
3951  {
3952  CollectIncidentElements(el.child[ch], tcoord, list);
3953  }
3954  }
3955 }
3956 
3957 void NCMesh::FindVertexCousins(int elem, int local, Array<int> &cousins) const
3958 {
3959  const Element &el = elements[elem];
3960 
3961  RefCoord coord[3];
3962  MFEM_ASSERT(geom_corners[el.Geom()], "unsupported geometry");
3963  std::memcpy(coord, geom_corners[el.Geom()][local], sizeof(coord));
3964 
3965  int root = GetVertexRootCoord(elem, coord);
3966 
3967  cousins.SetSize(0);
3968  CollectIncidentElements(root, coord, cousins);
3969 }
3970 
3971 
3972 //// Coarse/fine transformations ///////////////////////////////////////////////
3973 
3975 {
3976  MFEM_ASSERT(np == pm.np, "");
3977  for (int i = 0; i < np; i++)
3978  {
3979  MFEM_ASSERT(points[i].dim == pm.points[i].dim, "");
3980  for (int j = 0; j < points[i].dim; j++)
3981  {
3982  if (points[i].coord[j] != pm.points[i].coord[j]) { return false; }
3983  }
3984  }
3985  return true;
3986 }
3987 
3989 {
3990  point_matrix.SetSize(points[0].dim, np);
3991  for (int i = 0; i < np; i++)
3992  {
3993  for (int j = 0; j < points[0].dim; j++)
3994  {
3995  point_matrix(j, i) = points[i].coord[j];
3996  }
3997  }
3998 }
3999 
4001  Point(0), Point(1)
4002 );
4004  Point(0, 0), Point(1, 0), Point(0, 1)
4005 );
4007  Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)
4008 );
4010  Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)
4011 );
4013  Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0),
4014  Point(0, 0, 1), Point(1, 0, 1), Point(0, 1, 1)
4015 );
4017  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0),
4018  Point(0, 1, 0), Point(0, 0, 1)
4019 );
4021  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
4022  Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
4023 );
4024 
4026 {
4027  switch (geom)
4028  {
4029  case Geometry::SEGMENT: return pm_seg_identity;
4030  case Geometry::TRIANGLE: return pm_tri_identity;
4031  case Geometry::SQUARE: return pm_quad_identity;
4033  case Geometry::PRISM: return pm_prism_identity;
4035  case Geometry::CUBE: return pm_hex_identity;
4036  default:
4037  MFEM_ABORT("unsupported geometry " << geom);
4038  return pm_tri_identity;
4039  }
4040 }
4041 
4042 void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
4043  DenseMatrix& matrix)
4044 {
4045  PointMatrix pm = GetGeomIdentity(geom);
4046 
4047  while (*ref_path)
4048  {
4049  int ref_type = *ref_path++;
4050  int child = *ref_path++;
4051 
4052  // TODO: do this with the new child transform tables
4053 
4054  if (geom == Geometry::CUBE)
4055  {
4056  if (ref_type == 1) // split along X axis
4057  {
4058  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4059  Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
4060 
4061  if (child == 0)
4062  {
4063  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
4064  pm(4), mid45, mid67, pm(7));
4065  }
4066  else if (child == 1)
4067  {
4068  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
4069  mid45, pm(5), pm(6), mid67);
4070  }
4071  }
4072  else if (ref_type == 2) // split along Y axis
4073  {
4074  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4075  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4076 
4077  if (child == 0)
4078  {
4079  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4080  pm(4), pm(5), mid56, mid74);
4081  }
4082  else if (child == 1)
4083  {
4084  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4085  mid74, mid56, pm(6), pm(7));
4086  }
4087  }
4088  else if (ref_type == 4) // split along Z axis
4089  {
4090  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4091  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4092 
4093  if (child == 0)
4094  {
4095  pm = PointMatrix(pm(0), pm(1), pm(2), pm(3),
4096  mid04, mid15, mid26, mid37);
4097  }
4098  else if (child == 1)
4099  {
4100  pm = PointMatrix(mid04, mid15, mid26, mid37,
4101  pm(4), pm(5), pm(6), pm(7));
4102  }
4103  }
4104  else if (ref_type == 3) // XY split
4105  {
4106  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4107  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4108  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4109  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4110 
4111  Point midf0(mid23, mid12, mid01, mid30);
4112  Point midf5(mid45, mid56, mid67, mid74);
4113 
4114  if (child == 0)
4115  {
4116  pm = PointMatrix(pm(0), mid01, midf0, mid30,
4117  pm(4), mid45, midf5, mid74);
4118  }
4119  else if (child == 1)
4120  {
4121  pm = PointMatrix(mid01, pm(1), mid12, midf0,
4122  mid45, pm(5), mid56, midf5);
4123  }
4124  else if (child == 2)
4125  {
4126  pm = PointMatrix(midf0, mid12, pm(2), mid23,
4127  midf5, mid56, pm(6), mid67);
4128  }
4129  else if (child == 3)
4130  {
4131  pm = PointMatrix(mid30, midf0, mid23, pm(3),
4132  mid74, midf5, mid67, pm(7));
4133  }
4134  }
4135  else if (ref_type == 5) // XZ split
4136  {
4137  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4138  Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
4139  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4140  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4141 
4142  Point midf1(mid01, mid15, mid45, mid04);
4143  Point midf3(mid23, mid37, mid67, mid26);
4144 
4145  if (child == 0)
4146  {
4147  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
4148  mid04, midf1, midf3, mid37);
4149  }
4150  else if (child == 1)
4151  {
4152  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
4153  midf1, mid15, mid26, midf3);
4154  }
4155  else if (child == 2)
4156  {
4157  pm = PointMatrix(midf1, mid15, mid26, midf3,
4158  mid45, pm(5), pm(6), mid67);
4159  }
4160  else if (child == 3)
4161  {
4162  pm = PointMatrix(mid04, midf1, midf3, mid37,
4163  pm(4), mid45, mid67, pm(7));
4164  }
4165  }
4166  else if (ref_type == 6) // YZ split
4167  {
4168  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4169  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4170  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4171  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4172 
4173  Point midf2(mid12, mid26, mid56, mid15);
4174  Point midf4(mid30, mid04, mid74, mid37);
4175 
4176  if (child == 0)
4177  {
4178  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4179  mid04, mid15, midf2, midf4);
4180  }
4181  else if (child == 1)
4182  {
4183  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4184  midf4, midf2, mid26, mid37);
4185  }
4186  else if (child == 2)
4187  {
4188  pm = PointMatrix(mid04, mid15, midf2, midf4,
4189  pm(4), pm(5), mid56, mid74);
4190  }
4191  else if (child == 3)
4192  {
4193  pm = PointMatrix(midf4, midf2, mid26, mid37,
4194  mid74, mid56, pm(6), pm(7));
4195  }
4196  }
4197  else if (ref_type == 7) // full isotropic refinement
4198  {
4199  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4200  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4201  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4202  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4203  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4204  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4205 
4206  Point midf0(mid23, mid12, mid01, mid30);
4207  Point midf1(mid01, mid15, mid45, mid04);
4208  Point midf2(mid12, mid26, mid56, mid15);
4209  Point midf3(mid23, mid37, mid67, mid26);
4210  Point midf4(mid30, mid04, mid74, mid37);
4211  Point midf5(mid45, mid56, mid67, mid74);
4212 
4213  Point midel(midf1, midf3);
4214 
4215  if (child == 0)
4216  {
4217  pm = PointMatrix(pm(0), mid01, midf0, mid30,
4218  mid04, midf1, midel, midf4);
4219  }
4220  else if (child == 1)
4221  {
4222  pm = PointMatrix(mid01, pm(1), mid12, midf0,
4223  midf1, mid15, midf2, midel);
4224  }
4225  else if (child == 2)
4226  {
4227  pm = PointMatrix(midf0, mid12, pm(2), mid23,
4228  midel, midf2, mid26, midf3);
4229  }
4230  else if (child == 3)
4231  {
4232  pm = PointMatrix(mid30, midf0, mid23, pm(3),
4233  midf4, midel, midf3, mid37);
4234  }
4235  else if (child == 4)
4236  {
4237  pm = PointMatrix(mid04, midf1, midel, midf4,
4238  pm(4), mid45, midf5, mid74);
4239  }
4240  else if (child == 5)
4241  {
4242  pm = PointMatrix(midf1, mid15, midf2, midel,
4243  mid45, pm(5), mid56, midf5);
4244  }
4245  else if (child == 6)
4246  {
4247  pm = PointMatrix(midel, midf2, mid26, midf3,
4248  midf5, mid56, pm(6), mid67);
4249  }
4250  else if (child == 7)
4251  {
4252  pm = PointMatrix(midf4, midel, midf3, mid37,
4253  mid74, midf5, mid67, pm(7));
4254  }
4255  }
4256  }
4257  else if (geom == Geometry::PRISM)
4258  {
4259  if (ref_type < 4) // XY split
4260  {
4261  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4262  Point mid20(pm(2), pm(0)), mid34(pm(3), pm(4));
4263  Point mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4264 
4265  if (child == 0)
4266  {
4267  pm = PointMatrix(pm(0), mid01, mid20, pm(3), mid34, mid53);
4268  }
4269  else if (child == 1)
4270  {
4271  pm = PointMatrix(mid01, pm(1), mid12, mid34, pm(4), mid45);
4272  }
4273  else if (child == 2)
4274  {
4275  pm = PointMatrix(mid20, mid12, pm(2), mid53, mid45, pm(5));
4276  }
4277  else if (child == 3)
4278  {
4279  pm = PointMatrix(mid12, mid20, mid01, mid45, mid53, mid34);
4280  }
4281  }
4282  else if (ref_type == 4) // Z split
4283  {
4284  Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4285 
4286  if (child == 0)
4287  {
4288  pm = PointMatrix(pm(0), pm(1), pm(2), mid03, mid14, mid25);
4289  }
4290  else if (child == 1)
4291  {
4292  pm = PointMatrix(mid03, mid14, mid25, pm(3), pm(4), pm(5));
4293  }
4294  }
4295  else // ref_type > 4, iso split
4296  {
4297  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4298  Point mid34(pm(3), pm(4)), mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4299  Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4300 
4301  Point midf2(mid01, mid14, mid34, mid03);
4302  Point midf3(mid12, mid25, mid45, mid14);
4303  Point midf4(mid20, mid03, mid53, mid25);
4304 
4305  if (child == 0)
4306  {
4307  pm = PointMatrix(pm(0), mid01, mid20, mid03, midf2, midf4);
4308  }
4309  else if (child == 1)
4310  {
4311  pm = PointMatrix(mid01, pm(1), mid12, midf2, mid14, midf3);
4312  }
4313  else if (child == 2)
4314  {
4315  pm = PointMatrix(mid20, mid12, pm(2), midf4, midf3, mid25);
4316  }
4317  else if (child == 3)
4318  {
4319  pm = PointMatrix(mid12, mid20, mid01, midf3, midf4, midf2);
4320  }
4321  else if (child == 4)
4322  {
4323  pm = PointMatrix(mid03, midf2, midf4, pm(3), mid34, mid53);
4324  }
4325  else if (child == 5)
4326  {
4327  pm = PointMatrix(midf2, mid14, midf3, mid34, pm(4), mid45);
4328  }
4329  else if (child == 6)
4330  {
4331  pm = PointMatrix(midf4, midf3, mid25, mid53, mid45, pm(5));
4332  }
4333  else if (child == 7)
4334  {
4335  pm = PointMatrix(midf3, midf4, midf2, mid45, mid53, mid34);
4336  }
4337  }
4338  }
4339  else if (geom == Geometry::PYRAMID)
4340  {
4341  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4342  Point mid03(pm(0), pm(3)), mid12(pm(1), pm(2));
4343  Point mid04(pm(0), pm(4)), mid14(pm(1), pm(4));
4344  Point mid24(pm(2), pm(4)), mid34(pm(3), pm(4));
4345  Point midf0(mid23, mid12, mid01, mid03);
4346 
4347  if (child == 0) // Pyramid
4348  {
4349  pm = PointMatrix(pm(0), mid01, midf0, mid03, mid04);
4350  }
4351  else if (child == 1) // Pyramid
4352  {
4353  pm = PointMatrix(mid01, pm(1), mid12, midf0, mid14);
4354  }
4355  else if (child == 2) // Pyramid
4356  {
4357  pm = PointMatrix(midf0, mid12, pm(2), mid23, mid24);
4358  }
4359  else if (child == 3) // Pyramid
4360  {
4361  pm = PointMatrix(mid03, midf0, mid23, pm(3), mid34);
4362  }
4363  else if (child == 4) // Pyramid
4364  {
4365  pm = PointMatrix(mid24, mid14, mid04, mid34, midf0);
4366  }
4367  else if (child == 5) // Pyramid
4368  {
4369  pm = PointMatrix(mid04, mid14, mid24, mid34, pm(4));
4370  }
4371  else if (child == 6) // Tet
4372  {
4373  pm = PointMatrix(mid01, midf0, mid04, mid14);
4374  }
4375  else if (child == 7) // Tet
4376  {
4377  pm = PointMatrix(midf0, mid14, mid12, mid24);
4378  }
4379  else if (child == 8) // Tet
4380  {
4381  pm = PointMatrix(midf0, mid23, mid34, mid24);
4382  }
4383  else if (child == 9) // Tet
4384  {
4385  pm = PointMatrix(mid03, mid04, midf0, mid34);
4386  }
4387  }
4388  else if (geom == Geometry::TETRAHEDRON)
4389  {
4390  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid02(pm(2), pm(0));
4391  Point mid03(pm(0), pm(3)), mid13(pm(1), pm(3)), mid23(pm(2), pm(3));
4392 
4393  if (child == 0)
4394  {
4395  pm = PointMatrix(pm(0), mid01, mid02, mid03);
4396  }
4397  else if (child == 1)
4398  {
4399  pm = PointMatrix(mid01, pm(1), mid12, mid13);
4400  }
4401  else if (child == 2)
4402  {
4403  pm = PointMatrix(mid02, mid12, pm(2), mid23);
4404  }
4405  else if (child == 3)
4406  {
4407  pm = PointMatrix(mid03, mid13, mid23, pm(3));
4408  }
4409  else if (child == 4)
4410  {
4411  pm = PointMatrix(mid01, mid23, mid02, mid03);
4412  }
4413  else if (child == 5)
4414  {
4415  pm = PointMatrix(mid01, mid23, mid03, mid13);
4416  }
4417  else if (child == 6)
4418  {
4419  pm = PointMatrix(mid01, mid23, mid13, mid12);
4420  }
4421  else if (child == 7)
4422  {
4423  pm = PointMatrix(mid01, mid23, mid12, mid02);
4424  }
4425  }
4426  else if (geom == Geometry::SQUARE)
4427  {
4428  if (ref_type == 1) // X split
4429  {
4430  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4431 
4432  if (child == 0)
4433  {
4434  pm = PointMatrix(pm(0), mid01, mid23, pm(3));
4435  }
4436  else if (child == 1)
4437  {
4438  pm = PointMatrix(mid01, pm(1), pm(2), mid23);
4439  }
4440  }
4441  else if (ref_type == 2) // Y split
4442  {
4443  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4444 
4445  if (child == 0)
4446  {
4447  pm = PointMatrix(pm(0), pm(1), mid12, mid30);
4448  }
4449  else if (child == 1)
4450  {
4451  pm = PointMatrix(mid30, mid12, pm(2), pm(3));
4452  }
4453  }
4454  else if (ref_type == 3) // iso split
4455  {
4456  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4457  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4458  Point midel(mid01, mid23);
4459 
4460  if (child == 0)
4461  {
4462  pm = PointMatrix(pm(0), mid01, midel, mid30);
4463  }
4464  else if (child == 1)
4465  {
4466  pm = PointMatrix(mid01, pm(1), mid12, midel);
4467  }
4468  else if (child == 2)
4469  {
4470  pm = PointMatrix(midel, mid12, pm(2), mid23);
4471  }
4472  else if (child == 3)
4473  {
4474  pm = PointMatrix(mid30, midel, mid23, pm(3));
4475  }
4476  }
4477  }
4478  else if (geom == Geometry::TRIANGLE)
4479  {
4480  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4481 
4482  if (child == 0)
4483  {
4484  pm = PointMatrix(pm(0), mid01, mid20);
4485  }
4486  else if (child == 1)
4487  {
4488  pm = PointMatrix(mid01, pm(1), mid12);
4489  }
4490  else if (child == 2)
4491  {
4492  pm = PointMatrix(mid20, mid12, pm(2));
4493  }
4494  else if (child == 3)
4495  {
4496  pm = PointMatrix(mid12, mid20, mid01);
4497  }
4498  }
4499  else if (geom == Geometry::SEGMENT)
4500  {
4501  Point mid01(pm(0), pm(1));
4502 
4503  if (child == 0)
4504  {
4505  pm = PointMatrix(pm(0), mid01);
4506  }
4507  else if (child == 1)
4508  {
4509  pm = PointMatrix(mid01, pm(1));
4510  }
4511  }
4512  }
4513 
4514  // write the points to the matrix
4515  for (int i = 0; i < pm.np; i++)
4516  {
4517  for (int j = 0; j < pm(i).dim; j++)
4518  {
4519  matrix(j, i) = pm(i).coord[j];
4520  }
4521  }
4522 }
4523 
4525 {
4528 
4529  for (int i = 0; i < leaf_elements.Size(); i++)
4530  {
4531  int elem = leaf_elements[i];
4532  if (!IsGhost(elements[elem])) { coarse_elements.Append(elem); }
4533  }
4534 
4535  transforms.embeddings.DeleteAll();
4536 }
4537 
4538 void NCMesh::TraverseRefinements(int elem, int coarse_index,
4539  std::string &ref_path, RefPathMap &map)
4540 {
4541  Element &el = elements[elem];
4542  if (!el.ref_type)
4543  {
4544  int &matrix = map[ref_path];
4545  if (!matrix) { matrix = map.size(); }
4546 
4547  Embedding &emb = transforms.embeddings[el.index];
4548  emb.parent = coarse_index;
4549  emb.matrix = matrix - 1;
4550  emb.geom = el.Geom();
4551  emb.ghost = IsGhost(el);
4552  }
4553  else
4554  {
4555  MFEM_ASSERT(el.tet_type == 0, "not implemented");
4556 
4557  ref_path.push_back(el.ref_type);
4558  ref_path.push_back(0);
4559 
4560  for (int i = 0; i < MaxElemChildren; i++)
4561  {
4562  if (el.child[i] >= 0)
4563  {
4564  ref_path[ref_path.length()-1] = i;
4565  TraverseRefinements(el.child[i], coarse_index, ref_path, map);
4566  }
4567  }
4568  ref_path.resize(ref_path.length()-2);
4569  }
4570 }
4571 
4573 {
4574  MFEM_VERIFY(coarse_elements.Size() || !leaf_elements.Size(),
4575  "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
4576  " and Refine().");
4577 
4578  if (!transforms.embeddings.Size())
4579  {
4580  transforms.Clear();
4581  transforms.embeddings.SetSize(NElements);
4582 
4583  std::string ref_path;
4584  ref_path.reserve(100);
4585 
4586  RefPathMap path_map[Geometry::NumGeom];
4587  for (int g = 0; g < Geometry::NumGeom; g++)
4588  {
4589  path_map[g][ref_path] = 1; // empty path == identity
4590  }
4591 
4592  int used_geoms = 0;
4593  for (int i = 0; i < coarse_elements.Size(); i++)
4594  {
4595  int geom = elements[coarse_elements[i]].geom;
4596  TraverseRefinements(coarse_elements[i], i, ref_path, path_map[geom]);
4597  used_geoms |= (1 << geom);
4598  }
4599 
4600  for (int g = 0; g < Geometry::NumGeom; g++)
4601  {
4602  if (used_geoms & (1 << g))
4603  {
4604  Geometry::Type geom = Geometry::Type(g);
4605  const PointMatrix &identity = GetGeomIdentity(geom);
4606 
4608  .SetSize(Dim, identity.np, path_map[g].size());
4609 
4610  // calculate the point matrices
4611  RefPathMap::iterator it;
4612  for (it = path_map[g].begin(); it != path_map[g].end(); ++it)
4613  {
4614  GetPointMatrix(geom, it->first.c_str(),
4615  transforms.point_matrices[g](it->second-1));
4616  }
4617  }
4618  }
4619  }
4620  return transforms;
4621 }
4622 
4624 {
4625  MFEM_VERIFY(transforms.embeddings.Size() || !leaf_elements.Size(),
4626  "GetDerefinementTransforms() must be preceded by Derefine().");
4627 
4628  if (!transforms.IsInitialized())
4629  {
4630  std::map<int, int> mat_no[Geometry::NumGeom];
4631  for (int g = 0; g < Geometry::NumGeom; g++)
4632  {
4633  mat_no[g][0] = 1; // 0 == identity
4634  }
4635 
4636  // assign numbers to the different matrices used
4637  for (int i = 0; i < transforms.embeddings.Size(); i++)
4638  {
4639  Embedding &emb = transforms.embeddings[i];
4640  int code = emb.matrix; // see SetDerefMatrixCodes()
4641  if (code)
4642  {
4643  int &matrix = mat_no[emb.geom][code];
4644  if (!matrix) { matrix = mat_no[emb.geom].size(); }
4645 
4646  emb.matrix = matrix - 1;
4647  }
4648  }
4649 
4650  for (int g = 0; g < Geometry::NumGeom; g++)
4651  {
4652  if (Geoms & (1 << g))
4653  {
4654  Geometry::Type geom = Geometry::Type(g);
4655  const PointMatrix &identity = GetGeomIdentity(geom);
4656 
4658  .SetSize(Dim, identity.np, mat_no[geom].size());
4659 
4660  // calculate point matrices
4661  for (auto it = mat_no[geom].begin(); it != mat_no[geom].end(); ++it)
4662  {
4663  char path[3] = { 0, 0, 0 };
4664 
4665  int code = it->first;
4666  if (code)
4667  {
4668  path[0] = code >> 4; // ref_type (see SetDerefMatrixCodes())
4669  path[1] = code & 0xf; // child
4670  }
4671 
4672  GetPointMatrix(geom, path,
4673  transforms.point_matrices[geom](it->second-1));
4674  }
4675  }
4676  }
4677  }
4678  return transforms;
4679 }
4680 
4682  bool want_ghosts) const
4683 {
4684  Array<Connection> conn;
4685  conn.Reserve(embeddings.Size());
4686 
4687  int max_parent = -1;
4688  for (int i = 0; i < embeddings.Size(); i++)
4689  {
4690  const Embedding &emb = embeddings[i];
4691  if ((emb.parent >= 0) &&
4692  (!emb.ghost || want_ghosts))
4693  {
4694  conn.Append(Connection(emb.parent, i));
4695  max_parent = std::max(emb.parent, max_parent);
4696  }
4697  }
4698 
4699  conn.Sort(); // NOTE: unique is not necessary
4700  coarse_to_fine.MakeFromList(max_parent+1, conn);
4701 }
4702 
4704 {
4706  transforms.Clear();
4707 }
4708 
4710 {
4711  for (int i = 0; i < Geometry::NumGeom; i++)
4712  {
4713  point_matrices[i].SetSize(0, 0, 0);
4714  }
4715  embeddings.DeleteAll();
4716 }
4717 
4719 {
4720  // return true if point matrices are present for any geometry
4721  for (int i = 0; i < Geometry::NumGeom; i++)
4722  {
4723  if (point_matrices[i].SizeK()) { return true; }
4724  }
4725  return false;
4726 }
4727 
4729 {
4730  for (int g = 0; g < Geometry::NumGeom; ++g)
4731  {
4732  a.point_matrices[g].Swap(b.point_matrices[g]);
4733  }
4734  Swap(a.embeddings, b.embeddings);
4735 }
4736 
4737 
4738 //// SFC Ordering //////////////////////////////////////////////////////////////
4739 
4740 static int sgn(int x)
4741 {
4742  return (x < 0) ? -1 : (x > 0) ? 1 : 0;
4743 }
4744 
4745 static void HilbertSfc2D(int x, int y, int ax, int ay, int bx, int by,
4746  Array<int> &coords)
4747 {
4748  int w = std::abs(ax + ay);
4749  int h = std::abs(bx + by);
4750 
4751  int dax = sgn(ax), day = sgn(ay); // unit major direction ("right")
4752  int dbx = sgn(bx), dby = sgn(by); // unit orthogonal direction ("up")
4753 
4754  if (h == 1) // trivial row fill
4755  {
4756  for (int i = 0; i < w; i++, x += dax, y += day)
4757  {
4758  coords.Append(x);
4759  coords.Append(y);
4760  }
4761  return;
4762  }
4763  if (w == 1) // trivial column fill
4764  {
4765  for (int i = 0; i < h; i++, x += dbx, y += dby)
4766  {
4767  coords.Append(x);
4768  coords.Append(y);
4769  }
4770  return;
4771  }
4772 
4773  int ax2 = ax/2, ay2 = ay/2;
4774  int bx2 = bx/2, by2 = by/2;
4775 
4776  int w2 = std::abs(ax2 + ay2);
4777  int h2 = std::abs(bx2 + by2);
4778 
4779  if (2*w > 3*h) // long case: split in two parts only
4780  {
4781  if ((w2 & 0x1) && (w > 2))
4782  {
4783  ax2 += dax, ay2 += day; // prefer even steps
4784  }
4785 
4786  HilbertSfc2D(x, y, ax2, ay2, bx, by, coords);
4787  HilbertSfc2D(x+ax2, y+ay2, ax-ax2, ay-ay2, bx, by, coords);
4788  }
4789  else // standard case: one step up, one long horizontal step, one step down
4790  {
4791  if ((h2 & 0x1) && (h > 2))
4792  {
4793  bx2 += dbx, by2 += dby; // prefer even steps
4794  }
4795 
4796  HilbertSfc2D(x, y, bx2, by2, ax2, ay2, coords);
4797  HilbertSfc2D(x+bx2, y+by2, ax, ay, bx-bx2, by-by2, coords);
4798  HilbertSfc2D(x+(ax-dax)+(bx2-dbx), y+(ay-day)+(by2-dby),
4799  -bx2, -by2, -(ax-ax2), -(ay-ay2), coords);
4800  }
4801 }
4802 
4803 static void HilbertSfc3D(int x, int y, int z,
4804  int ax, int ay, int az,
4805  int bx, int by, int bz,
4806  int cx, int cy, int cz,
4807  Array<int> &coords)
4808 {
4809  int w = std::abs(ax + ay + az);
4810  int h = std::abs(bx + by + bz);
4811  int d = std::abs(cx + cy + cz);
4812 
4813  int dax = sgn(ax), day = sgn(ay), daz = sgn(az); // unit major dir ("right")
4814  int dbx = sgn(bx), dby = sgn(by), dbz = sgn(bz); // unit ortho dir ("forward")
4815  int dcx = sgn(cx), dcy = sgn(cy), dcz = sgn(cz); // unit ortho dir ("up")
4816 
4817  // trivial row/column fills
4818  if (h == 1 && d == 1)
4819  {
4820  for (int i = 0; i < w; i++, x += dax, y += day, z += daz)
4821  {
4822  coords.Append(x);
4823  coords.Append(y);
4824  coords.Append(z);
4825  }
4826  return;
4827  }
4828  if (w == 1 && d == 1)
4829  {
4830  for (int i = 0; i < h; i++, x += dbx, y += dby, z += dbz)
4831  {
4832  coords.Append(x);
4833  coords.Append(y);
4834  coords.Append(z);
4835  }
4836  return;
4837  }
4838  if (w == 1 && h == 1)
4839  {
4840  for (int i = 0; i < d; i++, x += dcx, y += dcy, z += dcz)
4841  {
4842  coords.Append(x);
4843  coords.Append(y);
4844  coords.Append(z);
4845  }
4846  return;
4847  }
4848 
4849  int ax2 = ax/2, ay2 = ay/2, az2 = az/2;
4850  int bx2 = bx/2, by2 = by/2, bz2 = bz/2;
4851  int cx2 = cx/2, cy2 = cy/2, cz2 = cz/2;
4852 
4853  int w2 = std::abs(ax2 + ay2 + az2);
4854  int h2 = std::abs(bx2 + by2 + bz2);
4855  int d2 = std::abs(cx2 + cy2 + cz2);
4856 
4857  // prefer even steps
4858  if ((w2 & 0x1) && (w > 2))
4859  {
4860  ax2 += dax, ay2 += day, az2 += daz;
4861  }
4862  if ((h2 & 0x1) && (h > 2))
4863  {
4864  bx2 += dbx, by2 += dby, bz2 += dbz;
4865  }
4866  if ((d2 & 0x1) && (d > 2))
4867  {
4868  cx2 += dcx, cy2 += dcy, cz2 += dcz;
4869  }
4870 
4871  // wide case, split in w only
4872  if ((2*w > 3*h) && (2*w > 3*d))
4873  {
4874  HilbertSfc3D(x, y, z,
4875  ax2, ay2, az2,
4876  bx, by, bz,
4877  cx, cy, cz, coords);
4878 
4879  HilbertSfc3D(x+ax2, y+ay2, z+az2,
4880  ax-ax2, ay-ay2, az-az2,
4881  bx, by, bz,
4882  cx, cy, cz, coords);
4883  }
4884  // do not split in d
4885  else if (3*h > 4*d)
4886  {
4887  HilbertSfc3D(x, y, z,
4888  bx2, by2, bz2,
4889  cx, cy, cz,
4890  ax2, ay2, az2, coords);
4891 
4892  HilbertSfc3D(x+bx2, y+by2, z+bz2,
4893  ax, ay, az,
4894  bx-bx2, by-by2, bz-bz2,
4895  cx, cy, cz, coords);
4896 
4897  HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
4898  y+(ay-day)+(by2-dby),
4899  z+(az-daz)+(bz2-dbz),
4900  -bx2, -by2, -bz2,
4901  cx, cy, cz,
4902  -(ax-ax2), -(ay-ay2), -(az-az2), coords);
4903  }
4904  // do not split in h
4905  else if (3*d > 4*h)
4906  {
4907  HilbertSfc3D(x, y, z,
4908  cx2, cy2, cz2,
4909  ax2, ay2, az2,
4910  bx, by, bz, coords);
4911 
4912  HilbertSfc3D(x+cx2, y+cy2, z+cz2,
4913  ax, ay, az,
4914  bx, by, bz,
4915  cx-cx2, cy-cy2, cz-cz2, coords);
4916 
4917  HilbertSfc3D(x+(ax-dax)+(cx2-dcx),
4918  y+(ay-day)+(cy2-dcy),
4919  z+(az-daz)+(cz2-dcz),
4920  -cx2, -cy2, -cz2,
4921  -(ax-ax2), -(ay-ay2), -(az-az2),
4922  bx, by, bz, coords);
4923  }
4924  // regular case, split in all w/h/d
4925  else
4926  {
4927  HilbertSfc3D(x, y, z,
4928  bx2, by2, bz2,
4929  cx2, cy2, cz2,
4930  ax2, ay2, az2, coords);
4931 
4932  HilbertSfc3D(x+bx2, y+by2, z+bz2,
4933  cx, cy, cz,
4934  ax2, ay2, az2,
4935  bx-bx2, by-by2, bz-bz2, coords);
4936 
4937  HilbertSfc3D(x+(bx2-dbx)+(cx-dcx),
4938  y+(by2-dby)+(cy-dcy),
4939  z+(bz2-dbz)+(cz-dcz),
4940  ax, ay, az,
4941  -bx2, -by2, -bz2,
4942  -(cx-cx2), -(cy-cy2), -(cz-cz2), coords);
4943 
4944  HilbertSfc3D(x+(ax-dax)+bx2+(cx-dcx),
4945  y+(ay-day)+by2+(cy-dcy),
4946  z+(az-daz)+bz2+(cz-dcz),
4947  -cx, -cy, -cz,
4948  -(ax-ax2), -(ay-ay2), -(az-az2),
4949  bx-bx2, by-by2, bz-bz2, coords);
4950 
4951  HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
4952  y+(ay-day)+(by2-dby),
4953  z+(az-daz)+(bz2-dbz),
4954  -bx2, -by2, -bz2,
4955  cx2, cy2, cz2,
4956  -(ax-ax2), -(ay-ay2), -(az-az2), coords);
4957  }
4958 }
4959 
4960 void NCMesh::GridSfcOrdering2D(int width, int height, Array<int> &coords)
4961 {
4962  coords.SetSize(0);
4963  coords.Reserve(2*width*height);
4964 
4965  if (width >= height)
4966  {
4967  HilbertSfc2D(0, 0, width, 0, 0, height, coords);
4968  }
4969  else
4970  {
4971  HilbertSfc2D(0, 0, 0, height, width, 0, coords);
4972  }
4973 }
4974 
4975 void NCMesh::GridSfcOrdering3D(int width, int height, int depth,
4976  Array<int> &coords)
4977 {
4978  coords.SetSize(0);
4979  coords.Reserve(3*width*height*depth);
4980 
4981  if (width >= height && width >= depth)
4982  {
4983  HilbertSfc3D(0, 0, 0,
4984  width, 0, 0,
4985  0, height, 0,
4986  0, 0, depth, coords);
4987  }
4988  else if (height >= width && height >= depth)
4989  {
4990  HilbertSfc3D(0, 0, 0,
4991  0, height, 0,
4992  width, 0, 0,
4993  0, 0, depth, coords);
4994  }
4995  else // depth >= width && depth >= height
4996  {
4997  HilbertSfc3D(0, 0, 0,
4998  0, 0, depth,
4999  width, 0, 0,
5000  0, height, 0, coords);
5001  }
5002 }
5003 
5004 
5005 //// Utility ///////////////////////////////////////////////////////////////////
5006 
5007 void NCMesh::GetEdgeVertices(const MeshId &edge_id, int vert_index[2],
5008  bool oriented) const
5009 {
5010  const Element &el = elements[edge_id.element];
5011  const GeomInfo& gi = GI[el.Geom()];
5012  const int* ev = gi.edges[(int) edge_id.local];
5013 
5014  int n0 = el.node[ev[0]], n1 = el.node[ev[1]];
5015  if (n0 > n1) { std::swap(n0, n1); }
5016 
5017  vert_index[0] = nodes[n0].vert_index;
5018  vert_index[1] = nodes[n1].vert_index;
5019 
5020  if (oriented && vert_index[0] > vert_index[1])
5021  {
5022  std::swap(vert_index[0], vert_index[1]);
5023  }
5024 }
5025 
5027 {
5028  const Element &el = elements[edge_id.element];
5029  const GeomInfo& gi = GI[el.Geom()];
5030  const int* ev = gi.edges[(int) edge_id.local];
5031 
5032  int v0 = nodes[el.node[ev[0]]].vert_index;
5033  int v1 = nodes[el.node[ev[1]]].vert_index;
5034 
5035  return ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1])) ? -1 : 1;
5036 }
5037 
5039  int vert_index[4], int edge_index[4],
5040  int edge_orientation[4]) const
5041 {
5042  MFEM_ASSERT(Dim >= 3, "");
5043 
5044  const Element &el = elements[face_id.element];
5045  const GeomInfo& gi = GI[el.Geom()];
5046 
5047  const int *fv = gi.faces[(int) face_id.local];
5048  const int nfv = gi.nfv[(int) face_id.local];
5049 
5050  vert_index[3] = edge_index[3] = -1;
5051  edge_orientation[3] = 0;
5052 
5053  for (int i = 0; i < nfv; i++)
5054  {
5055  vert_index[i] = nodes[el.node[fv[i]]].vert_index;
5056  }
5057 
5058  for (int i = 0; i < nfv; i++)
5059  {
5060  int j = i+1;
5061  if (j >= nfv) { j = 0; }
5062 
5063  int n1 = el.node[fv[i]];
5064  int n2 = el.node[fv[j]];
5065 
5066  const Node* en = nodes.Find(n1, n2);
5067  MFEM_ASSERT(en != NULL, "edge not found.");
5068 
5069  edge_index[i] = en->edge_index;
5070  edge_orientation[i] = (vert_index[i] < vert_index[j]) ? 1 : -1;
5071  }
5072 
5073  return nfv;
5074 }
5075 
5076 int NCMesh::GetEdgeMaster(int node) const
5077 {
5078  MFEM_ASSERT(node >= 0, "edge node not found.");
5079  const Node &nd = nodes[node];
5080 
5081  int p1 = nd.p1, p2 = nd.p2;
5082  MFEM_ASSERT(p1 != p2, "invalid edge node.");
5083 
5084  const Node &n1 = nodes[p1], &n2 = nodes[p2];
5085 
5086  int n1p1 = n1.p1, n1p2 = n1.p2;
5087  int n2p1 = n2.p1, n2p2 = n2.p2;
5088 
5089  if ((n2p1 != n2p2) && (p1 == n2p1 || p1 == n2p2))
5090  {
5091  // n1 is parent of n2:
5092  // (n1)--(nd)--(n2)------(*)
5093  if (n2.HasEdge()) { return p2; }
5094  else { return GetEdgeMaster(p2); }
5095  }
5096 
5097  if ((n1p1 != n1p2) && (p2 == n1p1 || p2 == n1p2))
5098  {
5099  // n2 is parent of n1:
5100  // (n2)--(nd)--(n1)------(*)
5101  if (n1.HasEdge()) { return p1; }
5102  else { return GetEdgeMaster(p1); }
5103  }
5104 
5105  return -1;
5106 }
5107 
5108 int NCMesh::GetEdgeMaster(int v1, int v2) const
5109 {
5110  int node = nodes.FindId(vertex_nodeId[v1], vertex_nodeId[v2]);
5111  MFEM_ASSERT(node >= 0 && nodes[node].HasEdge(), "(v1, v2) is not an edge.");
5112 
5113  int master = GetEdgeMaster(node);
5114  return (master >= 0) ? nodes[master].edge_index : -1;
5115 }
5116 
5117 int NCMesh::GetElementDepth(int i) const
5118 {
5119  int elem = leaf_elements[i];
5120  int depth = 0, parent;
5121  while ((parent = elements[elem].parent) != -1)
5122  {
5123  elem = parent;
5124  depth++;
5125  }
5126  return depth;
5127 }
5128 
5130 {
5131  int elem = leaf_elements[i];
5132  int parent, reduction = 1;
5133  while ((parent = elements[elem].parent) != -1)
5134  {
5135  if (elements[parent].ref_type & 1) { reduction *= 2; }
5136  if (elements[parent].ref_type & 2) { reduction *= 2; }
5137  if (elements[parent].ref_type & 4) { reduction *= 2; }
5138  elem = parent;
5139  }
5140  return reduction;
5141 }
5142 
5144  Array<int> &face_indices,
5145  Array<int> &face_attribs) const
5146 {
5147  const Element &el = elements[leaf_elements[leaf_elem]];
5148  const GeomInfo& gi = GI[el.Geom()];
5149 
5150  face_indices.SetSize(gi.nf);
5151  face_attribs.SetSize(gi.nf);
5152 
5153  for (int i = 0; i < gi.nf; i++)
5154  {
5155  const int* fv = gi.faces[i];
5156  const Face *face = faces.Find(el.node[fv[0]], el.node[fv[1]],
5157  el.node[fv[2]], el.node[fv[3]]);
5158  MFEM_ASSERT(face, "face not found");
5159  face_indices[i] = face->index;
5160  face_attribs[i] = face->attribute;
5161  }
5162 }
5163 
5164 void NCMesh::FindFaceNodes(int face, int node[4])
5165 {
5166  // Obtain face nodes from one of its elements (note that face->p1, p2, p3
5167  // cannot be used directly since they are not in order and p4 is missing).
5168 
5169  Face &fa = faces[face];
5170 
5171  int elem = fa.elem[0];
5172  if (elem < 0) { elem = fa.elem[1]; }
5173  MFEM_ASSERT(elem >= 0, "Face has no elements?");
5174 
5175  Element &el = elements[elem];
5176  int f = find_local_face(el.Geom(),
5177  find_node(el, fa.p1),
5178  find_node(el, fa.p2),
5179  find_node(el, fa.p3));
5180 
5181  const int* fv = GI[el.Geom()].faces[f];
5182  for (int i = 0; i < 4; i++)
5183  {
5184  node[i] = el.node[fv[i]];
5185  }
5186 }
5187 
5188 void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
5189  Array<int> &bdr_vertices, Array<int> &bdr_edges)
5190 {
5191  bdr_vertices.SetSize(0);
5192  bdr_edges.SetSize(0);
5193 
5194  if (Dim == 3)
5195  {
5196  GetFaceList(); // make sure 'boundary_faces' is up to date
5197 
5198  for (int i = 0; i < boundary_faces.Size(); i++)
5199  {
5200  int face = boundary_faces[i];
5201  if (bdr_attr_is_ess[faces[face].attribute - 1])
5202  {
5203  int node[4];
5204  FindFaceNodes(face, node);
5205  int nfv = (node[3] < 0) ? 3 : 4;
5206 
5207  for (int j = 0; j < nfv; j++)
5208  {
5209  bdr_vertices.Append(nodes[node[j]].vert_index);
5210 
5211  int enode = nodes.FindId(node[j], node[(j+1) % nfv]);
5212  MFEM_ASSERT(enode >= 0 && nodes[enode].HasEdge(), "Edge not found.");
5213  bdr_edges.Append(nodes[enode].edge_index);
5214 
5215  while ((enode = GetEdgeMaster(enode)) >= 0)
5216  {
5217  // append master edges that may not be accessible from any
5218  // boundary element, this happens in 3D in re-entrant corners
5219  bdr_edges.Append(nodes[enode].edge_index);
5220  }
5221  }
5222  }
5223  }
5224  }
5225  else if (Dim == 2)
5226  {
5227  GetEdgeList(); // make sure 'boundary_faces' is up to date
5228 
5229  for (int i = 0; i < boundary_faces.Size(); i++)
5230  {
5231  int face = boundary_faces[i];
5232  Face &fc = faces[face];
5233  if (bdr_attr_is_ess[fc.attribute - 1])
5234  {
5235  bdr_vertices.Append(nodes[fc.p1].vert_index);
5236  bdr_vertices.Append(nodes[fc.p3].vert_index);
5237  }
5238  }
5239  }
5240 
5241  bdr_vertices.Sort();
5242  bdr_vertices.Unique();
5243 
5244  bdr_edges.Sort();
5245  bdr_edges.Unique();
5246 }
5247 
5248 static int max4(int a, int b, int c, int d)
5249 {
5250  return std::max(std::max(a, b), std::max(c, d));
5251 }
5252 static int max6(int a, int b, int c, int d, int e, int f)
5253 {
5254  return std::max(max4(a, b, c, d), std::max(e, f));
5255 }
5256 static int max8(int a, int b, int c, int d, int e, int f, int g, int h)
5257 {
5258  return std::max(max4(a, b, c, d), max4(e, f, g, h));
5259 }
5260 
5261 int NCMesh::EdgeSplitLevel(int vn1, int vn2) const
5262 {
5263  int mid = nodes.FindId(vn1, vn2);
5264  if (mid < 0 || !nodes[mid].HasVertex()) { return 0; }
5265  return 1 + std::max(EdgeSplitLevel(vn1, mid), EdgeSplitLevel(mid, vn2));
5266 }
5267 
5268 int NCMesh::TriFaceSplitLevel(int vn1, int vn2, int vn3) const
5269 {
5270  int mid[3];
5271  if (TriFaceSplit(vn1, vn2, vn3, mid) &&
5272  faces.FindId(vn1, vn2, vn3) < 0)
5273  {
5274  return 1 + max4(TriFaceSplitLevel(vn1, mid[0], mid[2]),
5275  TriFaceSplitLevel(mid[0], vn2, mid[1]),
5276  TriFaceSplitLevel(mid[2], mid[1], vn3),
5277  TriFaceSplitLevel(mid[0], mid[1], mid[2]));
5278  }
5279  else // not split
5280  {
5281  return 0;
5282  }
5283 }
5284 
5285 void NCMesh::QuadFaceSplitLevel(int vn1, int vn2, int vn3, int vn4,
5286  int& h_level, int& v_level) const
5287 {
5288  int hl1, hl2, vl1, vl2;
5289  int mid[5];
5290 
5291  switch (QuadFaceSplitType(vn1, vn2, vn3, vn4, mid))
5292  {
5293  case 0: // not split
5294  h_level = v_level = 0;
5295  break;
5296 
5297  case 1: // vertical
5298  QuadFaceSplitLevel(vn1, mid[0], mid[2], vn4, hl1, vl1);
5299  QuadFaceSplitLevel(mid[0], vn2, vn3, mid[2], hl2, vl2);
5300  h_level = std::max(hl1, hl2);
5301  v_level = std::max(vl1, vl2) + 1;
5302  break;
5303 
5304  default: // horizontal
5305  QuadFaceSplitLevel(vn1, vn2, mid[1], mid[3], hl1, vl1);
5306  QuadFaceSplitLevel(mid[3], mid[1], vn3, vn4, hl2, vl2);
5307  h_level = std::max(hl1, hl2) + 1;
5308  v_level = std::max(vl1, vl2);
5309  }
5310 }
5311 
5312 void NCMesh::CountSplits(int elem, int splits[3]) const
5313 {
5314  const Element &el = elements[elem];
5315  const int* node = el.node;
5316  GeomInfo& gi = GI[el.Geom()];
5317 
5318  int elevel[MaxElemEdges];
5319  for (int i = 0; i < gi.ne; i++)
5320  {
5321  const int* ev = gi.edges[i];
5322  elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
5323  }
5324 
5325  int flevel[MaxElemFaces][2];
5326  if (Dim >= 3)
5327  {
5328  for (int i = 0; i < gi.nf; i++)
5329  {
5330  const int* fv = gi.faces[i];
5331  if (gi.nfv[i] == 4)
5332  {
5333  QuadFaceSplitLevel(node[fv[0]], node[fv[1]],
5334  node[fv[2]], node[fv[3]],
5335  flevel[i][1], flevel[i][0]);
5336  }
5337  else
5338  {
5339  flevel[i][1] = 0;
5340  flevel[i][0] =
5341  TriFaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]]);
5342  }
5343  }
5344  }
5345 
5346  if (el.Geom() == Geometry::CUBE)
5347  {
5348  splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
5349  elevel[0], elevel[2], elevel[4], elevel[6]);
5350 
5351  splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
5352  elevel[1], elevel[3], elevel[5], elevel[7]);
5353 
5354  splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
5355  elevel[8], elevel[9], elevel[10], elevel[11]);
5356  }
5357  else if (el.Geom() == Geometry::PRISM)
5358  {
5359  splits[0] = splits[1] =
5360  std::max(
5361  max6(flevel[0][0], flevel[1][0], 0,
5362  flevel[2][0], flevel[3][0], flevel[4][0]),
5363  max6(elevel[0], elevel[1], elevel[2],
5364  elevel[3], elevel[4], elevel[5]));
5365 
5366  splits[2] = max6(flevel[2][1], flevel[3][1], flevel[4][1],
5367  elevel[6], elevel[7], elevel[8]);
5368  }
5369  else if (el.Geom() == Geometry::PYRAMID)
5370  {
5371  splits[0] = std::max(
5372  max6(flevel[0][0], flevel[1][0], 0,
5373  flevel[2][0], flevel[3][0], flevel[4][0]),
5374  max8(elevel[0], elevel[1], elevel[2],
5375  elevel[3], elevel[4], elevel[5],
5376  elevel[6], elevel[7]));
5377 
5378  splits[1] = splits[0];
5379  splits[2] = splits[0];
5380  }
5381  else if (el.Geom() == Geometry::TETRAHEDRON)
5382  {
5383  splits[0] = std::max(
5384  max4(flevel[0][0], flevel[1][0], flevel[2][0], flevel[3][0]),
5385  max6(elevel[0], elevel[1], elevel[2],
5386  elevel[3], elevel[4], elevel[5]));
5387 
5388  splits[1] = splits[0];
5389  splits[2] = splits[0];
5390  }
5391  else if (el.Geom() == Geometry::SQUARE)
5392  {
5393  splits[0] = std::max(elevel[0], elevel[2]);
5394  splits[1] = std::max(elevel[1], elevel[3]);
5395  }
5396  else if (el.Geom() == Geometry::TRIANGLE)
5397  {
5398  splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
5399  splits[1] = splits[0];
5400  }
5401  else
5402  {
5403  MFEM_ABORT("Unsupported element geometry.");
5404  }
5405 }
5406 
5407 void NCMesh::GetLimitRefinements(Array<Refinement> &refinements, int max_level)
5408 {
5409  for (int i = 0; i < leaf_elements.Size(); i++)
5410  {
5411  if (IsGhost(elements[leaf_elements[i]])) { break; } // TODO: NElements
5412 
5413  int splits[3];
5414  CountSplits(leaf_elements[i], splits);
5415 
5416  char ref_type = 0;
5417  for (int k = 0; k < Dim; k++)
5418  {
5419  if (splits[k] > max_level)
5420  {
5421  ref_type |= (1 << k);
5422  }
5423  }
5424 
5425  if (ref_type)
5426  {
5427  if (Iso)
5428  {
5429  // iso meshes should only be modified by iso refinements
5430  ref_type = 7;
5431  }
5432  refinements.Append(Refinement(i, ref_type));
5433  }
5434  }
5435 }
5436 
5437 void NCMesh::LimitNCLevel(int max_nc_level)
5438 {
5439  MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
5440 
5441  while (1)
5442  {
5443  Array<Refinement> refinements;
5444  GetLimitRefinements(refinements, max_nc_level);
5445 
5446  if (!refinements.Size()) { break; }
5447 
5448  Refine(refinements);
5449  }
5450 }
5451 
5452 
5453 //// I/O ////////////////////////////////////////////////////////////////////////
5454 
5455 int NCMesh::PrintVertexParents(std::ostream *os) const
5456 {
5457  if (!os)
5458  {
5459  // count vertex nodes with parents
5460  int nv = 0;
5461  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5462  {
5463  if (node->HasVertex() && node->p1 != node->p2) { nv++; }
5464  }
5465  return nv;
5466  }
5467  else
5468  {
5469  // print the relations
5470  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5471  {
5472  if (node->HasVertex() && node->p1 != node->p2)
5473  {
5474  MFEM_ASSERT(nodes[node->p1].HasVertex(), "");
5475  MFEM_ASSERT(nodes[node->p2].HasVertex(), "");
5476 
5477  (*os) << node.index() << " " << node->p1 << " " << node->p2 << "\n";
5478  }
5479  }
5480  return 0;
5481  }
5482 }
5483 
5484 void NCMesh::LoadVertexParents(std::istream &input)
5485 {
5486  int nv;
5487  input >> nv;
5488  while (nv--)
5489  {
5490  int id, p1, p2;
5491  input >> id >> p1 >> p2;
5492  MFEM_VERIFY(input, "problem reading vertex parents.");
5493 
5494  MFEM_VERIFY(nodes.IdExists(id), "vertex " << id << " not found.");
5495  MFEM_VERIFY(nodes.IdExists(p1), "parent " << p1 << " not found.");
5496  MFEM_VERIFY(nodes.IdExists(p2), "parent " << p2 << " not found.");
5497 
5498  int check = nodes.FindId(p1, p2);
5499  MFEM_VERIFY(check < 0, "parents (" << p1 << ", " << p2 << ") already "
5500  "assigned to node " << check);
5501 
5502  // assign new parents for the node
5503  nodes.Reparent(id, p1, p2);
5504  }
5505 }
5506 
5507 int NCMesh::PrintBoundary(std::ostream *os) const
5508 {
5509  static const int nfv2geom[5] =
5510  {
5513  };
5514  int deg = (Dim == 2) ? 2 : 1; // for degenerate faces in 2D
5515 
5516  int count = 0;
5517  for (int i = 0; i < elements.Size(); i++)
5518  {
5519  const Element &el = elements[i];
5520  if (!el.IsLeaf()) { continue; }
5521 
5522  GeomInfo& gi = GI[el.Geom()];
5523  for (int k = 0; k < gi.nf; k++)
5524  {
5525  const int* fv = gi.faces[k];
5526  const int nfv = gi.nfv[k];
5527  const Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
5528  el.node[fv[2]], el.node[fv[3]]);
5529  MFEM_ASSERT(face != NULL, "face not found");
5530  if (face->Boundary())
5531  {
5532  if (!os) { count++; continue; }
5533 
5534  (*os) << face->attribute << " " << nfv2geom[nfv];
5535  for (int j = 0; j < nfv; j++)
5536  {
5537  (*os) << " " << el.node[fv[j*deg]];
5538  }
5539  (*os) << "\n";
5540  }
5541  }
5542  }
5543  return count;
5544 }
5545 
5546 void NCMesh::LoadBoundary(std::istream &input)
5547 {
5548  int nb, attr, geom;
5549  input >> nb;
5550  for (int i = 0; i < nb; i++)
5551  {
5552  input >> attr >> geom;
5553 
5554  int v1, v2, v3, v4;
5555  if (geom == Geometry::SQUARE)
5556  {
5557  input >> v1 >> v2 >> v3 >> v4;
5558  Face* face = faces.Get(v1, v2, v3, v4);
5559  face->attribute = attr;
5560  }
5561  else if (geom == Geometry::TRIANGLE)
5562  {
5563  input >> v1 >> v2 >> v3;
5564  Face* face = faces.Get(v1, v2, v3);
5565  face->attribute = attr;
5566  }
5567  else if (geom == Geometry::SEGMENT)
5568  {
5569  input >> v1 >> v2;
5570  Face* face = faces.Get(v1, v1, v2, v2);
5571  face->attribute = attr;
5572  }
5573  else if (geom == Geometry::POINT)
5574  {
5575  input >> v1;
5576  Face* face = faces.Get(v1, v1, v1, v1);
5577  face->attribute = attr;
5578  }
5579  else
5580  {
5581  MFEM_ABORT("unsupported boundary element geometry: " << geom);
5582  }
5583  }
5584 }
5585 
5586 void NCMesh::PrintCoordinates(std::ostream &os) const
5587 {
5588  int nv = coordinates.Size()/3;
5589  os << nv << "\n";
5590  if (!nv) { return; }
5591 
5592  os << spaceDim << "\n";
5593  for (int i = 0; i < nv; i++)
5594  {
5595  os << coordinates[3*i];
5596  for (int j = 1; j < spaceDim; j++)
5597  {
5598  os << " " << coordinates[3*i + j];
5599  }
5600  os << "\n";
5601  }
5602 }
5603 
5604 void NCMesh::LoadCoordinates(std::istream &input)
5605 {
5606  int nv;
5607  input >> nv;
5608  if (!nv) { return; }
5609 
5610  input >> spaceDim;
5611 
5612  coordinates.SetSize(nv * 3);
5613  coordinates = 0.0;
5614 
5615  for (int i = 0; i < nv; i++)
5616  {
5617  for (int j = 0; j < spaceDim; j++)
5618  {
5619  input >> coordinates[3*i + j];
5620  MFEM_VERIFY(input.good(), "unexpected EOF");
5621  }
5622  }
5623 }
5624 
5626 {
5627  for (int i = 0; i < root_state.Size(); i++)
5628  {
5629  if (root_state[i]) { return false; }
5630  }
5631  return true;
5632 }
5633 
5634 void NCMesh::Print(std::ostream &os) const
5635 {
5636  os << "MFEM NC mesh v1.0\n\n"
5637  "# NCMesh supported geometry types:\n"
5638  "# SEGMENT = 1\n"
5639  "# TRIANGLE = 2\n"
5640  "# SQUARE = 3\n"
5641  "# TETRAHEDRON = 4\n"
5642  "# CUBE = 5\n"
5643  "# PRISM = 6\n"
5644  "# PYRAMID = 7\n";
5645 
5646  os << "\ndimension\n" << Dim << "\n";
5647 
5648 #ifndef MFEM_USE_MPI
5649  if (MyRank != 0) // don't print this section in serial: default rank is 0
5650 #endif
5651  {
5652  os << "\nrank\n" << MyRank << "\n";
5653  }
5654 
5655  os << "\n# rank attr geom ref_type nodes/children";
5656  os << "\nelements\n" << elements.Size() << "\n";
5657 
5658  for (int i = 0; i < elements.Size(); i++)
5659  {
5660  const Element &el = elements[i];
5661  os << el.rank << " " << el.attribute << " ";
5662  if (el.parent == -2) { os << "-1\n"; continue; } // unused element
5663 
5664  os << int(el.geom) << " " << int(el.ref_type);
5665  for (int j = 0; j < MaxElemNodes && el.node[j] >= 0; j++)
5666  {
5667  os << " " << el.node[j];
5668  }
5669  os << "\n";
5670  }
5671 
5672  int nb = PrintBoundary(NULL);
5673  if (nb)
5674  {
5675  os << "\n# attr geom nodes";
5676  os << "\nboundary\n" << nb << "\n";
5677 
5678  PrintBoundary(&os);
5679  }
5680 
5681  int nvp = PrintVertexParents(NULL);
5682  if (nvp)
5683  {
5684  os << "\n# vert_id p1 p2";
5685  os << "\nvertex_parents\n" << nvp << "\n";
5686 
5687  PrintVertexParents(&os);
5688  }
5689 
5690  if (!ZeroRootStates()) // root_state section is optional
5691  {
5692  os << "\n# root element orientation";
5693  os << "\nroot_state\n" << root_state.Size() << "\n";
5694 
5695  for (int i = 0; i < root_state.Size(); i++)
5696  {
5697  os << root_state[i] << "\n";
5698  }
5699  }
5700 
5701  if (coordinates.Size())
5702  {
5703  os << "\n# top-level node coordinates";
5704  os << "\ncoordinates\n";
5705 
5706  PrintCoordinates(os);
5707  }
5708  else
5709  {
5710  // 'nodes' will be printed one level up by Mesh::Printer()
5711  }
5712 }
5713 
5715 {
5716  // initialize Element::parent
5717  for (int i = 0; i < elements.Size(); i++)
5718  {
5719  Element &el = elements[i];
5720  if (el.ref_type)
5721  {
5722  for (int j = 0; j < MaxElemChildren && el.child[j] >= 0; j++)
5723  {
5724  int child = el.child[j];
5725  MFEM_VERIFY(child < elements.Size(), "invalid mesh file: "
5726  "element " << i << " references invalid child " << child);
5727  elements[child].parent = i;
5728  }
5729  }
5730  }
5731 
5732  // count the root elements
5733  int nroots = 0;
5734  while (nroots < elements.Size() &&
5735  elements[nroots].parent == -1)
5736  {
5737  nroots++;
5738  }
5739  MFEM_VERIFY(nroots, "invalid mesh file: no root elements found.");
5740 
5741  // check that only the first 'nroot' elements are roots (have no parent)
5742  for (int i = nroots; i < elements.Size(); i++)
5743  {
5744  MFEM_VERIFY(elements[i].parent != -1,
5745  "invalid mesh file: only the first M elements can be roots. "
5746  "Found element " << i << " with no parent.");
5747  }
5748 
5749  // default root state is zero
5750  root_state.SetSize(nroots);
5751  root_state = 0;
5752 }
5753 
5755 {
5756  int ntop = 0;
5757  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5758  {
5759  if (node->p1 == node->p2) { ntop = node.index() + 1; }
5760  }
5761  return ntop;
5762 }
5763 
5764 NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
5765  : spaceDim(0), MyRank(0), Iso(true), Legacy(false)
5766 {
5767  is_nc = 1;
5768  if (version == 1) // old MFEM mesh v1.1 format
5769  {
5770  LoadLegacyFormat(input, curved, is_nc);
5771  Legacy = true;
5772  return;
5773  }
5774 
5775  MFEM_ASSERT(version == 10, "");
5776  std::string ident;
5777  int count;
5778 
5779  // load dimension
5780  skip_comment_lines(input, '#');
5781  input >> ident;
5782  MFEM_VERIFY(ident == "dimension", "Invalid mesh file: " << ident);
5783  input >> Dim;
5784 
5785  // load rank, if present
5786  skip_comment_lines(input, '#');
5787  input >> ident;
5788  if (ident == "rank")
5789  {
5790  input >> MyRank;
5791  MFEM_VERIFY(MyRank >= 0, "Invalid rank");
5792 
5793  skip_comment_lines(input, '#');
5794  input >> ident;
5795  }
5796 
5797  // load file SFC version, if present (for future changes to SFC ordering)
5798  if (ident == "sfc_version")
5799  {
5800  int sfc_version; // TODO future: store as class member
5801  input >> sfc_version;
5802  MFEM_VERIFY(sfc_version == 0,
5803  "Unsupported mesh file SFC version (" << sfc_version << "). "
5804  "Please update MFEM.");
5805 
5806  skip_comment_lines(input, '#');
5807  input >> ident;
5808  }
5809 
5810  // load elements
5811  MFEM_VERIFY(ident == "elements", "Invalid mesh file: " << ident);
5812  input >> count;
5813  for (int i = 0; i < count; i++)
5814  {
5815  int rank, attr, geom, ref_type;
5816  input >> rank >> attr >> geom;
5817 
5818  Geometry::Type type = Geometry::Type(geom);
5819  elements.Append(Element(type, attr));
5820 
5821  MFEM_ASSERT(elements.Size() == i+1, "");
5822  Element &el = elements[i];
5823  el.rank = rank;
5824 
5825  if (geom >= 0)
5826  {
5827  CheckSupportedGeom(type);
5828  GI[geom].InitGeom(type);
5829 
5830  input >> ref_type;
5831  MFEM_VERIFY(ref_type >= 0 && ref_type < 8, "");
5832  el.ref_type = ref_type;
5833 
5834  if (ref_type) // refined element
5835  {
5836  for (int j = 0; j < ref_type_num_children[ref_type]; j++)
5837  {
5838  input >> el.child[j];
5839  }
5840  if (Dim == 3 && ref_type != 7) { Iso = false; }
5841  }
5842  else // leaf element
5843  {
5844  for (int j = 0; j < GI<