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