MFEM  v3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ncmesh.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../fem/fem.hpp"
13 #include "ncmesh.hpp"
14 
15 #include <string>
16 #include <algorithm>
17 #include <cmath>
18 #include <climits> // INT_MAX
19 
20 namespace mfem
21 {
22 
23 NCMesh::GeomInfo NCMesh::GI[Geometry::NumGeom];
24 
25 static NCMesh::GeomInfo& gi_hex = NCMesh::GI[Geometry::CUBE];
26 static NCMesh::GeomInfo& gi_quad = NCMesh::GI[Geometry::SQUARE];
27 static NCMesh::GeomInfo& gi_tri = NCMesh::GI[Geometry::TRIANGLE];
28 
30 {
31  if (initialized) { return; }
32 
33  nv = elem->GetNVertices();
34  ne = elem->GetNEdges();
35  nf = elem->GetNFaces(nfv);
36 
37  for (int i = 0; i < ne; i++)
38  {
39  for (int j = 0; j < 2; j++)
40  {
41  edges[i][j] = elem->GetEdgeVertices(i)[j];
42  }
43  }
44  for (int i = 0; i < nf; i++)
45  {
46  for (int j = 0; j < nfv; j++)
47  {
48  faces[i][j] = elem->GetFaceVertices(i)[j];
49  }
50  }
51 
52  // in 2D we pretend to have faces too, so we can use Face::elem[2]
53  if (!nf)
54  {
55  for (int i = 0; i < ne; i++)
56  {
57  // make a degenerate face
58  faces[i][0] = faces[i][1] = edges[i][0];
59  faces[i][2] = faces[i][3] = edges[i][1];
60  }
61  nf = ne;
62  }
63 
64  initialized = true;
65 }
66 
67 
68 NCMesh::NCMesh(const Mesh *mesh, std::istream *vertex_parents)
69 {
70  Dim = mesh->Dimension();
71  spaceDim = mesh->SpaceDimension();
72 
73  // assume the mesh is anisotropic if we're loading a file
74  Iso = vertex_parents ? false : true;
75 
76  // examine elements and reserve the first node IDs for vertices
77  // (note: 'mesh' may not have vertices defined yet, e.g., on load)
78  int max_id = -1;
79  for (int i = 0; i < mesh->GetNE(); i++)
80  {
81  const mfem::Element *elem = mesh->GetElement(i);
82  const int *v = elem->GetVertices(), nv = elem->GetNVertices();
83  for (int j = 0; j < nv; j++)
84  {
85  max_id = std::max(max_id, v[j]);
86  }
87  }
88  for (int i = 0; i <= max_id; i++)
89  {
90  // top-level nodes are special: id == p1 == p2 == orig. vertex id
91  Node* node = nodes.Get(i, i);
92  MFEM_CONTRACT_VAR(node);
93  MFEM_ASSERT(node->id == i, "");
94  }
95 
96  // if a mesh file is being read, load the vertex hierarchy now;
97  // 'vertex_parents' must be at the appropriate section in the mesh file
98  if (vertex_parents)
99  {
100  LoadVertexParents(*vertex_parents);
101  }
102 
103  // create the NCMesh::Element struct for each Mesh element
104  for (int i = 0; i < mesh->GetNE(); i++)
105  {
106  const mfem::Element *elem = mesh->GetElement(i);
107  const int *v = elem->GetVertices();
108 
109  int geom = elem->GetGeometryType();
110  if (geom != Geometry::TRIANGLE &&
111  geom != Geometry::SQUARE &&
112  geom != Geometry::CUBE)
113  {
114  MFEM_ABORT("only triangles, quads and hexes are supported by NCMesh.");
115  }
116 
117  // initialize edge/face tables for this type of element
118  GI[geom].Initialize(elem);
119 
120  // create our Element struct for this element
121  Element* nc_elem = new Element(geom, elem->GetAttribute());
122  root_elements.Append(nc_elem);
123 
124  for (int j = 0; j < GI[geom].nv; j++)
125  {
126  Node* node = nodes.Peek(v[j]);
127  if (!node->vertex)
128  {
129  if (v[j] < mesh->GetNV())
130  {
131  // create a vertex in the node and initialize its position
132  const double* pos = mesh->GetVertex(v[j]);
133  node->vertex = new Vertex(pos[0], pos[1], pos[2]);
134  }
135  else
136  {
137  // the mesh may not have vertex positions defined yet
138  node->vertex = new Vertex(0.0, 0.0, 0.0);
139  }
140  }
141  nc_elem->node[j] = node;
142  }
143 
144  // increase reference count of all nodes the element is using
145  // (NOTE: this will also create and reference all edge and face nodes)
146  RefElementNodes(nc_elem);
147 
148  // make links from faces back to the element
149  RegisterFaces(nc_elem);
150  }
151 
152  // store boundary element attributes
153  for (int i = 0; i < mesh->GetNBE(); i++)
154  {
155  const mfem::Element *be = mesh->GetBdrElement(i);
156  const int *v = be->GetVertices();
157 
158  Node* node[4];
159  for (int i = 0; i < be->GetNVertices(); i++)
160  {
161  node[i] = nodes.Peek(v[i]);
162  MFEM_VERIFY(node[i], "boundary elements inconsistent.");
163  }
164 
166  {
167  Face* face = faces.Peek(node[0], node[1], node[2], node[3]);
168  MFEM_VERIFY(face, "boundary face not found.");
169  face->attribute = be->GetAttribute();
170  }
171  else if (be->GetType() == mfem::Element::SEGMENT)
172  {
173  Edge* edge = nodes.Peek(node[0], node[1])->edge;
174  MFEM_VERIFY(edge, "boundary edge not found.");
175  edge->attribute = be->GetAttribute();
176  }
177  else
178  {
179  MFEM_ABORT("only segment and quadrilateral boundary "
180  "elements are supported by NCMesh.");
181  }
182  }
183 
184  Update();
185 }
186 
188 {
190  UpdateVertices();
191 
192  face_list.Clear();
193  edge_list.Clear();
194 
196 }
197 
199 {
200  Element* new_elem = new Element(*elem);
201  if (elem->ref_type)
202  {
203  for (int i = 0; i < 8; i++)
204  {
205  if (elem->child[i])
206  {
207  new_elem->child[i] = CopyHierarchy(elem->child[i]);
208  new_elem->child[i]->parent = new_elem;
209  }
210  }
211  }
212  else
213  {
214  GeomInfo& gi = GI[(int) elem->geom];
215  for (int i = 0; i < gi.nv; i++)
216  {
217  new_elem->node[i] = nodes.Peek(elem->node[i]->id);
218  }
219  RegisterFaces(new_elem);
220  }
221  return new_elem;
222 }
223 
225 {
226  if (elem->ref_type)
227  {
228  for (int i = 0; i < 8; i++)
229  {
230  if (elem->child[i]) { DeleteHierarchy(elem->child[i]); }
231  }
232  }
233  else
234  {
235  UnrefElementNodes(elem);
236  }
237  delete elem;
238 }
239 
240 NCMesh::NCMesh(const NCMesh &other)
241  : Dim(other.Dim), spaceDim(other.spaceDim), Iso(other.Iso)
242  , nodes(other.nodes), faces(other.faces)
243 {
244  // NOTE: this copy constructor is used by ParNCMesh
245  root_elements.SetSize(other.root_elements.Size());
246  for (int i = 0; i < root_elements.Size(); i++)
247  {
249  }
250 
251  Update();
252 }
253 
255 {
256  for (int i = 0; i < root_elements.Size(); i++)
257  {
259  }
260 }
261 
262 
264 
266 {
267  MFEM_ASSERT(vertex, "can't create vertex here.");
268  vertex->Ref();
269 }
270 
272 {
273  if (!edge) { edge = new Edge; }
274  edge->Ref();
275 }
276 
278 {
279  MFEM_ASSERT(vertex, "cannot unref a nonexistent vertex.");
280  if (!vertex->Unref()) { vertex = NULL; }
281  if (!vertex && !edge) { nodes.Delete(this); }
282 }
283 
285 {
286  MFEM_ASSERT(node, "node not found.");
287  MFEM_ASSERT(node->edge, "cannot unref a nonexistent edge.");
288  if (!node->edge->Unref()) { node->edge = NULL; }
289  if (!node->vertex && !node->edge) { nodes.Delete(node); }
290 }
291 
293 {
294  std::memcpy(this, &other, sizeof(*this));
295  if (vertex) { vertex = new Vertex(*vertex); }
296  if (edge) { edge = new Edge(*edge); }
297 }
298 
300 {
301  MFEM_ASSERT(!vertex && !edge, "node was not unreffed properly.");
302  if (vertex) { delete vertex; }
303  if (edge) { delete edge; }
304 }
305 
307 {
308  Node** node = elem->node;
309  GeomInfo& gi = GI[(int) elem->geom];
310 
311  // ref all vertices
312  for (int i = 0; i < gi.nv; i++)
313  {
314  node[i]->RefVertex();
315  }
316 
317  // ref all edges (possibly creating them)
318  for (int i = 0; i < gi.ne; i++)
319  {
320  const int* ev = gi.edges[i];
321  nodes.Get(node[ev[0]], node[ev[1]])->RefEdge();
322  }
323 
324  // ref all faces (possibly creating them)
325  for (int i = 0; i < gi.nf; i++)
326  {
327  const int* fv = gi.faces[i];
328  faces.Get(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]])->Ref();
329  // NOTE: face->RegisterElement called elsewhere to avoid having
330  // to store 3 element pointers temporarily in the face when refining.
331  // See also NCMesh::RegisterFaces.
332  }
333 }
334 
336 {
337  Node** node = elem->node;
338  GeomInfo& gi = GI[(int) elem->geom];
339 
340  // unref all faces (possibly destroying them)
341  for (int i = 0; i < gi.nf; i++)
342  {
343  const int* fv = gi.faces[i];
344  Face* face = faces.Peek(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
345  face->ForgetElement(elem);
346  if (!face->Unref()) { faces.Delete(face); }
347  }
348 
349  // unref all edges (possibly destroying them)
350  for (int i = 0; i < gi.ne; i++)
351  {
352  const int* ev = gi.edges[i];
353  //nodes.Peek(node[ev[0]], node[ev[1]])->UnrefEdge(nodes); -- pre-aniso
354  Node::UnrefEdge(PeekAltParents(node[ev[0]], node[ev[1]]), nodes);
355  }
356 
357  // unref all vertices (possibly destroying them)
358  for (int i = 0; i < gi.nv; i++)
359  {
360  elem->node[i]->UnrefVertex(nodes);
361  }
362 }
363 
365  : Hashed4<Face>(id), attribute(-1), index(-1)
366 {
367  elem[0] = elem[1] = NULL;
368 }
369 
371 {
372  // (skip Hashed4 constructor)
373  std::memcpy(this, &other, sizeof(*this));
374  elem[0] = elem[1] = NULL;
375 }
376 
378 {
379  if (elem[0] == NULL) { elem[0] = e; }
380  else if (elem[1] == NULL) { elem[1] = e; }
381  else { MFEM_ABORT("can't have 3 elements in Face::elem[]."); }
382 }
383 
385 {
386  if (elem[0] == e) { elem[0] = NULL; }
387  else if (elem[1] == e) { elem[1] = NULL; }
388  else { MFEM_ABORT("element not found in Face::elem[]."); }
389 }
390 
392 {
393  GeomInfo& gi = GI[(int) elem->geom];
394  const int* fv = gi.faces[face_no];
395  Node** node = elem->node;
396  return faces.Peek(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
397 }
398 
399 void NCMesh::RegisterFaces(Element* elem, int* fattr)
400 {
401  GeomInfo& gi = GI[(int) elem->geom];
402  for (int i = 0; i < gi.nf; i++)
403  {
404  Face* face = GetFace(elem, i);
405  face->RegisterElement(elem);
406  if (fattr) { face->attribute = fattr[i]; }
407  }
408 }
409 
411 {
412  if (elem[0])
413  {
414  MFEM_ASSERT(!elem[1], "not a single element face.");
415  return elem[0];
416  }
417  else
418  {
419  MFEM_ASSERT(elem[1], "no elements in face.");
420  return elem[1];
421  }
422 }
423 
425 {
426  Node* mid = nodes.Peek(v1, v2);
427  if (!mid) // TODO: && !Iso ?
428  {
429  // In rare cases, a mid-face node exists under alternate parents w1, w2
430  // (see picture) instead of the requested parents v1, v2. This is an
431  // inconsistent situation that may exist temporarily as a result of
432  // "nodes.Reparent" while doing anisotropic splits, before forced
433  // refinements are all processed. This function attempts to retrieve such
434  // a node. An extra twist is that w1 and w2 may themselves need to be
435  // obtained using this very function.
436  //
437  // v1->p1 v1 v1->p2
438  // *------*------*
439  // | | |
440  // | |mid |
441  // w1 *------*------* w2
442  // | | |
443  // | | |
444  // *------*------*
445  // v2->p1 v2 v2->p2
446  //
447  // NOTE: this function would not be needed if the elements remembered
448  // pointers to their edge nodes. We have however opted to save memory
449  // at the cost of this computation, which is only necessary when forced
450  // refinements are being done.
451 
452  if ((v1->p1 != v1->p2) && (v2->p1 != v2->p2)) // non-top-level nodes?
453  {
454  Node *v1p1 = nodes.Peek(v1->p1), *v1p2 = nodes.Peek(v1->p2);
455  Node *v2p1 = nodes.Peek(v2->p1), *v2p2 = nodes.Peek(v2->p2);
456 
457  Node* w1 = PeekAltParents(v1p1, v2p1);
458  Node* w2 = w1 ? PeekAltParents(v1p2, v2p2) : NULL /* optimization */;
459 
460  if (!w1 || !w2) // one more try may be needed as p1, p2 are unordered
461  {
462  w1 = PeekAltParents(v1p1, v2p2);
463  w2 = w1 ? PeekAltParents(v1p2, v2p1) : NULL /* optimization */;
464  }
465 
466  if (w1 && w2) // got both alternate parents?
467  {
468  mid = nodes.Peek(w1, w2);
469  }
470  }
471  }
472  return mid;
473 }
474 
475 
477 
479  : geom(geom), ref_type(0), flag(0), index(-1), rank(0), attribute(attr)
480  , parent(NULL)
481 {
482  memset(node, 0, sizeof(node));
483 
484  // NOTE: in 2D the 8-element node/child arrays are not optimal, however,
485  // testing shows we would only save 17% of the total NCMesh memory if
486  // 4-element arrays were used (e.g. through templates); we thus prefer to
487  // keep the code as simple as possible.
488 }
489 
492  Node* n4, Node* n5, Node* n6, Node* n7,
493  int attr,
494  int fattr0, int fattr1, int fattr2,
495  int fattr3, int fattr4, int fattr5)
496 {
497  // create new unrefined element, initialize nodes
498  Element* e = new Element(Geometry::CUBE, attr);
499  e->node[0] = n0, e->node[1] = n1, e->node[2] = n2, e->node[3] = n3;
500  e->node[4] = n4, e->node[5] = n5, e->node[6] = n6, e->node[7] = n7;
501 
502  // get face nodes and assign face attributes
503  Face* f[6];
504  for (int i = 0; i < gi_hex.nf; i++)
505  {
506  const int* fv = gi_hex.faces[i];
507  f[i] = faces.Get(e->node[fv[0]], e->node[fv[1]],
508  e->node[fv[2]], e->node[fv[3]]);
509  }
510 
511  f[0]->attribute = fattr0, f[1]->attribute = fattr1;
512  f[2]->attribute = fattr2, f[3]->attribute = fattr3;
513  f[4]->attribute = fattr4, f[5]->attribute = fattr5;
514 
515  return e;
516 }
517 
520  int attr,
521  int eattr0, int eattr1, int eattr2, int eattr3)
522 {
523  // create new unrefined element, initialize nodes
524  Element* e = new Element(Geometry::SQUARE, attr);
525  e->node[0] = n0, e->node[1] = n1, e->node[2] = n2, e->node[3] = n3;
526 
527  // get edge nodes and assign edge attributes
528  Edge* edge[4];
529  for (int i = 0; i < gi_quad.ne; i++)
530  {
531  const int* ev = gi_quad.edges[i];
532  Node* node = nodes.Get(e->node[ev[0]], e->node[ev[1]]);
533  if (!node->edge) { node->edge = new Edge; }
534  edge[i] = node->edge;
535  }
536 
537  edge[0]->attribute = eattr0;
538  edge[1]->attribute = eattr1;
539  edge[2]->attribute = eattr2;
540  edge[3]->attribute = eattr3;
541 
542  return e;
543 }
544 
547  int attr, int eattr0, int eattr1, int eattr2)
548 {
549  // create new unrefined element, initialize nodes
550  Element* e = new Element(Geometry::TRIANGLE, attr);
551  e->node[0] = n0, e->node[1] = n1, e->node[2] = n2;
552 
553  // get edge nodes and assign edge attributes
554  Edge* edge[3];
555  for (int i = 0; i < gi_tri.ne; i++)
556  {
557  const int* ev = gi_tri.edges[i];
558  Node* node = nodes.Get(e->node[ev[0]], e->node[ev[1]]);
559  if (!node->edge) { node->edge = new Edge; }
560  edge[i] = node->edge;
561  }
562 
563  edge[0]->attribute = eattr0;
564  edge[1]->attribute = eattr1;
565  edge[2]->attribute = eattr2;
566 
567  return e;
568 }
569 
571 {
572  MFEM_ASSERT(v1->vertex && v2->vertex, "missing parent vertices.");
573 
574  // get the midpoint between v1 and v2
575  Vertex* v = new Vertex;
576  for (int i = 0; i < 3; i++)
577  {
578  v->pos[i] = (v1->vertex->pos[i] + v2->vertex->pos[i]) * 0.5;
579  }
580 
581  return v;
582 }
583 
585 {
586  // in 3D we must be careful about getting the mid-edge node
587  Node* mid = PeekAltParents(v1, v2);
588  if (!mid) { mid = nodes.Get(v1, v2); }
589  if (!mid->vertex) { mid->vertex = NewVertex(v1, v2); }
590  return mid;
591 }
592 
594 {
595  // simple version for 2D cases
596  Node* mid = nodes.Get(v1, v2);
597  if (!mid->vertex) { mid->vertex = NewVertex(v1, v2); }
598  return mid;
599 }
600 
603 {
604  // mid-face node can be created either from (e1, e3) or from (e2, e4)
605  Node* midf = nodes.Peek(e1, e3);
606  if (midf)
607  {
608  if (!midf->vertex) { midf->vertex = NewVertex(e1, e3); }
609  return midf;
610  }
611  else
612  {
613  midf = nodes.Get(e2, e4);
614  if (!midf->vertex) { midf->vertex = NewVertex(e2, e4); }
615  return midf;
616  }
617 }
618 
619 //
620 inline bool NCMesh::NodeSetX1(Node* node, Node** n)
621 { return node == n[0] || node == n[3] || node == n[4] || node == n[7]; }
622 
623 inline bool NCMesh::NodeSetX2(Node* node, Node** n)
624 { return node == n[1] || node == n[2] || node == n[5] || node == n[6]; }
625 
626 inline bool NCMesh::NodeSetY1(Node* node, Node** n)
627 { return node == n[0] || node == n[1] || node == n[4] || node == n[5]; }
628 
629 inline bool NCMesh::NodeSetY2(Node* node, Node** n)
630 { return node == n[2] || node == n[3] || node == n[6] || node == n[7]; }
631 
632 inline bool NCMesh::NodeSetZ1(Node* node, Node** n)
633 { return node == n[0] || node == n[1] || node == n[2] || node == n[3]; }
634 
635 inline bool NCMesh::NodeSetZ2(Node* node, Node** n)
636 { return node == n[4] || node == n[5] || node == n[6] || node == n[7]; }
637 
638 
639 void NCMesh::ForceRefinement(Node* v1, Node* v2, Node* v3, Node* v4)
640 {
641  // get the element this face belongs to
642  Face* face = faces.Peek(v1, v2, v3, v4);
643  if (!face) { return; }
644 
645  Element* elem = face->GetSingleElement();
646  MFEM_ASSERT(!elem->ref_type, "element already refined.");
647 
648  Node** nodes = elem->node;
649 
650  // schedule the right split depending on face orientation
651  if ((NodeSetX1(v1, nodes) && NodeSetX2(v2, nodes)) ||
652  (NodeSetX1(v2, nodes) && NodeSetX2(v1, nodes)))
653  {
654  ref_stack.Append(ElemRefType(elem, 1)); // X split
655  }
656  else if ((NodeSetY1(v1, nodes) && NodeSetY2(v2, nodes)) ||
657  (NodeSetY1(v2, nodes) && NodeSetY2(v1, nodes)))
658  {
659  ref_stack.Append(ElemRefType(elem, 2)); // Y split
660  }
661  else if ((NodeSetZ1(v1, nodes) && NodeSetZ2(v2, nodes)) ||
662  (NodeSetZ1(v2, nodes) && NodeSetZ2(v1, nodes)))
663  {
664  ref_stack.Append(ElemRefType(elem, 4)); // Z split
665  }
666  else
667  {
668  MFEM_ABORT("inconsistent element/face structure.");
669  }
670 }
671 
672 
673 void NCMesh::CheckAnisoFace(Node* v1, Node* v2, Node* v3, Node* v4,
674  Node* mid12, Node* mid34, int level)
675 {
676  // When a face is getting split anisotropically (without loss of generality
677  // we assume a "vertical" split here, see picture), it is important to make
678  // sure that the mid-face vertex (midf) has mid34 and mid12 as parents.
679  // This is necessary for the face traversal algorithm and at places like
680  // Refine() that assume the mid-edge nodes to be accessible through the right
681  // parents. However, midf may already exist under the parents mid41 and
682  // mid23. In that case we need to "reparent" midf, i.e., reinsert it to the
683  // hash-table under the correct parents. This doesn't affect other nodes as
684  // all IDs stay the same, only the face refinement "tree" is affected.
685  //
686  // v4 mid34 v3
687  // *------*------*
688  // | | |
689  // | |midf |
690  // mid41 *- - - *- - - * mid23
691  // | | |
692  // | | |
693  // *------*------*
694  // v1 mid12 v2
695  //
696  // This function is recursive, because the above applies to any node along
697  // the middle vertical edge. The function calls itself again for the bottom
698  // and upper half of the above picture.
699 
700  Node* mid23 = nodes.Peek(v2, v3);
701  Node* mid41 = nodes.Peek(v4, v1);
702  if (mid23 && mid41)
703  {
704  Node* midf = nodes.Peek(mid23, mid41);
705  if (midf)
706  {
707  nodes.Reparent(midf, mid12->id, mid34->id);
708 
709  CheckAnisoFace(v1, v2, mid23, mid41, mid12, midf, level+1);
710  CheckAnisoFace(mid41, mid23, v3, v4, midf, mid34, level+1);
711  return;
712  }
713  }
714 
715  // Also, this is the place where forced refinements begin. In the picture
716  // above, edges mid12-midf and midf-mid34 should actually exist in the
717  // neighboring elements, otherwise the mesh is inconsistent and needs to be
718  // fixed. Example: suppose an element is being refined isotropically (!)
719  // whose neighbors across some face look like this:
720  //
721  // *--------*--------*
722  // | d | e |
723  // *--------*--------*
724  // | c |
725  // *--------*--------*
726  // | | |
727  // | a | b |
728  // | | |
729  // *--------*--------*
730  //
731  // Element 'c' needs to be refined vertically for the mesh to remain valid.
732 
733  if (level > 0)
734  {
735  ForceRefinement(v1, v2, v3, v4);
736  }
737 }
738 
739 void NCMesh::CheckIsoFace(Node* v1, Node* v2, Node* v3, Node* v4,
740  Node* e1, Node* e2, Node* e3, Node* e4, Node* midf)
741 {
742  if (!Iso)
743  {
744  /* If anisotropic refinements are present in the mesh, we need to check
745  isotropically split faces as well, see second comment in
746  CheckAnisoFace above. */
747 
748  CheckAnisoFace(v1, v2, e2, e4, e1, midf);
749  CheckAnisoFace(e4, e2, v3, v4, midf, e3);
750  CheckAnisoFace(v4, v1, e1, e3, e4, midf);
751  CheckAnisoFace(e3, e1, v2, v3, midf, e2);
752  }
753 }
754 
755 
756 void NCMesh::RefineElement(Element* elem, char ref_type)
757 {
758  if (!ref_type) { return; }
759 
760  // handle elements that may have been (force-) refined already
761  if (elem->ref_type)
762  {
763  char remaining = ref_type & ~elem->ref_type;
764 
765  // do the remaining splits on the children
766  for (int i = 0; i < 8; i++)
767  {
768  if (elem->child[i]) { RefineElement(elem->child[i], remaining); }
769  }
770  return;
771  }
772 
773  Node** no = elem->node;
774  int attr = elem->attribute;
775 
776  Element* child[8];
777  memset(child, 0, sizeof(child));
778 
779  // create child elements
780  if (elem->geom == Geometry::CUBE)
781  {
782  // get parent's face attributes
783  int fa[6];
784  for (int i = 0; i < gi_hex.nf; i++)
785  {
786  const int* fv = gi_hex.faces[i];
787  Face* face = faces.Peek(no[fv[0]], no[fv[1]], no[fv[2]], no[fv[3]]);
788  fa[i] = face->attribute;
789  }
790 
791  // Vertex numbering is assumed to be as follows:
792  //
793  // 7 6
794  // +------------+ Faces: 0 bottom
795  // /| /| 1 front
796  // 4 / | 5 / | 2 right
797  // +------------+ | 3 back
798  // | | | | 4 left
799  // | +---------|--+ 5 top
800  // | / 3 | / 2 Z Y
801  // |/ |/ |/
802  // +------------+ *--X
803  // 0 1
804 
805  if (ref_type == 1) // split along X axis
806  {
807  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
808  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
809  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
810  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
811 
812  child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
813  no[4], mid45, mid67, no[7], attr,
814  fa[0], fa[1], -1, fa[3], fa[4], fa[5]);
815 
816  child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
817  mid45, no[5], no[6], mid67, attr,
818  fa[0], fa[1], fa[2], fa[3], -1, fa[5]);
819 
820  CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
821  CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
822  CheckAnisoFace(no[4], no[5], no[6], no[7], mid45, mid67);
823  CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
824  }
825  else if (ref_type == 2) // split along Y axis
826  {
827  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
828  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
829  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
830  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
831 
832  child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
833  no[4], no[5], mid56, mid74, attr,
834  fa[0], fa[1], fa[2], -1, fa[4], fa[5]);
835 
836  child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
837  mid74, mid56, no[6], no[7], attr,
838  fa[0], -1, fa[2], fa[3], fa[4], fa[5]);
839 
840  CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
841  CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
842  CheckAnisoFace(no[5], no[6], no[7], no[4], mid56, mid74);
843  CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
844  }
845  else if (ref_type == 4) // split along Z axis
846  {
847  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
848  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
849  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
850  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
851 
852  child[0] = NewHexahedron(no[0], no[1], no[2], no[3],
853  mid04, mid15, mid26, mid37, attr,
854  fa[0], fa[1], fa[2], fa[3], fa[4], -1);
855 
856  child[1] = NewHexahedron(mid04, mid15, mid26, mid37,
857  no[4], no[5], no[6], no[7], attr,
858  -1, fa[1], fa[2], fa[3], fa[4], fa[5]);
859 
860  CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
861  CheckAnisoFace(no[5], no[1], no[2], no[6], mid15, mid26);
862  CheckAnisoFace(no[6], no[2], no[3], no[7], mid26, mid37);
863  CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
864  }
865  else if (ref_type == 3) // XY split
866  {
867  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
868  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
869  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
870  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
871 
872  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
873  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
874  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
875  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
876 
877  Node* midf0 = GetMidFaceVertex(mid23, mid12, mid01, mid30);
878  Node* midf5 = GetMidFaceVertex(mid45, mid56, mid67, mid74);
879 
880  child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
881  no[4], mid45, midf5, mid74, attr,
882  fa[0], fa[1], -1, -1, fa[4], fa[5]);
883 
884  child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
885  mid45, no[5], mid56, midf5, attr,
886  fa[0], fa[1], fa[2], -1, -1, fa[5]);
887 
888  child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
889  midf5, mid56, no[6], mid67, attr,
890  fa[0], -1, fa[2], fa[3], -1, fa[5]);
891 
892  child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
893  mid74, midf5, mid67, no[7], attr,
894  fa[0], -1, -1, fa[3], fa[4], fa[5]);
895 
896  CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
897  CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
898  CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
899  CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
900 
901  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
902  CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
903  }
904  else if (ref_type == 5) // XZ split
905  {
906  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
907  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
908  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
909  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
910 
911  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
912  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
913  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
914  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
915 
916  Node* midf1 = GetMidFaceVertex(mid01, mid15, mid45, mid04);
917  Node* midf3 = GetMidFaceVertex(mid23, mid37, mid67, mid26);
918 
919  child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
920  mid04, midf1, midf3, mid37, attr,
921  fa[0], fa[1], -1, fa[3], fa[4], -1);
922 
923  child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
924  midf1, mid15, mid26, midf3, attr,
925  fa[0], fa[1], fa[2], fa[3], -1, -1);
926 
927  child[2] = NewHexahedron(midf1, mid15, mid26, midf3,
928  mid45, no[5], no[6], mid67, attr,
929  -1, fa[1], fa[2], fa[3], -1, fa[5]);
930 
931  child[3] = NewHexahedron(mid04, midf1, midf3, mid37,
932  no[4], mid45, mid67, no[7], attr,
933  -1, fa[1], -1, fa[3], fa[4], fa[5]);
934 
935  CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
936  CheckAnisoFace(no[2], no[6], no[5], no[1], mid26, mid15);
937  CheckAnisoFace(no[6], no[7], no[4], no[5], mid67, mid45);
938  CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
939 
940  CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
941  CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
942  }
943  else if (ref_type == 6) // YZ split
944  {
945  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
946  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
947  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
948  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
949 
950  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
951  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
952  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
953  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
954 
955  Node* midf2 = GetMidFaceVertex(mid12, mid26, mid56, mid15);
956  Node* midf4 = GetMidFaceVertex(mid30, mid04, mid74, mid37);
957 
958  child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
959  mid04, mid15, midf2, midf4, attr,
960  fa[0], fa[1], fa[2], -1, fa[4], -1);
961 
962  child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
963  midf4, midf2, mid26, mid37, attr,
964  fa[0], -1, fa[2], fa[3], fa[4], -1);
965 
966  child[2] = NewHexahedron(mid04, mid15, midf2, midf4,
967  no[4], no[5], mid56, mid74, attr,
968  -1, fa[1], fa[2], -1, fa[4], fa[5]);
969 
970  child[3] = NewHexahedron(midf4, midf2, mid26, mid37,
971  mid74, mid56, no[6], no[7], attr,
972  -1, -1, fa[2], fa[3], fa[4], fa[5]);
973 
974  CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
975  CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
976  CheckAnisoFace(no[3], no[7], no[6], no[2], mid37, mid26);
977  CheckAnisoFace(no[7], no[4], no[5], no[6], mid74, mid56);
978 
979  CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
980  CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
981  }
982  else if (ref_type == 7) // full isotropic refinement
983  {
984  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
985  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
986  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
987  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
988 
989  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
990  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
991  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
992  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
993 
994  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
995  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
996  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
997  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
998 
999  Node* midf0 = GetMidFaceVertex(mid23, mid12, mid01, mid30);
1000  Node* midf1 = GetMidFaceVertex(mid01, mid15, mid45, mid04);
1001  Node* midf2 = GetMidFaceVertex(mid12, mid26, mid56, mid15);
1002  Node* midf3 = GetMidFaceVertex(mid23, mid37, mid67, mid26);
1003  Node* midf4 = GetMidFaceVertex(mid30, mid04, mid74, mid37);
1004  Node* midf5 = GetMidFaceVertex(mid45, mid56, mid67, mid74);
1005 
1006  Node* midel = GetMidEdgeVertex(midf1, midf3);
1007 
1008  child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
1009  mid04, midf1, midel, midf4, attr,
1010  fa[0], fa[1], -1, -1, fa[4], -1);
1011 
1012  child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
1013  midf1, mid15, midf2, midel, attr,
1014  fa[0], fa[1], fa[2], -1, -1, -1);
1015 
1016  child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
1017  midel, midf2, mid26, midf3, attr,
1018  fa[0], -1, fa[2], fa[3], -1, -1);
1019 
1020  child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
1021  midf4, midel, midf3, mid37, attr,
1022  fa[0], -1, -1, fa[3], fa[4], -1);
1023 
1024  child[4] = NewHexahedron(mid04, midf1, midel, midf4,
1025  no[4], mid45, midf5, mid74, attr,
1026  -1, fa[1], -1, -1, fa[4], fa[5]);
1027 
1028  child[5] = NewHexahedron(midf1, mid15, midf2, midel,
1029  mid45, no[5], mid56, midf5, attr,
1030  -1, fa[1], fa[2], -1, -1, fa[5]);
1031 
1032  child[6] = NewHexahedron(midel, midf2, mid26, midf3,
1033  midf5, mid56, no[6], mid67, attr,
1034  -1, -1, fa[2], fa[3], -1, fa[5]);
1035 
1036  child[7] = NewHexahedron(midf4, midel, midf3, mid37,
1037  mid74, midf5, mid67, no[7], attr,
1038  -1, -1, -1, fa[3], fa[4], fa[5]);
1039 
1040  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1041  CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1042  CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1043  CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1044  CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1045  CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1046  }
1047  else
1048  {
1049  MFEM_ABORT("invalid refinement type.");
1050  }
1051 
1052  if (ref_type != 7) { Iso = false; }
1053  }
1054  else if (elem->geom == Geometry::SQUARE)
1055  {
1056  // get parent's edge attributes
1057  int ea0 = nodes.Peek(no[0], no[1])->edge->attribute;
1058  int ea1 = nodes.Peek(no[1], no[2])->edge->attribute;
1059  int ea2 = nodes.Peek(no[2], no[3])->edge->attribute;
1060  int ea3 = nodes.Peek(no[3], no[0])->edge->attribute;
1061 
1062  ref_type &= ~4; // ignore Z bit
1063 
1064  if (ref_type == 1) // X split
1065  {
1066  Node* mid01 = GetMidEdgeVertexSimple(no[0], no[1]);
1067  Node* mid23 = GetMidEdgeVertexSimple(no[2], no[3]);
1068 
1069  child[0] = NewQuadrilateral(no[0], mid01, mid23, no[3],
1070  attr, ea0, -1, ea2, ea3);
1071 
1072  child[1] = NewQuadrilateral(mid01, no[1], no[2], mid23,
1073  attr, ea0, ea1, ea2, -1);
1074  }
1075  else if (ref_type == 2) // Y split
1076  {
1077  Node* mid12 = GetMidEdgeVertexSimple(no[1], no[2]);
1078  Node* mid30 = GetMidEdgeVertexSimple(no[3], no[0]);
1079 
1080  child[0] = NewQuadrilateral(no[0], no[1], mid12, mid30,
1081  attr, ea0, ea1, -1, ea3);
1082 
1083  child[1] = NewQuadrilateral(mid30, mid12, no[2], no[3],
1084  attr, -1, ea1, ea2, ea3);
1085  }
1086  else if (ref_type == 3) // iso split
1087  {
1088  Node* mid01 = GetMidEdgeVertexSimple(no[0], no[1]);
1089  Node* mid12 = GetMidEdgeVertexSimple(no[1], no[2]);
1090  Node* mid23 = GetMidEdgeVertexSimple(no[2], no[3]);
1091  Node* mid30 = GetMidEdgeVertexSimple(no[3], no[0]);
1092 
1093  Node* midel = GetMidEdgeVertexSimple(mid01, mid23);
1094 
1095  child[0] = NewQuadrilateral(no[0], mid01, midel, mid30,
1096  attr, ea0, -1, -1, ea3);
1097 
1098  child[1] = NewQuadrilateral(mid01, no[1], mid12, midel,
1099  attr, ea0, ea1, -1, -1);
1100 
1101  child[2] = NewQuadrilateral(midel, mid12, no[2], mid23,
1102  attr, -1, ea1, ea2, -1);
1103 
1104  child[3] = NewQuadrilateral(mid30, midel, mid23, no[3],
1105  attr, -1, -1, ea2, ea3);
1106  }
1107  else
1108  {
1109  MFEM_ABORT("Invalid refinement type.");
1110  }
1111 
1112  if (ref_type != 3) { Iso = false; }
1113  }
1114  else if (elem->geom == Geometry::TRIANGLE)
1115  {
1116  // get parent's edge attributes
1117  int ea0 = nodes.Peek(no[0], no[1])->edge->attribute;
1118  int ea1 = nodes.Peek(no[1], no[2])->edge->attribute;
1119  int ea2 = nodes.Peek(no[2], no[0])->edge->attribute;
1120 
1121  ref_type = 3; // for consistence
1122 
1123  // isotropic split - the only ref_type available for triangles
1124  Node* mid01 = GetMidEdgeVertexSimple(no[0], no[1]);
1125  Node* mid12 = GetMidEdgeVertexSimple(no[1], no[2]);
1126  Node* mid20 = GetMidEdgeVertexSimple(no[2], no[0]);
1127 
1128  child[0] = NewTriangle(no[0], mid01, mid20, attr, ea0, -1, ea2);
1129  child[1] = NewTriangle(mid01, no[1], mid12, attr, ea0, ea1, -1);
1130  child[2] = NewTriangle(mid20, mid12, no[2], attr, -1, ea1, ea2);
1131  child[3] = NewTriangle(mid01, mid12, mid20, attr, -1, -1, -1);
1132  }
1133  else
1134  {
1135  MFEM_ABORT("Unsupported element geometry.");
1136  }
1137 
1138  // start using the nodes of the children, create edges & faces
1139  for (int i = 0; i < 8; i++)
1140  {
1141  if (child[i]) { RefElementNodes(child[i]); }
1142  }
1143 
1144  // sign off of all nodes of the parent, clean up unused nodes
1145  UnrefElementNodes(elem);
1146 
1147  // register the children in their faces once the parent is out of the way
1148  for (int i = 0; i < 8; i++)
1149  {
1150  if (child[i]) { RegisterFaces(child[i]); }
1151  }
1152 
1153  // make the children inherit our rank, set the parent pointer
1154  for (int i = 0; i < 8; i++)
1155  {
1156  if (child[i])
1157  {
1158  child[i]->rank = elem->rank;
1159  child[i]->parent = elem;
1160  }
1161  }
1162 
1163  // finish the refinement
1164  elem->ref_type = ref_type;
1165  memcpy(elem->child, child, sizeof(elem->child));
1166 }
1167 
1168 
1169 void NCMesh::Refine(const Array<Refinement>& refinements)
1170 {
1171  // push all refinements on the stack in reverse order
1172  for (int i = refinements.Size()-1; i >= 0; i--)
1173  {
1174  const Refinement& ref = refinements[i];
1175  ref_stack.Append(ElemRefType(leaf_elements[ref.index], ref.ref_type));
1176  }
1177 
1178  // keep refining as long as the stack contains something
1179  int nforced = 0;
1180  while (ref_stack.Size())
1181  {
1182  ElemRefType ref = ref_stack.Last();
1183  ref_stack.DeleteLast();
1184 
1185  int size = ref_stack.Size();
1186  RefineElement(ref.elem, ref.ref_type);
1187  nforced += ref_stack.Size() - size;
1188  }
1189 
1190  /* TODO: the current algorithm of forced refinements is not optimal. As
1191  forced refinements spread through the mesh, some may not be necessary
1192  in the end, since the affected elements may still be scheduled for
1193  refinement that could stop the propagation. We should introduce the
1194  member Element::ref_pending that would show the intended refinement in
1195  the batch. A forced refinement would be combined with ref_pending to
1196  (possibly) stop the propagation earlier.
1197 
1198  Update: what about a FIFO instead of ref_stack? */
1199 
1200 #ifdef MFEM_DEBUG
1201  std::cout << "Refined " << refinements.Size() << " + " << nforced
1202  << " elements" << std::endl;
1203 #endif
1204 
1205  Update();
1206 }
1207 
1208 
1210 
1212 {
1213  if (!elem->ref_type) { return; }
1214 
1215  Element* child[8];
1216  std::memcpy(child, elem->child, sizeof(child));
1217 
1218  // first make sure that all children are leaves, derefine them if not
1219  for (int i = 0; i < 8; i++)
1220  {
1221  if (child[i] && child[i]->ref_type)
1222  {
1223  DerefineElement(child[i]);
1224  }
1225  }
1226 
1227  // retrieve original corner nodes and face attributes from the children
1228  int fa[6];
1229  if (elem->geom == Geometry::CUBE)
1230  {
1231  const int table[7][8 + 6] =
1232  {
1233  { 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0 }, // 1 - X
1234  { 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1 }, // 2 - Y
1235  { 0, 1, 2, 3, 0, 1, 2, 3, 1, 1, 1, 3, 3, 3 }, // 3 - XY
1236  { 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1 }, // 4 - Z
1237  { 0, 1, 1, 0, 3, 2, 2, 3, 1, 1, 1, 3, 3, 3 }, // 5 - XZ
1238  { 0, 0, 1, 1, 2, 2, 3, 3, 0, 0, 0, 3, 3, 3 }, // 6 - YZ
1239  { 0, 1, 2, 3, 4, 5, 6, 7, 1, 1, 1, 7, 7, 7 } // 7 - iso
1240  };
1241  for (int i = 0; i < 8; i++)
1242  {
1243  elem->node[i] = child[table[elem->ref_type - 1][i]]->node[i];
1244  }
1245  for (int i = 0; i < 6; i++)
1246  {
1247  Element* ch = child[table[elem->ref_type - 1][i + 8]];
1248  const int* fv = gi_hex.faces[i];
1249  fa[i] = faces.Peek(ch->node[fv[0]], ch->node[fv[1]],
1250  ch->node[fv[2]], ch->node[fv[3]])->attribute;
1251  }
1252  }
1253  else if (elem->geom == Geometry::SQUARE)
1254  {
1255  const int table[3][4 + 4] =
1256  {
1257  { 0, 1, 1, 0, 1, 1, 0, 0 }, // 1 - X
1258  { 0, 0, 1, 1, 0, 0, 1, 1 }, // 2 - Y
1259  { 0, 1, 2, 3, 1, 1, 3, 3 } // 3 - iso
1260  };
1261  for (int i = 0; i < 4; i++)
1262  {
1263  elem->node[i] = child[table[elem->ref_type - 1][i]]->node[i];
1264  }
1265  for (int i = 0; i < 4; i++)
1266  {
1267  Element* ch = child[table[elem->ref_type - 1][i + 4]];
1268  const int* ev = gi_quad.edges[i];
1269  fa[i] = nodes.Peek(ch->node[ev[0]], ch->node[ev[1]])->edge->attribute;
1270  }
1271  }
1272  else if (elem->geom == Geometry::TRIANGLE)
1273  {
1274  for (int i = 0; i < 3; i++)
1275  {
1276  Element* ch = child[i];
1277  elem->node[i] = child[i]->node[i];
1278  const int* ev = gi_tri.edges[i];
1279  fa[i] = nodes.Peek(ch->node[ev[0]], ch->node[ev[1]])->edge->attribute;
1280  }
1281  }
1282  else
1283  {
1284  MFEM_ABORT("Unsupported element geometry.");
1285  }
1286 
1287  // sign in to all nodes again
1288  RefElementNodes(elem);
1289 
1290  // delete children, determine rank
1291  elem->rank = INT_MAX;
1292  for (int i = 0; i < 8; i++)
1293  {
1294  if (child[i])
1295  {
1296  elem->rank = std::min(elem->rank, child[i]->rank);
1297  DeleteHierarchy(child[i]);
1298  }
1299  }
1300 
1301  RegisterFaces(elem, fa);
1302 
1303  // set edge attributes (2D)
1304  // TODO: Edge::attribute should be removed
1305  if (Dim < 3)
1306  {
1307  Node** node = elem->node;
1308  GeomInfo& gi = GI[(int) elem->geom];
1309  for (int i = 0; i < gi.ne; i++)
1310  {
1311  const int* ev = gi.edges[i];
1312  nodes.Peek(node[ev[0]], node[ev[1]])->edge->attribute = fa[i];
1313  }
1314  }
1315 
1316  elem->ref_type = 0;
1317 }
1318 
1319 
1321 {
1322  if (!elem->ref_type) { return; }
1323 
1324  int total = 0, ref = 0, ghost = 0;
1325  for (int i = 0; i < 8; i++)
1326  {
1327  Element* ch = elem->child[i];
1328  if (ch)
1329  {
1330  total++;
1331  if (ch->ref_type) { ref++; break; }
1332  if (IsGhost(ch)) { ghost++; }
1333  }
1334  }
1335 
1336  if (!ref && ghost < total)
1337  {
1338  // can be derefined, add to list
1339  int next_row = list.Size() ? (list.Last().from + 1) : 0;
1340  for (int i = 0; i < 8; i++)
1341  {
1342  Element* ch = elem->child[i];
1343  if (ch) { list.Append(Connection(next_row, ch->index)); }
1344  }
1345  }
1346  else
1347  {
1348  for (int i = 0; i < 8; i++)
1349  {
1350  Element* ch = elem->child[i];
1351  if (ch) { CollectDerefinements(ch, list); }
1352  }
1353  }
1354 }
1355 
1357 {
1358  Array<Connection> list;
1359  list.Reserve(leaf_elements.Size());
1360 
1361  for (int i = 0; i < root_elements.Size(); i++)
1362  {
1364  }
1365 
1366  int size = list.Size() ? (list.Last().from + 1) : 0;
1367  derefinements.MakeFromList(size, list);
1368  return derefinements;
1369 }
1370 
1371 void NCMesh::CheckDerefinementNCLevel(const Table &deref_table,
1372  Array<int> &level_ok, int max_nc_level)
1373 {
1374  level_ok.SetSize(deref_table.Size());
1375  for (int i = 0; i < deref_table.Size(); i++)
1376  {
1377  const int* fine = deref_table.GetRow(i), size = deref_table.RowSize(i);
1378  Element* parent = leaf_elements[fine[0]]->parent;
1379 
1380  int ok = 1;
1381  for (int j = 0; j < size; j++)
1382  {
1383  int splits[3];
1384  CountSplits(leaf_elements[fine[j]], splits);
1385 
1386  for (int k = 0; k < Dim; k++)
1387  {
1388  if ((parent->ref_type & (1 << k)) &&
1389  splits[k] >= max_nc_level)
1390  {
1391  ok = 0; break;
1392  }
1393  }
1394  if (!ok) { break; }
1395  }
1396  level_ok[i] = ok;
1397  }
1398 }
1399 
1400 void NCMesh::Derefine(const Array<int> &derefs)
1401 {
1402  MFEM_VERIFY(Dim < 3 || Iso,
1403  "derefinement of 3D anisotropic meshes not implemented yet.");
1404 
1406 
1407  Array<Element*> coarse;
1408  leaf_elements.Copy(coarse);
1409 
1410  // perform the derefinements
1411  for (int i = 0; i < derefs.Size(); i++)
1412  {
1413  int row = derefs[i];
1414  MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
1415  "invalid derefinement number.");
1416 
1417  const int* fine = derefinements.GetRow(row);
1418  Element* parent = leaf_elements[fine[0]]->parent;
1419 
1420  // record the relation of the fine elements to their parent
1421  SetDerefMatrixCodes(parent, coarse);
1422 
1423  DerefineElement(parent);
1424  }
1425 
1426  // update leaf_elements, Element::index etc.
1427  Update();
1428 
1429  // link old fine elements to the new coarse elements
1430  for (int i = 0; i < coarse.Size(); i++)
1431  {
1432  transforms.embeddings[i].parent = coarse[i]->index;
1433  }
1434 }
1435 
1437 {
1438  int nfine = leaf_elements.Size();
1439 
1440  transforms.embeddings.SetSize(nfine);
1441  for (int i = 0; i < nfine; i++)
1442  {
1443  transforms.embeddings[i].parent = -1;
1444  transforms.embeddings[i].matrix = 0;
1445  }
1446 
1447  // this will tell GetDerefinementTransforms that transforms are not finished
1449 }
1450 
1452 {
1453  // encode the ref_type and child number for GetDerefinementTransforms()
1454  for (int i = 0; i < 8; i++)
1455  {
1456  Element* ch = parent->child[i];
1457  if (ch && ch->index >= 0)
1458  {
1459  int code = (parent->ref_type << 3) + i;
1460  transforms.embeddings[ch->index].matrix = code;
1461  coarse[ch->index] = parent;
1462  }
1463  }
1464 }
1465 
1466 
1468 
1470 {
1471  // (overridden in ParNCMesh to assign special indices to ghost vertices)
1472  int num_vert = 0;
1473  for (HashTable<Node>::Iterator it(nodes); it; ++it)
1474  {
1475  if (it->vertex) { it->vertex->index = num_vert++; }
1476  }
1477 
1478  vertex_nodeId.SetSize(num_vert);
1479 
1480  num_vert = 0;
1481  for (HashTable<Node>::Iterator it(nodes); it; ++it)
1482  {
1483  if (it->vertex) { vertex_nodeId[num_vert++] = it->id; }
1484  }
1485 }
1486 
1487 static char quad_hilbert_child_order[8][4] =
1488 {
1489  {0,1,2,3}, {0,3,2,1}, {1,2,3,0}, {1,0,3,2},
1490  {2,3,0,1}, {2,1,0,3}, {3,0,1,2}, {3,2,1,0}
1491 };
1492 static char quad_hilbert_child_state[8][4] =
1493 {
1494  {1,0,0,5}, {0,1,1,4}, {3,2,2,7}, {2,3,3,6},
1495  {5,4,4,1}, {4,5,5,0}, {7,6,6,3}, {6,7,7,2}
1496 };
1497 static char hex_hilbert_child_order[24][8] =
1498 {
1499  {0,1,2,3,7,6,5,4}, {0,3,7,4,5,6,2,1}, {0,4,5,1,2,6,7,3},
1500  {1,0,3,2,6,7,4,5}, {1,2,6,5,4,7,3,0}, {1,5,4,0,3,7,6,2},
1501  {2,1,5,6,7,4,0,3}, {2,3,0,1,5,4,7,6}, {2,6,7,3,0,4,5,1},
1502  {3,0,4,7,6,5,1,2}, {3,2,1,0,4,5,6,7}, {3,7,6,2,1,5,4,0},
1503  {4,0,1,5,6,2,3,7}, {4,5,6,7,3,2,1,0}, {4,7,3,0,1,2,6,5},
1504  {5,1,0,4,7,3,2,6}, {5,4,7,6,2,3,0,1}, {5,6,2,1,0,3,7,4},
1505  {6,2,3,7,4,0,1,5}, {6,5,1,2,3,0,4,7}, {6,7,4,5,1,0,3,2},
1506  {7,3,2,6,5,1,0,4}, {7,4,0,3,2,1,5,6}, {7,6,5,4,0,1,2,3}
1507 };
1508 static char hex_hilbert_child_state[24][8] =
1509 {
1510  {1,2,2,7,7,21,21,17}, {2,0,0,22,22,16,16,8}, {0,1,1,15,15,6,6,23},
1511  {4,5,5,10,10,18,18,14}, {5,3,3,19,19,13,13,11}, {3,4,4,12,12,9,9,20},
1512  {8,7,7,17,17,23,23,2}, {6,8,8,0,0,15,15,22}, {7,6,6,21,21,1,1,16},
1513  {11,10,10,14,14,20,20,5}, {9,11,11,3,3,12,12,19}, {10,9,9,18,18,4,4,13},
1514  {13,14,14,5,5,19,19,10}, {14,12,12,20,20,11,11,4}, {12,13,13,9,9,3,3,18},
1515  {16,17,17,2,2,22,22,7}, {17,15,15,23,23,8,8,1}, {15,16,16,6,6,0,0,21},
1516  {20,19,19,11,11,14,14,3}, {18,20,20,4,4,10,10,12}, {19,18,18,13,13,5,5,9},
1517  {23,22,22,8,8,17,17,0}, {21,23,23,1,1,7,7,15}, {22,21,21,16,16,2,2,6}
1518 };
1519 
1520 void NCMesh::CollectLeafElements(Element* elem, int state)
1521 {
1522  if (!elem->ref_type)
1523  {
1524  if (elem->rank >= 0) // skip elements beyond ghost layer in parallel
1525  {
1526  leaf_elements.Append(elem);
1527  }
1528  }
1529  else
1530  {
1531  if (elem->geom == Geometry::SQUARE && elem->ref_type == 3)
1532  {
1533  for (int i = 0; i < 4; i++)
1534  {
1535  int ch = quad_hilbert_child_order[state][i];
1536  int st = quad_hilbert_child_state[state][i];
1537  CollectLeafElements(elem->child[ch], st);
1538  }
1539  }
1540  else if (elem->geom == Geometry::CUBE && elem->ref_type == 7)
1541  {
1542  for (int i = 0; i < 8; i++)
1543  {
1544  int ch = hex_hilbert_child_order[state][i];
1545  int st = hex_hilbert_child_state[state][i];
1546  CollectLeafElements(elem->child[ch], st);
1547  }
1548  }
1549  else
1550  {
1551  for (int i = 0; i < 8; i++)
1552  {
1553  if (elem->child[i]) { CollectLeafElements(elem->child[i], state); }
1554  }
1555  }
1556  }
1557  elem->index = -1;
1558 }
1559 
1561 {
1562  // collect leaf elements from all roots
1563  leaf_elements.SetSize(0);
1564  for (int i = 0; i < root_elements.Size(); i++)
1565  {
1567  // TODO: root state should not always be 0, we need a precomputed array
1568  // with root element states to ensure continuity where possible, also
1569  // optimized ordering of the root elements themselves (Gecko?)
1570  }
1572 }
1573 
1575 {
1576  // (overridden in ParNCMesh to handle ghost elements)
1577  for (int i = 0; i < leaf_elements.Size(); i++)
1578  {
1579  leaf_elements[i]->index = i;
1580  }
1581 }
1582 
1584 {
1585  switch (geom)
1586  {
1587  case Geometry::CUBE: return new mfem::Hexahedron;
1588  case Geometry::SQUARE: return new mfem::Quadrilateral;
1589  case Geometry::TRIANGLE: return new mfem::Triangle;
1590  }
1591  MFEM_ABORT("invalid geometry");
1592  return NULL;
1593 }
1594 
1596  Array<mfem::Element*>& elements,
1597  Array<mfem::Element*>& boundary) const
1598 {
1599  // copy vertex coordinates
1600  vertices.SetSize(vertex_nodeId.Size());
1601  for (int i = 0; i < vertices.Size(); i++)
1602  {
1603  Node* node = nodes.Peek(vertex_nodeId[i]);
1604  vertices[i].SetCoords(node->vertex->pos);
1605  }
1606 
1607  elements.SetSize(leaf_elements.Size() - GetNumGhosts());
1608  elements.SetSize(0);
1609 
1610  boundary.SetSize(0);
1611 
1612  // create an mfem::Element for each leaf Element
1613  for (int i = 0; i < leaf_elements.Size(); i++)
1614  {
1615  Element* nc_elem = leaf_elements[i];
1616  if (IsGhost(nc_elem)) { continue; } // ParNCMesh
1617 
1618  Node** node = nc_elem->node;
1619  GeomInfo& gi = GI[(int) nc_elem->geom];
1620 
1621  mfem::Element* elem = NewMeshElement(nc_elem->geom);
1622  elements.Append(elem);
1623 
1624  elem->SetAttribute(nc_elem->attribute);
1625  for (int j = 0; j < gi.nv; j++)
1626  {
1627  elem->GetVertices()[j] = node[j]->vertex->index;
1628  }
1629 
1630  // create boundary elements
1631  if (nc_elem->geom == Geometry::CUBE)
1632  {
1633  for (int k = 0; k < gi.nf; k++)
1634  {
1635  const int* fv = gi.faces[k];
1636  Face* face = faces.Peek(node[fv[0]], node[fv[1]],
1637  node[fv[2]], node[fv[3]]);
1638  if (face->Boundary())
1639  {
1640  Quadrilateral* quad = new Quadrilateral;
1641  quad->SetAttribute(face->attribute);
1642  for (int j = 0; j < 4; j++)
1643  {
1644  quad->GetVertices()[j] = node[fv[j]]->vertex->index;
1645  }
1646  boundary.Append(quad);
1647  }
1648  }
1649  }
1650  else // quad & triangle boundary elements
1651  {
1652  for (int k = 0; k < gi.ne; k++)
1653  {
1654  const int* ev = gi.edges[k];
1655  Edge* edge = nodes.Peek(node[ev[0]], node[ev[1]])->edge;
1656  if (edge->Boundary())
1657  {
1658  Segment* segment = new Segment;
1659  segment->SetAttribute(edge->attribute);
1660  for (int j = 0; j < 2; j++)
1661  {
1662  segment->GetVertices()[j] = node[ev[j]]->vertex->index;
1663  }
1664  boundary.Append(segment);
1665  }
1666  }
1667  }
1668  }
1669 }
1670 
1672 {
1673  Table *edge_vertex = mesh->GetEdgeVertexTable();
1674 
1675  // get edge enumeration from the Mesh
1676  for (int i = 0; i < edge_vertex->Size(); i++)
1677  {
1678  const int *ev = edge_vertex->GetRow(i);
1679  Node* node = nodes.Peek(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
1680 
1681  MFEM_ASSERT(node && node->edge, "edge not found.");
1682  node->edge->index = i;
1683  }
1684 
1685  // get face enumeration from the Mesh
1686  for (int i = 0; i < mesh->GetNumFaces(); i++)
1687  {
1688  const int* fv = mesh->GetFace(i)->GetVertices();
1689  Face* face;
1690  if (Dim == 3)
1691  {
1692  MFEM_ASSERT(mesh->GetFace(i)->GetNVertices() == 4, "");
1693  face = faces.Peek(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
1694  vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
1695  }
1696  else
1697  {
1698  MFEM_ASSERT(mesh->GetFace(i)->GetNVertices() == 2, "");
1699  int n0 = vertex_nodeId[fv[0]], n1 = vertex_nodeId[fv[1]];
1700  face = faces.Peek(n0, n0, n1, n1); // look up degenerate face
1701 
1702 #ifdef MFEM_DEBUG
1703  // (non-ghost) edge and face numbers must match in 2D
1704  const int *ev = edge_vertex->GetRow(i);
1705  MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
1706  (ev[1] == fv[0] && ev[0] == fv[1]), "");
1707 #endif
1708  }
1709  MFEM_ASSERT(face, "face not found.");
1710  face->index = i;
1711  }
1712 }
1713 
1714 
1716 
1717 int NCMesh::FaceSplitType(Node* v1, Node* v2, Node* v3, Node* v4,
1718  Node* mid[4]) const
1719 {
1720  MFEM_ASSERT(Dim >= 3, "");
1721 
1722  // find edge nodes
1723  Node* e1 = nodes.Peek(v1, v2);
1724  Node* e2 = nodes.Peek(v2, v3);
1725  Node* e3 = e1 ? nodes.Peek(v3, v4) : NULL; // TODO: e1 && e1->vertex ?
1726  Node* e4 = e2 ? nodes.Peek(v4, v1) : NULL;
1727 
1728  // optional: return the mid-edge nodes if requested
1729  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
1730 
1731  // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
1732  Node *midf1 = NULL, *midf2 = NULL;
1733  if (e1 && e3) { midf1 = nodes.Peek(e1, e3); }
1734  if (e2 && e4) { midf2 = nodes.Peek(e2, e4); }
1735 
1736  // only one way to access the mid-face node must always exist
1737  MFEM_ASSERT(!(midf1 && midf2), "incorrectly split face!");
1738 
1739  if (!midf1 && !midf2) { return 0; } // face not split
1740 
1741  if (midf1) { return 1; } // face split "vertically"
1742  else { return 2; } // face split "horizontally"
1743 }
1744 
1745 int NCMesh::find_node(Element* elem, Node* node)
1746 {
1747  for (int i = 0; i < 8; i++)
1748  if (elem->node[i] == node) { return i; }
1749 
1750  MFEM_ABORT("Node not found.");
1751  return -1;
1752 }
1753 
1754 int NCMesh::find_node(Element* elem, int node_id)
1755 {
1756  for (int i = 0; i < 8; i++)
1757  if (elem->node[i]->id == node_id) { return i; }
1758 
1759  MFEM_ABORT("Node not found.");
1760  return -1;
1761 }
1762 
1763 int NCMesh::find_element_edge(Element* elem, int v0, int v1)
1764 {
1765  MFEM_ASSERT(!elem->ref_type, "");
1766  GeomInfo &gi = GI[(int) elem->geom];
1767  for (int i = 0; i < gi.ne; i++)
1768  {
1769  const int* ev = gi.edges[i];
1770  int n0 = elem->node[ev[0]]->id;
1771  int n1 = elem->node[ev[1]]->id;
1772  if ((n0 == v0 && n1 == v1) ||
1773  (n0 == v1 && n1 == v0)) { return i; }
1774  }
1775  MFEM_ABORT("Edge not found");
1776  return -1;
1777 }
1778 
1779 int NCMesh::find_hex_face(int a, int b, int c)
1780 {
1781  for (int i = 0; i < 6; i++)
1782  {
1783  const int* fv = gi_hex.faces[i];
1784  if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
1785  (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
1786  (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
1787  {
1788  return i;
1789  }
1790  }
1791  MFEM_ABORT("Face not found.");
1792  return -1;
1793 }
1794 
1796  Element* elem, DenseMatrix& mat) const
1797 {
1798  int master[4] =
1799  {
1800  find_node(elem, v0), find_node(elem, v1),
1801  find_node(elem, v2), find_node(elem, v3)
1802  };
1803 
1804  int local = find_hex_face(master[0], master[1], master[2]);
1805  const int* fv = gi_hex.faces[local];
1806 
1807  DenseMatrix tmp(mat);
1808  for (int i = 0, j; i < 4; i++)
1809  {
1810  for (j = 0; j < 4; j++)
1811  {
1812  if (fv[i] == master[j])
1813  {
1814  // "pm.column(i) = tmp.column(j)"
1815  for (int k = 0; k < mat.Height(); k++)
1816  {
1817  mat(k,i) = tmp(k,j);
1818  }
1819  break;
1820  }
1821  }
1822  MFEM_ASSERT(j != 4, "node not found.");
1823  }
1824  return local;
1825 }
1826 
1827 void NCMesh::TraverseFace(Node* v0, Node* v1, Node* v2, Node* v3,
1828  const PointMatrix& pm, int level)
1829 {
1830  if (level > 0)
1831  {
1832  // check if we made it to a face that is not split further
1833  Face* face = faces.Peek(v0, v1, v2, v3);
1834  if (face)
1835  {
1836  // we have a slave face, add it to the list
1837  Element* elem = face->GetSingleElement();
1838  face_list.slaves.push_back(Slave(face->index, elem, -1));
1839  DenseMatrix &mat = face_list.slaves.back().point_matrix;
1840  pm.GetMatrix(mat);
1841 
1842  // reorder the point matrix according to slave face orientation
1843  int local = ReorderFacePointMat(v0, v1, v2, v3, elem, mat);
1844  face_list.slaves.back().local = local;
1845 
1846  return;
1847  }
1848  }
1849 
1850  // we need to recurse deeper
1851  Node* mid[4];
1852  int split = FaceSplitType(v0, v1, v2, v3, mid);
1853 
1854  if (split == 1) // "X" split face
1855  {
1856  Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
1857 
1858  TraverseFace(v0, mid[0], mid[2], v3,
1859  PointMatrix(pm(0), mid0, mid2, pm(3)), level+1);
1860 
1861  TraverseFace(mid[0], v1, v2, mid[2],
1862  PointMatrix(mid0, pm(1), pm(2), mid2), level+1);
1863  }
1864  else if (split == 2) // "Y" split face
1865  {
1866  Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
1867 
1868  TraverseFace(v0, v1, mid[1], mid[3],
1869  PointMatrix(pm(0), pm(1), mid1, mid3), level+1);
1870 
1871  TraverseFace(mid[3], mid[1], v2, v3,
1872  PointMatrix(mid3, mid1, pm(2), pm(3)), level+1);
1873  }
1874 }
1875 
1877 {
1878  face_list.Clear();
1879  boundary_faces.SetSize(0);
1880 
1881  if (Dim < 3) { return; }
1882 
1883  // visit faces of leaf elements
1884  Array<char> processed_faces; // TODO: size
1885  for (int i = 0; i < leaf_elements.Size(); i++)
1886  {
1887  Element* elem = leaf_elements[i];
1888  MFEM_ASSERT(!elem->ref_type, "not a leaf element.");
1889 
1890  GeomInfo& gi = GI[(int) elem->geom];
1891  for (int j = 0; j < gi.nf; j++)
1892  {
1893  // get nodes for this face
1894  Node* node[4];
1895  for (int k = 0; k < 4; k++)
1896  {
1897  node[k] = elem->node[gi.faces[j][k]];
1898  }
1899 
1900  Face* face = faces.Peek(node[0], node[1], node[2], node[3]);
1901  MFEM_ASSERT(face, "face not found!");
1902 
1903  // tell ParNCMesh about the face
1904  ElementSharesFace(elem, face);
1905 
1906  // have we already processed this face? skip if yes
1907  int index = face->index;
1908  if (index >= processed_faces.Size())
1909  {
1910  processed_faces.SetSize(index + 1, 0);
1911  }
1912  else if (processed_faces[index])
1913  {
1914  continue;
1915  }
1916  processed_faces[index] = 1;
1917 
1918  if (face->ref_count == 2)
1919  {
1920  // this is a conforming face, add it to the list
1921  face_list.conforming.push_back(MeshId(index, elem, j));
1922  }
1923  else
1924  {
1925  PointMatrix pm(Point(0,0), Point(1,0), Point(1,1), Point(0,1));
1926 
1927  // this is either a master face or a slave face, but we can't
1928  // tell until we traverse the face refinement 'tree'...
1929  int sb = face_list.slaves.size();
1930  TraverseFace(node[0], node[1], node[2], node[3], pm, 0);
1931 
1932  int se = face_list.slaves.size();
1933  if (sb < se)
1934  {
1935  // found slaves, so this is a master face; add it to the list
1936  face_list.masters.push_back(Master(index, elem, j, sb, se));
1937 
1938  // also, set the master index for the slaves
1939  for (int i = sb; i < se; i++)
1940  {
1941  face_list.slaves[i].master = index;
1942  }
1943  }
1944  }
1945 
1946  if (face->Boundary()) { boundary_faces.Append(face); }
1947  }
1948  }
1949 }
1950 
1951 void NCMesh::TraverseEdge(Node* v0, Node* v1, double t0, double t1, int flags,
1952  int level)
1953 {
1954  Node* mid = nodes.Peek(v0, v1);
1955  if (!mid) { return; }
1956 
1957  if (mid->edge && level > 0)
1958  {
1959  // we have a slave edge, add it to the list
1960  edge_list.slaves.push_back(Slave(mid->edge->index, NULL, -1));
1961  Slave &sl = edge_list.slaves.back();
1962 
1963  sl.point_matrix.SetSize(1, 2);
1964  sl.point_matrix(0,0) = t0;
1965  sl.point_matrix(0,1) = t1;
1966 
1967  // handle slave edge orientation
1968  sl.edge_flags = flags;
1969  if (v0->vertex->index > v1->vertex->index) { sl.edge_flags |= 2; }
1970 
1971  // in 2D, get the element/local info from the degenerate face
1972  if (Dim == 2)
1973  {
1974  Face* face = faces.Peek(v0, v0, v1, v1);
1975  MFEM_ASSERT(face != NULL, "");
1976  sl.element = face->GetSingleElement();
1977  sl.local = find_element_edge(sl.element, v0->id, v1->id);
1978  }
1979  }
1980 
1981  // recurse deeper
1982  double tmid = (t0 + t1) / 2;
1983  TraverseEdge(v0, mid, t0, tmid, flags, level+1);
1984  TraverseEdge(mid, v1, tmid, t1, flags, level+1);
1985 }
1986 
1988 {
1989  edge_list.Clear();
1990  boundary_edges.SetSize(0);
1991 
1992  // visit edges of leaf elements
1993  Array<char> processed_edges; // TODO: size
1994  for (int i = 0; i < leaf_elements.Size(); i++)
1995  {
1996  Element* elem = leaf_elements[i];
1997  MFEM_ASSERT(!elem->ref_type, "not a leaf element.");
1998 
1999  GeomInfo& gi = GI[(int) elem->geom];
2000  for (int j = 0; j < gi.ne; j++)
2001  {
2002  // get nodes for this edge
2003  const int* ev = gi.edges[j];
2004  Node* node[2] = { elem->node[ev[0]], elem->node[ev[1]] };
2005 
2006  Node* edge = nodes.Peek(node[0], node[1]);
2007  MFEM_ASSERT(edge && edge->edge, "edge not found!");
2008 
2009  // tell ParNCMesh about the edge
2010  ElementSharesEdge(elem, edge->edge);
2011 
2012  // (2D only, store boundary edges)
2013  if (edge->edge->Boundary()) { boundary_edges.Append(edge); }
2014 
2015  // skip slave edges here, they will be reached from their masters
2016  if (GetEdgeMaster(edge)) { continue; }
2017 
2018  // have we already processed this edge? skip if yes
2019  int index = edge->edge->index;
2020  if (index >= processed_edges.Size())
2021  {
2022  processed_edges.SetSize(index + 1, 0);
2023  }
2024  else if (processed_edges[index])
2025  {
2026  continue;
2027  }
2028  processed_edges[index] = 1;
2029 
2030  // prepare edge interval for slave traversal, handle orientation
2031  double t0 = 0.0, t1 = 1.0;
2032  int flags = (node[0]->vertex->index > node[1]->vertex->index) ? 1 : 0;
2033 
2034  // try traversing the edge to find slave edges
2035  int sb = edge_list.slaves.size();
2036  TraverseEdge(node[0], node[1], t0, t1, flags, 0);
2037 
2038  int se = edge_list.slaves.size();
2039  if (sb < se)
2040  {
2041  // found slaves, this is a master face; add it to the list
2042  edge_list.masters.push_back(Master(index, elem, j, sb, se));
2043 
2044  // also, set the master index for the slaves
2045  for (int i = sb; i < se; i++)
2046  {
2047  edge_list.slaves[i].master = index;
2048  }
2049  }
2050  else
2051  {
2052  // no slaves, this is a conforming edge
2053  edge_list.conforming.push_back(MeshId(index, elem, j));
2054  }
2055  }
2056  }
2057 }
2058 
2060 {
2061  oriented_matrix = point_matrix;
2062 
2063  if (edge_flags)
2064  {
2065  MFEM_ASSERT(oriented_matrix.Height() == 1 &&
2066  oriented_matrix.Width() == 2, "not an edge point matrix");
2067 
2068  if (edge_flags & 1) // master inverted
2069  {
2070  oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
2071  oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
2072  }
2073  if (edge_flags & 2) // slave inverted
2074  {
2075  std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
2076  }
2077  }
2078 }
2079 
2080 
2082 
2084 {
2085  Node* mid = nodes.Peek(v0, v1);
2086  if (mid && mid->vertex)
2087  {
2088  indices.Append(mid->vertex->index);
2089 
2090  CollectEdgeVertices(v0, mid, indices);
2091  CollectEdgeVertices(mid, v1, indices);
2092  }
2093 }
2094 
2096  Array<int> &indices)
2097 {
2098  Node* mid[4];
2099  switch (FaceSplitType(v0, v1, v2, v3, mid))
2100  {
2101  case 1:
2102  indices.Append(mid[0]->vertex->index);
2103  indices.Append(mid[2]->vertex->index);
2104 
2105  CollectFaceVertices(v0, mid[0], mid[2], v3, indices);
2106  CollectFaceVertices(mid[0], v1, v2, mid[2], indices);
2107  break;
2108 
2109  case 2:
2110  indices.Append(mid[1]->vertex->index);
2111  indices.Append(mid[3]->vertex->index);
2112 
2113  CollectFaceVertices(v0, v1, mid[1], mid[3], indices);
2114  CollectFaceVertices(mid[3], mid[1], v2, v3, indices);
2115  break;
2116  }
2117 }
2118 
2120 {
2121  int nrows = leaf_elements.Size();
2122  int* I = new int[nrows + 1];
2123  int** JJ = new int*[nrows];
2124 
2125  Array<int> indices;
2126  indices.Reserve(128);
2127 
2128  // collect vertices coinciding with each element, including hanging vertices
2129  for (int i = 0; i < leaf_elements.Size(); i++)
2130  {
2131  Element* elem = leaf_elements[i];
2132  MFEM_ASSERT(!elem->ref_type, "not a leaf element.");
2133 
2134  GeomInfo& gi = GI[(int) elem->geom];
2135  Node** node = elem->node;
2136 
2137  indices.SetSize(0);
2138  for (int j = 0; j < gi.nv; j++)
2139  {
2140  indices.Append(node[j]->vertex->index);
2141  }
2142 
2143  for (int j = 0; j < gi.ne; j++)
2144  {
2145  const int* ev = gi.edges[j];
2146  CollectEdgeVertices(node[ev[0]], node[ev[1]], indices);
2147  }
2148 
2149  if (Dim >= 3)
2150  {
2151  for (int j = 0; j < gi.nf; j++)
2152  {
2153  const int* fv = gi.faces[j];
2154  CollectFaceVertices(node[fv[0]], node[fv[1]],
2155  node[fv[2]], node[fv[3]], indices);
2156  }
2157  }
2158 
2159  // temporarily store one row of the table
2160  indices.Sort();
2161  indices.Unique();
2162  int size = indices.Size();
2163  I[i] = size;
2164  JJ[i] = new int[size];
2165  memcpy(JJ[i], indices.GetData(), size * sizeof(int));
2166  }
2167 
2168  // finalize the I array of the table
2169  int nnz = 0;
2170  for (int i = 0; i < nrows; i++)
2171  {
2172  int cnt = I[i];
2173  I[i] = nnz;
2174  nnz += cnt;
2175  }
2176  I[nrows] = nnz;
2177 
2178  // copy the temporarily stored rows into one J array
2179  int *J = new int[nnz];
2180  nnz = 0;
2181  for (int i = 0; i < nrows; i++)
2182  {
2183  int cnt = I[i+1] - I[i];
2184  memcpy(J+nnz, JJ[i], cnt * sizeof(int));
2185  delete [] JJ[i];
2186  nnz += cnt;
2187  }
2188 
2189  element_vertex.SetIJ(I, J, nrows);
2191 
2192  delete [] JJ;
2193 }
2194 
2195 
2197  Array<Element*> *neighbors,
2198  Array<char> *neighbor_set)
2199 {
2200  // If A is the element-to-vertex table (see 'element_vertex') listing all
2201  // vertices touching each element, including hanging vertices, then A*A^T is
2202  // the element-to-neighbor table. Multiplying the element set with A*A^T
2203  // gives the neighbor set. To save memory, this function only computes the
2204  // action of A*A^T, the product itself is not stored anywhere.
2205 
2206  // TODO: we should further optimize the 'element_vertex' table by not storing
2207  // the obvious corner vertices in it. The table would therefore be empty
2208  // for conforming meshes and very cheap for NC meshes.
2209 
2211 
2212  int nleaves = leaf_elements.Size();
2213  MFEM_VERIFY(elem_set.Size() == nleaves, "");
2214  MFEM_ASSERT(element_vertex.Size() == nleaves, "");
2215 
2216  // step 1: vertices = A^T * elem_set, i.e, find all vertices touching the
2217  // element set
2218 
2219  Array<char> vertices(num_vertices);
2220  vertices = 0;
2221 
2222  for (int i = 0; i < nleaves; i++)
2223  {
2224  if (elem_set[i])
2225  {
2226  int *v = element_vertex.GetRow(i);
2227  int nv = element_vertex.RowSize(i);
2228  for (int j = 0; j < nv; j++)
2229  {
2230  vertices[v[j]] = 1;
2231  }
2232  }
2233  }
2234 
2235  // step 2: neighbors = A * vertices, i.e., find all elements coinciding with
2236  // vertices from step 1; NOTE: in the result we don't include elements from
2237  // the original set
2238 
2239  if (neighbor_set)
2240  {
2241  neighbor_set->SetSize(nleaves);
2242  *neighbor_set = 0;
2243  }
2244 
2245  for (int i = 0; i < nleaves; i++)
2246  {
2247  if (!elem_set[i])
2248  {
2249  int *v = element_vertex.GetRow(i);
2250  int nv = element_vertex.RowSize(i);
2251  for (int j = 0; j < nv; j++)
2252  {
2253  if (vertices[v[j]])
2254  {
2255  if (neighbors) { neighbors->Append(leaf_elements[i]); }
2256  if (neighbor_set) { (*neighbor_set)[i] = 1; }
2257  break;
2258  }
2259  }
2260  }
2261  }
2262 }
2263 
2264 static bool sorted_lists_intersect(const int* a, const int* b, int na, int nb)
2265 {
2266  if (!na || !nb) { return false; }
2267  int a_last = a[na-1], b_last = b[nb-1];
2268  if (*b < *a) { goto l2; } // woo-hoo! I always wanted to use a goto! :)
2269 l1:
2270  if (a_last < *b) { return false; }
2271  while (*a < *b) { a++; }
2272  if (*a == *b) { return true; }
2273 l2:
2274  if (b_last < *a) { return false; }
2275  while (*b < *a) { b++; }
2276  if (*a == *b) { return true; }
2277  goto l1;
2278 }
2279 
2281  Array<Element*> &neighbors,
2282  const Array<Element*> *search_set)
2283 {
2285 
2286  Array<int> vert, tmp;
2287 
2288  int *v1, nv1;
2289  if (!elem->ref_type)
2290  {
2291  v1 = element_vertex.GetRow(elem->index);
2292  nv1 = element_vertex.RowSize(elem->index);
2293  }
2294  else // support for non-leaf 'elem', collect vertices of all children
2295  {
2296  Array<const Element*> stack;
2297  stack.Reserve(32);
2298  stack.Append(elem);
2299 
2300  while (stack.Size())
2301  {
2302  const Element* e = stack.Last();
2303  stack.DeleteLast();
2304  if (!e->ref_type)
2305  {
2306  element_vertex.GetRow(e->index, tmp);
2307  vert.Append(tmp);
2308  }
2309  else
2310  {
2311  for (int i = 0; i < 8; i++)
2312  {
2313  if (e->child[i]) { stack.Append(e->child[i]); }
2314  }
2315  }
2316  }
2317  vert.Sort();
2318  vert.Unique();
2319 
2320  v1 = vert.GetData();
2321  nv1 = vert.Size();
2322  }
2323 
2324  if (!search_set) { search_set = &leaf_elements; }
2325 
2326  for (int i = 0; i < search_set->Size(); i++)
2327  {
2328  Element* testme = (*search_set)[i];
2329  if (testme != elem)
2330  {
2331  int *v2 = element_vertex.GetRow(testme->index);
2332  int nv2 = element_vertex.RowSize(testme->index);
2333 
2334  if (sorted_lists_intersect(v1, v2, nv1, nv2))
2335  {
2336  neighbors.Append(testme);
2337  }
2338  }
2339  }
2340 }
2341 
2343  Array<Element*> &expanded,
2344  const Array<Element*> *search_set)
2345 {
2347 
2348  Array<char> vertices(num_vertices);
2349  vertices = 0;
2350 
2351  for (int i = 0; i < elements.Size(); i++)
2352  {
2353  int index = elements[i]->index;
2354  int *v = element_vertex.GetRow(index);
2355  int nv = element_vertex.RowSize(index);
2356  for (int j = 0; j < nv; j++)
2357  {
2358  vertices[v[j]] = 1;
2359  }
2360  }
2361 
2362  if (!search_set)
2363  {
2364  search_set = &leaf_elements;
2365  }
2366 
2367  expanded.SetSize(0);
2368  for (int i = 0; i < search_set->Size(); i++)
2369  {
2370  Element* testme = (*search_set)[i];
2371  int *v = element_vertex.GetRow(testme->index);
2372  int nv = element_vertex.RowSize(testme->index);
2373  for (int j = 0; j < nv; j++)
2374  {
2375  if (vertices[v[j]])
2376  {
2377  expanded.Append(testme);
2378  break;
2379  }
2380  }
2381  }
2382 }
2383 
2384 #ifdef MFEM_DEBUG
2386 {
2387  Array<Element*> neighbors;
2388  FindSetNeighbors(elem_set, &neighbors);
2389 
2390  for (int i = 0; i < neighbors.Size(); i++)
2391  {
2392  elem_set[neighbors[i]->index] = 2;
2393  }
2394 }
2395 #endif
2396 
2397 
2399 
2401 {
2402  point_matrix.SetSize(points[0].dim, np);
2403  for (int i = 0; i < np; i++)
2404  {
2405  for (int j = 0; j < points[0].dim; j++)
2406  {
2407  point_matrix(j, i) = points[i].coord[j];
2408  }
2409  }
2410 }
2411 
2413  Point(0, 0), Point(1, 0), Point(0, 1)
2414 );
2416  Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)
2417 );
2419  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
2420  Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
2421 );
2422 
2424 {
2425  switch (geom)
2426  {
2427  case Geometry::TRIANGLE: return pm_tri_identity;
2428  case Geometry::SQUARE: return pm_quad_identity;
2429  case Geometry::CUBE: return pm_hex_identity;
2430  default:
2431  MFEM_ABORT("unsupported geometry."); throw;
2432  }
2433 }
2434 
2435 void NCMesh::GetPointMatrix(int geom, const char* ref_path, DenseMatrix& matrix)
2436 {
2437  PointMatrix pm = GetGeomIdentity(geom);
2438 
2439  while (*ref_path)
2440  {
2441  int ref_type = *ref_path++;
2442  int child = *ref_path++;
2443 
2444  if (geom == Geometry::CUBE)
2445  {
2446  if (ref_type == 1) // split along X axis
2447  {
2448  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2449  Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
2450 
2451  if (child == 0)
2452  {
2453  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
2454  pm(4), mid45, mid67, pm(7));
2455  }
2456  else if (child == 1)
2457  {
2458  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
2459  mid45, pm(5), pm(6), mid67);
2460  }
2461  }
2462  else if (ref_type == 2) // split along Y axis
2463  {
2464  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2465  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2466 
2467  if (child == 0)
2468  {
2469  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
2470  pm(4), pm(5), mid56, mid74);
2471  }
2472  else if (child == 1)
2473  {
2474  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
2475  mid74, mid56, pm(6), pm(7));
2476  }
2477  }
2478  else if (ref_type == 4) // split along Z axis
2479  {
2480  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2481  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2482 
2483  if (child == 0)
2484  {
2485  pm = PointMatrix(pm(0), pm(1), pm(2), pm(3),
2486  mid04, mid15, mid26, mid37);
2487  }
2488  else if (child == 1)
2489  {
2490  pm = PointMatrix(mid04, mid15, mid26, mid37,
2491  pm(4), pm(5), pm(6), pm(7));
2492  }
2493  }
2494  else if (ref_type == 3) // XY split
2495  {
2496  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2497  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2498  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2499  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2500 
2501  Point midf0(mid23, mid12, mid01, mid30);
2502  Point midf5(mid45, mid56, mid67, mid74);
2503 
2504  if (child == 0)
2505  {
2506  pm = PointMatrix(pm(0), mid01, midf0, mid30,
2507  pm(4), mid45, midf5, mid74);
2508  }
2509  else if (child == 1)
2510  {
2511  pm = PointMatrix(mid01, pm(1), mid12, midf0,
2512  mid45, pm(5), mid56, midf5);
2513  }
2514  else if (child == 2)
2515  {
2516  pm = PointMatrix(midf0, mid12, pm(2), mid23,
2517  midf5, mid56, pm(6), mid67);
2518  }
2519  else if (child == 3)
2520  {
2521  pm = PointMatrix(mid30, midf0, mid23, pm(3),
2522  mid74, midf5, mid67, pm(7));
2523  }
2524  }
2525  else if (ref_type == 5) // XZ split
2526  {
2527  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2528  Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
2529  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2530  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2531 
2532  Point midf1(mid01, mid15, mid45, mid04);
2533  Point midf3(mid23, mid37, mid67, mid26);
2534 
2535  if (child == 0)
2536  {
2537  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
2538  mid04, midf1, midf3, mid37);
2539  }
2540  else if (child == 1)
2541  {
2542  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
2543  midf1, mid15, mid26, midf3);
2544  }
2545  else if (child == 2)
2546  {
2547  pm = PointMatrix(midf1, mid15, mid26, midf3,
2548  mid45, pm(5), pm(6), mid67);
2549  }
2550  else if (child == 3)
2551  {
2552  pm = PointMatrix(mid04, midf1, midf3, mid37,
2553  pm(4), mid45, mid67, pm(7));
2554  }
2555  }
2556  else if (ref_type == 6) // YZ split
2557  {
2558  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2559  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2560  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2561  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2562 
2563  Point midf2(mid12, mid26, mid56, mid15);
2564  Point midf4(mid30, mid04, mid74, mid37);
2565 
2566  if (child == 0)
2567  {
2568  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
2569  mid04, mid15, midf2, midf4);
2570  }
2571  else if (child == 1)
2572  {
2573  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
2574  midf4, midf2, mid26, mid37);
2575  }
2576  else if (child == 2)
2577  {
2578  pm = PointMatrix(mid04, mid15, midf2, midf4,
2579  pm(4), pm(5), mid56, mid74);
2580  }
2581  else if (child == 3)
2582  {
2583  pm = PointMatrix(midf4, midf2, mid26, mid37,
2584  mid74, mid56, pm(6), pm(7));
2585  }
2586  }
2587  else if (ref_type == 7) // full isotropic refinement
2588  {
2589  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2590  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2591  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2592  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2593  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2594  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2595 
2596  Point midf0(mid23, mid12, mid01, mid30);
2597  Point midf1(mid01, mid15, mid45, mid04);
2598  Point midf2(mid12, mid26, mid56, mid15);
2599  Point midf3(mid23, mid37, mid67, mid26);
2600  Point midf4(mid30, mid04, mid74, mid37);
2601  Point midf5(mid45, mid56, mid67, mid74);
2602 
2603  Point midel(midf1, midf3);
2604 
2605  if (child == 0)
2606  {
2607  pm = PointMatrix(pm(0), mid01, midf0, mid30,
2608  mid04, midf1, midel, midf4);
2609  }
2610  else if (child == 1)
2611  {
2612  pm = PointMatrix(mid01, pm(1), mid12, midf0,
2613  midf1, mid15, midf2, midel);
2614  }
2615  else if (child == 2)
2616  {
2617  pm = PointMatrix(midf0, mid12, pm(2), mid23,
2618  midel, midf2, mid26, midf3);
2619  }
2620  else if (child == 3)
2621  {
2622  pm = PointMatrix(mid30, midf0, mid23, pm(3),
2623  midf4, midel, midf3, mid37);
2624  }
2625  else if (child == 4)
2626  {
2627  pm = PointMatrix(mid04, midf1, midel, midf4,
2628  pm(4), mid45, midf5, mid74);
2629  }
2630  else if (child == 5)
2631  {
2632  pm = PointMatrix(midf1, mid15, midf2, midel,
2633  mid45, pm(5), mid56, midf5);
2634  }
2635  else if (child == 6)
2636  {
2637  pm = PointMatrix(midel, midf2, mid26, midf3,
2638  midf5, mid56, pm(6), mid67);
2639  }
2640  else if (child == 7)
2641  {
2642  pm = PointMatrix(midf4, midel, midf3, mid37,
2643  mid74, midf5, mid67, pm(7));
2644  }
2645  }
2646  }
2647  else if (geom == Geometry::SQUARE)
2648  {
2649  if (ref_type == 1) // X split
2650  {
2651  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2652 
2653  if (child == 0)
2654  {
2655  pm = PointMatrix(pm(0), mid01, mid23, pm(3));
2656  }
2657  else if (child == 1)
2658  {
2659  pm = PointMatrix(mid01, pm(1), pm(2), mid23);
2660  }
2661  }
2662  else if (ref_type == 2) // Y split
2663  {
2664  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2665 
2666  if (child == 0)
2667  {
2668  pm = PointMatrix(pm(0), pm(1), mid12, mid30);
2669  }
2670  else if (child == 1)
2671  {
2672  pm = PointMatrix(mid30, mid12, pm(2), pm(3));
2673  }
2674  }
2675  else if (ref_type == 3) // iso split
2676  {
2677  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2678  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2679  Point midel(mid01, mid23);
2680 
2681  if (child == 0)
2682  {
2683  pm = PointMatrix(pm(0), mid01, midel, mid30);
2684  }
2685  else if (child == 1)
2686  {
2687  pm = PointMatrix(mid01, pm(1), mid12, midel);
2688  }
2689  else if (child == 2)
2690  {
2691  pm = PointMatrix(midel, mid12, pm(2), mid23);
2692  }
2693  else if (child == 3)
2694  {
2695  pm = PointMatrix(mid30, midel, mid23, pm(3));
2696  }
2697  }
2698  }
2699  else if (geom == Geometry::TRIANGLE)
2700  {
2701  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
2702 
2703  if (child == 0)
2704  {
2705  pm = PointMatrix(pm(0), mid01, mid20);
2706  }
2707  else if (child == 1)
2708  {
2709  pm = PointMatrix(mid01, pm(1), mid12);
2710  }
2711  else if (child == 2)
2712  {
2713  pm = PointMatrix(mid20, mid12, pm(2));
2714  }
2715  else if (child == 3)
2716  {
2717  pm = PointMatrix(mid01, mid12, mid20);
2718  }
2719  }
2720  }
2721 
2722  // write the points to the matrix
2723  for (int i = 0; i < pm.np; i++)
2724  {
2725  for (int j = 0; j < pm(i).dim; j++)
2726  {
2727  matrix(j, i) = pm(i).coord[j];
2728  }
2729  }
2730 }
2731 
2733 {
2734  coarse_elements.SetSize(leaf_elements.Size());
2735  coarse_elements.SetSize(0);
2736 
2737  for (int i = 0; i < leaf_elements.Size(); i++)
2738  {
2739  Element* e = leaf_elements[i];
2740  if (!IsGhost(e)) { coarse_elements.Append(e); }
2741  }
2742 
2743  transforms.embeddings.DeleteAll();
2744 }
2745 
2746 void NCMesh::TraverseRefinements(Element* elem, int coarse_index,
2747  std::string &ref_path, RefPathMap &map)
2748 {
2749  if (!elem->ref_type)
2750  {
2751  int &matrix = map[ref_path];
2752  if (!matrix) { matrix = map.size(); }
2753 
2754  Embedding &emb = transforms.embeddings[elem->index];
2755  emb.parent = coarse_index;
2756  emb.matrix = matrix - 1;
2757  }
2758  else
2759  {
2760  ref_path.push_back(elem->ref_type);
2761  ref_path.push_back(0);
2762 
2763  for (int i = 0; i < 8; i++)
2764  {
2765  if (elem->child[i])
2766  {
2767  ref_path[ref_path.length()-1] = i;
2768  TraverseRefinements(elem->child[i], coarse_index, ref_path, map);
2769  }
2770  }
2771  ref_path.resize(ref_path.length()-2);
2772  }
2773 }
2774 
2776 {
2777  MFEM_VERIFY(coarse_elements.Size() || !leaf_elements.Size(),
2778  "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
2779  " and Refine().");
2780 
2781  if (!transforms.embeddings.Size())
2782  {
2783  transforms.embeddings.SetSize(leaf_elements.Size());
2784 
2785  std::string ref_path;
2786  ref_path.reserve(100);
2787 
2788  RefPathMap map;
2789  map[ref_path] = 1; // identity
2790 
2791  for (int i = 0; i < coarse_elements.Size(); i++)
2792  {
2793  TraverseRefinements(coarse_elements[i], i, ref_path, map);
2794  }
2795 
2796  MFEM_ASSERT(root_elements.Size(), "");
2797  int geom = root_elements[0]->geom;
2798  const PointMatrix &identity = GetGeomIdentity(geom);
2799 
2800  transforms.point_matrices.SetSize(Dim, identity.np, map.size());
2801 
2802  // calculate the point matrices
2803  for (RefPathMap::iterator it = map.begin(); it != map.end(); ++it)
2804  {
2805  GetPointMatrix(geom, it->first.c_str(),
2806  transforms.point_matrices(it->second-1));
2807  }
2808  }
2809  return transforms;
2810 }
2811 
2813 {
2814  MFEM_VERIFY(transforms.embeddings.Size() || !leaf_elements.Size(),
2815  "GetDerefinementTransforms() must be preceded by Derefine().");
2816 
2818  {
2819  std::map<int, int> mat_no;
2820  mat_no[0] = 1; // identity
2821 
2822  // assign numbers to the different matrices used
2823  for (int i = 0; i < transforms.embeddings.Size(); i++)
2824  {
2825  int code = transforms.embeddings[i].matrix;
2826  if (code)
2827  {
2828  int &matrix = mat_no[code];
2829  if (!matrix) { matrix = mat_no.size(); }
2830  transforms.embeddings[i].matrix = matrix - 1;
2831  }
2832  }
2833 
2834  MFEM_ASSERT(root_elements.Size(), "");
2835  int geom = root_elements[0]->geom;
2836  const PointMatrix &identity = GetGeomIdentity(geom);
2837 
2838  transforms.point_matrices.SetSize(Dim, identity.np, mat_no.size());
2839 
2840  std::map<int, int>::iterator it;
2841  for (it = mat_no.begin(); it != mat_no.end(); ++it)
2842  {
2843  char path[3];
2844  int code = it->first;
2845  path[0] = code >> 3; // ref_type (see SetDerefMatrixCodes())
2846  path[1] = code & 7; // child
2847  path[2] = 0;
2848 
2849  GetPointMatrix(geom, path, transforms.point_matrices(it->second-1));
2850  }
2851  }
2852  return transforms;
2853 }
2854 
2856 {
2857  coarse_elements.DeleteAll();
2858  transforms.embeddings.DeleteAll();
2860 }
2861 
2862 
2864 
2866 {
2867  MFEM_ASSERT(node != NULL && node->p1 != node->p2, "Invalid edge node.");
2868 
2869  Node *n1 = nodes.Peek(node->p1);
2870  Node *n2 = nodes.Peek(node->p2);
2871 
2872  if ((n2->p1 != n2->p2) && (n1->id == n2->p1 || n1->id == n2->p2))
2873  {
2874  // n1 is parent of n2:
2875  // (n1)--(n)--(n2)------(*) or (*)------(n2)--(n)--(n1)
2876  if (n2->edge) { return n2; }
2877  else { return GetEdgeMaster(n2); }
2878  }
2879 
2880  if ((n1->p1 != n1->p2) && (n2->id == n1->p1 || n2->id == n1->p2))
2881  {
2882  // n2 is parent of n1:
2883  // (n2)--(n)--(n1)------(*) or (*)------(n1)--(n)--(n2)
2884  if (n1->edge) { return n1; }
2885  else { return GetEdgeMaster(n1); }
2886  }
2887 
2888  return NULL;
2889 }
2890 
2891 int NCMesh::GetEdgeMaster(int v1, int v2) const
2892 {
2893  Node* node = nodes.Peek(vertex_nodeId[v1], vertex_nodeId[v2]);
2894  MFEM_ASSERT(node->edge != NULL, "(v1, v2) is not an edge.");
2895 
2896  Node* master = GetEdgeMaster(node);
2897  return master ? master->edge->index : -1;
2898 }
2899 
2900 int NCMesh::GetElementDepth(int i) const
2901 {
2902  Element* elem = leaf_elements[i];
2903  int depth = 0;
2904  while (elem->parent)
2905  {
2906  elem = elem->parent;
2907  depth++;
2908  }
2909  return depth;
2910 }
2911 
2912 void NCMesh::find_face_nodes(const Face *face, Node* node[4])
2913 {
2914  // Obtain face nodes from one of its elements (note that face->p1, p2, p3
2915  // cannot be used directly since they are not in order and p4 is missing).
2916 
2917  Element* elem = face->elem[0];
2918  if (!elem) { elem = face->elem[1]; }
2919  MFEM_ASSERT(elem, "Face has no elements?");
2920 
2921  int f = find_hex_face(find_node(elem, face->p1),
2922  find_node(elem, face->p2),
2923  find_node(elem, face->p3));
2924 
2925  const int* fv = GI[Geometry::CUBE].faces[f];
2926  for (int i = 0; i < 4; i++)
2927  {
2928  node[i] = elem->node[fv[i]];
2929  }
2930 }
2931 
2932 void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
2933  Array<int> &bdr_vertices, Array<int> &bdr_edges)
2934 {
2935  bdr_vertices.SetSize(0);
2936  bdr_edges.SetSize(0);
2937 
2938  if (Dim == 3)
2939  {
2940  GetFaceList();
2941  for (int i = 0; i < boundary_faces.Size(); i++)
2942  {
2943  Face* face = boundary_faces[i];
2944  if (bdr_attr_is_ess[face->attribute - 1])
2945  {
2946  Node* node[4];
2947  find_face_nodes(face, node);
2948 
2949  for (int j = 0; j < 4; j++)
2950  {
2951  bdr_vertices.Append(node[j]->vertex->index);
2952 
2953  Node* edge = nodes.Peek(node[j], node[(j+1) % 4]);
2954  MFEM_ASSERT(edge && edge->edge, "Edge not found.");
2955  bdr_edges.Append(edge->edge->index);
2956 
2957  while ((edge = GetEdgeMaster(edge)) != NULL)
2958  {
2959  // append master edges that may not be accessible from any
2960  // boundary element, this happens in 3D in re-entrant corners
2961  bdr_edges.Append(edge->edge->index);
2962  }
2963  }
2964  }
2965  }
2966  }
2967  else if (Dim == 2)
2968  {
2969  GetEdgeList();
2970  for (int i = 0; i < boundary_edges.Size(); i++)
2971  {
2972  Node* edge = boundary_edges[i];
2973  if (bdr_attr_is_ess[edge->edge->attribute - 1])
2974  {
2975  bdr_vertices.Append(nodes.Peek(edge->p1)->vertex->index);
2976  bdr_vertices.Append(nodes.Peek(edge->p2)->vertex->index);
2977  }
2978  }
2979  }
2980 
2981  bdr_vertices.Sort();
2982  bdr_vertices.Unique();
2983 
2984  bdr_edges.Sort();
2985  bdr_edges.Unique();
2986 }
2987 
2988 int NCMesh::EdgeSplitLevel(Node *v1, Node *v2) const
2989 {
2990  Node* mid = nodes.Peek(v1, v2);
2991  if (!mid || !mid->vertex) { return 0; }
2992  return 1 + std::max(EdgeSplitLevel(v1, mid), EdgeSplitLevel(mid, v2));
2993 }
2994 
2995 void NCMesh::FaceSplitLevel(Node* v1, Node* v2, Node* v3, Node* v4,
2996  int& h_level, int& v_level) const
2997 {
2998  int hl1, hl2, vl1, vl2;
2999  Node* mid[4];
3000 
3001  switch (FaceSplitType(v1, v2, v3, v4, mid))
3002  {
3003  case 0: // not split
3004  h_level = v_level = 0;
3005  break;
3006 
3007  case 1: // vertical
3008  FaceSplitLevel(v1, mid[0], mid[2], v4, hl1, vl1);
3009  FaceSplitLevel(mid[0], v2, v3, mid[2], hl2, vl2);
3010  h_level = std::max(hl1, hl2);
3011  v_level = std::max(vl1, vl2) + 1;
3012  break;
3013 
3014  default: // horizontal
3015  FaceSplitLevel(v1, v2, mid[1], mid[3], hl1, vl1);
3016  FaceSplitLevel(mid[3], mid[1], v3, v4, hl2, vl2);
3017  h_level = std::max(hl1, hl2) + 1;
3018  v_level = std::max(vl1, vl2);
3019  }
3020 }
3021 
3022 static int max8(int a, int b, int c, int d, int e, int f, int g, int h)
3023 {
3024  return std::max(std::max(std::max(a, b), std::max(c, d)),
3025  std::max(std::max(e, f), std::max(g, h)));
3026 }
3027 
3028 void NCMesh::CountSplits(Element* elem, int splits[3]) const
3029 {
3030  Node** node = elem->node;
3031  GeomInfo& gi = GI[(int) elem->geom];
3032 
3033  int elevel[12];
3034  for (int i = 0; i < gi.ne; i++)
3035  {
3036  const int* ev = gi.edges[i];
3037  elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
3038  }
3039 
3040  if (elem->geom == Geometry::CUBE)
3041  {
3042  int flevel[6][2];
3043  for (int i = 0; i < gi.nf; i++)
3044  {
3045  const int* fv = gi.faces[i];
3046  FaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]],
3047  flevel[i][1], flevel[i][0]);
3048  }
3049 
3050  splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
3051  elevel[0], elevel[2], elevel[4], elevel[6]);
3052 
3053  splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
3054  elevel[1], elevel[3], elevel[5], elevel[7]);
3055 
3056  splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
3057  elevel[8], elevel[9], elevel[10], elevel[11]);
3058  }
3059  else if (elem->geom == Geometry::SQUARE)
3060  {
3061  splits[0] = std::max(elevel[0], elevel[2]);
3062  splits[1] = std::max(elevel[1], elevel[3]);
3063  }
3064  else if (elem->geom == Geometry::TRIANGLE)
3065  {
3066  splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
3067  splits[1] = splits[0];
3068  }
3069  else
3070  {
3071  MFEM_ABORT("Unsupported element geometry.");
3072  }
3073 }
3074 
3075 void NCMesh::GetLimitRefinements(Array<Refinement> &refinements, int max_level)
3076 {
3077  for (int i = 0; i < leaf_elements.Size(); i++)
3078  {
3079  if (IsGhost(leaf_elements[i])) { break; }
3080 
3081  int splits[3];
3082  CountSplits(leaf_elements[i], splits);
3083 
3084  char ref_type = 0;
3085  for (int k = 0; k < Dim; k++)
3086  {
3087  if (splits[k] > max_level)
3088  {
3089  ref_type |= (1 << k);
3090  }
3091  }
3092 
3093  if (ref_type)
3094  {
3095  if (Iso)
3096  {
3097  // iso meshes should only be modified by iso refinements
3098  ref_type = 7;
3099  }
3100  refinements.Append(Refinement(i, ref_type));
3101  }
3102  }
3103 }
3104 
3105 void NCMesh::LimitNCLevel(int max_nc_level)
3106 {
3107  MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
3108 
3109  while (1)
3110  {
3111  Array<Refinement> refinements;
3112  GetLimitRefinements(refinements, max_nc_level);
3113 
3114  if (!refinements.Size()) { break; }
3115 
3116  Refine(refinements);
3117  }
3118 }
3119 
3120 void NCMesh::PrintVertexParents(std::ostream &out) const
3121 {
3122  // count vertices with parents
3123  int nv = 0;
3124  for (HashTable<Node>::Iterator it(nodes); it; ++it)
3125  {
3126  if (it->vertex && it->p1 != it->p2) { nv++; }
3127  }
3128  out << nv << "\n";
3129 
3130  // print the relations
3131  for (HashTable<Node>::Iterator it(nodes); it; ++it)
3132  {
3133  if (it->vertex && it->p1 != it->p2)
3134  {
3135  Node *p1 = nodes.Peek(it->p1);
3136  Node *p2 = nodes.Peek(it->p2);
3137 
3138  MFEM_ASSERT(p1 && p1->vertex, "");
3139  MFEM_ASSERT(p2 && p2->vertex, "");
3140 
3141  out << it->vertex->index << " "
3142  << p1->vertex->index << " "
3143  << p2->vertex->index << "\n";
3144  }
3145  }
3146 }
3147 
3148 void NCMesh::LoadVertexParents(std::istream &input)
3149 {
3150  int nv;
3151  input >> nv;
3152  while (nv--)
3153  {
3154  int id, p1, p2;
3155  input >> id >> p1 >> p2;
3156  MFEM_VERIFY(input, "problem reading vertex parents.");
3157 
3158  Node* node = nodes.Peek(id);
3159  MFEM_VERIFY(node, "vertex " << id << " not found.");
3160  MFEM_VERIFY(nodes.Peek(p1), "parent " << p1 << " not found.");
3161  MFEM_VERIFY(nodes.Peek(p2), "parent " << p2 << " not found.");
3162 
3163  // assign new parents for the node
3164  nodes.Reparent(node, p1, p2);
3165  }
3166 }
3167 
3169 {
3170  for (int i = 0; i < vertices.Size(); i++)
3171  {
3172  Node* node = nodes.Peek(i);
3173  MFEM_ASSERT(node && node->vertex, "");
3174 
3175  const double* pos = vertices[i]();
3176  memcpy(node->vertex->pos, pos, sizeof(node->vertex->pos));
3177  }
3178 }
3179 
3180 static int ref_type_num_children[8] = {0, 2, 2, 4, 2, 4, 4, 8 };
3181 
3182 int NCMesh::PrintElements(std::ostream &out, Element* elem,
3183  int &coarse_id) const
3184 {
3185  if (elem->ref_type)
3186  {
3187  int child_id[8], nch = 0;
3188  for (int i = 0; i < 8; i++)
3189  {
3190  if (elem->child[i])
3191  {
3192  child_id[nch++] = PrintElements(out, elem->child[i], coarse_id);
3193  }
3194  }
3195  MFEM_ASSERT(nch == ref_type_num_children[(int) elem->ref_type], "");
3196 
3197  out << (int) elem->ref_type;
3198  for (int i = 0; i < nch; i++)
3199  {
3200  out << " " << child_id[i];
3201  }
3202  out << "\n";
3203  return coarse_id++; // return new id for this coarse element
3204  }
3205  else
3206  {
3207  return elem->index;
3208  }
3209 }
3210 
3211 void NCMesh::PrintCoarseElements(std::ostream &out) const
3212 {
3213  // print the number of non-leaf elements
3214  int ne = 0;
3215  for (int i = 0; i < root_elements.Size(); i++)
3216  {
3217  ne += CountElements(root_elements[i]);
3218  }
3219  out << (ne - leaf_elements.Size()) << "\n";
3220 
3221  // print the hierarchy recursively
3222  int coarse_id = leaf_elements.Size();
3223  for (int i = 0; i < root_elements.Size(); i++)
3224  {
3225  PrintElements(out, root_elements[i], coarse_id);
3226  }
3227 }
3228 
3229 void NCMesh::LoadCoarseElements(std::istream &input)
3230 {
3231  int ne;
3232  input >> ne;
3233 
3234  Array<Element*> coarse, leaves;
3235  coarse.Reserve(ne);
3236  leaf_elements.Copy(leaves);
3237  int nleaf = leaves.Size();
3238  bool iso = true;
3239 
3240  // load the coarse elements
3241  while (ne--)
3242  {
3243  int ref_type;
3244  input >> ref_type;
3245 
3246  Element* elem = new Element(0, 0);
3247  elem->ref_type = ref_type;
3248 
3249  if (Dim == 3 && ref_type != 7) { iso = false; }
3250 
3251  // load child IDs and convert to Element*
3252  int nch = ref_type_num_children[ref_type];
3253  for (int i = 0, id; i < nch; i++)
3254  {
3255  input >> id;
3256  MFEM_VERIFY(id >= 0, "");
3257  MFEM_VERIFY(id < nleaf || id - nleaf < coarse.Size(),
3258  "coarse element cannot be referenced before it is "
3259  "defined (id=" << id << ").");
3260 
3261  Element* &child = (id < nleaf) ? leaves[id] : coarse[id - nleaf];
3262 
3263  MFEM_VERIFY(child, "element " << id << " cannot have two parents.");
3264  elem->child[i] = child;
3265  child->parent = elem;
3266  child = NULL; // make sure the child can't be used again
3267 
3268  if (!i) // copy geom and attribute from first child
3269  {
3270  elem->geom = elem->child[i]->geom;
3271  elem->attribute = elem->child[i]->attribute;
3272  }
3273  }
3274 
3275  // keep a list of coarse elements (and their IDs, implicitly)
3276  coarse.Append(elem);
3277  }
3278 
3279  // elements that have no parents are the original 'root_elements'
3280  root_elements.SetSize(0);
3281  for (int i = 0; i < coarse.Size(); i++)
3282  {
3283  if (coarse[i]) { root_elements.Append(coarse[i]); }
3284  }
3285  for (int i = 0; i < leaves.Size(); i++)
3286  {
3287  if (leaves[i]) { root_elements.Append(leaves[i]); }
3288  }
3289 
3290  // set the Iso flag (must be false if there are 3D aniso refinements)
3291  Iso = iso;
3292 }
3293 
3295 {
3296  int n = 1;
3297  if (elem->ref_type)
3298  {
3299  for (int i = 0; i < 8; i++)
3300  {
3301  if (elem->child[i]) { n += CountElements(elem->child[i]); }
3302  }
3303  }
3304  return n;
3305 }
3306 
3308 {
3309  int pmsize = 0;
3310  if (slaves.size())
3311  {
3312  const DenseMatrix &pm = slaves[0].point_matrix;
3313  pmsize = pm.Width() * pm.Height() * sizeof(double);
3314  }
3315 
3316  return conforming.capacity() * sizeof(MeshId) +
3317  masters.capacity() * sizeof(Master) +
3318  slaves.capacity() * sizeof(Slave) +
3319  slaves.size() * pmsize;
3320 }
3321 
3322 void NCMesh::CountObjects(int &nelem, int &nvert, int &nedges) const
3323 {
3324  nelem = nvert = nedges = 0;
3325  for (int i = 0; i < root_elements.Size(); i++)
3326  {
3327  nelem += CountElements(root_elements[i]);
3328  }
3329  for (HashTable<Node>::Iterator it(nodes); it; ++it)
3330  {
3331  if (it->vertex) { nvert++; }
3332  if (it->edge) { nedges++; }
3333  }
3334 }
3335 
3337 {
3338  int nelem, nvert, nedges;
3339  CountObjects(nelem, nvert, nedges);
3340 
3341  return nelem * sizeof(Element) +
3342  nvert * sizeof(Vertex) +
3343  nedges * sizeof(Edge) +
3344  nodes.MemoryUsage() +
3345  faces.MemoryUsage() +
3346  root_elements.MemoryUsage() +
3347  leaf_elements.MemoryUsage() +
3351  boundary_faces.MemoryUsage() +
3352  boundary_edges.MemoryUsage() +
3354  ref_stack.MemoryUsage() +
3355  coarse_elements.MemoryUsage() +
3356  sizeof(*this);
3357 }
3358 
3360 {
3361  int nelem, nvert, nedges;
3362  CountObjects(nelem, nvert, nedges);
3363 
3364  std::cout << nelem * sizeof(Element) << " elements\n"
3365  << nvert * sizeof(Vertex) << " vertices\n"
3366  << nedges * sizeof(Edge) << " edges\n";
3367 
3368  nodes.PrintMemoryDetail(); std::cout << " nodes\n";
3369  faces.PrintMemoryDetail(); std::cout << " faces\n";
3370 
3371  std::cout << root_elements.MemoryUsage() << " root_elements\n"
3372  << leaf_elements.MemoryUsage() << " leaf_elements\n"
3373  << vertex_nodeId.MemoryUsage() << " vertex_nodeId\n"
3374  << face_list.MemoryUsage() << " face_list\n"
3375  << edge_list.MemoryUsage() << " edge_list\n"
3376  << boundary_faces.MemoryUsage() << " boundary_faces\n"
3377  << boundary_edges.MemoryUsage() << " boundary_edges\n"
3378  << element_vertex.MemoryUsage() << " element_vertex\n"
3379  << ref_stack.MemoryUsage() << " ref_stack\n"
3380  << coarse_elements.MemoryUsage() << " coarse_elements\n"
3381  << sizeof(*this) << " NCMesh" << std::endl;
3382 }
3383 
3384 #ifdef MFEM_DEBUG
3386 {
3387  for (int i = 0; i < leaf_elements.Size(); i++)
3388  {
3389  Element* elem = leaf_elements[i];
3390  for (int j = 0; j < Dim; j++)
3391  {
3392  double sum = 0.0;
3393  int count = 0;
3394  for (int k = 0; k < 8; k++)
3395  {
3396  if (elem->node[k])
3397  {
3398  sum += elem->node[k]->vertex->pos[j];
3399  count++;
3400  }
3401  }
3402  std::cout << sum / count << " ";
3403  }
3404  std::cout << "\n";
3405  }
3406 }
3407 #endif
3408 
3409 } // namespace mfem
NCList face_list
lazy-initialized list of faces, see GetFaceList
Definition: ncmesh.hpp:434
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
Definition: ncmesh.hpp:435
NCMesh(const Mesh *mesh, std::istream *vertex_parents=NULL)
Definition: ncmesh.cpp:68
int Size() const
Logical size of the array.
Definition: array.hpp:109
Element * NewHexahedron(Node *n0, Node *n1, Node *n2, Node *n3, Node *n4, Node *n5, Node *n6, Node *n7, int attr, int fattr0, int fattr1, int fattr2, int fattr3, int fattr4, int fattr5)
Definition: ncmesh.cpp:491
CoarseFineTransformations transforms
storage for data returned by Get[De]RefinementTransforms()
Definition: ncmesh.hpp:668
void FindNeighbors(const Element *elem, Array< Element * > &neighbors, const Array< Element * > *search_set=NULL)
Definition: ncmesh.cpp:2280
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:529
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:36
void Unique()
Definition: array.hpp:197
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
Definition: mesh.cpp:3470
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:2812
int index
edge number in the Mesh
Definition: ncmesh.hpp:322
Array< Element * > coarse_elements
state of leaf_elements before Refine(), set by MarkCoarseLevel()
Definition: ncmesh.hpp:671
virtual int GetNEdges() const =0
Iterator over items contained in the HashTable.
Definition: hash.hpp:160
int PrintElements(std::ostream &out, Element *elem, int &coarse_id) const
Definition: ncmesh.cpp:3182
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
static const int NumGeom
Definition: geom.hpp:34
void PrintVertexParents(std::ostream &out) const
I/O: Print the &quot;vertex_parents&quot; section of the mesh file (ver. &gt;= 1.1).
Definition: ncmesh.cpp:3120
void GetPointMatrix(int geom, const char *ref_path, DenseMatrix &matrix)
Definition: ncmesh.cpp:2435
Node * PeekAltParents(Node *v1, Node *v2)
Definition: ncmesh.cpp:424
virtual void LimitNCLevel(int max_nc_level)
Definition: ncmesh.cpp:3105
void PrintMemoryDetail() const
Definition: ncmesh.cpp:3359
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:499
int index
element number in the Mesh, -1 if refined
Definition: ncmesh.hpp:396
const Geometry::Type geom
virtual int GetNumGhosts() const
Definition: ncmesh.hpp:451
void DebugLeafOrder() const
Print the space-filling curve formed by the sequence of leaf elements.
Definition: ncmesh.cpp:3385
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
int SizeK() const
Definition: densemat.hpp:595
T * GetData()
Returns the data.
Definition: array.hpp:91
int attribute
boundary element attribute, -1 if internal edge (2D)
Definition: ncmesh.hpp:321
void ClearTransforms()
Free all internal data created by the above three functions.
Definition: ncmesh.cpp:2855
bool Boundary() const
Definition: ncmesh.hpp:327
Table derefinements
possible derefinements, see GetDerefinementTable
Definition: ncmesh.hpp:467
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
void InitDerefTransforms()
Definition: ncmesh.cpp:1436
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
void RefElementNodes(Element *elem)
Definition: ncmesh.cpp:306
char geom
Geometry::Type of the element.
Definition: ncmesh.hpp:393
static int find_node(Element *elem, Node *node)
Definition: ncmesh.cpp:1745
void SetDerefMatrixCodes(Element *parent, Array< Element * > &coarse)
Definition: ncmesh.cpp:1451
int FaceSplitType(Node *v1, Node *v2, Node *v3, Node *v4, Node *mid[4]=NULL) const
Definition: ncmesh.cpp:1717
void SetSize(int i, int j, int k)
Definition: densemat.hpp:597
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:496
virtual void ElementSharesFace(Element *elem, Face *face)
Definition: ncmesh.hpp:542
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
Definition: ncmesh.cpp:2059
const Element * GetFace(int i) const
Definition: mesh.hpp:547
static PointMatrix pm_tri_identity
Definition: ncmesh.hpp:654
int index
vertex number in the Mesh
Definition: ncmesh.hpp:309
void NeighborExpand(const Array< Element * > &elements, Array< Element * > &expanded, const Array< Element * > *search_set=NULL)
Definition: ncmesh.cpp:2342
void CountObjects(int &nelem, int &nvert, int &nedges) const
Definition: ncmesh.cpp:3322
std::vector< MeshId > conforming
Definition: ncmesh.hpp:169
Data type quadrilateral element.
void GetLimitRefinements(Array< Refinement > &refinements, int max_level)
Definition: ncmesh.cpp:3075
void Delete(ItemT *item)
Remove an item from the hash table and also delete the item itself.
Definition: hash.hpp:438
Element * parent
parent element, NULL if this is a root element
Definition: ncmesh.hpp:404
Vertex * NewVertex(Node *v1, Node *v2)
Definition: ncmesh.cpp:570
bool NodeSetY1(Node *node, Node **n)
Definition: ncmesh.cpp:626
int spaceDim
dimensions of the elements and the vertex coordinates
Definition: ncmesh.hpp:276
Array< int > vertex_nodeId
Definition: ncmesh.hpp:432
int attribute
boundary element attribute, -1 if internal face
Definition: ncmesh.hpp:368
int GetGeometryType() const
Definition: element.hpp:47
int index
Mesh element number.
Definition: ncmesh.hpp:35
Element * element
NCMesh::Element containing this vertex/edge/face.
Definition: ncmesh.hpp:135
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Definition: ncmesh.hpp:181
virtual void OnMeshUpdated(Mesh *mesh)
Definition: ncmesh.cpp:1671
int Width() const
Returns the number of TYPE II elements (after Finalize() is called).
Definition: table.cpp:285
Element * NewQuadrilateral(Node *n0, Node *n1, Node *n2, Node *n3, int attr, int eattr0, int eattr1, int eattr2, int eattr3)
Definition: ncmesh.cpp:519
void DebugNeighbors(Array< char > &elem_set)
Definition: ncmesh.cpp:2385
void CollectEdgeVertices(Node *v0, Node *v1, Array< int > &indices)
Definition: ncmesh.cpp:2083
virtual bool IsGhost(const Element *elem) const
Definition: ncmesh.hpp:450
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3039
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
Definition: segment.cpp:40
Data type hexahedron element.
Definition: hexahedron.hpp:22
int dim
Definition: ex3.cpp:47
virtual const int * GetEdgeVertices(int) const =0
int num_vertices
width of the table
Definition: ncmesh.hpp:441
int index
face number in the Mesh
Definition: ncmesh.hpp:369
const Table & GetDerefinementTable()
Definition: ncmesh.cpp:1356
void MakeFromList(int nrows, const Array< Connection > &list)
Definition: table.cpp:264
Table element_vertex
leaf-element to vertex table, see FindSetNeighbors
Definition: ncmesh.hpp:440
DenseTensor point_matrices
matrices for IsoparametricTransformation
Definition: ncmesh.hpp:53
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
Definition: operator.hpp:35
int local
local number within &#39;element&#39;
Definition: ncmesh.hpp:134
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:376
Element * elem[2]
up to 2 elements sharing the face
Definition: ncmesh.hpp:370
static void UnrefEdge(Node *node, HashTable< Node > &nodes)
Definition: ncmesh.cpp:284
void GetMatrix(DenseMatrix &point_matrix) const
Definition: ncmesh.cpp:2400
void Clear()
Definition: table.cpp:346
void CheckAnisoFace(Node *v1, Node *v2, Node *v3, Node *v4, Node *mid12, Node *mid34, int level=0)
Definition: ncmesh.cpp:673
virtual void Derefine(const Array< int > &derefs)
Definition: ncmesh.cpp:1400
void RegisterFaces(Element *elem, int *fattr=NULL)
Definition: ncmesh.cpp:399
std::vector< Master > masters
Definition: ncmesh.hpp:170
virtual void UpdateVertices()
update Vertex::index and vertex_nodeId
Definition: ncmesh.cpp:1469
void RegisterElement(Element *e)
Definition: ncmesh.cpp:377
int ReorderFacePointMat(Node *v0, Node *v1, Node *v2, Node *v3, Element *elem, DenseMatrix &mat) const
Definition: ncmesh.cpp:1795
long MemoryUsage() const
Definition: table.cpp:379
Element * CopyHierarchy(Element *elem)
Definition: ncmesh.cpp:198
void RefineElement(Element *elem, char ref_type)
Definition: ncmesh.cpp:756
void FaceSplitLevel(Node *v1, Node *v2, Node *v3, Node *v4, int &h_level, int &v_level) const
Definition: ncmesh.cpp:2995
Element(int geom, int attr)
Definition: ncmesh.cpp:478
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:122
Array< Embedding > embeddings
fine element positions in their parents
Definition: ncmesh.hpp:54
virtual void BuildEdgeList()
Definition: ncmesh.cpp:1987
const CoarseFineTransformations & GetRefinementTransforms()
Definition: ncmesh.cpp:2775
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:50
void LoadCoarseElements(std::istream &input)
I/O: Load the element refinement hierarchy from a mesh file.
Definition: ncmesh.cpp:3229
Identifies a vertex/edge/face in both Mesh and NCMesh.
Definition: ncmesh.hpp:131
Data type triangle element.
Definition: triangle.hpp:23
void MarkCoarseLevel()
Definition: ncmesh.cpp:2732
const Element * GetElement(int i) const
Definition: mesh.hpp:539
virtual void ElementSharesEdge(Element *elem, Edge *edge)
Definition: ncmesh.hpp:541
long MemoryUsage() const
Definition: array.hpp:221
void Sort()
Sorts the array. This requires operator&lt; to be defined for T.
Definition: array.hpp:189
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:86
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
Definition: ncmesh.hpp:79
bool NodeSetX1(Node *node, Node **n)
Definition: ncmesh.cpp:620
int Dimension() const
Definition: mesh.hpp:523
virtual ~NCMesh()
Definition: ncmesh.cpp:254
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:2932
static GeomInfo GI[Geometry::NumGeom]
Definition: ncmesh.hpp:712
void DeleteHierarchy(Element *elem)
Definition: ncmesh.cpp:224
Vertex * vertex
Definition: ncmesh.hpp:342
int SpaceDimension() const
Definition: mesh.hpp:524
bool Iso
true if the mesh only contains isotropic refinements
Definition: ncmesh.hpp:277
void TraverseEdge(Node *v0, Node *v1, double t0, double t1, int flags, int level)
Definition: ncmesh.cpp:1951
void DerefineElement(Element *elem)
Definition: ncmesh.cpp:1211
std::map< std::string, int > RefPathMap
Definition: ncmesh.hpp:662
Array< Node * > boundary_edges
subset of all edges, set by BuildEdgeList
Definition: ncmesh.hpp:438
void FindSetNeighbors(const Array< char > &elem_set, Array< Element * > *neighbors, Array< char > *neighbor_set=NULL)
Definition: ncmesh.cpp:2196
std::vector< Slave > slaves
Definition: ncmesh.hpp:171
virtual void BuildFaceList()
Definition: ncmesh.cpp:1876
void ForceRefinement(Node *v1, Node *v2, Node *v3, Node *v4)
Definition: ncmesh.cpp:639
void SetVertexPositions(const Array< mfem::Vertex > &vertices)
I/O: Set positions of all vertices (used by mesh loader).
Definition: ncmesh.cpp:3168
int edge_flags
edge orientation flags
Definition: ncmesh.hpp:156
bool Boundary() const
Definition: ncmesh.hpp:375
double pos[3]
3D position
Definition: ncmesh.hpp:310
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:331
static const PointMatrix & GetGeomIdentity(int geom)
Definition: ncmesh.cpp:2423
Element * child[8]
2-8 children (if ref_type != 0)
Definition: ncmesh.hpp:402
Array< Element * > leaf_elements
Definition: ncmesh.hpp:430
Face * GetFace(Element *elem, int face_no)
Definition: ncmesh.cpp:391
mfem::Element * NewMeshElement(int geom) const
Definition: ncmesh.cpp:1583
void DeleteLast()
Delete the last entry.
Definition: array.hpp:151
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition: table.hpp:24
void UpdateLeafElements()
Definition: ncmesh.cpp:1560
bool NodeSetZ1(Node *node, Node **n)
Definition: ncmesh.cpp:632
Node * GetMidEdgeVertexSimple(Node *v1, Node *v2)
Definition: ncmesh.cpp:593
long MemoryUsage() const
Definition: ncmesh.cpp:3307
void BuildElementToVertexTable()
Definition: ncmesh.cpp:2119
HashTable< Node > nodes
Definition: ncmesh.hpp:416
int GetNV() const
Definition: mesh.hpp:493
void UnrefElementNodes(Element *elem)
Definition: ncmesh.cpp:335
virtual int GetNFaces(int &nFaceVertices) const =0
static PointMatrix pm_quad_identity
Definition: ncmesh.hpp:655
void CheckIsoFace(Node *v1, Node *v2, Node *v3, Node *v4, Node *e1, Node *e2, Node *e3, Node *e4, Node *midf)
Definition: ncmesh.cpp:739
void SetIJ(int *newI, int *newJ, int newsize=-1)
Definition: table.cpp:196
Node * GetMidFaceVertex(Node *e1, Node *e2, Node *e3, Node *e4)
Definition: ncmesh.cpp:602
void CollectFaceVertices(Node *v0, Node *v1, Node *v2, Node *v3, Array< int > &indices)
Definition: ncmesh.cpp:2095
bool NodeSetZ2(Node *node, Node **n)
Definition: ncmesh.cpp:635
void Initialize(const mfem::Element *elem)
Definition: ncmesh.cpp:29
bool NodeSetX2(Node *node, Node **n)
Definition: ncmesh.cpp:623
void GetMeshComponents(Array< mfem::Vertex > &vertices, Array< mfem::Element * > &elements, Array< mfem::Element * > &boundary) const
Return the basic Mesh arrays for the current finest level.
Definition: ncmesh.cpp:1595
Node * node[8]
element corners (if ref_type == 0)
Definition: ncmesh.hpp:401
void CollectLeafElements(Element *elem, int state)
Definition: ncmesh.cpp:1520
Element * NewTriangle(Node *n0, Node *n1, Node *n2, int attr, int eattr0, int eattr1, int eattr2)
Definition: ncmesh.cpp:546
virtual int GetNVertices() const =0
int parent
element index in the coarse mesh
Definition: ncmesh.hpp:44
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:53
virtual void AssignLeafIndices()
Definition: ncmesh.cpp:1574
static PointMatrix pm_hex_identity
Definition: ncmesh.hpp:656
static int find_hex_face(int a, int b, int c)
Definition: ncmesh.cpp:1779
virtual const int * GetFaceVertices(int fi) const =0
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
Definition: ncmesh.cpp:1371
Array< ElemRefType > ref_stack
stack of scheduled refinements (temporary)
Definition: ncmesh.hpp:465
T & Last()
Return the last element in the array.
Definition: array.hpp:411
static void find_face_nodes(const Face *face, Node *node[4])
Definition: ncmesh.cpp:2912
Array< Element * > root_elements
Definition: ncmesh.hpp:414
virtual void Update()
Definition: ncmesh.cpp:187
long MemoryUsage() const
Return total number of bytes allocated.
Definition: ncmesh.cpp:3336
Node * GetMidEdgeVertex(Node *v1, Node *v2)
Definition: ncmesh.cpp:584
int EdgeSplitLevel(Node *v1, Node *v2) const
Definition: ncmesh.cpp:2988
Element * GetSingleElement() const
Return one of elem[0] or elem[1] and make sure the other is NULL.
Definition: ncmesh.cpp:410
void PrintCoarseElements(std::ostream &out) const
I/O: Print the &quot;coarse_elements&quot; section of the mesh file (ver. &gt;= 1.1).
Definition: ncmesh.cpp:3211
static int find_element_edge(Element *elem, int v0, int v1)
Definition: ncmesh.cpp:1763
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:51
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition: ncmesh.hpp:188
void TraverseFace(Node *v0, Node *v1, Node *v2, Node *v3, const PointMatrix &pm, int level)
Definition: ncmesh.cpp:1827
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
void UpdateElementToVertexTable()
Definition: ncmesh.hpp:579
int RowSize(int i) const
Definition: table.hpp:102
int GetElementDepth(int i) const
Return the distance of leaf &#39;i&#39; from the root.
Definition: ncmesh.cpp:2900
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:71
void CountSplits(Element *elem, int splits[3]) const
Definition: ncmesh.cpp:3028
void ForgetElement(Element *e)
Definition: ncmesh.cpp:384
char ref_type
bit mask of X,Y,Z refinements (bits 0,1,2 respectively)
Definition: ncmesh.hpp:394
bool NodeSetY2(Node *node, Node **n)
Definition: ncmesh.cpp:629
void LoadVertexParents(std::istream &input)
Definition: ncmesh.cpp:3148
int CountElements(Element *elem) const
Definition: ncmesh.cpp:3294
Abstract data type element.
Definition: element.hpp:27
Data type line segment element.
Definition: segment.hpp:22
Array< Face * > boundary_faces
subset of all faces, set by BuildFaceList
Definition: ncmesh.hpp:437
void TraverseRefinements(Element *elem, int coarse_index, std::string &ref_path, RefPathMap &map)
Definition: ncmesh.cpp:2746
int rank
processor number (ParNCMesh)
Definition: ncmesh.hpp:397
void CollectDerefinements(Element *elem, Array< Connection > &list)
Definition: ncmesh.cpp:1320
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:543
virtual int GetType() const =0
Returns element&#39;s type.
DenseMatrix point_matrix
position within the master edge/face
Definition: ncmesh.hpp:157
HashTable< Face > faces
Definition: ncmesh.hpp:417
Defines the position of a fine element within a coarse element.
Definition: ncmesh.hpp:42
virtual void Refine(const Array< Refinement > &refinements)
Definition: ncmesh.cpp:1169
int matrix
index into CoarseFineTransformations::point_matrices
Definition: ncmesh.hpp:45
int GetEdgeMaster(int v1, int v2) const
Definition: ncmesh.cpp:2891
void UnrefVertex(HashTable< Node > &nodes)
Definition: ncmesh.cpp:277