MFEM  v4.5.1
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 #else
1560  MFEM_CONTRACT_VAR(nforced);
1561 #endif
1562 
1563  ref_stack.DeleteAll();
1564  shadow.DeleteAll();
1565 
1566  Update();
1567 }
1568 
1569 
1570 //// Derefinement //////////////////////////////////////////////////////////////
1571 
1573 {
1574  if (!el.ref_type) { return el.node[index]; }
1575 
1576  // need to retrieve node from a child element (there is always a child
1577  // that inherited the parent's corner under the same index)
1578  int ch;
1579  switch (el.Geom())
1580  {
1581  case Geometry::CUBE:
1582  ch = el.child[hex_deref_table[el.ref_type - 1][index]];
1583  break;
1584 
1585  case Geometry::PRISM:
1586  ch = prism_deref_table[el.ref_type - 1][index];
1587  MFEM_ASSERT(ch != -1, "");
1588  ch = el.child[ch];
1589  break;
1590 
1591  case Geometry::SQUARE:
1592  ch = el.child[quad_deref_table[el.ref_type - 1][index]];
1593  break;
1594 
1595  case Geometry::TETRAHEDRON:
1596  case Geometry::TRIANGLE:
1597  ch = el.child[index];
1598  break;
1599 
1600  default:
1601  ch = 0; // suppress compiler warning
1602  MFEM_ABORT("Unsupported element geometry.");
1603  }
1604  return RetrieveNode(elements[ch], index);
1605 }
1606 
1607 
1609 {
1610  Element &el = elements[elem];
1611  if (!el.ref_type) { return; }
1612 
1613  int child[8];
1614  std::memcpy(child, el.child, sizeof(child));
1615 
1616  // first make sure that all children are leaves, derefine them if not
1617  for (int i = 0; i < 8 && child[i] >= 0; i++)
1618  {
1619  if (elements[child[i]].ref_type)
1620  {
1621  DerefineElement(child[i]);
1622  }
1623  }
1624 
1625  int faces_attribute[6];
1626  int ref_type_key = el.ref_type - 1;
1627 
1628  for (int i = 0; i < 8; i++) { el.node[i] = -1; }
1629 
1630  // retrieve original corner nodes and face attributes from the children
1631  if (el.Geom() == Geometry::CUBE)
1632  {
1633  // Sets corner nodes from childs
1634  constexpr int nb_cube_childs = 8;
1635  for (int i = 0; i < nb_cube_childs; i++)
1636  {
1637  const int child_local_index = hex_deref_table[ref_type_key][i];
1638  const int child_global_index = child[child_local_index];
1639  Element &ch = elements[child_global_index];
1640  el.node[i] = ch.node[i];
1641  }
1642  // Sets faces attributes from childs' faces
1643  constexpr int nb_cube_faces = 6;
1644  for (int i = 0; i < nb_cube_faces; i++)
1645  {
1646  const int child_local_index = hex_deref_table[ref_type_key]
1647  [i + nb_cube_childs];
1648  const int child_global_index = child[child_local_index];
1649  Element &ch = elements[child_global_index];
1650  const int* fv = GI[el.Geom()].faces[i];
1651  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1652  ch.node[fv[2]], ch.node[fv[3]])
1653  ->attribute;
1654  }
1655  }
1656  else if (el.Geom() == Geometry::PRISM)
1657  {
1658  MFEM_ASSERT(prism_deref_table[ref_type_key][0] != -1,
1659  "invalid prism refinement");
1660  constexpr int nb_prism_childs = 6;
1661  for (int i = 0; i < nb_prism_childs; i++)
1662  {
1663  const int child_local_index = prism_deref_table[ref_type_key][i];
1664  const int child_global_index = child[child_local_index];
1665  Element &ch = elements[child_global_index];
1666  el.node[i] = ch.node[i];
1667  }
1668  el.node[6] = el.node[7] = -1;
1669 
1670  constexpr int nb_prism_faces = 5;
1671  for (int i = 0; i < nb_prism_faces; i++)
1672  {
1673  const int child_local_index = prism_deref_table[ref_type_key]
1674  [i + nb_prism_childs];
1675  const int child_global_index = child[child_local_index];
1676  Element &ch = elements[child_global_index];
1677  const int* fv = GI[el.Geom()].faces[i];
1678  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1679  ch.node[fv[2]], ch.node[fv[3]])
1680  ->attribute;
1681  }
1682  }
1683  else if (el.Geom() == Geometry::TETRAHEDRON)
1684  {
1685  for (int i = 0; i < 4; i++)
1686  {
1687  Element& ch1 = elements[child[i]];
1688  Element& ch2 = elements[child[(i+1) & 0x3]];
1689  el.node[i] = ch1.node[i];
1690  const int* fv = GI[el.Geom()].faces[i];
1691  faces_attribute[i] = faces.Find(ch2.node[fv[0]], ch2.node[fv[1]],
1692  ch2.node[fv[2]], ch2.node[fv[3]])
1693  ->attribute;
1694  }
1695  }
1696  else if (el.Geom() == Geometry::SQUARE)
1697  {
1698  constexpr int nb_square_childs = 4;
1699  for (int i = 0; i < nb_square_childs; i++)
1700  {
1701  const int child_local_index = quad_deref_table[ref_type_key][i];
1702  const int child_global_index = child[child_local_index];
1703  Element &ch = elements[child_global_index];
1704  el.node[i] = ch.node[i];
1705  }
1706  constexpr int nb_square_faces = 4;
1707  for (int i = 0; i < nb_square_faces; i++)
1708  {
1709  const int child_local_index = quad_deref_table[ref_type_key]
1710  [i + nb_square_childs];
1711  const int child_global_index = child[child_local_index];
1712  Element &ch = elements[child_global_index];
1713  const int* fv = GI[el.Geom()].faces[i];
1714  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1715  ch.node[fv[2]], ch.node[fv[3]])
1716  ->attribute;
1717  }
1718  }
1719  else if (el.Geom() == Geometry::TRIANGLE)
1720  {
1721  constexpr int nb_triangle_childs = 3;
1722  for (int i = 0; i < nb_triangle_childs; i++)
1723  {
1724  Element& ch = elements[child[i]];
1725  el.node[i] = ch.node[i];
1726  const int* fv = GI[el.Geom()].faces[i];
1727  faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
1728  ch.node[fv[2]], ch.node[fv[3]])
1729  ->attribute;
1730  }
1731  }
1732  else if (el.Geom() == Geometry::SEGMENT)
1733  {
1734  constexpr int nb_segment_childs = 2;
1735  for (int i = 0; i < nb_segment_childs; i++)
1736  {
1737  int ni = elements[child[i]].node[i];
1738  el.node[i] = ni;
1739  faces_attribute[i] = faces.Find(ni, ni, ni, ni)->attribute;
1740  }
1741  }
1742  else
1743  {
1744  MFEM_ABORT("Unsupported element geometry.");
1745  }
1746 
1747  // sign in to all nodes
1748  ReferenceElement(elem);
1749 
1750  int buf[8*6];
1751  Array<int> childFaces(buf, 8*6);
1752  childFaces.SetSize(0);
1753 
1754  // delete children, determine rank
1755  el.rank = std::numeric_limits<int>::max();
1756  for (int i = 0; i < 8 && child[i] >= 0; i++)
1757  {
1758  el.rank = std::min(el.rank, elements[child[i]].rank);
1759  UnreferenceElement(child[i], childFaces);
1760  FreeElement(child[i]);
1761  }
1762 
1763  RegisterFaces(elem, faces_attribute);
1764 
1765  // delete unused faces
1766  childFaces.Sort();
1767  childFaces.Unique();
1768  DeleteUnusedFaces(childFaces);
1769 
1770  el.ref_type = 0;
1771 }
1772 
1773 
1775 {
1776  Element &el = elements[elem];
1777  if (!el.ref_type) { return; }
1778 
1779  int total = 0, ref = 0, ghost = 0;
1780  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
1781  {
1782  total++;
1783  Element &ch = elements[el.child[i]];
1784  if (ch.ref_type) { ref++; break; }
1785  if (IsGhost(ch)) { ghost++; }
1786  }
1787 
1788  if (!ref && ghost < total)
1789  {
1790  // can be derefined, add to list
1791  int next_row = list.Size() ? (list.Last().from + 1) : 0;
1792  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
1793  {
1794  Element &ch = elements[el.child[i]];
1795  list.Append(Connection(next_row, ch.index));
1796  }
1797  }
1798  else
1799  {
1800  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
1801  {
1802  CollectDerefinements(el.child[i], list);
1803  }
1804  }
1805 }
1806 
1808 {
1809  Array<Connection> list;
1810  list.Reserve(leaf_elements.Size());
1811 
1812  for (int i = 0; i < root_state.Size(); i++)
1813  {
1814  CollectDerefinements(i, list);
1815  }
1816 
1817  int size = list.Size() ? (list.Last().from + 1) : 0;
1818  derefinements.MakeFromList(size, list);
1819  return derefinements;
1820 }
1821 
1822 void NCMesh::CheckDerefinementNCLevel(const Table &deref_table,
1823  Array<int> &level_ok, int max_nc_level)
1824 {
1825  level_ok.SetSize(deref_table.Size());
1826  for (int i = 0; i < deref_table.Size(); i++)
1827  {
1828  const int* fine = deref_table.GetRow(i), size = deref_table.RowSize(i);
1829  Element &parent = elements[elements[leaf_elements[fine[0]]].parent];
1830 
1831  int ok = 1;
1832  for (int j = 0; j < size; j++)
1833  {
1834  int splits[3];
1835  CountSplits(leaf_elements[fine[j]], splits);
1836 
1837  for (int k = 0; k < Dim; k++)
1838  {
1839  if ((parent.ref_type & (1 << k)) &&
1840  splits[k] >= max_nc_level)
1841  {
1842  ok = 0; break;
1843  }
1844  }
1845  if (!ok) { break; }
1846  }
1847  level_ok[i] = ok;
1848  }
1849 }
1850 
1851 void NCMesh::Derefine(const Array<int> &derefs)
1852 {
1853  MFEM_VERIFY(Dim < 3 || Iso,
1854  "derefinement of 3D anisotropic meshes not implemented yet.");
1855 
1857 
1858  Array<int> fine_coarse;
1859  leaf_elements.Copy(fine_coarse);
1860 
1861  // perform the derefinements
1862  for (int i = 0; i < derefs.Size(); i++)
1863  {
1864  int row = derefs[i];
1865  MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
1866  "invalid derefinement number.");
1867 
1868  const int* fine = derefinements.GetRow(row);
1869  int parent = elements[leaf_elements[fine[0]]].parent;
1870 
1871  // record the relation of the fine elements to their parent
1872  SetDerefMatrixCodes(parent, fine_coarse);
1873 
1874  DerefineElement(parent);
1875  }
1876 
1877  // update leaf_elements, Element::index etc.
1878  Update();
1879 
1880  // link old fine elements to the new coarse elements
1881  for (int i = 0; i < fine_coarse.Size(); i++)
1882  {
1883  transforms.embeddings[i].parent = elements[fine_coarse[i]].index;
1884  }
1885 }
1886 
1888 {
1889  int nfine = leaf_elements.Size();
1890 
1891  // this will tell GetDerefinementTransforms that transforms are not finished
1892  transforms.Clear();
1893 
1894  transforms.embeddings.SetSize(nfine);
1895  for (int i = 0; i < nfine; i++)
1896  {
1897  Embedding &emb = transforms.embeddings[i];
1898  emb.parent = -1;
1899  emb.matrix = 0;
1900  Element &el = elements[leaf_elements[i]];
1901  emb.geom = el.Geom();
1902  emb.ghost = IsGhost(el);
1903  }
1904 }
1905 
1906 void NCMesh::SetDerefMatrixCodes(int parent, Array<int> &fine_coarse)
1907 {
1908  // encode the ref_type and child number for GetDerefinementTransforms()
1909  Element &prn = elements[parent];
1910  for (int i = 0; i < 8 && prn.child[i] >= 0; i++)
1911  {
1912  Element &ch = elements[prn.child[i]];
1913  if (ch.index >= 0)
1914  {
1915  int code = (prn.ref_type << 4) | i;
1916  transforms.embeddings[ch.index].matrix = code;
1917  fine_coarse[ch.index] = parent;
1918  }
1919  }
1920 }
1921 
1922 
1923 //// Mesh Interface ////////////////////////////////////////////////////////////
1924 
1925 void NCMesh::CollectLeafElements(int elem, int state, Array<int> &ghosts,
1926  int &counter)
1927 {
1928  Element &el = elements[elem];
1929  if (!el.ref_type)
1930  {
1931  if (el.rank >= 0) // skip elements beyond the ghost layer in parallel
1932  {
1933  if (!IsGhost(el))
1934  {
1935  leaf_elements.Append(elem);
1936  }
1937  else
1938  {
1939  // in parallel (or in serial loading a parallel file), collect
1940  // elements of neighboring ranks in a separate array
1941  ghosts.Append(elem);
1942  }
1943 
1944  // assign the SFC index (temporarily, will be replaced by Mesh index)
1945  el.index = counter++;
1946  }
1947  else
1948  {
1949  // elements beyond the ghost layer are invalid and don't appear in
1950  // 'leaf_elements' (also for performance reasons)
1951  el.index = -1;
1952  }
1953  }
1954  else // Refined element
1955  {
1956  // in non-leaf elements, the 'rank' and 'index' members have no meaning
1957  el.rank = -1;
1958  el.index = -1;
1959 
1960  // recurse to subtrees; try to order leaf elements along a space-filling
1961  // curve by changing the order the children are visited at each level
1962  if (el.Geom() == Geometry::SQUARE && el.ref_type == Refinement::XY)
1963  {
1964  for (int i = 0; i < 4; i++)
1965  {
1966  int ch = quad_hilbert_child_order[state][i];
1967  int st = quad_hilbert_child_state[state][i];
1968  CollectLeafElements(el.child[ch], st, ghosts, counter);
1969  }
1970  }
1971  else if (el.Geom() == Geometry::CUBE && el.ref_type == Refinement::XYZ)
1972  {
1973  for (int i = 0; i < 8; i++)
1974  {
1975  int ch = hex_hilbert_child_order[state][i];
1976  int st = hex_hilbert_child_state[state][i];
1977  CollectLeafElements(el.child[ch], st, ghosts, counter);
1978  }
1979  }
1980  else // no space filling curve tables yet for remaining cases
1981  {
1982  for (int i = 0; i < 8; i++)
1983  {
1984  if (el.child[i] >= 0)
1985  {
1986  CollectLeafElements(el.child[i], state, ghosts, counter);
1987  }
1988  }
1989  }
1990  }
1991 }
1992 
1994 {
1995  Array<int> ghosts;
1996 
1997  // collect leaf elements in leaf_elements and ghosts elements in ghosts from
1998  // all roots
2000  for (int i = 0, counter = 0; i < root_state.Size(); i++)
2001  {
2002  CollectLeafElements(i, root_state[i], ghosts, counter);
2003  }
2004 
2006  NGhostElements = ghosts.Size();
2007 
2008  // append ghost elements at the end of 'leaf_element' (if any)
2009  // and assign the final (Mesh) indices of leaves
2010  leaf_elements.Append(ghosts);
2012 
2013  for (int i = 0; i < leaf_elements.Size(); i++)
2014  {
2015  Element &el = elements[leaf_elements[i]];
2016  leaf_sfc_index[i] = el.index;
2017  el.index = i;
2018  }
2019 }
2020 
2022 {
2023 #ifndef MFEM_NCMESH_OLD_VERTEX_ORDERING
2024  // This method assigns indices to vertices (Node::vert_index) that will
2025  // be seen by the Mesh class and the rest of MFEM. We must be careful to:
2026  //
2027  // 1. Stay compatible with the conforming code, which expects top-level
2028  // (original) vertices to be indexed first, otherwise GridFunctions
2029  // defined on a conforming mesh would no longer be valid when the
2030  // mesh is converted to an NC mesh.
2031  //
2032  // 2. Make sure serial NCMesh is compatible with the parallel ParNCMesh,
2033  // so it is possible to read parallel partial solutions in serial code
2034  // (e.g., serial GLVis). This means handling ghost elements, if present.
2035  //
2036  // 3. Assign vertices in a globally consistent order for parallel meshes:
2037  // if two vertices i,j are shared by two ranks r1,r2, and i<j on r1,
2038  // then i<j on r2 as well. This is true for top-level vertices but also
2039  // for the remaining shared vertices thanks to the globally consistent
2040  // SFC ordering of the leaf elements. This property reduces communication
2041  // and simplifies ParNCMesh.
2042 
2043  // STEP 1: begin by splitting vertices into 4 classes:
2044  // - local top-level vertices (code -1)
2045  // - local non-top level vertices (code -2)
2046  // - ghost (non-local) vertices (code -3)
2047  // - vertices beyond the ghost layer (code -4)
2048 
2049  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2050  {
2051  node->vert_index = -4; // assume beyond ghost layer
2052  }
2053 
2054  for (int i = 0; i < leaf_elements.Size(); i++)
2055  {
2056  Element &el = elements[leaf_elements[i]];
2057  for (int j = 0; j < GI[el.Geom()].nv; j++)
2058  {
2059  Node &nd = nodes[el.node[j]];
2060  if (el.rank == MyRank)
2061  {
2062  if (nd.p1 == nd.p2) // local top-level vertex
2063  {
2064  if (nd.vert_index < -1) { nd.vert_index = -1; }
2065  }
2066  else // local non-top-level vertex
2067  {
2068  if (nd.vert_index < -2) { nd.vert_index = -2; }
2069  }
2070  }
2071  else // ghost vertex
2072  {
2073  if (nd.vert_index < -3) { nd.vert_index = -3; }
2074  }
2075  }
2076  }
2077 
2078  // STEP 2: assign indices of top-level local vertices, in original order
2079 
2080  NVertices = 0;
2081  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2082  {
2083  if (node->vert_index == -1)
2084  {
2085  node->vert_index = NVertices++;
2086  }
2087  }
2088 
2089  // STEP 3: go over all elements (local and ghost) in SFC order and assign
2090  // remaining local vertices in that order.
2091 
2092  Array<int> sfc_order(leaf_elements.Size());
2093  for (int i = 0; i < sfc_order.Size(); i++)
2094  {
2095  sfc_order[leaf_sfc_index[i]] = leaf_elements[i];
2096  }
2097 
2098  for (int i = 0; i < sfc_order.Size(); i++)
2099  {
2100  const Element &el = elements[sfc_order[i]];
2101  for (int j = 0; j < GI[el.Geom()].nv; j++)
2102  {
2103  Node &nd = nodes[el.node[j]];
2104  if (nd.vert_index == -2) { nd.vert_index = NVertices++; }
2105  }
2106  }
2107 
2108  // STEP 4: create the mapping from Mesh vertex index to NCMesh node index
2109 
2111  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2112  {
2113  if (node->HasVertex() && node->vert_index >= 0)
2114  {
2115  MFEM_ASSERT(node->vert_index < vertex_nodeId.Size(), "");
2116  vertex_nodeId[node->vert_index] = node.index();
2117  }
2118  }
2119 
2120  // STEP 5: assign remaining ghost vertices, ignore vertices beyond the
2121  // ghost layer
2122 
2123  NGhostVertices = 0;
2124  for (int i = 0; i < sfc_order.Size(); i++)
2125  {
2126  const Element &el = elements[sfc_order[i]];
2127  for (int j = 0; j < GI[el.Geom()].nv; j++)
2128  {
2129  Node &nd = nodes[el.node[j]];
2130  if (nd.vert_index == -3)
2131  {
2133  }
2134  }
2135  }
2136 
2137 #else // old ordering for debugging/testing only
2138  bool parallel = false;
2139 #ifdef MFEM_USE_MPI
2140  if (dynamic_cast<ParNCMesh*>(this)) { parallel = true; }
2141 #endif
2142 
2143  if (!parallel)
2144  {
2145  NVertices = 0;
2146  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2147  {
2148  if (node->HasVertex()) { node->vert_index = NVertices++; }
2149  }
2150 
2152 
2153  NVertices = 0;
2154  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2155  {
2156  if (node->HasVertex()) { vertex_nodeId[NVertices++] = node.index(); }
2157  }
2158  }
2159  else
2160  {
2161  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2162  {
2163  if (node->HasVertex()) { node->vert_index = -1; }
2164  }
2165 
2166  NVertices = 0;
2167  for (int i = 0; i < leaf_elements.Size(); i++)
2168  {
2169  Element &el = elements[leaf_elements[i]];
2170  if (el.rank == MyRank)
2171  {
2172  for (int j = 0; j < GI[el.Geom()].nv; j++)
2173  {
2174  int &vindex = nodes[el.node[j]].vert_index;
2175  if (vindex < 0) { vindex = NVertices++; }
2176  }
2177  }
2178  }
2179 
2181  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2182  {
2183  if (node->HasVertex() && node->vert_index >= 0)
2184  {
2185  vertex_nodeId[node->vert_index] = node.index();
2186  }
2187  }
2188 
2189  NGhostVertices = 0;
2190  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2191  {
2192  if (node->HasVertex() && node->vert_index < 0)
2193  {
2194  node->vert_index = NVertices + (NGhostVertices++);
2195  }
2196  }
2197  }
2198 #endif
2199 }
2200 
2201 void NCMesh::InitRootState(int root_count)
2202 {
2203  root_state.SetSize(root_count);
2204  root_state = 0;
2205 
2206  char* node_order;
2207  int nch;
2208 
2209  switch (elements[0].Geom()) // TODO: mixed meshes
2210  {
2211  case Geometry::SQUARE:
2212  nch = 4;
2213  node_order = (char*) quad_hilbert_child_order;
2214  break;
2215 
2216  case Geometry::CUBE:
2217  nch = 8;
2218  node_order = (char*) hex_hilbert_child_order;
2219  break;
2220 
2221  default:
2222  return; // do nothing, all states stay zero
2223  }
2224 
2225  int entry_node = -2;
2226 
2227  // process the root element sequence
2228  for (int i = 0; i < root_count; i++)
2229  {
2230  Element &el = elements[i];
2231 
2232  int v_in = FindNodeExt(el, entry_node, false);
2233  if (v_in < 0) { v_in = 0; }
2234 
2235  // determine which nodes are shared with the next element
2236  bool shared[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
2237  if (i+1 < root_count)
2238  {
2239  Element &next = elements[i+1];
2240  for (int j = 0; j < nch; j++)
2241  {
2242  int node = FindNodeExt(el, RetrieveNode(next, j), false);
2243  if (node >= 0) { shared[node] = true; }
2244  }
2245  }
2246 
2247  // select orientation that starts in v_in and exits in shared node
2248  int state = Dim*v_in;
2249  for (int j = 0; j < Dim; j++)
2250  {
2251  if (shared[(int) node_order[nch*(state + j) + nch-1]])
2252  {
2253  state += j;
2254  break;
2255  }
2256  }
2257 
2258  root_state[i] = state;
2259 
2260  entry_node = RetrieveNode(el, node_order[nch*state + nch-1]);
2261  }
2262 }
2263 
2265 {
2266  switch (geom)
2267  {
2268  case Geometry::CUBE: return new mfem::Hexahedron;
2269  case Geometry::PRISM: return new mfem::Wedge;
2270  case Geometry::TETRAHEDRON: return new mfem::Tetrahedron;
2271  case Geometry::SQUARE: return new mfem::Quadrilateral;
2272  case Geometry::TRIANGLE: return new mfem::Triangle;
2273  }
2274  MFEM_ABORT("invalid geometry");
2275  return NULL;
2276 }
2277 
2278 const double* NCMesh::CalcVertexPos(int node) const
2279 {
2280  const Node &nd = nodes[node];
2281  if (nd.p1 == nd.p2) // top-level vertex
2282  {
2283  return &coordinates[3*nd.p1];
2284  }
2285 
2286  TmpVertex &tv = tmp_vertex[node];
2287  if (tv.valid) { return tv.pos; }
2288 
2289  MFEM_VERIFY(tv.visited == false, "cyclic vertex dependencies.");
2290  tv.visited = true;
2291 
2292  const double* pos1 = CalcVertexPos(nd.p1);
2293  const double* pos2 = CalcVertexPos(nd.p2);
2294 
2295  for (int i = 0; i < 3; i++)
2296  {
2297  tv.pos[i] = (pos1[i] + pos2[i]) * 0.5;
2298  }
2299  tv.valid = true;
2300  return tv.pos;
2301 }
2302 
2304 {
2305  mesh.vertices.SetSize(vertex_nodeId.Size());
2306  if (coordinates.Size())
2307  {
2308  // calculate vertex positions from stored top-level vertex coordinates
2309  tmp_vertex = new TmpVertex[nodes.NumIds()];
2310  for (int i = 0; i < mesh.vertices.Size(); i++)
2311  {
2312  mesh.vertices[i].SetCoords(spaceDim, CalcVertexPos(vertex_nodeId[i]));
2313  }
2314  delete [] tmp_vertex;
2315  }
2316  // NOTE: if the mesh is curved ('coordinates' is empty), mesh.vertices are
2317  // left uninitialized here; they will be initialized later by the Mesh from
2318  // Nodes -- here we just make sure mesh.vertices has the correct size.
2319 
2320  mesh.elements.SetSize(0);
2321 
2322  mesh.boundary.SetSize(0);
2323 
2324  // create an mfem::Element for each leaf Element
2325  for (int i = 0; i < NElements; i++)
2326  {
2327  const Element &nc_elem = elements[leaf_elements[i]];
2328 
2329  const int* node = nc_elem.node;
2330  GeomInfo& gi = GI[(int) nc_elem.geom];
2331 
2332  mfem::Element* elem = mesh.NewElement(nc_elem.geom);
2333  mesh.elements.Append(elem);
2334 
2335  elem->SetAttribute(nc_elem.attribute);
2336  for (int j = 0; j < gi.nv; j++)
2337  {
2338  elem->GetVertices()[j] = nodes[node[j]].vert_index;
2339  }
2340 
2341  // create boundary elements
2342  // TODO: use boundary_faces?
2343  for (int k = 0; k < gi.nf; k++)
2344  {
2345  const int* fv = gi.faces[k];
2346  const int nfv = gi.nfv[k];
2347  const Face* face = faces.Find(node[fv[0]], node[fv[1]],
2348  node[fv[2]], node[fv[3]]);
2349  if (face->Boundary())
2350  {
2351  if ((nc_elem.geom == Geometry::CUBE) ||
2352  (nc_elem.geom == Geometry::PRISM && nfv == 4))
2353  {
2354  auto* quad = (Quadrilateral*) mesh.NewElement(Geometry::SQUARE);
2355  quad->SetAttribute(face->attribute);
2356  for (int j = 0; j < 4; j++)
2357  {
2358  quad->GetVertices()[j] = nodes[node[fv[j]]].vert_index;
2359  }
2360  mesh.boundary.Append(quad);
2361  }
2362  else if (nc_elem.geom == Geometry::PRISM ||
2363  nc_elem.geom == Geometry::TETRAHEDRON)
2364  {
2365  MFEM_ASSERT(nfv == 3, "");
2366  auto* tri = (Triangle*) mesh.NewElement(Geometry::TRIANGLE);
2367  tri->SetAttribute(face->attribute);
2368  for (int j = 0; j < 3; j++)
2369  {
2370  tri->GetVertices()[j] = nodes[node[fv[j]]].vert_index;
2371  }
2372  mesh.boundary.Append(tri);
2373  }
2374  else if (nc_elem.geom == Geometry::SQUARE ||
2375  nc_elem.geom == Geometry::TRIANGLE)
2376  {
2377  auto* segment = (Segment*) mesh.NewElement(Geometry::SEGMENT);
2378  segment->SetAttribute(face->attribute);
2379  for (int j = 0; j < 2; j++)
2380  {
2381  segment->GetVertices()[j] = nodes[node[fv[2*j]]].vert_index;
2382  }
2383  mesh.boundary.Append(segment);
2384  }
2385  else
2386  {
2387  MFEM_ASSERT(nc_elem.geom == Geometry::SEGMENT, "");
2388  auto* point = (mfem::Point*) mesh.NewElement(Geometry::POINT);
2389  point->SetAttribute(face->attribute);
2390  point->GetVertices()[0] = nodes[node[fv[0]]].vert_index;
2391  mesh.boundary.Append(point);
2392  }
2393  }
2394  }
2395  }
2396 }
2397 
2399 {
2400  //// PART 1: pull indices of regular edges/faces from the Mesh
2401 
2402  NEdges = mesh->GetNEdges();
2403  NFaces = mesh->GetNumFaces();
2404  if (Dim < 2) { NFaces = 0; }
2405  // clear Node::edge_index and Face::index
2406  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2407  {
2408  if (node->HasEdge()) { node->edge_index = -1; }
2409  }
2410  for (auto face = faces.begin(); face != faces.end(); ++face)
2411  {
2412  face->index = -1;
2413  }
2414 
2415  // get edge enumeration from the Mesh
2416  Table *edge_vertex = mesh->GetEdgeVertexTable();
2417  for (int i = 0; i < edge_vertex->Size(); i++)
2418  {
2419  const int *ev = edge_vertex->GetRow(i);
2420  Node* node = nodes.Find(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
2421 
2422  MFEM_ASSERT(node && node->HasEdge(),
2423  "edge (" << ev[0] << "," << ev[1] << ") not found, "
2424  "node = " << node);
2425 
2426  node->edge_index = i;
2427  }
2428 
2429  // get face enumeration from the Mesh, initialize 'face_geom'
2431  for (int i = 0; i < NFaces; i++)
2432  {
2433  const int* fv = mesh->GetFace(i)->GetVertices();
2434  const int nfv = mesh->GetFace(i)->GetNVertices();
2435 
2436  Face* face;
2437  if (Dim == 3)
2438  {
2439  if (nfv == 4)
2440  {
2442  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2443  vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
2444  }
2445  else
2446  {
2447  MFEM_ASSERT(nfv == 3, "");
2449  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
2450  vertex_nodeId[fv[2]]);
2451  }
2452  }
2453  else
2454  {
2455  MFEM_ASSERT(nfv == 2, "");
2457  int n0 = vertex_nodeId[fv[0]], n1 = vertex_nodeId[fv[1]];
2458  face = faces.Find(n0, n0, n1, n1); // look up degenerate face
2459 
2460 #ifdef MFEM_DEBUG
2461  // (non-ghost) edge and face numbers must match in 2D
2462  const int *ev = edge_vertex->GetRow(i);
2463  MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
2464  (ev[1] == fv[0] && ev[0] == fv[1]), "");
2465 #endif
2466  }
2467 
2468  MFEM_VERIFY(face, "face not found.");
2469  face->index = i;
2470  }
2471 
2472  //// PART 2: assign indices of ghost edges/faces, if any
2473 
2474  // count ghost edges and assign their indices
2475  NGhostEdges = 0;
2476  for (auto node = nodes.begin(); node != nodes.end(); ++node)
2477  {
2478  if (node->HasEdge() && node->edge_index < 0)
2479  {
2480  node->edge_index = NEdges + (NGhostEdges++);
2481  }
2482  }
2483 
2484  // count ghost faces
2485  NGhostFaces = 0;
2486  for (auto face = faces.begin(); face != faces.end(); ++face)
2487  {
2488  if (face->index < 0) { NGhostFaces++; }
2489  }
2490 
2491  if (Dim == 2)
2492  {
2493  // in 2D we have fake faces because of DG
2494  MFEM_ASSERT(NFaces == NEdges, "");
2495  MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
2496  }
2497 
2498  // resize face_geom (default_geom is for slave faces beyond the ghost layer)
2499  Geometry::Type default_geom = Geometry::SQUARE;
2500  face_geom.SetSize(NFaces + NGhostFaces, default_geom);
2501 
2502  // update 'face_geom' for ghost faces, assign ghost face indices
2503  int nghosts = 0;
2504  for (int i = 0; i < NGhostElements; i++)
2505  {
2506  Element &el = elements[leaf_elements[NElements + i]]; // ghost element
2507  GeomInfo &gi = GI[el.Geom()];
2508 
2509  for (int j = 0; j < gi.nf; j++)
2510  {
2511  const int *fv = gi.faces[j];
2512  Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
2513  el.node[fv[2]], el.node[fv[3]]);
2514  MFEM_ASSERT(face, "face not found!");
2515 
2516  if (face->index < 0)
2517  {
2518  face->index = NFaces + (nghosts++);
2519 
2520  // store the face geometry
2521  static const Geometry::Type types[5] =
2522  {
2525  };
2526  face_geom[face->index] = types[gi.nfv[j]];
2527  }
2528  }
2529  }
2530 
2531  // assign valid indices also to faces beyond the ghost layer
2532  for (auto face = faces.begin(); face != faces.end(); ++face)
2533  {
2534  if (face->index < 0) { face->index = NFaces + (nghosts++); }
2535  }
2536  MFEM_ASSERT(nghosts == NGhostFaces, "");
2537 }
2538 
2539 
2540 //// Face/edge lists ///////////////////////////////////////////////////////////
2541 
2542 int NCMesh::QuadFaceSplitType(int v1, int v2, int v3, int v4,
2543  int mid[5]) const
2544 {
2545  MFEM_ASSERT(Dim >= 3, "");
2546 
2547  // find edge nodes
2548  int e1 = FindMidEdgeNode(v1, v2);
2549  int e2 = FindMidEdgeNode(v2, v3);
2550  int e3 = (e1 >= 0 && nodes[e1].HasVertex()) ? FindMidEdgeNode(v3, v4) : -1;
2551  int e4 = (e2 >= 0 && nodes[e2].HasVertex()) ? FindMidEdgeNode(v4, v1) : -1;
2552 
2553  // optional: return the mid-edge nodes if requested
2554  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
2555 
2556  // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
2557  int midf1 = -1, midf2 = -1;
2558  if (e1 >= 0 && e3 >= 0) { midf1 = FindMidEdgeNode(e1, e3); }
2559  if (e2 >= 0 && e4 >= 0) { midf2 = FindMidEdgeNode(e2, e4); }
2560 
2561  // get proper node if shadow node exists
2562  if (midf1 >= 0 && midf1 == midf2)
2563  {
2564  const Node &nd = nodes[midf1];
2565  if (nd.p1 != e1 && nd.p2 != e1) { midf1 = -1; }
2566  if (nd.p1 != e2 && nd.p2 != e2) { midf2 = -1; }
2567  }
2568 
2569  // only one way to access the mid-face node must always exist
2570  MFEM_ASSERT(!(midf1 >= 0 && midf2 >= 0), "incorrectly split face!");
2571 
2572  if (midf1 < 0 && midf2 < 0) // face not split
2573  {
2574  if (mid) { mid[4] = -1; }
2575  return 0;
2576  }
2577  else if (midf1 >= 0) // face split "vertically"
2578  {
2579  if (mid) { mid[4] = midf1; }
2580  return 1;
2581  }
2582  else // face split "horizontally"
2583  {
2584  if (mid) { mid[4] = midf2; }
2585  return 2;
2586  }
2587 }
2588 
2589 bool NCMesh::TriFaceSplit(int v1, int v2, int v3, int mid[3]) const
2590 {
2591  int e1 = nodes.FindId(v1, v2);
2592  if (e1 < 0 || !nodes[e1].HasVertex()) { return false; }
2593 
2594  int e2 = nodes.FindId(v2, v3);
2595  if (e2 < 0 || !nodes[e2].HasVertex()) { return false; }
2596 
2597  int e3 = nodes.FindId(v3, v1);
2598  if (e3 < 0 || !nodes[e3].HasVertex()) { return false; }
2599 
2600  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3; }
2601 
2602  // NOTE: face (v1, v2, v3) still needs to be checked
2603  return true;
2604 }
2605 
2606 int NCMesh::find_node(const Element &el, int node)
2607 {
2608  for (int i = 0; i < 8; i++)
2609  {
2610  if (el.node[i] == node) { return i; }
2611  }
2612  MFEM_ABORT("Node not found.");
2613  return -1;
2614 }
2615 
2616 int NCMesh::FindNodeExt(const Element &el, int node, bool abort)
2617 {
2618  for (int i = 0; i < GI[el.Geom()].nv; i++)
2619  {
2620  if (RetrieveNode(el, i) == node) { return i; }
2621  }
2622  if (abort) { MFEM_ABORT("Node not found."); }
2623  return -1;
2624 }
2625 
2626 int NCMesh::find_element_edge(const Element &el, int vn0, int vn1, bool abort)
2627 {
2628  MFEM_ASSERT(!el.ref_type, "");
2629 
2630  GeomInfo &gi = GI[el.Geom()];
2631  for (int i = 0; i < gi.ne; i++)
2632  {
2633  const int* ev = gi.edges[i];
2634  int n0 = el.node[ev[0]];
2635  int n1 = el.node[ev[1]];
2636  if ((n0 == vn0 && n1 == vn1) ||
2637  (n0 == vn1 && n1 == vn0)) { return i; }
2638  }
2639 
2640  if (abort) { MFEM_ABORT("Edge (" << vn0 << ", " << vn1 << ") not found"); }
2641  return -1;
2642 }
2643 
2644 int NCMesh::find_local_face(int geom, int a, int b, int c)
2645 {
2646  GeomInfo &gi = GI[geom];
2647  for (int i = 0; i < gi.nf; i++)
2648  {
2649  const int* fv = gi.faces[i];
2650  if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
2651  (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
2652  (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
2653  {
2654  return i;
2655  }
2656  }
2657  MFEM_ABORT("Face not found.");
2658  return -1;
2659 }
2660 
2661 
2662 /// Hash function for a PointMatrix, used in MatrixMap::map.
2663 struct PointMatrixHash
2664 {
2665  std::size_t operator()(const NCMesh::PointMatrix &pm) const
2666  {
2667  MFEM_ASSERT(sizeof(double) == sizeof(std::uint64_t), "");
2668 
2669  // This is a variation on "Hashing an array of floats" from here:
2670  // https://cs.stackexchange.com/questions/37952
2671  std::uint64_t hash = 0xf9ca9ba106acbba9; // random initial value
2672  for (int i = 0; i < pm.np; i++)
2673  {
2674  for (int j = 0; j < pm.points[i].dim; j++)
2675  {
2676  // mix the doubles by adding their binary representations
2677  // many times over (note: 31 is 11111 in binary)
2678  double coord = pm.points[i].coord[j];
2679  hash = 31*hash + *((std::uint64_t*) &coord);
2680  }
2681  }
2682  return hash; // return the lowest bits of the huge sum
2683  }
2684 };
2685 
2686 /** Helper container to keep track of point matrices encountered during
2687  * face/edge traversal and to assign unique indices to them.
2688  */
2689 struct MatrixMap
2690 {
2691  int GetIndex(const NCMesh::PointMatrix &pm)
2692  {
2693  int &index = map[pm];
2694  if (!index) { index = map.size(); }
2695  return index - 1;
2696  }
2697 
2698  void ExportMatrices(Array<DenseMatrix*> &point_matrices) const
2699  {
2700  point_matrices.SetSize(map.size());
2701  for (const auto &pair : map)
2702  {
2703  DenseMatrix* mat = new DenseMatrix();
2704  pair.first.GetMatrix(*mat);
2705  point_matrices[pair.second - 1] = mat;
2706  }
2707  }
2708 
2709  void DumpBucketSizes() const
2710  {
2711  for (unsigned i = 0; i < map.bucket_count(); i++)
2712  {
2713  mfem::out << map.bucket_size(i) << " ";
2714  }
2715  }
2716 
2717 private:
2718  std::unordered_map<NCMesh::PointMatrix, int, PointMatrixHash> map;
2719 };
2720 
2721 
2722 int NCMesh::ReorderFacePointMat(int v0, int v1, int v2, int v3,
2723  int elem, const PointMatrix &pm,
2724  PointMatrix &reordered) const
2725 {
2726  const Element &el = elements[elem];
2727  int master[4] =
2728  {
2729  find_node(el, v0), find_node(el, v1), find_node(el, v2),
2730  (v3 >= 0) ? find_node(el, v3) : -1
2731  };
2732  int nfv = (v3 >= 0) ? 4 : 3;
2733 
2734  int local = find_local_face(el.Geom(), master[0], master[1], master[2]);
2735  const int* fv = GI[el.Geom()].faces[local];
2736 
2737  reordered.np = pm.np;
2738  for (int i = 0, j; i < nfv; i++)
2739  {
2740  for (j = 0; j < nfv; j++)
2741  {
2742  if (fv[i] == master[j])
2743  {
2744  reordered.points[i] = pm.points[j];
2745  break;
2746  }
2747  }
2748  MFEM_ASSERT(j != nfv, "node not found.");
2749  }
2750  return local;
2751 }
2752 
2753 void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
2754  const PointMatrix& pm, int level,
2755  Face* eface[4], MatrixMap &matrix_map)
2756 {
2757  if (level > 0)
2758  {
2759  // check if we made it to a face that is not split further
2760  Face* fa = faces.Find(vn0, vn1, vn2, vn3);
2761  if (fa)
2762  {
2763  // we have a slave face, add it to the list
2764  int elem = fa->GetSingleElement();
2765  face_list.slaves.Append(
2766  Slave(fa->index, elem, -1, Geometry::SQUARE));
2767  Slave &sl = face_list.slaves.Last();
2768 
2769  // reorder the point matrix according to slave face orientation
2770  PointMatrix pm_r;
2771  sl.local = ReorderFacePointMat(vn0, vn1, vn2, vn3, elem, pm, pm_r);
2772  sl.matrix = matrix_map.GetIndex(pm_r);
2773 
2774  eface[0] = eface[2] = fa;
2775  eface[1] = eface[3] = fa;
2776 
2777  return;
2778  }
2779  }
2780 
2781  // we need to recurse deeper
2782  int mid[5];
2783  int split = QuadFaceSplitType(vn0, vn1, vn2, vn3, mid);
2784 
2785  Face *ef[2][4];
2786  if (split == 1) // "X" split face
2787  {
2788  Point pmid0(pm(0), pm(1)), pmid2(pm(2), pm(3));
2789 
2790  TraverseQuadFace(vn0, mid[0], mid[2], vn3,
2791  PointMatrix(pm(0), pmid0, pmid2, pm(3)),
2792  level+1, ef[0], matrix_map);
2793 
2794  TraverseQuadFace(mid[0], vn1, vn2, mid[2],
2795  PointMatrix(pmid0, pm(1), pm(2), pmid2),
2796  level+1, ef[1], matrix_map);
2797 
2798  eface[1] = ef[1][1];
2799  eface[3] = ef[0][3];
2800  eface[0] = eface[2] = NULL;
2801  }
2802  else if (split == 2) // "Y" split face
2803  {
2804  Point pmid1(pm(1), pm(2)), pmid3(pm(3), pm(0));
2805 
2806  TraverseQuadFace(vn0, vn1, mid[1], mid[3],
2807  PointMatrix(pm(0), pm(1), pmid1, pmid3),
2808  level+1, ef[0], matrix_map);
2809 
2810  TraverseQuadFace(mid[3], mid[1], vn2, vn3,
2811  PointMatrix(pmid3, pmid1, pm(2), pm(3)),
2812  level+1, ef[1], matrix_map);
2813 
2814  eface[0] = ef[0][0];
2815  eface[2] = ef[1][2];
2816  eface[1] = eface[3] = NULL;
2817  }
2818 
2819  // check for a prism edge constrained by the master face
2820  if (HavePrisms() && mid[4] >= 0)
2821  {
2822  Node& enode = nodes[mid[4]];
2823  if (enode.HasEdge())
2824  {
2825  // process the edge only if it's not shared by slave faces
2826  // within this master face (i.e. the edge is "hidden")
2827  const int fi[3][2] = {{0, 0}, {1, 3}, {2, 0}};
2828  if (!ef[0][fi[split][0]] && !ef[1][fi[split][1]])
2829  {
2830  MFEM_ASSERT(enode.edge_refc == 1, "");
2831 
2832  MeshId buf[4];
2833  Array<MeshId> eid(buf, 4);
2834 
2835  (split == 1) ? FindEdgeElements(mid[0], vn1, vn2, mid[2], eid)
2836  /* */ : FindEdgeElements(mid[3], vn0, vn1, mid[1], eid);
2837 
2838  MFEM_ASSERT(eid.Size() > 0, "edge prism not found");
2839  MFEM_ASSERT(eid.Size() < 2, "non-unique edge prism");
2840 
2841  // create a slave face record with a degenerate point matrix
2842  face_list.slaves.Append(
2843  Slave(-1 - enode.edge_index,
2844  eid[0].element, eid[0].local, Geometry::SQUARE));
2845  Slave &sl = face_list.slaves.Last();
2846 
2847  if (split == 1)
2848  {
2849  Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
2850  int v1 = nodes[mid[0]].vert_index;
2851  int v2 = nodes[mid[2]].vert_index;
2852  sl.matrix =
2853  matrix_map.GetIndex(
2854  (v1 < v2) ? PointMatrix(mid0, mid2, mid2, mid0) :
2855  /* */ PointMatrix(mid2, mid0, mid0, mid2));
2856  }
2857  else
2858  {
2859  Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
2860  int v1 = nodes[mid[1]].vert_index;
2861  int v2 = nodes[mid[3]].vert_index;
2862  sl.matrix =
2863  matrix_map.GetIndex(
2864  (v1 < v2) ? PointMatrix(mid1, mid3, mid3, mid1) :
2865  /* */ PointMatrix(mid3, mid1, mid1, mid3));
2866  }
2867  }
2868  }
2869  }
2870 }
2871 
2872 void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
2873  MatrixMap &matrix_map)
2874 {
2875  int mid = nodes.FindId(vn0, vn1);
2876  if (mid < 0) { return; }
2877 
2878  const Node &nd = nodes[mid];
2879  if (nd.HasEdge())
2880  {
2881  // check if the edge is already a master in 'edge_list'
2882  int type;
2883  const MeshId &eid = edge_list.LookUp(nd.edge_index, &type);
2884  if (type == 1)
2885  {
2886  // in this case we need to add an edge-face constraint, because the
2887  // master edge is really a (face-)slave itself
2888 
2889  face_list.slaves.Append(
2890  Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
2891 
2892  int v0index = nodes[vn0].vert_index;
2893  int v1index = nodes[vn1].vert_index;
2894 
2895  face_list.slaves.Last().matrix =
2896  matrix_map.GetIndex((v0index < v1index) ? PointMatrix(p0, p1, p0)
2897  /* */ : PointMatrix(p1, p0, p1));
2898 
2899  return; // no need to continue deeper
2900  }
2901  }
2902 
2903  // recurse deeper
2904  Point pmid(p0, p1);
2905  TraverseTetEdge(vn0, mid, p0, pmid, matrix_map);
2906  TraverseTetEdge(mid, vn1, pmid, p1, matrix_map);
2907 }
2908 
2909 bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
2910  const PointMatrix& pm, int level,
2911  MatrixMap &matrix_map)
2912 {
2913  if (level > 0)
2914  {
2915  // check if we made it to a face that is not split further
2916  Face* fa = faces.Find(vn0, vn1, vn2);
2917  if (fa)
2918  {
2919  // we have a slave face, add it to the list
2920  int elem = fa->GetSingleElement();
2921  face_list.slaves.Append(
2922  Slave(fa->index, elem, -1, Geometry::TRIANGLE));
2923  Slave &sl = face_list.slaves.Last();
2924 
2925  // reorder the point matrix according to slave face orientation
2926  PointMatrix pm_r;
2927  sl.local = ReorderFacePointMat(vn0, vn1, vn2, -1, elem, pm, pm_r);
2928  sl.matrix = matrix_map.GetIndex(pm_r);
2929 
2930  return true;
2931  }
2932  }
2933 
2934  int mid[3];
2935  if (TriFaceSplit(vn0, vn1, vn2, mid))
2936  {
2937  Point pmid0(pm(0), pm(1)), pmid1(pm(1), pm(2)), pmid2(pm(2), pm(0));
2938  bool b[4];
2939 
2940  b[0] = TraverseTriFace(vn0, mid[0], mid[2],
2941  PointMatrix(pm(0), pmid0, pmid2),
2942  level+1, matrix_map);
2943 
2944  b[1] = TraverseTriFace(mid[0], vn1, mid[1],
2945  PointMatrix(pmid0, pm(1), pmid1),
2946  level+1, matrix_map);
2947 
2948  b[2] = TraverseTriFace(mid[2], mid[1], vn2,
2949  PointMatrix(pmid2, pmid1, pm(2)),
2950  level+1, matrix_map);
2951 
2952  b[3] = TraverseTriFace(mid[1], mid[2], mid[0],
2953  PointMatrix(pmid1, pmid2, pmid0),
2954  level+1, matrix_map);
2955 
2956  // traverse possible tet edges constrained by the master face
2957  if (HaveTets() && !b[3])
2958  {
2959  if (!b[1]) { TraverseTetEdge(mid[0],mid[1], pmid0,pmid1, matrix_map); }
2960  if (!b[2]) { TraverseTetEdge(mid[1],mid[2], pmid1,pmid2, matrix_map); }
2961  if (!b[0]) { TraverseTetEdge(mid[2],mid[0], pmid2,pmid0, matrix_map); }
2962  }
2963  }
2964 
2965  return false;
2966 }
2967 
2969 {
2970  face_list.Clear();
2971  if (Dim < 3) { return; }
2972 
2973  if (HaveTets()) { GetEdgeList(); } // needed by TraverseTetEdge()
2974 
2976 
2977  Array<char> processed_faces(faces.NumIds());
2978  processed_faces = 0;
2979 
2980  MatrixMap matrix_maps[Geometry::NumGeom];
2981 
2982  // visit faces of leaf elements
2983  for (int i = 0; i < leaf_elements.Size(); i++)
2984  {
2985  int elem = leaf_elements[i];
2986  Element &el = elements[elem];
2987  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
2988 
2989  GeomInfo& gi = GI[el.Geom()];
2990  for (int j = 0; j < gi.nf; j++)
2991  {
2992  // get nodes for this face
2993  int node[4];
2994  for (int k = 0; k < 4; k++)
2995  {
2996  node[k] = el.node[gi.faces[j][k]];
2997  }
2998 
2999  int face = faces.FindId(node[0], node[1], node[2], node[3]);
3000  MFEM_ASSERT(face >= 0, "face not found!");
3001 
3002  // tell ParNCMesh about the face
3003  ElementSharesFace(elem, j, face);
3004 
3005  // have we already processed this face? skip if yes
3006  if (processed_faces[face]) { continue; }
3007  processed_faces[face] = 1;
3008 
3009  int fgeom = (node[3] >= 0) ? Geometry::SQUARE : Geometry::TRIANGLE;
3010 
3011  Face &fa = faces[face];
3012  if (fa.elem[0] >= 0 && fa.elem[1] >= 0)
3013  {
3014  // this is a conforming face, add it to the list
3015  face_list.conforming.Append(MeshId(fa.index, elem, j, fgeom));
3016  }
3017  else
3018  {
3019  // this is either a master face or a slave face, but we can't
3020  // tell until we traverse the face refinement 'tree'...
3021  int sb = face_list.slaves.Size();
3022  if (fgeom == Geometry::SQUARE)
3023  {
3024  Face* dummy[4];
3025  TraverseQuadFace(node[0], node[1], node[2], node[3],
3026  pm_quad_identity, 0, dummy, matrix_maps[fgeom]);
3027  }
3028  else
3029  {
3030  TraverseTriFace(node[0], node[1], node[2],
3031  pm_tri_identity, 0, matrix_maps[fgeom]);
3032  }
3033 
3034  int se = face_list.slaves.Size();
3035  if (sb < se)
3036  {
3037  // found slaves, so this is a master face; add it to the list
3038  face_list.masters.Append(
3039  Master(fa.index, elem, j, fgeom, sb, se));
3040 
3041  // also, set the master index for the slaves
3042  for (int ii = sb; ii < se; ii++)
3043  {
3044  face_list.slaves[ii].master = fa.index;
3045  }
3046  }
3047  }
3048 
3049  if (fa.Boundary()) { boundary_faces.Append(face); }
3050  }
3051  }
3052 
3053  // export unique point matrices
3054  for (int i = 0; i < Geometry::NumGeom; i++)
3055  {
3056  matrix_maps[i].ExportMatrices(face_list.point_matrices[i]);
3057  }
3058 }
3059 
3060 void NCMesh::TraverseEdge(int vn0, int vn1, double t0, double t1, int flags,
3061  int level, MatrixMap &matrix_map)
3062 {
3063  int mid = nodes.FindId(vn0, vn1);
3064  if (mid < 0) { return; }
3065 
3066  Node &nd = nodes[mid];
3067  if (nd.HasEdge() && level > 0)
3068  {
3069  // we have a slave edge, add it to the list
3070  edge_list.slaves.Append(Slave(nd.edge_index, -1, -1, Geometry::SEGMENT));
3071 
3072  Slave &sl = edge_list.slaves.Last();
3073  sl.matrix = matrix_map.GetIndex(PointMatrix(Point(t0), Point(t1)));
3074 
3075  // handle slave edge orientation
3076  sl.edge_flags = flags;
3077  int v0index = nodes[vn0].vert_index;
3078  int v1index = nodes[vn1].vert_index;
3079  if (v0index > v1index) { sl.edge_flags |= 2; }
3080  }
3081 
3082  // recurse deeper
3083  double tmid = (t0 + t1) / 2;
3084  TraverseEdge(vn0, mid, t0, tmid, flags, level+1, matrix_map);
3085  TraverseEdge(mid, vn1, tmid, t1, flags, level+1, matrix_map);
3086 }
3087 
3089 {
3090  edge_list.Clear();
3091  if (Dim < 3) { boundary_faces.SetSize(0); }
3092 
3093  Array<char> processed_edges(nodes.NumIds());
3094  processed_edges = 0;
3095 
3096  Array<int> edge_element(nodes.NumIds());
3097  Array<signed char> edge_local(nodes.NumIds());
3098  edge_local = -1;
3099 
3100  MatrixMap matrix_map;
3101 
3102  // visit edges of leaf elements
3103  for (int i = 0; i < leaf_elements.Size(); i++)
3104  {
3105  int elem = leaf_elements[i];
3106  Element &el = elements[elem];
3107  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3108 
3109  GeomInfo& gi = GI[el.Geom()];
3110  for (int j = 0; j < gi.ne; j++)
3111  {
3112  // get nodes for this edge
3113  const int* ev = gi.edges[j];
3114  int node[2] = { el.node[ev[0]], el.node[ev[1]] };
3115 
3116  int enode = nodes.FindId(node[0], node[1]);
3117  MFEM_ASSERT(enode >= 0, "edge node not found!");
3118 
3119  Node &nd = nodes[enode];
3120  MFEM_ASSERT(nd.HasEdge(), "edge not found!");
3121 
3122  // tell ParNCMesh about the edge
3123  ElementSharesEdge(elem, j, enode);
3124 
3125  // (2D only, store boundary faces)
3126  if (Dim <= 2)
3127  {
3128  int face = faces.FindId(node[0], node[0], node[1], node[1]);
3129  MFEM_ASSERT(face >= 0, "face not found!");
3130  if (faces[face].Boundary()) { boundary_faces.Append(face); }
3131  }
3132 
3133  // store element/local for later
3134  edge_element[nd.edge_index] = elem;
3135  edge_local[nd.edge_index] = j;
3136 
3137  // skip slave edges here, they will be reached from their masters
3138  if (GetEdgeMaster(enode) >= 0) { continue; }
3139 
3140  // have we already processed this edge? skip if yes
3141  if (processed_edges[enode]) { continue; }
3142  processed_edges[enode] = 1;
3143 
3144  // prepare edge interval for slave traversal, handle orientation
3145  double t0 = 0.0, t1 = 1.0;
3146  int v0index = nodes[node[0]].vert_index;
3147  int v1index = nodes[node[1]].vert_index;
3148  int flags = (v0index > v1index) ? 1 : 0;
3149 
3150  // try traversing the edge to find slave edges
3151  int sb = edge_list.slaves.Size();
3152  TraverseEdge(node[0], node[1], t0, t1, flags, 0, matrix_map);
3153 
3154  int se = edge_list.slaves.Size();
3155  if (sb < se)
3156  {
3157  // found slaves, this is a master face; add it to the list
3158  edge_list.masters.Append(
3159  Master(nd.edge_index, elem, j, Geometry::SEGMENT, sb, se));
3160 
3161  // also, set the master index for the slaves
3162  for (int ii = sb; ii < se; ii++)
3163  {
3164  edge_list.slaves[ii].master = nd.edge_index;
3165  }
3166  }
3167  else
3168  {
3169  // no slaves, this is a conforming edge
3170  edge_list.conforming.Append(MeshId(nd.edge_index, elem, j));
3171  }
3172  }
3173  }
3174 
3175  // fix up slave edge element/local
3176  for (int i = 0; i < edge_list.slaves.Size(); i++)
3177  {
3178  Slave &sl = edge_list.slaves[i];
3179  int local = edge_local[sl.index];
3180  if (local >= 0)
3181  {
3182  sl.local = local;
3183  sl.element = edge_element[sl.index];
3184  }
3185  }
3186 
3187  // export unique point matrices
3188  matrix_map.ExportMatrices(edge_list.point_matrices[Geometry::SEGMENT]);
3189 }
3190 
3192 {
3193  int total = NVertices + NGhostVertices;
3194 
3195  vertex_list.Clear();
3196  vertex_list.conforming.Reserve(total);
3197 
3198  Array<char> processed_vertices(total);
3199  processed_vertices = 0;
3200 
3201  // analogously to above, visit vertices of leaf elements
3202  for (int i = 0; i < leaf_elements.Size(); i++)
3203  {
3204  int elem = leaf_elements[i];
3205  Element &el = elements[elem];
3206 
3207  for (int j = 0; j < GI[el.Geom()].nv; j++)
3208  {
3209  int node = el.node[j];
3210  Node &nd = nodes[node];
3211 
3212  int index = nd.vert_index;
3213  if (index >= 0)
3214  {
3215  ElementSharesVertex(elem, j, node);
3216 
3217  if (processed_vertices[index]) { continue; }
3218  processed_vertices[index] = 1;
3219 
3220  vertex_list.conforming.Append(MeshId(index, elem, j));
3221  }
3222  }
3223  }
3224 }
3225 
3227  DenseMatrix &oriented_matrix) const
3228 {
3229  oriented_matrix = *(point_matrices[slave.Geom()][slave.matrix]);
3230 
3231  if (slave.edge_flags)
3232  {
3233  MFEM_ASSERT(oriented_matrix.Height() == 1 &&
3234  oriented_matrix.Width() == 2, "not an edge point matrix");
3235 
3236  if (slave.edge_flags & 1) // master inverted
3237  {
3238  oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
3239  oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
3240  }
3241  if (slave.edge_flags & 2) // slave inverted
3242  {
3243  std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
3244  }
3245  }
3246 }
3247 
3249 {
3250  conforming.DeleteAll();
3251  masters.DeleteAll();
3252  slaves.DeleteAll();
3253 
3254  for (int i = 0; i < Geometry::NumGeom; i++)
3255  {
3256  for (int j = 0; j < point_matrices[i].Size(); j++)
3257  {
3258  delete point_matrices[i][j];
3259  }
3260  point_matrices[i].DeleteAll();
3261  }
3262 
3263  inv_index.DeleteAll();
3264 }
3265 
3267 {
3268  return conforming.Size() + masters.Size() + slaves.Size();
3269 }
3270 
3271 const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
3272 {
3273  if (!inv_index.Size())
3274  {
3275  int max_index = -1;
3276  for (int i = 0; i < conforming.Size(); i++)
3277  {
3278  max_index = std::max(conforming[i].index, max_index);
3279  }
3280  for (int i = 0; i < masters.Size(); i++)
3281  {
3282  max_index = std::max(masters[i].index, max_index);
3283  }
3284  for (int i = 0; i < slaves.Size(); i++)
3285  {
3286  if (slaves[i].index < 0) { continue; }
3287  max_index = std::max(slaves[i].index, max_index);
3288  }
3289 
3290  inv_index.SetSize(max_index + 1);
3291  inv_index = -1;
3292 
3293  for (int i = 0; i < conforming.Size(); i++)
3294  {
3295  inv_index[conforming[i].index] = (i << 2);
3296  }
3297  for (int i = 0; i < masters.Size(); i++)
3298  {
3299  inv_index[masters[i].index] = (i << 2) + 1;
3300  }
3301  for (int i = 0; i < slaves.Size(); i++)
3302  {
3303  if (slaves[i].index < 0) { continue; }
3304  inv_index[slaves[i].index] = (i << 2) + 2;
3305  }
3306  }
3307 
3308  MFEM_ASSERT(index >= 0 && index < inv_index.Size(), "");
3309  int key = inv_index[index];
3310 
3311  if (!type)
3312  {
3313  MFEM_VERIFY(key >= 0, "entity not found.");
3314  }
3315  else // return entity type if requested, don't abort when not found
3316  {
3317  *type = (key >= 0) ? (key & 0x3) : -1;
3318 
3319  static MeshId invalid;
3320  if (*type < 0) { return invalid; } // not found
3321  }
3322 
3323  // return found entity MeshId
3324  switch (key & 0x3)
3325  {
3326  case 0: return conforming[key >> 2];
3327  case 1: return masters[key >> 2];
3328  case 2: return slaves[key >> 2];
3329  default: MFEM_ABORT("internal error"); return conforming[0];
3330  }
3331 }
3332 
3333 
3334 //// Neighbors /////////////////////////////////////////////////////////////////
3335 
3336 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
3337 {
3338  int mid = nodes.FindId(v0, v1);
3339  if (mid >= 0 && nodes[mid].HasVertex())
3340  {
3341  indices.Append(mid);
3342 
3343  CollectEdgeVertices(v0, mid, indices);
3344  CollectEdgeVertices(mid, v1, indices);
3345  }
3346 }
3347 
3348 void NCMesh::CollectTriFaceVertices(int v0, int v1, int v2, Array<int> &indices)
3349 {
3350  int mid[3];
3351  if (TriFaceSplit(v0, v1, v2, mid))
3352  {
3353  for (int i = 0; i < 3; i++)
3354  {
3355  indices.Append(mid[i]);
3356  }
3357 
3358  CollectTriFaceVertices(v0, mid[0], mid[2], indices);
3359  CollectTriFaceVertices(mid[0], v1, mid[1], indices);
3360  CollectTriFaceVertices(mid[2], mid[1], v2, indices);
3361  CollectTriFaceVertices(mid[0], mid[1], mid[2], indices);
3362 
3363  if (HaveTets()) // possible edge-face contact
3364  {
3365  CollectEdgeVertices(mid[0], mid[1], indices);
3366  CollectEdgeVertices(mid[1], mid[2], indices);
3367  CollectEdgeVertices(mid[2], mid[0], indices);
3368  }
3369  }
3370 }
3371 
3372 void NCMesh::CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
3373  Array<int> &indices)
3374 {
3375  int mid[5];
3376  switch (QuadFaceSplitType(v0, v1, v2, v3, mid))
3377  {
3378  case 1:
3379  indices.Append(mid[0]);
3380  indices.Append(mid[2]);
3381 
3382  CollectQuadFaceVertices(v0, mid[0], mid[2], v3, indices);
3383  CollectQuadFaceVertices(mid[0], v1, v2, mid[2], indices);
3384 
3385  if (HavePrisms()) // possible edge-face contact
3386  {
3387  CollectEdgeVertices(mid[0], mid[2], indices);
3388  }
3389  break;
3390 
3391  case 2:
3392  indices.Append(mid[1]);
3393  indices.Append(mid[3]);
3394 
3395  CollectQuadFaceVertices(v0, v1, mid[1], mid[3], indices);
3396  CollectQuadFaceVertices(mid[3], mid[1], v2, v3, indices);
3397 
3398  if (HavePrisms()) // possible edge-face contact
3399  {
3400  CollectEdgeVertices(mid[1], mid[3], indices);
3401  }
3402  break;
3403  }
3404 }
3405 
3407 {
3408  int nrows = leaf_elements.Size();
3409  int* I = Memory<int>(nrows + 1);
3410  int** JJ = new int*[nrows];
3411 
3412  Array<int> indices;
3413  indices.Reserve(128);
3414 
3415  // collect vertices coinciding with each element, including hanging vertices
3416  for (int i = 0; i < leaf_elements.Size(); i++)
3417  {
3418  int elem = leaf_elements[i];
3419  Element &el = elements[elem];
3420  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
3421 
3422  GeomInfo& gi = GI[el.Geom()];
3423  int* node = el.node;
3424 
3425  indices.SetSize(0);
3426  for (int j = 0; j < gi.ne; j++)
3427  {
3428  const int* ev = gi.edges[j];
3429  CollectEdgeVertices(node[ev[0]], node[ev[1]], indices);
3430  }
3431 
3432  if (Dim >= 3)
3433  {
3434  for (int j = 0; j < gi.nf; j++)
3435  {
3436  const int* fv = gi.faces[j];
3437  if (gi.nfv[j] == 4)
3438  {
3439  CollectQuadFaceVertices(node[fv[0]], node[fv[1]],
3440  node[fv[2]], node[fv[3]], indices);
3441  }
3442  else
3443  {
3444  CollectTriFaceVertices(node[fv[0]], node[fv[1]], node[fv[2]],
3445  indices);
3446  }
3447  }
3448  }
3449 
3450  // temporarily store one row of the table
3451  indices.Sort();
3452  indices.Unique();
3453  int size = indices.Size();
3454  I[i] = size;
3455  JJ[i] = new int[size];
3456  std::memcpy(JJ[i], indices.GetData(), size * sizeof(int));
3457  }
3458 
3459  // finalize the I array of the table
3460  int nnz = 0;
3461  for (int i = 0; i < nrows; i++)
3462  {
3463  int cnt = I[i];
3464  I[i] = nnz;
3465  nnz += cnt;
3466  }
3467  I[nrows] = nnz;
3468 
3469  // copy the temporarily stored rows into one J array
3470  int *J = Memory<int>(nnz);
3471  nnz = 0;
3472  for (int i = 0; i < nrows; i++)
3473  {
3474  int cnt = I[i+1] - I[i];
3475  std::memcpy(J+nnz, JJ[i], cnt * sizeof(int));
3476  delete [] JJ[i];
3477  nnz += cnt;
3478  }
3479 
3480  element_vertex.SetIJ(I, J, nrows);
3481 
3482  delete [] JJ;
3483 }
3484 
3485 
3487  Array<int> *neighbors,
3488  Array<char> *neighbor_set)
3489 {
3490  // If A is the element-to-vertex table (see 'element_vertex') listing all
3491  // vertices touching each element, including hanging vertices, then A*A^T is
3492  // the element-to-neighbor table. Multiplying the element set with A*A^T
3493  // gives the neighbor set. To save memory, this function only computes the
3494  // action of A*A^T, the product itself is not stored anywhere.
3495 
3496  // Optimization: the 'element_vertex' table does not store the obvious
3497  // corner nodes in it. The table is therefore empty for conforming meshes.
3498 
3500 
3501  int nleaves = leaf_elements.Size();
3502  MFEM_VERIFY(elem_set.Size() == nleaves, "");
3503  MFEM_ASSERT(element_vertex.Size() == nleaves, "");
3504 
3505  // step 1: vertices = A^T * elem_set, i.e, find all vertices touching the
3506  // element set
3507 
3508  Array<char> vmark(nodes.NumIds());
3509  vmark = 0;
3510 
3511  for (int i = 0; i < nleaves; i++)
3512  {
3513  if (elem_set[i])
3514  {
3515  int *v = element_vertex.GetRow(i);
3516  int nv = element_vertex.RowSize(i);
3517  for (int j = 0; j < nv; j++)
3518  {
3519  vmark[v[j]] = 1;
3520  }
3521 
3522  Element &el = elements[leaf_elements[i]];
3523  nv = GI[el.Geom()].nv;
3524  for (int j = 0; j < nv; j++)
3525  {
3526  vmark[el.node[j]] = 1;
3527  }
3528  }
3529  }
3530 
3531  // step 2: neighbors = A * vertices, i.e., find all elements coinciding with
3532  // vertices from step 1; NOTE: in the result we don't include elements from
3533  // the original set
3534 
3535  if (neighbor_set)
3536  {
3537  neighbor_set->SetSize(nleaves);
3538  *neighbor_set = 0;
3539  }
3540 
3541  for (int i = 0; i < nleaves; i++)
3542  {
3543  if (!elem_set[i])
3544  {
3545  bool hit = false;
3546 
3547  int *v = element_vertex.GetRow(i);
3548  int nv = element_vertex.RowSize(i);
3549  for (int j = 0; j < nv; j++)
3550  {
3551  if (vmark[v[j]]) { hit = true; break; }
3552  }
3553 
3554  if (!hit)
3555  {
3556  Element &el = elements[leaf_elements[i]];
3557  nv = GI[el.Geom()].nv;
3558  for (int j = 0; j < nv; j++)
3559  {
3560  if (vmark[el.node[j]]) { hit = true; break; }
3561  }
3562  }
3563 
3564  if (hit)
3565  {
3566  if (neighbors) { neighbors->Append(leaf_elements[i]); }
3567  if (neighbor_set) { (*neighbor_set)[i] = 1; }
3568  }
3569  }
3570  }
3571 }
3572 
3573 static bool sorted_lists_intersect(const int* a, const int* b, int na, int nb)
3574 {
3575  if (!na || !nb) { return false; }
3576  int a_last = a[na-1], b_last = b[nb-1];
3577  if (*b < *a) { goto l2; } // woo-hoo! I always wanted to use a goto! :)
3578 l1:
3579  if (a_last < *b) { return false; }
3580  while (*a < *b) { a++; }
3581  if (*a == *b) { return true; }
3582 l2:
3583  if (b_last < *a) { return false; }
3584  while (*b < *a) { b++; }
3585  if (*a == *b) { return true; }
3586  goto l1;
3587 }
3588 
3589 void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
3590  const Array<int> *search_set)
3591 {
3592  // TODO future: this function is inefficient. For a single element, an
3593  // octree neighbor search algorithm would be better. However, the octree
3594  // neighbor algorithm is hard to get right in the multi-octree case due to
3595  // the different orientations of the octrees (i.e., the root elements).
3596 
3598 
3599  // sorted list of all vertex nodes touching 'elem'
3600  Array<int> vert;
3601  vert.Reserve(128);
3602 
3603  // support for non-leaf 'elem', collect vertices of all children
3604  Array<int> stack;
3605  stack.Reserve(64);
3606  stack.Append(elem);
3607 
3608  while (stack.Size())
3609  {
3610  Element &el = elements[stack.Last()];
3611  stack.DeleteLast();
3612 
3613  if (!el.ref_type)
3614  {
3615  int *v = element_vertex.GetRow(el.index);
3616  int nv = element_vertex.RowSize(el.index);
3617  for (int i = 0; i < nv; i++)
3618  {
3619  vert.Append(v[i]);
3620  }
3621 
3622  nv = GI[el.Geom()].nv;
3623  for (int i = 0; i < nv; i++)
3624  {
3625  vert.Append(el.node[i]);
3626  }
3627  }
3628  else
3629  {
3630  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
3631  {
3632  stack.Append(el.child[i]);
3633  }
3634  }
3635  }
3636 
3637  vert.Sort();
3638  vert.Unique();
3639 
3640  int *v1 = vert.GetData();
3641  int nv1 = vert.Size();
3642 
3643  if (!search_set) { search_set = &leaf_elements; }
3644 
3645  // test *all* potential neighbors from the search set
3646  for (int i = 0; i < search_set->Size(); i++)
3647  {
3648  int testme = (*search_set)[i];
3649  if (testme != elem)
3650  {
3651  Element &el = elements[testme];
3652  int *v2 = element_vertex.GetRow(el.index);
3653  int nv2 = element_vertex.RowSize(el.index);
3654 
3655  bool hit = sorted_lists_intersect(v1, v2, nv1, nv2);
3656 
3657  if (!hit)
3658  {
3659  int nv = GI[el.Geom()].nv;
3660  for (int j = 0; j < nv; j++)
3661  {
3662  hit = sorted_lists_intersect(&el.node[j], v1, 1, nv1);
3663  if (hit) { break; }
3664  }
3665  }
3666 
3667  if (hit) { neighbors.Append(testme); }
3668  }
3669  }
3670 }
3671 
3673  Array<int> &expanded,
3674  const Array<int> *search_set)
3675 {
3677 
3678  Array<char> vmark(nodes.NumIds());
3679  vmark = 0;
3680 
3681  for (int i = 0; i < elems.Size(); i++)
3682  {
3683  Element &el = elements[elems[i]];
3684 
3685  int *v = element_vertex.GetRow(el.index);
3686  int nv = element_vertex.RowSize(el.index);
3687  for (int j = 0; j < nv; j++)
3688  {
3689  vmark[v[j]] = 1;
3690  }
3691 
3692  nv = GI[el.Geom()].nv;
3693  for (int j = 0; j < nv; j++)
3694  {
3695  vmark[el.node[j]] = 1;
3696  }
3697  }
3698 
3699  if (!search_set)
3700  {
3701  search_set = &leaf_elements;
3702  }
3703 
3704  expanded.SetSize(0);
3705  for (int i = 0; i < search_set->Size(); i++)
3706  {
3707  int testme = (*search_set)[i];
3708  Element &el = elements[testme];
3709  bool hit = false;
3710 
3711  int *v = element_vertex.GetRow(el.index);
3712  int nv = element_vertex.RowSize(el.index);
3713  for (int j = 0; j < nv; j++)
3714  {
3715  if (vmark[v[j]]) { hit = true; break; }
3716  }
3717 
3718  if (!hit)
3719  {
3720  nv = GI[el.Geom()].nv;
3721  for (int j = 0; j < nv; j++)
3722  {
3723  if (vmark[el.node[j]]) { hit = true; break; }
3724  }
3725  }
3726 
3727  if (hit) { expanded.Append(testme); }
3728  }
3729 }
3730 
3731 int NCMesh::GetVertexRootCoord(int elem, RefCoord coord[3]) const
3732 {
3733  while (1)
3734  {
3735  const Element &el = elements[elem];
3736  if (el.parent < 0) { return elem; }
3737 
3738  const Element &pa = elements[el.parent];
3739  MFEM_ASSERT(pa.ref_type, "internal error");
3740 
3741  int ch = 0;
3742  while (ch < 8 && pa.child[ch] != elem) { ch++; }
3743  MFEM_ASSERT(ch < 8, "internal error");
3744 
3745  MFEM_ASSERT(geom_parent[el.Geom()], "unsupported geometry");
3746  const RefTrf &tr = geom_parent[el.Geom()][(int) pa.ref_type][ch];
3747  tr.Apply(coord, coord);
3748 
3749  elem = el.parent;
3750  }
3751 }
3752 
3753 static bool RefPointInside(Geometry::Type geom, const RefCoord pt[3])
3754 {
3755  switch (geom)
3756  {
3757  case Geometry::SQUARE:
3758  return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
3759  (pt[1] >= 0) && (pt[1] <= T_ONE);
3760 
3761  case Geometry::CUBE:
3762  return (pt[0] >= 0) && (pt[0] <= T_ONE) &&
3763  (pt[1] >= 0) && (pt[1] <= T_ONE) &&
3764  (pt[2] >= 0) && (pt[2] <= T_ONE);
3765 
3766  case Geometry::TRIANGLE:
3767  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE);
3768 
3769  case Geometry::PRISM:
3770  return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE) &&
3771  (pt[2] >= 0) && (pt[2] <= T_ONE);
3772 
3773  default:
3774  MFEM_ABORT("unsupported geometry");
3775  return false;
3776  }
3777 }
3778 
3779 void NCMesh::CollectIncidentElements(int elem, const RefCoord coord[3],
3780  Array<int> &list) const
3781 {
3782  const Element &el = elements[elem];
3783  if (!el.ref_type)
3784  {
3785  list.Append(elem);
3786  return;
3787  }
3788 
3789  RefCoord tcoord[3];
3790  for (int ch = 0; ch < 8 && el.child[ch] >= 0; ch++)
3791  {
3792  const RefTrf &tr = geom_child[el.Geom()][(int) el.ref_type][ch];
3793  tr.Apply(coord, tcoord);
3794 
3795  if (RefPointInside(el.Geom(), tcoord))
3796  {
3797  CollectIncidentElements(el.child[ch], tcoord, list);
3798  }
3799  }
3800 }
3801 
3802 void NCMesh::FindVertexCousins(int elem, int local, Array<int> &cousins) const
3803 {
3804  const Element &el = elements[elem];
3805 
3806  RefCoord coord[3];
3807  MFEM_ASSERT(geom_corners[el.Geom()], "unsupported geometry");
3808  std::memcpy(coord, geom_corners[el.Geom()][local], sizeof(coord));
3809 
3810  int root = GetVertexRootCoord(elem, coord);
3811 
3812  cousins.SetSize(0);
3813  CollectIncidentElements(root, coord, cousins);
3814 }
3815 
3816 
3817 //// Coarse/fine transformations ///////////////////////////////////////////////
3818 
3820 {
3821  MFEM_ASSERT(np == pm.np, "");
3822  for (int i = 0; i < np; i++)
3823  {
3824  MFEM_ASSERT(points[i].dim == pm.points[i].dim, "");
3825  for (int j = 0; j < points[i].dim; j++)
3826  {
3827  if (points[i].coord[j] != pm.points[i].coord[j]) { return false; }
3828  }
3829  }
3830  return true;
3831 }
3832 
3834 {
3835  point_matrix.SetSize(points[0].dim, np);
3836  for (int i = 0; i < np; i++)
3837  {
3838  for (int j = 0; j < points[0].dim; j++)
3839  {
3840  point_matrix(j, i) = points[i].coord[j];
3841  }
3842  }
3843 }
3844 
3846  Point(0), Point(1)
3847 );
3849  Point(0, 0), Point(1, 0), Point(0, 1)
3850 );
3852  Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)
3853 );
3855  Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)
3856 );
3858  Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0),
3859  Point(0, 0, 1), Point(1, 0, 1), Point(0, 1, 1)
3860 );
3862  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
3863  Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
3864 );
3865 
3867 {
3868  switch (geom)
3869  {
3870  case Geometry::SEGMENT: return pm_seg_identity;
3871  case Geometry::TRIANGLE: return pm_tri_identity;
3872  case Geometry::SQUARE: return pm_quad_identity;
3874  case Geometry::PRISM: return pm_prism_identity;
3875  case Geometry::CUBE: return pm_hex_identity;
3876  default:
3877  MFEM_ABORT("unsupported geometry " << geom);
3878  return pm_tri_identity;
3879  }
3880 }
3881 
3882 void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
3883  DenseMatrix& matrix)
3884 {
3885  PointMatrix pm = GetGeomIdentity(geom);
3886 
3887  while (*ref_path)
3888  {
3889  int ref_type = *ref_path++;
3890  int child = *ref_path++;
3891 
3892  // TODO: do this with the new child transform tables
3893 
3894  if (geom == Geometry::CUBE)
3895  {
3896  if (ref_type == 1) // split along X axis
3897  {
3898  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
3899  Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
3900 
3901  if (child == 0)
3902  {
3903  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
3904  pm(4), mid45, mid67, pm(7));
3905  }
3906  else if (child == 1)
3907  {
3908  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
3909  mid45, pm(5), pm(6), mid67);
3910  }
3911  }
3912  else if (ref_type == 2) // split along Y axis
3913  {
3914  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
3915  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
3916 
3917  if (child == 0)
3918  {
3919  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
3920  pm(4), pm(5), mid56, mid74);
3921  }
3922  else if (child == 1)
3923  {
3924  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
3925  mid74, mid56, pm(6), pm(7));
3926  }
3927  }
3928  else if (ref_type == 4) // split along Z axis
3929  {
3930  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
3931  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
3932 
3933  if (child == 0)
3934  {
3935  pm = PointMatrix(pm(0), pm(1), pm(2), pm(3),
3936  mid04, mid15, mid26, mid37);
3937  }
3938  else if (child == 1)
3939  {
3940  pm = PointMatrix(mid04, mid15, mid26, mid37,
3941  pm(4), pm(5), pm(6), pm(7));
3942  }
3943  }
3944  else if (ref_type == 3) // XY split
3945  {
3946  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
3947  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
3948  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
3949  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
3950 
3951  Point midf0(mid23, mid12, mid01, mid30);
3952  Point midf5(mid45, mid56, mid67, mid74);
3953 
3954  if (child == 0)
3955  {
3956  pm = PointMatrix(pm(0), mid01, midf0, mid30,
3957  pm(4), mid45, midf5, mid74);
3958  }
3959  else if (child == 1)
3960  {
3961  pm = PointMatrix(mid01, pm(1), mid12, midf0,
3962  mid45, pm(5), mid56, midf5);
3963  }
3964  else if (child == 2)
3965  {
3966  pm = PointMatrix(midf0, mid12, pm(2), mid23,
3967  midf5, mid56, pm(6), mid67);
3968  }
3969  else if (child == 3)
3970  {
3971  pm = PointMatrix(mid30, midf0, mid23, pm(3),
3972  mid74, midf5, mid67, pm(7));
3973  }
3974  }
3975  else if (ref_type == 5) // XZ split
3976  {
3977  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
3978  Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
3979  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
3980  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
3981 
3982  Point midf1(mid01, mid15, mid45, mid04);
3983  Point midf3(mid23, mid37, mid67, mid26);
3984 
3985  if (child == 0)
3986  {
3987  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
3988  mid04, midf1, midf3, mid37);
3989  }
3990  else if (child == 1)
3991  {
3992  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
3993  midf1, mid15, mid26, midf3);
3994  }
3995  else if (child == 2)
3996  {
3997  pm = PointMatrix(midf1, mid15, mid26, midf3,
3998  mid45, pm(5), pm(6), mid67);
3999  }
4000  else if (child == 3)
4001  {
4002  pm = PointMatrix(mid04, midf1, midf3, mid37,
4003  pm(4), mid45, mid67, pm(7));
4004  }
4005  }
4006  else if (ref_type == 6) // YZ split
4007  {
4008  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4009  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
4010  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4011  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4012 
4013  Point midf2(mid12, mid26, mid56, mid15);
4014  Point midf4(mid30, mid04, mid74, mid37);
4015 
4016  if (child == 0)
4017  {
4018  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
4019  mid04, mid15, midf2, midf4);
4020  }
4021  else if (child == 1)
4022  {
4023  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
4024  midf4, midf2, mid26, mid37);
4025  }
4026  else if (child == 2)
4027  {
4028  pm = PointMatrix(mid04, mid15, midf2, midf4,
4029  pm(4), pm(5), mid56, mid74);
4030  }
4031  else if (child == 3)
4032  {
4033  pm = PointMatrix(midf4, midf2, mid26, mid37,
4034  mid74, mid56, pm(6), pm(7));
4035  }
4036  }
4037  else if (ref_type == 7) // full isotropic refinement
4038  {
4039  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4040  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4041  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
4042  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
4043  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
4044  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
4045 
4046  Point midf0(mid23, mid12, mid01, mid30);
4047  Point midf1(mid01, mid15, mid45, mid04);
4048  Point midf2(mid12, mid26, mid56, mid15);
4049  Point midf3(mid23, mid37, mid67, mid26);
4050  Point midf4(mid30, mid04, mid74, mid37);
4051  Point midf5(mid45, mid56, mid67, mid74);
4052 
4053  Point midel(midf1, midf3);
4054 
4055  if (child == 0)
4056  {
4057  pm = PointMatrix(pm(0), mid01, midf0, mid30,
4058  mid04, midf1, midel, midf4);
4059  }
4060  else if (child == 1)
4061  {
4062  pm = PointMatrix(mid01, pm(1), mid12, midf0,
4063  midf1, mid15, midf2, midel);
4064  }
4065  else if (child == 2)
4066  {
4067  pm = PointMatrix(midf0, mid12, pm(2), mid23,
4068  midel, midf2, mid26, midf3);
4069  }
4070  else if (child == 3)
4071  {
4072  pm = PointMatrix(mid30, midf0, mid23, pm(3),
4073  midf4, midel, midf3, mid37);
4074  }
4075  else if (child == 4)
4076  {
4077  pm = PointMatrix(mid04, midf1, midel, midf4,
4078  pm(4), mid45, midf5, mid74);
4079  }
4080  else if (child == 5)
4081  {
4082  pm = PointMatrix(midf1, mid15, midf2, midel,
4083  mid45, pm(5), mid56, midf5);
4084  }
4085  else if (child == 6)
4086  {
4087  pm = PointMatrix(midel, midf2, mid26, midf3,
4088  midf5, mid56, pm(6), mid67);
4089  }
4090  else if (child == 7)
4091  {
4092  pm = PointMatrix(midf4, midel, midf3, mid37,
4093  mid74, midf5, mid67, pm(7));
4094  }
4095  }
4096  }
4097  else if (geom == Geometry::PRISM)
4098  {
4099  if (ref_type < 4) // XY split
4100  {
4101  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4102  Point mid20(pm(2), pm(0)), mid34(pm(3), pm(4));
4103  Point mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4104 
4105  if (child == 0)
4106  {
4107  pm = PointMatrix(pm(0), mid01, mid20, pm(3), mid34, mid53);
4108  }
4109  else if (child == 1)
4110  {
4111  pm = PointMatrix(mid01, pm(1), mid12, mid34, pm(4), mid45);
4112  }
4113  else if (child == 2)
4114  {
4115  pm = PointMatrix(mid20, mid12, pm(2), mid53, mid45, pm(5));
4116  }
4117  else if (child == 3)
4118  {
4119  pm = PointMatrix(mid12, mid20, mid01, mid45, mid53, mid34);
4120  }
4121  }
4122  else if (ref_type == 4) // Z split
4123  {
4124  Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4125 
4126  if (child == 0)
4127  {
4128  pm = PointMatrix(pm(0), pm(1), pm(2), mid03, mid14, mid25);
4129  }
4130  else if (child == 1)
4131  {
4132  pm = PointMatrix(mid03, mid14, mid25, pm(3), pm(4), pm(5));
4133  }
4134  }
4135  else // ref_type > 4, iso split
4136  {
4137  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4138  Point mid34(pm(3), pm(4)), mid45(pm(4), pm(5)), mid53(pm(5), pm(3));
4139  Point mid03(pm(0), pm(3)), mid14(pm(1), pm(4)), mid25(pm(2), pm(5));
4140 
4141  Point midf2(mid01, mid14, mid34, mid03);
4142  Point midf3(mid12, mid25, mid45, mid14);
4143  Point midf4(mid20, mid03, mid53, mid25);
4144 
4145  if (child == 0)
4146  {
4147  pm = PointMatrix(pm(0), mid01, mid20, mid03, midf2, midf4);
4148  }
4149  else if (child == 1)
4150  {
4151  pm = PointMatrix(mid01, pm(1), mid12, midf2, mid14, midf3);
4152  }
4153  else if (child == 2)
4154  {
4155  pm = PointMatrix(mid20, mid12, pm(2), midf4, midf3, mid25);
4156  }
4157  else if (child == 3)
4158  {
4159  pm = PointMatrix(mid12, mid20, mid01, midf3, midf4, midf2);
4160  }
4161  else if (child == 4)
4162  {
4163  pm = PointMatrix(mid03, midf2, midf4, pm(3), mid34, mid53);
4164  }
4165  else if (child == 5)
4166  {
4167  pm = PointMatrix(midf2, mid14, midf3, mid34, pm(4), mid45);
4168  }
4169  else if (child == 6)
4170  {
4171  pm = PointMatrix(midf4, midf3, mid25, mid53, mid45, pm(5));
4172  }
4173  else if (child == 7)
4174  {
4175  pm = PointMatrix(midf3, midf4, midf2, mid45, mid53, mid34);
4176  }
4177  }
4178  }
4179  else if (geom == Geometry::TETRAHEDRON)
4180  {
4181  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid02(pm(2), pm(0));
4182  Point mid03(pm(0), pm(3)), mid13(pm(1), pm(3)), mid23(pm(2), pm(3));
4183 
4184  if (child == 0)
4185  {
4186  pm = PointMatrix(pm(0), mid01, mid02, mid03);
4187  }
4188  else if (child == 1)
4189  {
4190  pm = PointMatrix(mid01, pm(1), mid12, mid13);
4191  }
4192  else if (child == 2)
4193  {
4194  pm = PointMatrix(mid02, mid12, pm(2), mid23);
4195  }
4196  else if (child == 3)
4197  {
4198  pm = PointMatrix(mid03, mid13, mid23, pm(3));
4199  }
4200  else if (child == 4)
4201  {
4202  pm = PointMatrix(mid01, mid23, mid02, mid03);
4203  }
4204  else if (child == 5)
4205  {
4206  pm = PointMatrix(mid01, mid23, mid03, mid13);
4207  }
4208  else if (child == 6)
4209  {
4210  pm = PointMatrix(mid01, mid23, mid13, mid12);
4211  }
4212  else if (child == 7)
4213  {
4214  pm = PointMatrix(mid01, mid23, mid12, mid02);
4215  }
4216  }
4217  else if (geom == Geometry::SQUARE)
4218  {
4219  if (ref_type == 1) // X split
4220  {
4221  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
4222 
4223  if (child == 0)
4224  {
4225  pm = PointMatrix(pm(0), mid01, mid23, pm(3));
4226  }
4227  else if (child == 1)
4228  {
4229  pm = PointMatrix(mid01, pm(1), pm(2), mid23);
4230  }
4231  }
4232  else if (ref_type == 2) // Y split
4233  {
4234  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
4235 
4236  if (child == 0)
4237  {
4238  pm = PointMatrix(pm(0), pm(1), mid12, mid30);
4239  }
4240  else if (child == 1)
4241  {
4242  pm = PointMatrix(mid30, mid12, pm(2), pm(3));
4243  }
4244  }
4245  else if (ref_type == 3) // iso split
4246  {
4247  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
4248  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
4249  Point midel(mid01, mid23);
4250 
4251  if (child == 0)
4252  {
4253  pm = PointMatrix(pm(0), mid01, midel, mid30);
4254  }
4255  else if (child == 1)
4256  {
4257  pm = PointMatrix(mid01, pm(1), mid12, midel);
4258  }
4259  else if (child == 2)
4260  {
4261  pm = PointMatrix(midel, mid12, pm(2), mid23);
4262  }
4263  else if (child == 3)
4264  {
4265  pm = PointMatrix(mid30, midel, mid23, pm(3));
4266  }
4267  }
4268  }
4269  else if (geom == Geometry::TRIANGLE)
4270  {
4271  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
4272 
4273  if (child == 0)
4274  {
4275  pm = PointMatrix(pm(0), mid01, mid20);
4276  }
4277  else if (child == 1)
4278  {
4279  pm = PointMatrix(mid01, pm(1), mid12);
4280  }
4281  else if (child == 2)
4282  {
4283  pm = PointMatrix(mid20, mid12, pm(2));
4284  }
4285  else if (child == 3)
4286  {
4287  pm = PointMatrix(mid12, mid20, mid01);
4288  }
4289  }
4290  }
4291 
4292  // write the points to the matrix
4293  for (int i = 0; i < pm.np; i++)
4294  {
4295  for (int j = 0; j < pm(i).dim; j++)
4296  {
4297  matrix(j, i) = pm(i).coord[j];
4298  }
4299  }
4300 }
4301 
4303 {
4306 
4307  for (int i = 0; i < leaf_elements.Size(); i++)
4308  {
4309  int elem = leaf_elements[i];
4310  if (!IsGhost(elements[elem])) { coarse_elements.Append(elem); }
4311  }
4312 
4313  transforms.embeddings.DeleteAll();
4314 }
4315 
4316 void NCMesh::TraverseRefinements(int elem, int coarse_index,
4317  std::string &ref_path, RefPathMap &map)
4318 {
4319  Element &el = elements[elem];
4320  if (!el.ref_type)
4321  {
4322  int &matrix = map[ref_path];
4323  if (!matrix) { matrix = map.size(); }
4324 
4325  Embedding &emb = transforms.embeddings[el.index];
4326  emb.parent = coarse_index;
4327  emb.matrix = matrix - 1;
4328  emb.geom = el.Geom();
4329  emb.ghost = IsGhost(el);
4330  }
4331  else
4332  {
4333  MFEM_ASSERT(el.tet_type == 0, "not implemented");
4334 
4335  ref_path.push_back(el.ref_type);
4336  ref_path.push_back(0);
4337 
4338  for (int i = 0; i < 8; i++)
4339  {
4340  if (el.child[i] >= 0)
4341  {
4342  ref_path[ref_path.length()-1] = i;
4343  TraverseRefinements(el.child[i], coarse_index, ref_path, map);
4344  }
4345  }
4346  ref_path.resize(ref_path.length()-2);
4347  }
4348 }
4349 
4351 {
4352  MFEM_VERIFY(coarse_elements.Size() || !leaf_elements.Size(),
4353  "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
4354  " and Refine().");
4355 
4356  if (!transforms.embeddings.Size())
4357  {
4358  transforms.Clear();
4359  transforms.embeddings.SetSize(NElements);
4360 
4361  std::string ref_path;
4362  ref_path.reserve(100);
4363 
4364  RefPathMap path_map[Geometry::NumGeom];
4365  for (int g = 0; g < Geometry::NumGeom; g++)
4366  {
4367  path_map[g][ref_path] = 1; // empty path == identity
4368  }
4369 
4370  int used_geoms = 0;
4371  for (int i = 0; i < coarse_elements.Size(); i++)
4372  {
4373  int geom = elements[coarse_elements[i]].geom;
4374  TraverseRefinements(coarse_elements[i], i, ref_path, path_map[geom]);
4375  used_geoms |= (1 << geom);
4376  }
4377 
4378  for (int g = 0; g < Geometry::NumGeom; g++)
4379  {
4380  if (used_geoms & (1 << g))
4381  {
4382  Geometry::Type geom = Geometry::Type(g);
4383  const PointMatrix &identity = GetGeomIdentity(geom);
4384 
4386  .SetSize(Dim, identity.np, path_map[g].size());
4387 
4388  // calculate the point matrices
4389  RefPathMap::iterator it;
4390  for (it = path_map[g].begin(); it != path_map[g].end(); ++it)
4391  {
4392  GetPointMatrix(geom, it->first.c_str(),
4393  transforms.point_matrices[g](it->second-1));
4394  }
4395  }
4396  }
4397  }
4398  return transforms;
4399 }
4400 
4402 {
4403  MFEM_VERIFY(transforms.embeddings.Size() || !leaf_elements.Size(),
4404  "GetDerefinementTransforms() must be preceded by Derefine().");
4405 
4406  if (!transforms.IsInitialized())
4407  {
4408  std::map<int, int> mat_no[Geometry::NumGeom];
4409  for (int g = 0; g < Geometry::NumGeom; g++)
4410  {
4411  mat_no[g][0] = 1; // 0 == identity
4412  }
4413 
4414  // assign numbers to the different matrices used
4415  for (int i = 0; i < transforms.embeddings.Size(); i++)
4416  {
4417  Embedding &emb = transforms.embeddings[i];
4418  int code = emb.matrix; // see SetDerefMatrixCodes()
4419  if (code)
4420  {
4421  int &matrix = mat_no[emb.geom][code];
4422  if (!matrix) { matrix = mat_no[emb.geom].size(); }
4423 
4424  emb.matrix = matrix - 1;
4425  }
4426  }
4427 
4428  for (int g = 0; g < Geometry::NumGeom; g++)
4429  {
4430  if (Geoms & (1 << g))
4431  {
4432  Geometry::Type geom = Geometry::Type(g);
4433  const PointMatrix &identity = GetGeomIdentity(geom);
4434 
4436  .SetSize(Dim, identity.np, mat_no[geom].size());
4437 
4438  // calculate point matrices
4439  for (auto it = mat_no[geom].begin(); it != mat_no[geom].end(); ++it)
4440  {
4441  char path[3] = { 0, 0, 0 };
4442 
4443  int code = it->first;
4444  if (code)
4445  {
4446  path[0] = code >> 4; // ref_type (see SetDerefMatrixCodes())
4447  path[1] = code & 0xf; // child
4448  }
4449 
4450  GetPointMatrix(geom, path,
4451  transforms.point_matrices[geom](it->second-1));
4452  }
4453  }
4454  }
4455  }
4456  return transforms;
4457 }
4458 
4460  bool want_ghosts) const
4461 {
4462  Array<Connection> conn;
4463  conn.Reserve(embeddings.Size());
4464 
4465  int max_parent = -1;
4466  for (int i = 0; i < embeddings.Size(); i++)
4467  {
4468  const Embedding &emb = embeddings[i];
4469  if ((emb.parent >= 0) &&
4470  (!emb.ghost || want_ghosts))
4471  {
4472  conn.Append(Connection(emb.parent, i));
4473  max_parent = std::max(emb.parent, max_parent);
4474  }
4475  }
4476 
4477  conn.Sort(); // NOTE: unique is not necessary
4478  coarse_to_fine.MakeFromList(max_parent+1, conn);
4479 }
4480 
4482 {
4484  transforms.Clear();
4485 }
4486 
4488 {
4489  for (int i = 0; i < Geometry::NumGeom; i++)
4490  {
4491  point_matrices[i].SetSize(0, 0, 0);
4492  }
4493  embeddings.DeleteAll();
4494 }
4495 
4497 {
4498  // return true if point matrices are present for any geometry
4499  for (int i = 0; i < Geometry::NumGeom; i++)
4500  {
4501  if (point_matrices[i].SizeK()) { return true; }
4502  }
4503  return false;
4504 }
4505 
4507 {
4508  for (int g = 0; g < Geometry::NumGeom; ++g)
4509  {
4510  a.point_matrices[g].Swap(b.point_matrices[g]);
4511  }
4512  Swap(a.embeddings, b.embeddings);
4513 }
4514 
4515 
4516 //// SFC Ordering //////////////////////////////////////////////////////////////
4517 
4518 static int sgn(int x)
4519 {
4520  return (x < 0) ? -1 : (x > 0) ? 1 : 0;
4521 }
4522 
4523 static void HilbertSfc2D(int x, int y, int ax, int ay, int bx, int by,
4524  Array<int> &coords)
4525 {
4526  int w = std::abs(ax + ay);
4527  int h = std::abs(bx + by);
4528 
4529  int dax = sgn(ax), day = sgn(ay); // unit major direction ("right")
4530  int dbx = sgn(bx), dby = sgn(by); // unit orthogonal direction ("up")
4531 
4532  if (h == 1) // trivial row fill
4533  {
4534  for (int i = 0; i < w; i++, x += dax, y += day)
4535  {
4536  coords.Append(x);
4537  coords.Append(y);
4538  }
4539  return;
4540  }
4541  if (w == 1) // trivial column fill
4542  {
4543  for (int i = 0; i < h; i++, x += dbx, y += dby)
4544  {
4545  coords.Append(x);
4546  coords.Append(y);
4547  }
4548  return;
4549  }
4550 
4551  int ax2 = ax/2, ay2 = ay/2;
4552  int bx2 = bx/2, by2 = by/2;
4553 
4554  int w2 = std::abs(ax2 + ay2);
4555  int h2 = std::abs(bx2 + by2);
4556 
4557  if (2*w > 3*h) // long case: split in two parts only
4558  {
4559  if ((w2 & 0x1) && (w > 2))
4560  {
4561  ax2 += dax, ay2 += day; // prefer even steps
4562  }
4563 
4564  HilbertSfc2D(x, y, ax2, ay2, bx, by, coords);
4565  HilbertSfc2D(x+ax2, y+ay2, ax-ax2, ay-ay2, bx, by, coords);
4566  }
4567  else // standard case: one step up, one long horizontal step, one step down
4568  {
4569  if ((h2 & 0x1) && (h > 2))
4570  {
4571  bx2 += dbx, by2 += dby; // prefer even steps
4572  }
4573 
4574  HilbertSfc2D(x, y, bx2, by2, ax2, ay2, coords);
4575  HilbertSfc2D(x+bx2, y+by2, ax, ay, bx-bx2, by-by2, coords);
4576  HilbertSfc2D(x+(ax-dax)+(bx2-dbx), y+(ay-day)+(by2-dby),
4577  -bx2, -by2, -(ax-ax2), -(ay-ay2), coords);
4578  }
4579 }
4580 
4581 static void HilbertSfc3D(int x, int y, int z,
4582  int ax, int ay, int az,
4583  int bx, int by, int bz,
4584  int cx, int cy, int cz,
4585  Array<int> &coords)
4586 {
4587  int w = std::abs(ax + ay + az);
4588  int h = std::abs(bx + by + bz);
4589  int d = std::abs(cx + cy + cz);
4590 
4591  int dax = sgn(ax), day = sgn(ay), daz = sgn(az); // unit major dir ("right")
4592  int dbx = sgn(bx), dby = sgn(by), dbz = sgn(bz); // unit ortho dir ("forward")
4593  int dcx = sgn(cx), dcy = sgn(cy), dcz = sgn(cz); // unit ortho dir ("up")
4594 
4595  // trivial row/column fills
4596  if (h == 1 && d == 1)
4597  {
4598  for (int i = 0; i < w; i++, x += dax, y += day, z += daz)
4599  {
4600  coords.Append(x);
4601  coords.Append(y);
4602  coords.Append(z);
4603  }
4604  return;
4605  }
4606  if (w == 1 && d == 1)
4607  {
4608  for (int i = 0; i < h; i++, x += dbx, y += dby, z += dbz)
4609  {
4610  coords.Append(x);
4611  coords.Append(y);
4612  coords.Append(z);
4613  }
4614  return;
4615  }
4616  if (w == 1 && h == 1)
4617  {
4618  for (int i = 0; i < d; i++, x += dcx, y += dcy, z += dcz)
4619  {
4620  coords.Append(x);
4621  coords.Append(y);
4622  coords.Append(z);
4623  }
4624  return;
4625  }
4626 
4627  int ax2 = ax/2, ay2 = ay/2, az2 = az/2;
4628  int bx2 = bx/2, by2 = by/2, bz2 = bz/2;
4629  int cx2 = cx/2, cy2 = cy/2, cz2 = cz/2;
4630 
4631  int w2 = std::abs(ax2 + ay2 + az2);
4632  int h2 = std::abs(bx2 + by2 + bz2);
4633  int d2 = std::abs(cx2 + cy2 + cz2);
4634 
4635  // prefer even steps
4636  if ((w2 & 0x1) && (w > 2))
4637  {
4638  ax2 += dax, ay2 += day, az2 += daz;
4639  }
4640  if ((h2 & 0x1) && (h > 2))
4641  {
4642  bx2 += dbx, by2 += dby, bz2 += dbz;
4643  }
4644  if ((d2 & 0x1) && (d > 2))
4645  {
4646  cx2 += dcx, cy2 += dcy, cz2 += dcz;
4647  }
4648 
4649  // wide case, split in w only
4650  if ((2*w > 3*h) && (2*w > 3*d))
4651  {
4652  HilbertSfc3D(x, y, z,
4653  ax2, ay2, az2,
4654  bx, by, bz,
4655  cx, cy, cz, coords);
4656 
4657  HilbertSfc3D(x+ax2, y+ay2, z+az2,
4658  ax-ax2, ay-ay2, az-az2,
4659  bx, by, bz,
4660  cx, cy, cz, coords);
4661  }
4662  // do not split in d
4663  else if (3*h > 4*d)
4664  {
4665  HilbertSfc3D(x, y, z,
4666  bx2, by2, bz2,
4667  cx, cy, cz,
4668  ax2, ay2, az2, coords);
4669 
4670  HilbertSfc3D(x+bx2, y+by2, z+bz2,
4671  ax, ay, az,
4672  bx-bx2, by-by2, bz-bz2,
4673  cx, cy, cz, coords);
4674 
4675  HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
4676  y+(ay-day)+(by2-dby),
4677  z+(az-daz)+(bz2-dbz),
4678  -bx2, -by2, -bz2,
4679  cx, cy, cz,
4680  -(ax-ax2), -(ay-ay2), -(az-az2), coords);
4681  }
4682  // do not split in h
4683  else if (3*d > 4*h)
4684  {
4685  HilbertSfc3D(x, y, z,
4686  cx2, cy2, cz2,
4687  ax2, ay2, az2,
4688  bx, by, bz, coords);
4689 
4690  HilbertSfc3D(x+cx2, y+cy2, z+cz2,
4691  ax, ay, az,
4692  bx, by, bz,
4693  cx-cx2, cy-cy2, cz-cz2, coords);
4694 
4695  HilbertSfc3D(x+(ax-dax)+(cx2-dcx),
4696  y+(ay-day)+(cy2-dcy),
4697  z+(az-daz)+(cz2-dcz),
4698  -cx2, -cy2, -cz2,
4699  -(ax-ax2), -(ay-ay2), -(az-az2),
4700  bx, by, bz, coords);
4701  }
4702  // regular case, split in all w/h/d
4703  else
4704  {
4705  HilbertSfc3D(x, y, z,
4706  bx2, by2, bz2,
4707  cx2, cy2, cz2,
4708  ax2, ay2, az2, coords);
4709 
4710  HilbertSfc3D(x+bx2, y+by2, z+bz2,
4711  cx, cy, cz,
4712  ax2, ay2, az2,
4713  bx-bx2, by-by2, bz-bz2, coords);
4714 
4715  HilbertSfc3D(x+(bx2-dbx)+(cx-dcx),
4716  y+(by2-dby)+(cy-dcy),
4717  z+(bz2-dbz)+(cz-dcz),
4718  ax, ay, az,
4719  -bx2, -by2, -bz2,
4720  -(cx-cx2), -(cy-cy2), -(cz-cz2), coords);
4721 
4722  HilbertSfc3D(x+(ax-dax)+bx2+(cx-dcx),
4723  y+(ay-day)+by2+(cy-dcy),
4724  z+(az-daz)+bz2+(cz-dcz),
4725  -cx, -cy, -cz,
4726  -(ax-ax2), -(ay-ay2), -(az-az2),
4727  bx-bx2, by-by2, bz-bz2, coords);
4728 
4729  HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
4730  y+(ay-day)+(by2-dby),
4731  z+(az-daz)+(bz2-dbz),
4732  -bx2, -by2, -bz2,
4733  cx2, cy2, cz2,
4734  -(ax-ax2), -(ay-ay2), -(az-az2), coords);
4735  }
4736 }
4737 
4738 void NCMesh::GridSfcOrdering2D(int width, int height, Array<int> &coords)
4739 {
4740  coords.SetSize(0);
4741  coords.Reserve(2*width*height);
4742 
4743  if (width >= height)
4744  {
4745  HilbertSfc2D(0, 0, width, 0, 0, height, coords);
4746  }
4747  else
4748  {
4749  HilbertSfc2D(0, 0, 0, height, width, 0, coords);
4750  }
4751 }
4752 
4753 void NCMesh::GridSfcOrdering3D(int width, int height, int depth,
4754  Array<int> &coords)
4755 {
4756  coords.SetSize(0);
4757  coords.Reserve(3*width*height*depth);
4758 
4759  if (width >= height && width >= depth)
4760  {
4761  HilbertSfc3D(0, 0, 0,
4762  width, 0, 0,
4763  0, height, 0,
4764  0, 0, depth, coords);
4765  }
4766  else if (height >= width && height >= depth)
4767  {
4768  HilbertSfc3D(0, 0, 0,
4769  0, height, 0,
4770  width, 0, 0,
4771  0, 0, depth, coords);
4772  }
4773  else // depth >= width && depth >= height
4774  {
4775  HilbertSfc3D(0, 0, 0,
4776  0, 0, depth,
4777  width, 0, 0,
4778  0, height, 0, coords);
4779  }
4780 }
4781 
4782 
4783 //// Utility ///////////////////////////////////////////////////////////////////
4784 
4785 void NCMesh::GetEdgeVertices(const MeshId &edge_id, int vert_index[2],
4786  bool oriented) const
4787 {
4788  const Element &el = elements[edge_id.element];
4789  const GeomInfo& gi = GI[el.Geom()];
4790  const int* ev = gi.edges[(int) edge_id.local];
4791 
4792  int n0 = el.node[ev[0]], n1 = el.node[ev[1]];
4793  if (n0 > n1) { std::swap(n0, n1); }
4794 
4795  vert_index[0] = nodes[n0].vert_index;
4796  vert_index[1] = nodes[n1].vert_index;
4797 
4798  if (oriented && vert_index[0] > vert_index[1])
4799  {
4800  std::swap(vert_index[0], vert_index[1]);
4801  }
4802 }
4803 
4805 {
4806  const Element &el = elements[edge_id.element];
4807  const GeomInfo& gi = GI[el.Geom()];
4808  const int* ev = gi.edges[(int) edge_id.local];
4809 
4810  int v0 = nodes[el.node[ev[0]]].vert_index;
4811  int v1 = nodes[el.node[ev[1]]].vert_index;
4812 
4813  return ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1])) ? -1 : 1;
4814 }
4815 
4817  int vert_index[4], int edge_index[4],
4818  int edge_orientation[4]) const
4819 {
4820  MFEM_ASSERT(Dim >= 3, "");
4821 
4822  const Element &el = elements[face_id.element];
4823  const GeomInfo& gi = GI[el.Geom()];
4824 
4825  const int *fv = gi.faces[(int) face_id.local];
4826  const int nfv = gi.nfv[(int) face_id.local];
4827 
4828  vert_index[3] = edge_index[3] = -1;
4829  edge_orientation[3] = 0;
4830 
4831  for (int i = 0; i < nfv; i++)
4832  {
4833  vert_index[i] = nodes[el.node[fv[i]]].vert_index;
4834  }
4835 
4836  for (int i = 0; i < nfv; i++)
4837  {
4838  int j = i+1;
4839  if (j >= nfv) { j = 0; }
4840 
4841  int n1 = el.node[fv[i]];
4842  int n2 = el.node[fv[j]];
4843 
4844  const Node* en = nodes.Find(n1, n2);
4845  MFEM_ASSERT(en != NULL, "edge not found.");
4846 
4847  edge_index[i] = en->edge_index;
4848  edge_orientation[i] = (vert_index[i] < vert_index[j]) ? 1 : -1;
4849  }
4850 
4851  return nfv;
4852 }
4853 
4854 int NCMesh::GetEdgeMaster(int node) const
4855 {
4856  MFEM_ASSERT(node >= 0, "edge node not found.");
4857  const Node &nd = nodes[node];
4858 
4859  int p1 = nd.p1, p2 = nd.p2;
4860  MFEM_ASSERT(p1 != p2, "invalid edge node.");
4861 
4862  const Node &n1 = nodes[p1], &n2 = nodes[p2];
4863 
4864  int n1p1 = n1.p1, n1p2 = n1.p2;
4865  int n2p1 = n2.p1, n2p2 = n2.p2;
4866 
4867  if ((n2p1 != n2p2) && (p1 == n2p1 || p1 == n2p2))
4868  {
4869  // n1 is parent of n2:
4870  // (n1)--(nd)--(n2)------(*)
4871  if (n2.HasEdge()) { return p2; }
4872  else { return GetEdgeMaster(p2); }
4873  }
4874 
4875  if ((n1p1 != n1p2) && (p2 == n1p1 || p2 == n1p2))
4876  {
4877  // n2 is parent of n1:
4878  // (n2)--(nd)--(n1)------(*)
4879  if (n1.HasEdge()) { return p1; }
4880  else { return GetEdgeMaster(p1); }
4881  }
4882 
4883  return -1;
4884 }
4885 
4886 int NCMesh::GetEdgeMaster(int v1, int v2) const
4887 {
4888  int node = nodes.FindId(vertex_nodeId[v1], vertex_nodeId[v2]);
4889  MFEM_ASSERT(node >= 0 && nodes[node].HasEdge(), "(v1, v2) is not an edge.");
4890 
4891  int master = GetEdgeMaster(node);
4892  return (master >= 0) ? nodes[master].edge_index : -1;
4893 }
4894 
4895 int NCMesh::GetElementDepth(int i) const
4896 {
4897  int elem = leaf_elements[i];
4898  int depth = 0, parent;
4899  while ((parent = elements[elem].parent) != -1)
4900  {
4901  elem = parent;
4902  depth++;
4903  }
4904  return depth;
4905 }
4906 
4908 {
4909  int elem = leaf_elements[i];
4910  int parent, reduction = 1;
4911  while ((parent = elements[elem].parent) != -1)
4912  {
4913  if (elements[parent].ref_type & 1) { reduction *= 2; }
4914  if (elements[parent].ref_type & 2) { reduction *= 2; }
4915  if (elements[parent].ref_type & 4) { reduction *= 2; }
4916  elem = parent;
4917  }
4918  return reduction;
4919 }
4920 
4922  Array<int> &face_indices,
4923  Array<int> &face_attribs) const
4924 {
4925  const Element &el = elements[leaf_elements[leaf_elem]];
4926  const GeomInfo& gi = GI[el.Geom()];
4927 
4928  face_indices.SetSize(gi.nf);
4929  face_attribs.SetSize(gi.nf);
4930 
4931  for (int i = 0; i < gi.nf; i++)
4932  {
4933  const int* fv = gi.faces[i];
4934  const Face *face = faces.Find(el.node[fv[0]], el.node[fv[1]],
4935  el.node[fv[2]], el.node[fv[3]]);
4936  MFEM_ASSERT(face, "face not found");
4937  face_indices[i] = face->index;
4938  face_attribs[i] = face->attribute;
4939  }
4940 }
4941 
4942 void NCMesh::FindFaceNodes(int face, int node[4])
4943 {
4944  // Obtain face nodes from one of its elements (note that face->p1, p2, p3
4945  // cannot be used directly since they are not in order and p4 is missing).
4946 
4947  Face &fa = faces[face];
4948 
4949  int elem = fa.elem[0];
4950  if (elem < 0) { elem = fa.elem[1]; }
4951  MFEM_ASSERT(elem >= 0, "Face has no elements?");
4952 
4953  Element &el = elements[elem];
4954  int f = find_local_face(el.Geom(),
4955  find_node(el, fa.p1),
4956  find_node(el, fa.p2),
4957  find_node(el, fa.p3));
4958 
4959  const int* fv = GI[el.Geom()].faces[f];
4960  for (int i = 0; i < 4; i++)
4961  {
4962  node[i] = el.node[fv[i]];
4963  }
4964 }
4965 
4966 void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
4967  Array<int> &bdr_vertices, Array<int> &bdr_edges)
4968 {
4969  bdr_vertices.SetSize(0);
4970  bdr_edges.SetSize(0);
4971 
4972  if (Dim == 3)
4973  {
4974  GetFaceList(); // make sure 'boundary_faces' is up to date
4975 
4976  for (int i = 0; i < boundary_faces.Size(); i++)
4977  {
4978  int face = boundary_faces[i];
4979  if (bdr_attr_is_ess[faces[face].attribute - 1])
4980  {
4981  int node[4];
4982  FindFaceNodes(face, node);
4983  int nfv = (node[3] < 0) ? 3 : 4;
4984 
4985  for (int j = 0; j < nfv; j++)
4986  {
4987  bdr_vertices.Append(nodes[node[j]].vert_index);
4988 
4989  int enode = nodes.FindId(node[j], node[(j+1) % nfv]);
4990  MFEM_ASSERT(enode >= 0 && nodes[enode].HasEdge(), "Edge not found.");
4991  bdr_edges.Append(nodes[enode].edge_index);
4992 
4993  while ((enode = GetEdgeMaster(enode)) >= 0)
4994  {
4995  // append master edges that may not be accessible from any
4996  // boundary element, this happens in 3D in re-entrant corners
4997  bdr_edges.Append(nodes[enode].edge_index);
4998  }
4999  }
5000  }
5001  }
5002  }
5003  else if (Dim == 2)
5004  {
5005  GetEdgeList(); // make sure 'boundary_faces' is up to date
5006 
5007  for (int i = 0; i < boundary_faces.Size(); i++)
5008  {
5009  int face = boundary_faces[i];
5010  Face &fc = faces[face];
5011  if (bdr_attr_is_ess[fc.attribute - 1])
5012  {
5013  bdr_vertices.Append(nodes[fc.p1].vert_index);
5014  bdr_vertices.Append(nodes[fc.p3].vert_index);
5015  }
5016  }
5017  }
5018 
5019  bdr_vertices.Sort();
5020  bdr_vertices.Unique();
5021 
5022  bdr_edges.Sort();
5023  bdr_edges.Unique();
5024 }
5025 
5026 static int max4(int a, int b, int c, int d)
5027 {
5028  return std::max(std::max(a, b), std::max(c, d));
5029 }
5030 static int max6(int a, int b, int c, int d, int e, int f)
5031 {
5032  return std::max(max4(a, b, c, d), std::max(e, f));
5033 }
5034 static int max8(int a, int b, int c, int d, int e, int f, int g, int h)
5035 {
5036  return std::max(max4(a, b, c, d), max4(e, f, g, h));
5037 }
5038 
5039 int NCMesh::EdgeSplitLevel(int vn1, int vn2) const
5040 {
5041  int mid = nodes.FindId(vn1, vn2);
5042  if (mid < 0 || !nodes[mid].HasVertex()) { return 0; }
5043  return 1 + std::max(EdgeSplitLevel(vn1, mid), EdgeSplitLevel(mid, vn2));
5044 }
5045 
5046 int NCMesh::TriFaceSplitLevel(int vn1, int vn2, int vn3) const
5047 {
5048  int mid[3];
5049  if (TriFaceSplit(vn1, vn2, vn3, mid) &&
5050  faces.FindId(vn1, vn2, vn3) < 0)
5051  {
5052  return 1 + max4(TriFaceSplitLevel(vn1, mid[0], mid[2]),
5053  TriFaceSplitLevel(mid[0], vn2, mid[1]),
5054  TriFaceSplitLevel(mid[2], mid[1], vn3),
5055  TriFaceSplitLevel(mid[0], mid[1], mid[2]));
5056  }
5057  else // not split
5058  {
5059  return 0;
5060  }
5061 }
5062 
5063 void NCMesh::QuadFaceSplitLevel(int vn1, int vn2, int vn3, int vn4,
5064  int& h_level, int& v_level) const
5065 {
5066  int hl1, hl2, vl1, vl2;
5067  int mid[5];
5068 
5069  switch (QuadFaceSplitType(vn1, vn2, vn3, vn4, mid))
5070  {
5071  case 0: // not split
5072  h_level = v_level = 0;
5073  break;
5074 
5075  case 1: // vertical
5076  QuadFaceSplitLevel(vn1, mid[0], mid[2], vn4, hl1, vl1);
5077  QuadFaceSplitLevel(mid[0], vn2, vn3, mid[2], hl2, vl2);
5078  h_level = std::max(hl1, hl2);
5079  v_level = std::max(vl1, vl2) + 1;
5080  break;
5081 
5082  default: // horizontal
5083  QuadFaceSplitLevel(vn1, vn2, mid[1], mid[3], hl1, vl1);
5084  QuadFaceSplitLevel(mid[3], mid[1], vn3, vn4, hl2, vl2);
5085  h_level = std::max(hl1, hl2) + 1;
5086  v_level = std::max(vl1, vl2);
5087  }
5088 }
5089 
5090 void NCMesh::CountSplits(int elem, int splits[3]) const
5091 {
5092  const Element &el = elements[elem];
5093  const int* node = el.node;
5094  GeomInfo& gi = GI[el.Geom()];
5095 
5096  int elevel[12];
5097  for (int i = 0; i < gi.ne; i++)
5098  {
5099  const int* ev = gi.edges[i];
5100  elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
5101  }
5102 
5103  int flevel[6][2];
5104  if (Dim >= 3)
5105  {
5106  for (int i = 0; i < gi.nf; i++)
5107  {
5108  const int* fv = gi.faces[i];
5109  if (gi.nfv[i] == 4)
5110  {
5111  QuadFaceSplitLevel(node[fv[0]], node[fv[1]],
5112  node[fv[2]], node[fv[3]],
5113  flevel[i][1], flevel[i][0]);
5114  }
5115  else
5116  {
5117  flevel[i][1] = 0;
5118  flevel[i][0] =
5119  TriFaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]]);
5120  }
5121  }
5122  }
5123 
5124  if (el.Geom() == Geometry::CUBE)
5125  {
5126  splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
5127  elevel[0], elevel[2], elevel[4], elevel[6]);
5128 
5129  splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
5130  elevel[1], elevel[3], elevel[5], elevel[7]);
5131 
5132  splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
5133  elevel[8], elevel[9], elevel[10], elevel[11]);
5134  }
5135  else if (el.Geom() == Geometry::PRISM)
5136  {
5137  splits[0] = splits[1] =
5138  std::max(
5139  max6(flevel[0][0], flevel[1][0], 0,
5140  flevel[2][0], flevel[3][0], flevel[4][0]),
5141  max6(elevel[0], elevel[1], elevel[2],
5142  elevel[3], elevel[4], elevel[5]));
5143 
5144  splits[2] = max6(flevel[2][1], flevel[3][1], flevel[4][1],
5145  elevel[6], elevel[7], elevel[8]);
5146  }
5147  else if (el.Geom() == Geometry::TETRAHEDRON)
5148  {
5149  splits[0] = std::max(
5150  max4(flevel[0][0], flevel[1][0], flevel[2][0], flevel[3][0]),
5151  max6(elevel[0], elevel[1], elevel[2],
5152  elevel[3], elevel[4], elevel[5]));
5153 
5154  splits[1] = splits[0];
5155  splits[2] = splits[0];
5156  }
5157  else if (el.Geom() == Geometry::SQUARE)
5158  {
5159  splits[0] = std::max(elevel[0], elevel[2]);
5160  splits[1] = std::max(elevel[1], elevel[3]);
5161  }
5162  else if (el.Geom() == Geometry::TRIANGLE)
5163  {
5164  splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
5165  splits[1] = splits[0];
5166  }
5167  else
5168  {
5169  MFEM_ABORT("Unsupported element geometry.");
5170  }
5171 }
5172 
5173 void NCMesh::GetLimitRefinements(Array<Refinement> &refinements, int max_level)
5174 {
5175  for (int i = 0; i < leaf_elements.Size(); i++)
5176  {
5177  if (IsGhost(elements[leaf_elements[i]])) { break; } // TODO: NElements
5178 
5179  int splits[3];
5180  CountSplits(leaf_elements[i], splits);
5181 
5182  char ref_type = 0;
5183  for (int k = 0; k < Dim; k++)
5184  {
5185  if (splits[k] > max_level)
5186  {
5187  ref_type |= (1 << k);
5188  }
5189  }
5190 
5191  if (ref_type)
5192  {
5193  if (Iso)
5194  {
5195  // iso meshes should only be modified by iso refinements
5196  ref_type = 7;
5197  }
5198  refinements.Append(Refinement(i, ref_type));
5199  }
5200  }
5201 }
5202 
5203 void NCMesh::LimitNCLevel(int max_nc_level)
5204 {
5205  MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
5206 
5207  while (1)
5208  {
5209  Array<Refinement> refinements;
5210  GetLimitRefinements(refinements, max_nc_level);
5211 
5212  if (!refinements.Size()) { break; }
5213 
5214  Refine(refinements);
5215  }
5216 }
5217 
5218 
5219 //// I/O ////////////////////////////////////////////////////////////////////////
5220 
5221 int NCMesh::PrintVertexParents(std::ostream *os) const
5222 {
5223  if (!os)
5224  {
5225  // count vertex nodes with parents
5226  int nv = 0;
5227  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5228  {
5229  if (node->HasVertex() && node->p1 != node->p2) { nv++; }
5230  }
5231  return nv;
5232  }
5233  else
5234  {
5235  // print the relations
5236  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5237  {
5238  if (node->HasVertex() && node->p1 != node->p2)
5239  {
5240  MFEM_ASSERT(nodes[node->p1].HasVertex(), "");
5241  MFEM_ASSERT(nodes[node->p2].HasVertex(), "");
5242 
5243  (*os) << node.index() << " " << node->p1 << " " << node->p2 << "\n";
5244  }
5245  }
5246  return 0;
5247  }
5248 }
5249 
5250 void NCMesh::LoadVertexParents(std::istream &input)
5251 {
5252  int nv;
5253  input >> nv;
5254  while (nv--)
5255  {
5256  int id, p1, p2;
5257  input >> id >> p1 >> p2;
5258  MFEM_VERIFY(input, "problem reading vertex parents.");
5259 
5260  MFEM_VERIFY(nodes.IdExists(id), "vertex " << id << " not found.");
5261  MFEM_VERIFY(nodes.IdExists(p1), "parent " << p1 << " not found.");
5262  MFEM_VERIFY(nodes.IdExists(p2), "parent " << p2 << " not found.");
5263 
5264  int check = nodes.FindId(p1, p2);
5265  MFEM_VERIFY(check < 0, "parents (" << p1 << ", " << p2 << ") already "
5266  "assigned to node " << check);
5267 
5268  // assign new parents for the node
5269  nodes.Reparent(id, p1, p2);
5270  }
5271 }
5272 
5273 int NCMesh::PrintBoundary(std::ostream *os) const
5274 {
5275  static const int nfv2geom[5] =
5276  {
5279  };
5280  int deg = (Dim == 2) ? 2 : 1; // for degenerate faces in 2D
5281 
5282  int count = 0;
5283  for (int i = 0; i < elements.Size(); i++)
5284  {
5285  const Element &el = elements[i];
5286  if (!el.IsLeaf()) { continue; }
5287 
5288  GeomInfo& gi = GI[el.Geom()];
5289  for (int k = 0; k < gi.nf; k++)
5290  {
5291  const int* fv = gi.faces[k];
5292  const int nfv = gi.nfv[k];
5293  const Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
5294  el.node[fv[2]], el.node[fv[3]]);
5295  MFEM_ASSERT(face != NULL, "face not found");
5296  if (face->Boundary())
5297  {
5298  if (!os) { count++; continue; }
5299 
5300  (*os) << face->attribute << " " << nfv2geom[nfv];
5301  for (int j = 0; j < nfv; j++)
5302  {
5303  (*os) << " " << el.node[fv[j*deg]];
5304  }
5305  (*os) << "\n";
5306  }
5307  }
5308  }
5309  return count;
5310 }
5311 
5312 void NCMesh::LoadBoundary(std::istream &input)
5313 {
5314  int nb, attr, geom;
5315  input >> nb;
5316  for (int i = 0; i < nb; i++)
5317  {
5318  input >> attr >> geom;
5319 
5320  int v1, v2, v3, v4;
5321  if (geom == Geometry::SQUARE)
5322  {
5323  input >> v1 >> v2 >> v3 >> v4;
5324  Face* face = faces.Get(v1, v2, v3, v4);
5325  face->attribute = attr;
5326  }
5327  else if (geom == Geometry::TRIANGLE)
5328  {
5329  input >> v1 >> v2 >> v3;
5330  Face* face = faces.Get(v1, v2, v3);
5331  face->attribute = attr;
5332  }
5333  else if (geom == Geometry::SEGMENT)
5334  {
5335  input >> v1 >> v2;
5336  Face* face = faces.Get(v1, v1, v2, v2);
5337  face->attribute = attr;
5338  }
5339  else if (geom == Geometry::POINT)
5340  {
5341  input >> v1;
5342  Face* face = faces.Get(v1, v1, v1, v1);
5343  face->attribute = attr;
5344  }
5345  else
5346  {
5347  MFEM_ABORT("unsupported boundary element geometry: " << geom);
5348  }
5349  }
5350 }
5351 
5352 void NCMesh::PrintCoordinates(std::ostream &os) const
5353 {
5354  int nv = coordinates.Size()/3;
5355  os << nv << "\n";
5356  if (!nv) { return; }
5357 
5358  os << spaceDim << "\n";
5359  for (int i = 0; i < nv; i++)
5360  {
5361  os << coordinates[3*i];
5362  for (int j = 1; j < spaceDim; j++)
5363  {
5364  os << " " << coordinates[3*i + j];
5365  }
5366  os << "\n";
5367  }
5368 }
5369 
5370 void NCMesh::LoadCoordinates(std::istream &input)
5371 {
5372  int nv;
5373  input >> nv;
5374  if (!nv) { return; }
5375 
5376  input >> spaceDim;
5377 
5378  coordinates.SetSize(nv * 3);
5379  coordinates = 0.0;
5380 
5381  for (int i = 0; i < nv; i++)
5382  {
5383  for (int j = 0; j < spaceDim; j++)
5384  {
5385  input >> coordinates[3*i + j];
5386  MFEM_VERIFY(input.good(), "unexpected EOF");
5387  }
5388  }
5389 }
5390 
5392 {
5393  for (int i = 0; i < root_state.Size(); i++)
5394  {
5395  if (root_state[i]) { return false; }
5396  }
5397  return true;
5398 }
5399 
5400 void NCMesh::Print(std::ostream &os) const
5401 {
5402  os << "MFEM NC mesh v1.0\n\n"
5403  "# NCMesh supported geometry types:\n"
5404  "# SEGMENT = 1\n"
5405  "# TRIANGLE = 2\n"
5406  "# SQUARE = 3\n"
5407  "# TETRAHEDRON = 4\n"
5408  "# CUBE = 5\n"
5409  "# PRISM = 6\n";
5410 
5411  os << "\ndimension\n" << Dim << "\n";
5412 
5413 #ifndef MFEM_USE_MPI
5414  if (MyRank != 0) // don't print this section in serial: default rank is 0
5415 #endif
5416  {
5417  os << "\nrank\n" << MyRank << "\n";
5418  }
5419 
5420  os << "\n# rank attr geom ref_type nodes/children";
5421  os << "\nelements\n" << elements.Size() << "\n";
5422 
5423  for (int i = 0; i < elements.Size(); i++)
5424  {
5425  const Element &el = elements[i];
5426  os << el.rank << " " << el.attribute << " ";
5427  if (el.parent == -2) { os << "-1\n"; continue; } // unused element
5428 
5429  os << int(el.geom) << " " << int(el.ref_type);
5430  for (int j = 0; j < 8 && el.node[j] >= 0; j++)
5431  {
5432  os << " " << el.node[j];
5433  }
5434  os << "\n";
5435  }
5436 
5437  int nb = PrintBoundary(NULL);
5438  if (nb)
5439  {
5440  os << "\n# attr geom nodes";
5441  os << "\nboundary\n" << nb << "\n";
5442 
5443  PrintBoundary(&os);
5444  }
5445 
5446  int nvp = PrintVertexParents(NULL);
5447  if (nvp)
5448  {
5449  os << "\n# vert_id p1 p2";
5450  os << "\nvertex_parents\n" << nvp << "\n";
5451 
5452  PrintVertexParents(&os);
5453  }
5454 
5455  if (!ZeroRootStates()) // root_state section is optional
5456  {
5457  os << "\n# root element orientation";
5458  os << "\nroot_state\n" << root_state.Size() << "\n";
5459 
5460  for (int i = 0; i < root_state.Size(); i++)
5461  {
5462  os << root_state[i] << "\n";
5463  }
5464  }
5465 
5466  if (coordinates.Size())
5467  {
5468  os << "\n# top-level node coordinates";
5469  os << "\ncoordinates\n";
5470 
5471  PrintCoordinates(os);
5472  }
5473  else
5474  {
5475  // 'nodes' will be printed one level up by Mesh::Printer()
5476  }
5477 }
5478 
5480 {
5481  // initialize Element::parent
5482  for (int i = 0; i < elements.Size(); i++)
5483  {
5484  Element &el = elements[i];
5485  if (el.ref_type)
5486  {
5487  for (int j = 0; j < 8 && el.child[j] >= 0; j++)
5488  {
5489  int child = el.child[j];
5490  MFEM_VERIFY(child < elements.Size(), "invalid mesh file: "
5491  "element " << i << " references invalid child " << child);
5492  elements[child].parent = i;
5493  }
5494  }
5495  }
5496 
5497  // count the root elements
5498  int nroots = 0;
5499  while (nroots < elements.Size() &&
5500  elements[nroots].parent == -1)
5501  {
5502  nroots++;
5503  }
5504  MFEM_VERIFY(nroots, "invalid mesh file: no root elements found.");
5505 
5506  // check that only the first 'nroot' elements are roots (have no parent)
5507  for (int i = nroots; i < elements.Size(); i++)
5508  {
5509  MFEM_VERIFY(elements[i].parent != -1,
5510  "invalid mesh file: only the first M elements can be roots. "
5511  "Found element " << i << " with no parent.");
5512  }
5513 
5514  // default root state is zero
5515  root_state.SetSize(nroots);
5516  root_state = 0;
5517 }
5518 
5520 {
5521  int ntop = 0;
5522  for (auto node = nodes.cbegin(); node != nodes.cend(); ++node)
5523  {
5524  if (node->p1 == node->p2) { ntop = node.index() + 1; }
5525  }
5526  return ntop;
5527 }
5528 
5529 NCMesh::NCMesh(std::istream &input, int version, int &curved, int &is_nc)
5530  : spaceDim(0), MyRank(0), Iso(true), Legacy(false)
5531 {
5532  is_nc = 1;
5533  if (version == 1) // old MFEM mesh v1.1 format
5534  {
5535  LoadLegacyFormat(input, curved, is_nc);
5536  Legacy = true;
5537  return;
5538  }
5539 
5540  MFEM_ASSERT(version == 10, "");
5541  std::string ident;
5542  int count;
5543 
5544  // load dimension
5545  skip_comment_lines(input, '#');
5546  input >> ident;
5547  MFEM_VERIFY(ident == "dimension", "Invalid mesh file: " << ident);
5548  input >> Dim;
5549 
5550  // load rank, if present
5551  skip_comment_lines(input, '#');
5552  input >> ident;
5553  if (ident == "rank")
5554  {
5555  input >> MyRank;
5556  MFEM_VERIFY(MyRank >= 0, "Invalid rank");
5557 
5558  skip_comment_lines(input, '#');
5559  input >> ident;
5560  }
5561 
5562  // load file SFC version, if present (for future changes to SFC ordering)
5563  if (ident == "sfc_version")
5564  {
5565  int sfc_version; // TODO future: store as class member
5566  input >> sfc_version;
5567  MFEM_VERIFY(sfc_version == 0,
5568  "Unsupported mesh file SFC version (" << sfc_version << "). "
5569  "Please update MFEM.");
5570 
5571  skip_comment_lines(input, '#');
5572  input >> ident;
5573  }
5574 
5575  // load elements
5576  MFEM_VERIFY(ident == "elements", "Invalid mesh file: " << ident);
5577  input >> count;
5578  for (int i = 0; i < count; i++)
5579  {
5580  int rank, attr, geom, ref_type;
5581  input >> rank >> attr >> geom;
5582 
5583  Geometry::Type type = Geometry::Type(geom);
5584  elements.Append(Element(type, attr));
5585 
5586  MFEM_ASSERT(elements.Size() == i+1, "");
5587  Element &el = elements[i];
5588  el.rank = rank;
5589 
5590  if (geom >= 0)
5591  {
5592  CheckSupportedGeom(type);
5593  GI[geom].InitGeom(type);
5594 
5595  input >> ref_type;
5596  MFEM_VERIFY(ref_type >= 0 && ref_type < 8, "");
5597  el.ref_type = ref_type;
5598 
5599  if (ref_type) // refined element
5600  {
5601  for (int j = 0; j < ref_type_num_children[ref_type]; j++)
5602  {
5603  input >> el.child[j];
5604  }
5605  if (Dim == 3 && ref_type != 7) { Iso = false; }
5606  }
5607  else // leaf element
5608  {
5609  for (int j = 0; j < GI[geom].nv; j++)
5610  {
5611  int id;
5612  input >> id;
5613  el.node[j] = id;
5614  nodes.Alloc(id, id, id);
5615  // NOTE: nodes that won't get parents assigned will
5616  // stay hashed with p1 == p2 == id (top-level nodes)
5617  }
5618  }
5619  }
5620  else
5621  {
5622  el.parent = -2; // mark as unused
5624  }
5625  }
5626 
5627  InitRootElements();
5628  InitGeomFlags();
5629 
5630  skip_comment_lines(input, '#');
5631  input >> ident;
5632 
5633  // load boundary
5634  if (ident == "boundary")
5635  {
5636  LoadBoundary(input);
5637 
5638  skip_comment_lines(input, '#');
5639  input >> ident;
5640  }
5641 
5642  // load vertex hierarchy
5643  if (ident == "vertex_parents")
5644  {
5645  LoadVertexParents(input);
5646 
5647  skip_comment_lines(input, '#');
5648  input >> ident;
5649  }
5650 
5651  // load root states
5652  if (ident == "root_state")
5653  {
5654  input >> count;
5655  MFEM_VERIFY(count <= root_state.Size(), "Too many root states");
5656  for (int i = 0; i < count; i++)
5657  {
5658  input >> root_state[i];
5659  }
5660 
5661  skip_comment_lines(input, '#');
5662  input >> ident;
5663  }
5664 
5665  // load coordinates or nodes
5666  if (ident == "coordinates")
5667  {
5668  LoadCoordinates(input);
5669 
5670  MFEM_VERIFY(coordinates.Size()/3 >= CountTopLevelNodes(),
5671  "Invalid mesh file: not all top-level nodes are covered by "
5672  "the 'coordinates' section of the mesh file.");
5673  curved = 0;
5674  }
5675  else if (ident == "nodes")
5676  {
5677  coordinates.SetSize(0); // this means the mesh is curved
5678 
5679  // prepare to read the nodes
5680  input >> std::ws;
5681  curved = 1;
5682  }
5683  else
5684  {
5685  MFEM_ABORT("Invalid mesh file: either 'coordinates' or "
5686  "'nodes' must be present");
5687  }
5688 
5689  // create edge nodes and faces
5690  nodes.UpdateUnused();
5691  for (int i = 0; i < elements.Size(); i++)
5692  {
5693  if (elements[i].IsLeaf())
5694  {
5695  ReferenceElement(i);
5696  RegisterFaces(i);
5697  }
5698  }
5699 
5700  Update();
5701 }
5702 
5703 void NCMesh::CopyElements(int elem,
5704  const BlockArray<Element> &tmp_elements)
5705 {
5706  Element &el = elements[elem];
5707  if (el.ref_type)
5708  {
5709  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
5710  {
5711  int old_id = el.child[i];
5712  // here we know 'free_element_ids' is empty
5713  int new_id = elements.Append(tmp_elements[old_id]);
5714  el.child[i] = new_id;
5715  elements[new_id].parent = elem;
5716  CopyElements(new_id, tmp_elements);
5717  }
5718  }
5719 }
5720 
5721 void NCMesh::LoadCoarseElements(std::istream &input)
5722 {
5723  int ne;
5724  input >> ne;
5725 
5726  bool iso = true;
5727 
5728  // load the coarse elements
5729  while (ne--)
5730  {
5731  int ref_type;
5732  input >> ref_type;
5733 
5734  int elem = AddElement(Element(Geometry::INVALID, 0));
5735  Element &el = elements[elem];
5736  el.ref_type = ref_type;
5737 
5738  if (Dim == 3 && ref_type != 7) { iso = false; }
5739 
5740  // load child IDs and make parent-child links
5741  int nch = ref_type_num_children[ref_type];
5742  for (int i = 0, id; i < nch; i++)
5743  {
5744  input >> id;
5745  MFEM_VERIFY(id >= 0, "");
5746  MFEM_VERIFY(id < elements.Size(),
5747  "coarse element cannot be referenced before it is "
5748  "defined (id=" << id << ").");
5749 
5750  Element &child = elements[id];
5751  MFEM_VERIFY(child.parent == -1,
5752  "element " << id << " cannot have two parents.");
5753 
5754  el.child[i] = id;
5755  child.parent = elem;
5756 
5757  if (!i) // copy geom and attribute from first child
5758  {
5759  el.geom = child.geom;
5760  el.attribute = child.attribute;
5761  }
5762  }
5763  }
5764 
5765  // prepare for reordering the elements
5766  BlockArray<Element> tmp_elements;
5767  elements.Swap(tmp_elements);
5768 
5769  // copy roots, they need to be at the beginning of 'elements'
5770  int root_count = 0;
5771  for (auto el = tmp_elements.begin(); el != tmp_elements.end(); ++el)
5772  {
5773  if (el->parent == -1)
5774  {
5775  elements.Append(*el); // same as AddElement()
5776  root_count++;
5777  }
5778  }
5779 
5780  // copy the rest of the hierarchy
5781  for (int i = 0; i < root_count; i++)
5782  {
5783  CopyElements(i, tmp_elements);
5784  }
5785 
5786  // set the Iso flag (must be false if there are 3D aniso refinements)
5787  Iso = iso;
5788 
5789  InitRootState(root_count);
5790 }
5791 
5792 void NCMesh::LoadLegacyFormat(std::istream &input, int &curved, int &is_nc)
5793 {
5794  MFEM_ASSERT(elements.Size() == 0, "");
5795  MFEM_ASSERT(nodes.Size() == 0, "");
5796  MFEM_ASSERT(free_element_ids.Size() == 0, "");
5797 
5798  std::string ident;
5799  int count, attr, geom;
5800 
5801  // load dimension
5802  skip_comment_lines(input, '#');
5803  input >> ident;
5804  MFEM_VERIFY(ident == "dimension", "invalid mesh file");
5805  input >> Dim;
5806 
5807  // load elements
5808  skip_comment_lines(input, '#');
5809  input >> ident;
5810  MFEM_VERIFY(ident == "elements", "invalid mesh file");
5811 
5812  input >> count;
5813  for (int i = 0; i < count; i++)
5814  {
5815  input >> attr >> geom;
5816 
5817  Geometry::Type type = Geometry::Type(geom);
5818  CheckSupportedGeom(type);
5819  GI[geom].InitGeom(type);
5820 
5821  int eid = AddElement(Element(type, attr));
5822  MFEM_ASSERT(eid == i, "");
5823 
5824  Element &el = elements[eid];
5825  for (int j = 0; j < GI[geom].nv; j++)
5826  {
5827  int id;
5828  input >> id;
5829  el.node[j] = id;
5830  nodes.Alloc(id, id, id); // see comment in NCMesh::NCMesh
5831  }
5832  el.index = i; // needed for file leaf order below
5833  }
5834 
5835  // load boundary
5836  skip_comment_lines(input, '#');
5837  input >> ident;
5838  MFEM_VERIFY(ident == "boundary", "invalid mesh file");
5839 
5840  LoadBoundary(input);
5841 
5842  // load vertex hierarchy
5843  skip_comment_lines(input, '#');
5844  input >> ident;
5845  if (ident == "vertex_parents")
5846  {
5847  LoadVertexParents(input);
5848  is_nc = 1;
5849 
5850  skip_comment_lines(input, '#');
5851  input >> ident;
5852  }
5853  else
5854  {
5855  // no "vertex_parents" section: this file needs to be treated as a
5856  // conforming mesh for complete backward compatibility with MFEM 4.2
5857  is_nc = 0;
5858  }
5859 
5860  // load element hierarchy
5861  if (ident == "coarse_elements")
5862  {
5863  LoadCoarseElements(input);
5864 
5865  skip_comment_lines(input, '#');
5866  input >> ident;
5867  }
5868  else
5869  {
5870  // no element hierarchy -> all elements are roots
5871  InitRootState(elements.Size());
5872  }
5873  InitGeomFlags();
5874 
5875  // load vertices
5876  MFEM_VERIFY(ident == "vertices", "invalid mesh file");
5877  int nvert;
5878  input >> nvert;
5879  input >> std::ws >> ident;
5880  if (ident != "nodes")
5881  {