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