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