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