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