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