MFEM  v3.4 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 <cmath>
17 #include <climits> // INT_MAX
18
19 namespace mfem
20 {
21
22 NCMesh::GeomInfo NCMesh::GI[Geometry::NumGeom];
23
24 NCMesh::GeomInfo& NCMesh::gi_hex = NCMesh::GI[Geometry::CUBE];
26 NCMesh::GeomInfo& NCMesh::gi_tri = NCMesh::GI[Geometry::TRIANGLE];
27
29 {
30  if (initialized) { return; }
31
32  nv = elem->GetNVertices();
33  ne = elem->GetNEdges();
34  nf = elem->GetNFaces(nfv);
35
36  for (int i = 0; i < ne; i++)
37  {
38  for (int j = 0; j < 2; j++)
39  {
40  edges[i][j] = elem->GetEdgeVertices(i)[j];
41  }
42  }
43  for (int i = 0; i < nf; i++)
44  {
45  for (int j = 0; j < nfv; j++)
46  {
47  faces[i][j] = elem->GetFaceVertices(i)[j];
48  }
49  }
50
51  // in 2D we pretend to have faces too, so we can use Face::elem[2]
52  if (!nf)
53  {
54  for (int i = 0; i < ne; i++)
55  {
56  // make a degenerate face
57  faces[i][0] = faces[i][1] = edges[i][0];
58  faces[i][2] = faces[i][3] = edges[i][1];
59  }
60  nf = ne;
61  }
62
63  initialized = true;
64 }
65
66
67 NCMesh::NCMesh(const Mesh *mesh, std::istream *vertex_parents)
68 {
69  Dim = mesh->Dimension();
70  spaceDim = mesh->SpaceDimension();
71
72  // assume the mesh is anisotropic if we're loading a file
73  Iso = vertex_parents ? false : true;
74
75  // examine elements and reserve the first node IDs for vertices
76  // (note: 'mesh' may not have vertices defined yet, e.g., on load)
77  int max_id = -1;
78  for (int i = 0; i < mesh->GetNE(); i++)
79  {
80  const mfem::Element *elem = mesh->GetElement(i);
81  const int *v = elem->GetVertices(), nv = elem->GetNVertices();
82  for (int j = 0; j < nv; j++)
83  {
84  max_id = std::max(max_id, v[j]);
85  }
86  }
87  for (int id = 0; id <= max_id; id++)
88  {
89  // top-level nodes are special: id == p1 == p2 == orig. vertex id
90  int node = nodes.GetId(id, id);
91  MFEM_CONTRACT_VAR(node);
92  MFEM_ASSERT(node == id, "");
93  }
94
95  // if a mesh file is being read, load the vertex hierarchy now;
96  // 'vertex_parents' must be at the appropriate section in the mesh file
97  if (vertex_parents)
98  {
100  }
101  else
102  {
103  top_vertex_pos.SetSize(3*mesh->GetNV());
104  for (int i = 0; i < mesh->GetNV(); i++)
105  {
106  memcpy(&top_vertex_pos[3*i], mesh->GetVertex(i), 3*sizeof(double));
107  }
108  }
109
110  // create the NCMesh::Element struct for each Mesh element
111  root_count = mesh->GetNE();
112  for (int i = 0; i < root_count; i++)
113  {
114  const mfem::Element *elem = mesh->GetElement(i);
115
116  int geom = elem->GetGeometryType();
117  if (geom != Geometry::TRIANGLE &&
118  geom != Geometry::SQUARE &&
119  geom != Geometry::CUBE)
120  {
121  MFEM_ABORT("only triangles, quads and hexes are supported by NCMesh.");
122  }
123
124  // initialize edge/face tables for this type of element
125  GI[geom].Initialize(elem);
126
127  // create our Element struct for this element
128  int root_id = AddElement(Element(geom, elem->GetAttribute()));
129  MFEM_ASSERT(root_id == i, "");
130  Element &root_elem = elements[root_id];
131
132  const int *v = elem->GetVertices();
133  for (int j = 0; j < GI[geom].nv; j++)
134  {
135  root_elem.node[j] = v[j];
136  }
137
138  // increase reference count of all nodes the element is using
139  // (NOTE: this will also create and reference all edge nodes and faces)
140  RefElement(root_id);
141
142  // make links from faces back to the element
143  RegisterFaces(root_id);
144  }
145
146  // store boundary element attributes
147  for (int i = 0; i < mesh->GetNBE(); i++)
148  {
149  const mfem::Element *be = mesh->GetBdrElement(i);
150  const int *v = be->GetVertices();
151
153  {
154  Face* face = faces.Find(v[0], v[1], v[2], v[3]);
156  face->attribute = be->GetAttribute();
157  }
158  else if (be->GetType() == mfem::Element::SEGMENT)
159  {
160  Face* face = faces.Find(v[0], v[0], v[1], v[1]);
162  face->attribute = be->GetAttribute();
163  }
164  else
165  {
166  MFEM_ABORT("only segment and quadrilateral boundary "
167  "elements are supported by NCMesh.");
168  }
169  }
170
171  Update();
172 }
173
174 NCMesh::NCMesh(const NCMesh &other)
175  : Dim(other.Dim)
176  , spaceDim(other.spaceDim)
177  , Iso(other.Iso)
178  , nodes(other.nodes)
179  , faces(other.faces)
180  , elements(other.elements)
181  , root_count(other.root_count)
182 {
185  Update();
186 }
187
189 {
191  UpdateVertices();
192
193  vertex_list.Clear();
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  NVertices = 0;
1376  for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
1377  {
1378  if (node->HasVertex()) { node->vert_index = NVertices++; }
1379  }
1380
1382
1383  NVertices = 0;
1384  for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
1385  {
1386  if (node->HasVertex()) { vertex_nodeId[NVertices++] = 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 #ifdef MFEM_DEBUG
1508  TmpVertex &tv = tmp_vertex[node]; // to make DebugDump work
1509 #else
1510  TmpVertex &tv = tmp_vertex[nd.vert_index];
1511 #endif
1512  if (tv.valid) { return tv.pos; }
1513
1514  MFEM_VERIFY(tv.visited == false, "cyclic vertex dependencies.");
1515  tv.visited = true;
1516
1517  const double* pos1 = CalcVertexPos(nd.p1);
1518  const double* pos2 = CalcVertexPos(nd.p2);
1519
1520  for (int i = 0; i < 3; i++)
1521  {
1522  tv.pos[i] = (pos1[i] + pos2[i]) * 0.5;
1523  }
1524  tv.valid = true;
1525  return tv.pos;
1526 }
1527
1529  Array<mfem::Element*>& melements,
1530  Array<mfem::Element*>& mboundary) const
1531 {
1532  mvertices.SetSize(vertex_nodeId.Size());
1533  if (top_vertex_pos.Size())
1534  {
1535  // calculate vertex positions from stored top-level vertex coordinates
1536  tmp_vertex = new TmpVertex[nodes.NumIds()];
1537  for (int i = 0; i < mvertices.Size(); i++)
1538  {
1539  mvertices[i].SetCoords(spaceDim, CalcVertexPos(vertex_nodeId[i]));
1540  }
1541  delete [] tmp_vertex;
1542  }
1543  // NOTE: if the mesh is curved (top_vertex_pos is empty), mvertices are left
1544  // uninitialized here; they will be initialized later by the Mesh from Nodes
1545  // - here we just make sure mvertices has the correct size.
1546
1547  melements.SetSize(leaf_elements.Size() - GetNumGhostElements());
1548  melements.SetSize(0);
1549
1550  mboundary.SetSize(0);
1551
1552  // create an mfem::Element for each leaf Element
1553  for (int i = 0; i < leaf_elements.Size(); i++)
1554  {
1555  const Element &nc_elem = elements[leaf_elements[i]];
1556  if (IsGhost(nc_elem)) { continue; } // ParNCMesh
1557
1558  const int* node = nc_elem.node;
1559  GeomInfo& gi = GI[(int) nc_elem.geom];
1560
1561  mfem::Element* elem = NewMeshElement(nc_elem.geom);
1562  melements.Append(elem);
1563
1564  elem->SetAttribute(nc_elem.attribute);
1565  for (int j = 0; j < gi.nv; j++)
1566  {
1567  elem->GetVertices()[j] = nodes[node[j]].vert_index;
1568  }
1569
1570  // create boundary elements
1571  for (int k = 0; k < gi.nf; k++)
1572  {
1573  const int* fv = gi.faces[k];
1574  const Face* face = faces.Find(node[fv[0]], node[fv[1]],
1575  node[fv[2]], node[fv[3]]);
1576  if (face->Boundary())
1577  {
1578  if (nc_elem.geom == Geometry::CUBE)
1579  {
1582  for (int j = 0; j < 4; j++)
1583  {
1585  }
1587  }
1588  else
1589  {
1590  Segment* segment = new Segment;
1591  segment->SetAttribute(face->attribute);
1592  for (int j = 0; j < 2; j++)
1593  {
1594  segment->GetVertices()[j] = nodes[node[fv[2*j]]].vert_index;
1595  }
1596  mboundary.Append(segment);
1597  }
1598  }
1599  }
1600  }
1601 }
1602
1604 {
1605  Table *edge_vertex = mesh->GetEdgeVertexTable();
1606
1607  // get edge enumeration from the Mesh
1608  for (int i = 0; i < edge_vertex->Size(); i++)
1609  {
1610  const int *ev = edge_vertex->GetRow(i);
1611  Node* node = nodes.Find(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
1612
1614  node->edge_index = i;
1615  }
1616
1617  // get face enumeration from the Mesh
1618  for (int i = 0; i < mesh->GetNumFaces(); i++)
1619  {
1620  const int* fv = mesh->GetFace(i)->GetVertices();
1621  Face* face;
1622  if (Dim == 3)
1623  {
1624  MFEM_ASSERT(mesh->GetFace(i)->GetNVertices() == 4, "");
1625  face = faces.Find(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
1626  vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
1627  }
1628  else
1629  {
1630  MFEM_ASSERT(mesh->GetFace(i)->GetNVertices() == 2, "");
1631  int n0 = vertex_nodeId[fv[0]], n1 = vertex_nodeId[fv[1]];
1632  face = faces.Find(n0, n0, n1, n1); // look up degenerate face
1633
1634 #ifdef MFEM_DEBUG
1635  // (non-ghost) edge and face numbers must match in 2D
1636  const int *ev = edge_vertex->GetRow(i);
1637  MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
1638  (ev[1] == fv[0] && ev[0] == fv[1]), "");
1639 #endif
1640  }
1642  face->index = i;
1643  }
1644
1645  NEdges = mesh->GetNEdges();
1646  NFaces = mesh->GetNumFaces();
1647 }
1648
1649
1650 //// Face/edge lists ///////////////////////////////////////////////////////////
1651
1652 int NCMesh::FaceSplitType(int v1, int v2, int v3, int v4,
1653  int mid[4]) const
1654 {
1655  MFEM_ASSERT(Dim >= 3, "");
1656
1657  // find edge nodes
1658  int e1 = nodes.FindId(v1, v2);
1659  int e2 = nodes.FindId(v2, v3);
1660  int e3 = (e1 >= 0 && nodes[e1].HasVertex()) ? nodes.FindId(v3, v4) : -1;
1661  int e4 = (e2 >= 0 && nodes[e2].HasVertex()) ? nodes.FindId(v4, v1) : -1;
1662
1663  // optional: return the mid-edge nodes if requested
1664  if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
1665
1666  // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
1667  int midf1 = -1, midf2 = -1;
1668  if (e1 >= 0 && e3 >= 0) { midf1 = nodes.FindId(e1, e3); }
1669  if (e2 >= 0 && e4 >= 0) { midf2 = nodes.FindId(e2, e4); }
1670
1671  // only one way to access the mid-face node must always exist
1672  MFEM_ASSERT(!(midf1 >= 0 && midf2 >= 0), "incorrectly split face!");
1673
1674  if (midf1 < 0 && midf2 < 0) { return 0; } // face not split
1675  else if (midf1 >= 0) { return 1; } // face split "vertically"
1676  else { return 2; } // face split "horizontally"
1677 }
1678
1679 int NCMesh::find_node(const Element &el, int node)
1680 {
1681  for (int i = 0; i < 8; i++)
1682  {
1683  if (el.node[i] == node) { return i; }
1684  }
1686  return -1;
1687 }
1688
1689 int NCMesh::find_element_edge(const Element &el, int vn0, int vn1)
1690 {
1691  MFEM_ASSERT(!el.ref_type, "");
1692  GeomInfo &gi = GI[(int) el.geom];
1693  for (int i = 0; i < gi.ne; i++)
1694  {
1695  const int* ev = gi.edges[i];
1696  int n0 = el.node[ev[0]];
1697  int n1 = el.node[ev[1]];
1698  if ((n0 == vn0 && n1 == vn1) ||
1699  (n0 == vn1 && n1 == vn0)) { return i; }
1700  }
1702  return -1;
1703 }
1704
1705 int NCMesh::find_hex_face(int a, int b, int c)
1706 {
1707  for (int i = 0; i < 6; i++)
1708  {
1709  const int* fv = gi_hex.faces[i];
1710  if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
1711  (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
1712  (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
1713  {
1714  return i;
1715  }
1716  }
1718  return -1;
1719 }
1720
1721 int NCMesh::ReorderFacePointMat(int v0, int v1, int v2, int v3,
1722  int elem, DenseMatrix& mat) const
1723 {
1724  const Element &el = elements[elem];
1725  int master[4] =
1726  {
1727  find_node(el, v0), find_node(el, v1),
1728  find_node(el, v2), find_node(el, v3)
1729  };
1730
1731  int local = find_hex_face(master[0], master[1], master[2]);
1732  const int* fv = gi_hex.faces[local];
1733
1734  DenseMatrix tmp(mat);
1735  for (int i = 0, j; i < 4; i++)
1736  {
1737  for (j = 0; j < 4; j++)
1738  {
1739  if (fv[i] == master[j])
1740  {
1741  // "pm.column(i) = tmp.column(j)"
1742  for (int k = 0; k < mat.Height(); k++)
1743  {
1744  mat(k,i) = tmp(k,j);
1745  }
1746  break;
1747  }
1748  }
1750  }
1751  return local;
1752 }
1753
1754 void NCMesh::TraverseFace(int vn0, int vn1, int vn2, int vn3,
1755  const PointMatrix& pm, int level)
1756 {
1757  if (level > 0)
1758  {
1759  // check if we made it to a face that is not split further
1760  Face* fa = faces.Find(vn0, vn1, vn2, vn3);
1761  if (fa)
1762  {
1763  // we have a slave face, add it to the list
1764  int elem = fa->GetSingleElement();
1765  face_list.slaves.push_back(Slave(fa->index, elem, -1));
1766  DenseMatrix &mat = face_list.slaves.back().point_matrix;
1767  pm.GetMatrix(mat);
1768
1769  // reorder the point matrix according to slave face orientation
1770  int local = ReorderFacePointMat(vn0, vn1, vn2, vn3, elem, mat);
1771  face_list.slaves.back().local = local;
1772
1773  return;
1774  }
1775  }
1776
1777  // we need to recurse deeper
1778  int mid[4];
1779  int split = FaceSplitType(vn0, vn1, vn2, vn3, mid);
1780
1781  if (split == 1) // "X" split face
1782  {
1783  Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
1784
1785  TraverseFace(vn0, mid[0], mid[2], vn3,
1786  PointMatrix(pm(0), mid0, mid2, pm(3)), level+1);
1787
1788  TraverseFace(mid[0], vn1, vn2, mid[2],
1789  PointMatrix(mid0, pm(1), pm(2), mid2), level+1);
1790  }
1791  else if (split == 2) // "Y" split face
1792  {
1793  Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
1794
1795  TraverseFace(vn0, vn1, mid[1], mid[3],
1796  PointMatrix(pm(0), pm(1), mid1, mid3), level+1);
1797
1798  TraverseFace(mid[3], mid[1], vn2, vn3,
1799  PointMatrix(mid3, mid1, pm(2), pm(3)), level+1);
1800  }
1801 }
1802
1804 {
1805  face_list.Clear();
1807
1808  if (Dim < 3) { return; }
1809
1810  Array<char> processed_faces(faces.NumIds());
1811  processed_faces = 0;
1812
1813  // visit faces of leaf elements
1814  for (int i = 0; i < leaf_elements.Size(); i++)
1815  {
1816  int elem = leaf_elements[i];
1817  Element &el = elements[elem];
1818  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
1819
1820  GeomInfo& gi = GI[(int) el.geom];
1821  for (int j = 0; j < gi.nf; j++)
1822  {
1823  // get nodes for this face
1824  int node[4];
1825  for (int k = 0; k < 4; k++)
1826  {
1827  node[k] = el.node[gi.faces[j][k]];
1828  }
1829
1830  int face = faces.FindId(node[0], node[1], node[2], node[3]);
1832
1833  // tell ParNCMesh about the face
1834  ElementSharesFace(elem, face);
1835
1836  // have we already processed this face? skip if yes
1837  if (processed_faces[face]) { continue; }
1838  processed_faces[face] = 1;
1839
1840  Face &fa = faces[face];
1841  if (fa.elem[0] >= 0 && fa.elem[1] >= 0)
1842  {
1843  // this is a conforming face, add it to the list
1844  face_list.conforming.push_back(MeshId(fa.index, elem, j));
1845  }
1846  else
1847  {
1848  PointMatrix pm(Point(0,0), Point(1,0), Point(1,1), Point(0,1));
1849
1850  // this is either a master face or a slave face, but we can't
1851  // tell until we traverse the face refinement 'tree'...
1852  int sb = face_list.slaves.size();
1853  TraverseFace(node[0], node[1], node[2], node[3], pm, 0);
1854
1855  int se = face_list.slaves.size();
1856  if (sb < se)
1857  {
1858  // found slaves, so this is a master face; add it to the list
1859  face_list.masters.push_back(Master(fa.index, elem, j, sb, se));
1860
1861  // also, set the master index for the slaves
1862  for (int i = sb; i < se; i++)
1863  {
1864  face_list.slaves[i].master = fa.index;
1865  }
1866  }
1867  }
1868
1869  if (fa.Boundary()) { boundary_faces.Append(face); }
1870  }
1871  }
1872 }
1873
1874 void NCMesh::TraverseEdge(int vn0, int vn1, double t0, double t1, int flags,
1875  int level)
1876 {
1877  int mid = nodes.FindId(vn0, vn1);
1878  if (mid < 0) { return; }
1879
1880  Node &nd = nodes[mid];
1881  if (nd.HasEdge() && level > 0)
1882  {
1883  // we have a slave edge, add it to the list
1884  edge_list.slaves.push_back(Slave(nd.edge_index, -1, -1));
1885  Slave &sl = edge_list.slaves.back();
1886
1887  sl.point_matrix.SetSize(1, 2);
1888  sl.point_matrix(0,0) = t0;
1889  sl.point_matrix(0,1) = t1;
1890
1891  // handle slave edge orientation
1892  sl.edge_flags = flags;
1893  int v0index = nodes[vn0].vert_index;
1894  int v1index = nodes[vn1].vert_index;
1895  if (v0index > v1index) { sl.edge_flags |= 2; }
1896
1897  // in 2D, get the element/local info from the degenerate face
1898  if (Dim == 2)
1899  {
1900  Face* fa = faces.Find(vn0, vn0, vn1, vn1);
1901  MFEM_ASSERT(fa != NULL, "");
1902  sl.element = fa->GetSingleElement();
1903  sl.local = find_element_edge(elements[sl.element], vn0, vn1);
1904  }
1905  }
1906
1907  // recurse deeper
1908  double tmid = (t0 + t1) / 2;
1909  TraverseEdge(vn0, mid, t0, tmid, flags, level+1);
1910  TraverseEdge(mid, vn1, tmid, t1, flags, level+1);
1911 }
1912
1914 {
1915  edge_list.Clear();
1916  if (Dim <= 2)
1917  {
1919  }
1920
1921  Array<char> processed_edges(nodes.NumIds());
1922  processed_edges = 0;
1923
1924  Array<int> edge_element(nodes.NumIds());
1925  Array<char> edge_local(nodes.NumIds());
1926  edge_local = -1;
1927
1928  // visit edges of leaf elements
1929  for (int i = 0; i < leaf_elements.Size(); i++)
1930  {
1931  int elem = leaf_elements[i];
1932  Element &el = elements[elem];
1933  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
1934
1935  GeomInfo& gi = GI[(int) el.geom];
1936  for (int j = 0; j < gi.ne; j++)
1937  {
1938  // get nodes for this edge
1939  const int* ev = gi.edges[j];
1940  int node[2] = { el.node[ev[0]], el.node[ev[1]] };
1941
1942  int enode = nodes.FindId(node[0], node[1]);
1944
1945  Node &nd = nodes[enode];
1947
1948  // tell ParNCMesh about the edge
1949  ElementSharesEdge(elem, enode);
1950
1951  // (2D only, store boundary faces)
1952  if (Dim <= 2)
1953  {
1954  int face = faces.FindId(node[0], node[0], node[1], node[1]);
1956  if (faces[face].Boundary()) { boundary_faces.Append(face); }
1957  }
1958
1959  // store element/local for later
1960  edge_element[nd.edge_index] = elem;
1961  edge_local[nd.edge_index] = j;
1962
1963  // skip slave edges here, they will be reached from their masters
1964  if (GetEdgeMaster(enode) >= 0) { continue; }
1965
1966  // have we already processed this edge? skip if yes
1967  if (processed_edges[enode]) { continue; }
1968  processed_edges[enode] = 1;
1969
1970  // prepare edge interval for slave traversal, handle orientation
1971  double t0 = 0.0, t1 = 1.0;
1972  int v0index = nodes[node[0]].vert_index;
1973  int v1index = nodes[node[1]].vert_index;
1974  int flags = (v0index > v1index) ? 1 : 0;
1975
1976  // try traversing the edge to find slave edges
1977  int sb = edge_list.slaves.size();
1978  TraverseEdge(node[0], node[1], t0, t1, flags, 0);
1979
1980  int se = edge_list.slaves.size();
1981  if (sb < se)
1982  {
1983  // found slaves, this is a master face; add it to the list
1984  edge_list.masters.push_back(Master(nd.edge_index, elem, j, sb, se));
1985
1986  // also, set the master index for the slaves
1987  for (int i = sb; i < se; i++)
1988  {
1989  edge_list.slaves[i].master = nd.edge_index;
1990  }
1991  }
1992  else
1993  {
1994  // no slaves, this is a conforming edge
1995  edge_list.conforming.push_back(MeshId(nd.edge_index, elem, j));
1996  }
1997  }
1998  }
1999
2000  // fix up slave edge element/local
2001  for (unsigned i = 0; i < edge_list.slaves.size(); i++)
2002  {
2003  Slave &sl = edge_list.slaves[i];
2004  int local = edge_local[sl.index];
2005  if (local >= 0)
2006  {
2007  sl.local = local;
2008  sl.element = edge_element[sl.index];
2009  }
2010  }
2011 }
2012
2014 {
2015  int total = NVertices + GetNumGhostVertices();
2016
2017  vertex_list.Clear();
2018  vertex_list.conforming.reserve(total);
2019
2020  Array<char> processed_vertices(total);
2021  processed_vertices = 0;
2022
2023  // analogously to above, visit vertices of leaf elements
2024  for (int i = 0; i < leaf_elements.Size(); i++)
2025  {
2026  int elem = leaf_elements[i];
2027  Element &el = elements[elem];
2028
2029  for (int j = 0; j < GI[(int) el.geom].nv; j++)
2030  {
2031  int node = el.node[j];
2032  Node &nd = nodes[node];
2033
2034  int index = nd.vert_index;
2035  if (index >= 0)
2036  {
2037  ElementSharesVertex(elem, node);
2038
2039  if (processed_vertices[index]) { continue; }
2040  processed_vertices[index] = 1;
2041
2042  vertex_list.conforming.push_back(MeshId(index, elem, j));
2043  }
2044  }
2045  }
2046 }
2047
2049 {
2050  oriented_matrix = point_matrix;
2051
2052  if (edge_flags)
2053  {
2054  MFEM_ASSERT(oriented_matrix.Height() == 1 &&
2055  oriented_matrix.Width() == 2, "not an edge point matrix");
2056
2057  if (edge_flags & 1) // master inverted
2058  {
2059  oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
2060  oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
2061  }
2062  if (edge_flags & 2) // slave inverted
2063  {
2064  std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
2065  }
2066  }
2067 }
2068
2069 void NCMesh::NCList::Clear(bool hard)
2070 {
2071  if (!hard)
2072  {
2073  conforming.clear();
2074  masters.clear();
2075  slaves.clear();
2076  }
2077  else
2078  {
2079  NCList empty;
2080  conforming.swap(empty.conforming);
2081  masters.swap(empty.masters);
2082  slaves.swap(empty.slaves);
2083  }
2084  inv_index.DeleteAll();
2085 }
2086
2088 {
2089  return conforming.size() + masters.size() + slaves.size();
2090 }
2091
2092 const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
2093 {
2094  if (!inv_index.Size())
2095  {
2096  int max_index = -1;
2097  for (unsigned i = 0; i < conforming.size(); i++)
2098  {
2099  max_index = std::max(conforming[i].index, max_index);
2100  }
2101  for (unsigned i = 0; i < masters.size(); i++)
2102  {
2103  max_index = std::max(masters[i].index, max_index);
2104  }
2105  for (unsigned i = 0; i < slaves.size(); i++)
2106  {
2107  max_index = std::max(slaves[i].index, max_index);
2108  }
2109
2110  inv_index.SetSize(max_index + 1);
2111  inv_index = -1;
2112
2113  for (unsigned i = 0; i < conforming.size(); i++)
2114  {
2115  inv_index[conforming[i].index] = (i << 2);
2116  }
2117  for (unsigned i = 0; i < masters.size(); i++)
2118  {
2119  inv_index[masters[i].index] = (i << 2) + 1;
2120  }
2121  for (unsigned i = 0; i < slaves.size(); i++)
2122  {
2123  inv_index[slaves[i].index] = (i << 2) + 2;
2124  }
2125  }
2126
2127  MFEM_ASSERT(index >= 0 && index < inv_index.Size(), "");
2128  int key = inv_index[index];
2130
2131  if (type) { *type = key & 0x3; }
2132
2133  switch (key & 0x3)
2134  {
2135  case 0: return conforming[key >> 2];
2136  case 1: return masters[key >> 2];
2137  case 2: return slaves[key >> 2];
2138  default: MFEM_ABORT("internal error"); return conforming[0];
2139  }
2140 }
2141
2142
2143 //// Neighbors /////////////////////////////////////////////////////////////////
2144
2145 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
2146 {
2147  int mid = nodes.FindId(v0, v1);
2148  if (mid >= 0 && nodes[mid].HasVertex())
2149  {
2150  indices.Append(mid);
2151
2152  CollectEdgeVertices(v0, mid, indices);
2153  CollectEdgeVertices(mid, v1, indices);
2154  }
2155 }
2156
2157 void NCMesh::CollectFaceVertices(int v0, int v1, int v2, int v3,
2158  Array<int> &indices)
2159 {
2160  int mid[4];
2161  switch (FaceSplitType(v0, v1, v2, v3, mid))
2162  {
2163  case 1:
2164  indices.Append(mid[0]);
2165  indices.Append(mid[2]);
2166
2167  CollectFaceVertices(v0, mid[0], mid[2], v3, indices);
2168  CollectFaceVertices(mid[0], v1, v2, mid[2], indices);
2169  break;
2170
2171  case 2:
2172  indices.Append(mid[1]);
2173  indices.Append(mid[3]);
2174
2175  CollectFaceVertices(v0, v1, mid[1], mid[3], indices);
2176  CollectFaceVertices(mid[3], mid[1], v2, v3, indices);
2177  break;
2178  }
2179 }
2180
2182 {
2183  int nrows = leaf_elements.Size();
2184  int* I = new int[nrows + 1];
2185  int** JJ = new int*[nrows];
2186
2187  Array<int> indices;
2188  indices.Reserve(128);
2189
2190  // collect vertices coinciding with each element, including hanging vertices
2191  for (int i = 0; i < leaf_elements.Size(); i++)
2192  {
2193  int elem = leaf_elements[i];
2194  Element &el = elements[elem];
2195  MFEM_ASSERT(!el.ref_type, "not a leaf element.");
2196
2197  GeomInfo& gi = GI[(int) el.geom];
2198  int* node = el.node;
2199
2200  indices.SetSize(0);
2201  for (int j = 0; j < gi.ne; j++)
2202  {
2203  const int* ev = gi.edges[j];
2204  CollectEdgeVertices(node[ev[0]], node[ev[1]], indices);
2205  }
2206
2207  if (Dim >= 3)
2208  {
2209  for (int j = 0; j < gi.nf; j++)
2210  {
2211  const int* fv = gi.faces[j];
2212  CollectFaceVertices(node[fv[0]], node[fv[1]],
2213  node[fv[2]], node[fv[3]], indices);
2214  }
2215  }
2216
2217  // temporarily store one row of the table
2218  indices.Sort();
2219  indices.Unique();
2220  int size = indices.Size();
2221  I[i] = size;
2222  JJ[i] = new int[size];
2223  memcpy(JJ[i], indices.GetData(), size * sizeof(int));
2224  }
2225
2226  // finalize the I array of the table
2227  int nnz = 0;
2228  for (int i = 0; i < nrows; i++)
2229  {
2230  int cnt = I[i];
2231  I[i] = nnz;
2232  nnz += cnt;
2233  }
2234  I[nrows] = nnz;
2235
2236  // copy the temporarily stored rows into one J array
2237  int *J = new int[nnz];
2238  nnz = 0;
2239  for (int i = 0; i < nrows; i++)
2240  {
2241  int cnt = I[i+1] - I[i];
2242  memcpy(J+nnz, JJ[i], cnt * sizeof(int));
2243  delete [] JJ[i];
2244  nnz += cnt;
2245  }
2246
2247  element_vertex.SetIJ(I, J, nrows);
2248
2249  delete [] JJ;
2250 }
2251
2252
2254  Array<int> *neighbors,
2255  Array<char> *neighbor_set)
2256 {
2257  // If A is the element-to-vertex table (see 'element_vertex') listing all
2258  // vertices touching each element, including hanging vertices, then A*A^T is
2259  // the element-to-neighbor table. Multiplying the element set with A*A^T
2260  // gives the neighbor set. To save memory, this function only computes the
2261  // action of A*A^T, the product itself is not stored anywhere.
2262
2263  // Optimization: the 'element_vertex' table does not store the obvious
2264  // corner nodes in it. The table is therefore empty for conforming meshes.
2265
2267
2268  int nleaves = leaf_elements.Size();
2269  MFEM_VERIFY(elem_set.Size() == nleaves, "");
2270  MFEM_ASSERT(element_vertex.Size() == nleaves, "");
2271
2272  // step 1: vertices = A^T * elem_set, i.e, find all vertices touching the
2273  // element set
2274
2275  Array<char> vmark(nodes.NumIds());
2276  vmark = 0;
2277
2278  for (int i = 0; i < nleaves; i++)
2279  {
2280  if (elem_set[i])
2281  {
2282  int *v = element_vertex.GetRow(i);
2283  int nv = element_vertex.RowSize(i);
2284  for (int j = 0; j < nv; j++)
2285  {
2286  vmark[v[j]] = 1;
2287  }
2288
2289  Element &el = elements[leaf_elements[i]];
2290  nv = GI[(int) el.geom].nv;
2291  for (int j = 0; j < nv; j++)
2292  {
2293  vmark[el.node[j]] = 1;
2294  }
2295  }
2296  }
2297
2298  // step 2: neighbors = A * vertices, i.e., find all elements coinciding with
2299  // vertices from step 1; NOTE: in the result we don't include elements from
2300  // the original set
2301
2302  if (neighbor_set)
2303  {
2304  neighbor_set->SetSize(nleaves);
2305  *neighbor_set = 0;
2306  }
2307
2308  for (int i = 0; i < nleaves; i++)
2309  {
2310  if (!elem_set[i])
2311  {
2312  bool hit = false;
2313
2314  int *v = element_vertex.GetRow(i);
2315  int nv = element_vertex.RowSize(i);
2316  for (int j = 0; j < nv; j++)
2317  {
2318  if (vmark[v[j]]) { hit = true; break; }
2319  }
2320
2321  if (!hit)
2322  {
2323  Element &el = elements[leaf_elements[i]];
2324  nv = GI[(int) el.geom].nv;
2325  for (int j = 0; j < nv; j++)
2326  {
2327  if (vmark[el.node[j]]) { hit = true; break; }
2328  }
2329  }
2330
2331  if (hit)
2332  {
2333  if (neighbors) { neighbors->Append(leaf_elements[i]); }
2334  if (neighbor_set) { (*neighbor_set)[i] = 1; }
2335  }
2336  }
2337  }
2338 }
2339
2340 static bool sorted_lists_intersect(const int* a, const int* b, int na, int nb)
2341 {
2342  if (!na || !nb) { return false; }
2343  int a_last = a[na-1], b_last = b[nb-1];
2344  if (*b < *a) { goto l2; } // woo-hoo! I always wanted to use a goto! :)
2345 l1:
2346  if (a_last < *b) { return false; }
2347  while (*a < *b) { a++; }
2348  if (*a == *b) { return true; }
2349 l2:
2350  if (b_last < *a) { return false; }
2351  while (*b < *a) { b++; }
2352  if (*a == *b) { return true; }
2353  goto l1;
2354 }
2355
2356 void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
2357  const Array<int> *search_set)
2358 {
2359  // TODO future: this function is inefficient. For a single element, an
2360  // octree neighbor search algorithm would be better. However, the octree
2361  // neighbor algorithm is hard to get right in the multi-octree case due to
2362  // the different orientations of the octrees (i.e., the root elements).
2363
2365
2366  // sorted list of all vertex nodes touching 'elem'
2367  Array<int> vert;
2368  vert.Reserve(128);
2369
2370  // support for non-leaf 'elem', collect vertices of all children
2371  Array<int> stack;
2372  stack.Reserve(64);
2373  stack.Append(elem);
2374
2375  while (stack.Size())
2376  {
2377  Element &el = elements[stack.Last()];
2378  stack.DeleteLast();
2379
2380  if (!el.ref_type)
2381  {
2382  int *v = element_vertex.GetRow(el.index);
2383  int nv = element_vertex.RowSize(el.index);
2384  for (int i = 0; i < nv; i++)
2385  {
2386  vert.Append(v[i]);
2387  }
2388
2389  nv = GI[(int) el.geom].nv;
2390  for (int i = 0; i < nv; i++)
2391  {
2392  vert.Append(el.node[i]);
2393  }
2394  }
2395  else
2396  {
2397  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
2398  {
2399  stack.Append(el.child[i]);
2400  }
2401  }
2402  }
2403
2404  vert.Sort();
2405  vert.Unique();
2406
2407  int *v1 = vert.GetData();
2408  int nv1 = vert.Size();
2409
2410  if (!search_set) { search_set = &leaf_elements; }
2411
2412  // test *all* potential neighbors from the search set
2413  for (int i = 0; i < search_set->Size(); i++)
2414  {
2415  int testme = (*search_set)[i];
2416  if (testme != elem)
2417  {
2418  Element &el = elements[testme];
2419  int *v2 = element_vertex.GetRow(el.index);
2420  int nv2 = element_vertex.RowSize(el.index);
2421
2422  bool hit = sorted_lists_intersect(v1, v2, nv1, nv2);
2423
2424  if (!hit)
2425  {
2426  int nv = GI[(int) el.geom].nv;
2427  for (int j = 0; j < nv; j++)
2428  {
2429  hit = sorted_lists_intersect(&el.node[j], v1, 1, nv1);
2430  if (hit) { break; }
2431  }
2432  }
2433
2434  if (hit) { neighbors.Append(testme); }
2435  }
2436  }
2437 }
2438
2440  Array<int> &expanded,
2441  const Array<int> *search_set)
2442 {
2444
2445  Array<char> vmark(nodes.NumIds());
2446  vmark = 0;
2447
2448  for (int i = 0; i < elems.Size(); i++)
2449  {
2450  Element &el = elements[elems[i]];
2451
2452  int *v = element_vertex.GetRow(el.index);
2453  int nv = element_vertex.RowSize(el.index);
2454  for (int j = 0; j < nv; j++)
2455  {
2456  vmark[v[j]] = 1;
2457  }
2458
2459  nv = GI[(int) el.geom].nv;
2460  for (int j = 0; j < nv; j++)
2461  {
2462  vmark[el.node[j]] = 1;
2463  }
2464  }
2465
2466  if (!search_set)
2467  {
2468  search_set = &leaf_elements;
2469  }
2470
2471  expanded.SetSize(0);
2472  for (int i = 0; i < search_set->Size(); i++)
2473  {
2474  int testme = (*search_set)[i];
2475  Element &el = elements[testme];
2476  bool hit = false;
2477
2478  int *v = element_vertex.GetRow(el.index);
2479  int nv = element_vertex.RowSize(el.index);
2480  for (int j = 0; j < nv; j++)
2481  {
2482  if (vmark[v[j]]) { hit = true; break; }
2483  }
2484
2485  if (!hit)
2486  {
2487  nv = GI[(int) el.geom].nv;
2488  for (int j = 0; j < nv; j++)
2489  {
2490  if (vmark[el.node[j]]) { hit = true; break; }
2491  }
2492  }
2493
2494  if (hit) { expanded.Append(testme); }
2495  }
2496 }
2497
2498
2499 //// Coarse/fine transformations ///////////////////////////////////////////////
2500
2502 {
2503  point_matrix.SetSize(points[0].dim, np);
2504  for (int i = 0; i < np; i++)
2505  {
2506  for (int j = 0; j < points[0].dim; j++)
2507  {
2508  point_matrix(j, i) = points[i].coord[j];
2509  }
2510  }
2511 }
2512
2514  Point(0, 0), Point(1, 0), Point(0, 1)
2515 );
2517  Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)
2518 );
2520  Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
2521  Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
2522 );
2523
2525 {
2526  switch (geom)
2527  {
2528  case Geometry::TRIANGLE: return pm_tri_identity;
2530  case Geometry::CUBE: return pm_hex_identity;
2531  default:
2532  MFEM_ABORT("unsupported geometry.");
2533  return pm_tri_identity;
2534  }
2535 }
2536
2537 void NCMesh::GetPointMatrix(int geom, const char* ref_path, DenseMatrix& matrix)
2538 {
2539  PointMatrix pm = GetGeomIdentity(geom);
2540
2541  while (*ref_path)
2542  {
2543  int ref_type = *ref_path++;
2544  int child = *ref_path++;
2545
2546  if (geom == Geometry::CUBE)
2547  {
2548  if (ref_type == 1) // split along X axis
2549  {
2550  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2551  Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
2552
2553  if (child == 0)
2554  {
2555  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
2556  pm(4), mid45, mid67, pm(7));
2557  }
2558  else if (child == 1)
2559  {
2560  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
2561  mid45, pm(5), pm(6), mid67);
2562  }
2563  }
2564  else if (ref_type == 2) // split along Y axis
2565  {
2566  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2567  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2568
2569  if (child == 0)
2570  {
2571  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
2572  pm(4), pm(5), mid56, mid74);
2573  }
2574  else if (child == 1)
2575  {
2576  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
2577  mid74, mid56, pm(6), pm(7));
2578  }
2579  }
2580  else if (ref_type == 4) // split along Z axis
2581  {
2582  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2583  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2584
2585  if (child == 0)
2586  {
2587  pm = PointMatrix(pm(0), pm(1), pm(2), pm(3),
2588  mid04, mid15, mid26, mid37);
2589  }
2590  else if (child == 1)
2591  {
2592  pm = PointMatrix(mid04, mid15, mid26, mid37,
2593  pm(4), pm(5), pm(6), pm(7));
2594  }
2595  }
2596  else if (ref_type == 3) // XY split
2597  {
2598  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2599  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2600  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2601  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2602
2603  Point midf0(mid23, mid12, mid01, mid30);
2604  Point midf5(mid45, mid56, mid67, mid74);
2605
2606  if (child == 0)
2607  {
2608  pm = PointMatrix(pm(0), mid01, midf0, mid30,
2609  pm(4), mid45, midf5, mid74);
2610  }
2611  else if (child == 1)
2612  {
2613  pm = PointMatrix(mid01, pm(1), mid12, midf0,
2614  mid45, pm(5), mid56, midf5);
2615  }
2616  else if (child == 2)
2617  {
2618  pm = PointMatrix(midf0, mid12, pm(2), mid23,
2619  midf5, mid56, pm(6), mid67);
2620  }
2621  else if (child == 3)
2622  {
2623  pm = PointMatrix(mid30, midf0, mid23, pm(3),
2624  mid74, midf5, mid67, pm(7));
2625  }
2626  }
2627  else if (ref_type == 5) // XZ split
2628  {
2629  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2630  Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
2631  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2632  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2633
2634  Point midf1(mid01, mid15, mid45, mid04);
2635  Point midf3(mid23, mid37, mid67, mid26);
2636
2637  if (child == 0)
2638  {
2639  pm = PointMatrix(pm(0), mid01, mid23, pm(3),
2640  mid04, midf1, midf3, mid37);
2641  }
2642  else if (child == 1)
2643  {
2644  pm = PointMatrix(mid01, pm(1), pm(2), mid23,
2645  midf1, mid15, mid26, midf3);
2646  }
2647  else if (child == 2)
2648  {
2649  pm = PointMatrix(midf1, mid15, mid26, midf3,
2650  mid45, pm(5), pm(6), mid67);
2651  }
2652  else if (child == 3)
2653  {
2654  pm = PointMatrix(mid04, midf1, midf3, mid37,
2655  pm(4), mid45, mid67, pm(7));
2656  }
2657  }
2658  else if (ref_type == 6) // YZ split
2659  {
2660  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2661  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2662  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2663  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2664
2665  Point midf2(mid12, mid26, mid56, mid15);
2666  Point midf4(mid30, mid04, mid74, mid37);
2667
2668  if (child == 0)
2669  {
2670  pm = PointMatrix(pm(0), pm(1), mid12, mid30,
2671  mid04, mid15, midf2, midf4);
2672  }
2673  else if (child == 1)
2674  {
2675  pm = PointMatrix(mid30, mid12, pm(2), pm(3),
2676  midf4, midf2, mid26, mid37);
2677  }
2678  else if (child == 2)
2679  {
2680  pm = PointMatrix(mid04, mid15, midf2, midf4,
2681  pm(4), pm(5), mid56, mid74);
2682  }
2683  else if (child == 3)
2684  {
2685  pm = PointMatrix(midf4, midf2, mid26, mid37,
2686  mid74, mid56, pm(6), pm(7));
2687  }
2688  }
2689  else if (ref_type == 7) // full isotropic refinement
2690  {
2691  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2692  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2693  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2694  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2695  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2696  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2697
2698  Point midf0(mid23, mid12, mid01, mid30);
2699  Point midf1(mid01, mid15, mid45, mid04);
2700  Point midf2(mid12, mid26, mid56, mid15);
2701  Point midf3(mid23, mid37, mid67, mid26);
2702  Point midf4(mid30, mid04, mid74, mid37);
2703  Point midf5(mid45, mid56, mid67, mid74);
2704
2705  Point midel(midf1, midf3);
2706
2707  if (child == 0)
2708  {
2709  pm = PointMatrix(pm(0), mid01, midf0, mid30,
2710  mid04, midf1, midel, midf4);
2711  }
2712  else if (child == 1)
2713  {
2714  pm = PointMatrix(mid01, pm(1), mid12, midf0,
2715  midf1, mid15, midf2, midel);
2716  }
2717  else if (child == 2)
2718  {
2719  pm = PointMatrix(midf0, mid12, pm(2), mid23,
2720  midel, midf2, mid26, midf3);
2721  }
2722  else if (child == 3)
2723  {
2724  pm = PointMatrix(mid30, midf0, mid23, pm(3),
2725  midf4, midel, midf3, mid37);
2726  }
2727  else if (child == 4)
2728  {
2729  pm = PointMatrix(mid04, midf1, midel, midf4,
2730  pm(4), mid45, midf5, mid74);
2731  }
2732  else if (child == 5)
2733  {
2734  pm = PointMatrix(midf1, mid15, midf2, midel,
2735  mid45, pm(5), mid56, midf5);
2736  }
2737  else if (child == 6)
2738  {
2739  pm = PointMatrix(midel, midf2, mid26, midf3,
2740  midf5, mid56, pm(6), mid67);
2741  }
2742  else if (child == 7)
2743  {
2744  pm = PointMatrix(midf4, midel, midf3, mid37,
2745  mid74, midf5, mid67, pm(7));
2746  }
2747  }
2748  }
2749  else if (geom == Geometry::SQUARE)
2750  {
2751  if (ref_type == 1) // X split
2752  {
2753  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2754
2755  if (child == 0)
2756  {
2757  pm = PointMatrix(pm(0), mid01, mid23, pm(3));
2758  }
2759  else if (child == 1)
2760  {
2761  pm = PointMatrix(mid01, pm(1), pm(2), mid23);
2762  }
2763  }
2764  else if (ref_type == 2) // Y split
2765  {
2766  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2767
2768  if (child == 0)
2769  {
2770  pm = PointMatrix(pm(0), pm(1), mid12, mid30);
2771  }
2772  else if (child == 1)
2773  {
2774  pm = PointMatrix(mid30, mid12, pm(2), pm(3));
2775  }
2776  }
2777  else if (ref_type == 3) // iso split
2778  {
2779  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2780  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2781  Point midel(mid01, mid23);
2782
2783  if (child == 0)
2784  {
2785  pm = PointMatrix(pm(0), mid01, midel, mid30);
2786  }
2787  else if (child == 1)
2788  {
2789  pm = PointMatrix(mid01, pm(1), mid12, midel);
2790  }
2791  else if (child == 2)
2792  {
2793  pm = PointMatrix(midel, mid12, pm(2), mid23);
2794  }
2795  else if (child == 3)
2796  {
2797  pm = PointMatrix(mid30, midel, mid23, pm(3));
2798  }
2799  }
2800  }
2801  else if (geom == Geometry::TRIANGLE)
2802  {
2803  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
2804
2805  if (child == 0)
2806  {
2807  pm = PointMatrix(pm(0), mid01, mid20);
2808  }
2809  else if (child == 1)
2810  {
2811  pm = PointMatrix(mid01, pm(1), mid12);
2812  }
2813  else if (child == 2)
2814  {
2815  pm = PointMatrix(mid20, mid12, pm(2));
2816  }
2817  else if (child == 3)
2818  {
2819  pm = PointMatrix(mid01, mid12, mid20);
2820  }
2821  }
2822  }
2823
2824  // write the points to the matrix
2825  for (int i = 0; i < pm.np; i++)
2826  {
2827  for (int j = 0; j < pm(i).dim; j++)
2828  {
2829  matrix(j, i) = pm(i).coord[j];
2830  }
2831  }
2832 }
2833
2835 {
2838
2839  for (int i = 0; i < leaf_elements.Size(); i++)
2840  {
2841  int elem = leaf_elements[i];
2842  if (!IsGhost(elements[elem])) { coarse_elements.Append(elem); }
2843  }
2844
2845  transforms.embeddings.DeleteAll();
2846 }
2847
2848 void NCMesh::TraverseRefinements(int elem, int coarse_index,
2849  std::string &ref_path, RefPathMap &map)
2850 {
2851  Element &el = elements[elem];
2852  if (!el.ref_type)
2853  {
2854  int &matrix = map[ref_path];
2855  if (!matrix) { matrix = map.size(); }
2856
2857  Embedding &emb = transforms.embeddings[el.index];
2858  emb.parent = coarse_index;
2859  emb.matrix = matrix - 1;
2860  }
2861  else
2862  {
2863  ref_path.push_back(el.ref_type);
2864  ref_path.push_back(0);
2865
2866  for (int i = 0; i < 8; i++)
2867  {
2868  if (el.child[i] >= 0)
2869  {
2870  ref_path[ref_path.length()-1] = i;
2871  TraverseRefinements(el.child[i], coarse_index, ref_path, map);
2872  }
2873  }
2874  ref_path.resize(ref_path.length()-2);
2875  }
2876 }
2877
2879 {
2880  MFEM_VERIFY(coarse_elements.Size() || !leaf_elements.Size(),
2881  "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
2882  " and Refine().");
2883
2884  if (!transforms.embeddings.Size())
2885  {
2887
2888  std::string ref_path;
2889  ref_path.reserve(100);
2890
2891  RefPathMap map;
2892  map[ref_path] = 1; // identity
2893
2894  for (int i = 0; i < coarse_elements.Size(); i++)
2895  {
2896  TraverseRefinements(coarse_elements[i], i, ref_path, map);
2897  }
2898
2899  MFEM_ASSERT(elements.Size() > free_element_ids.Size(), "");
2900  int geom = elements[0].geom;
2901  const PointMatrix &identity = GetGeomIdentity(geom);
2902
2903  transforms.point_matrices.SetSize(Dim, identity.np, map.size());
2904
2905  // calculate the point matrices
2906  for (RefPathMap::iterator it = map.begin(); it != map.end(); ++it)
2907  {
2908  GetPointMatrix(geom, it->first.c_str(),
2909  transforms.point_matrices(it->second-1));
2910  }
2911  }
2912  return transforms;
2913 }
2914
2916 {
2917  MFEM_VERIFY(transforms.embeddings.Size() || !leaf_elements.Size(),
2918  "GetDerefinementTransforms() must be preceded by Derefine().");
2919
2921  {
2922  std::map<int, int> mat_no;
2923  mat_no[0] = 1; // identity
2924
2925  // assign numbers to the different matrices used
2926  for (int i = 0; i < transforms.embeddings.Size(); i++)
2927  {
2928  int code = transforms.embeddings[i].matrix;
2929  if (code)
2930  {
2931  int &matrix = mat_no[code];
2932  if (!matrix) { matrix = mat_no.size(); }
2933  transforms.embeddings[i].matrix = matrix - 1;
2934  }
2935  }
2936
2937  MFEM_ASSERT(elements.Size() > free_element_ids.Size(), "");
2938  int geom = elements[0].geom;
2939  const PointMatrix &identity = GetGeomIdentity(geom);
2940
2941  transforms.point_matrices.SetSize(Dim, identity.np, mat_no.size());
2942
2943  std::map<int, int>::iterator it;
2944  for (it = mat_no.begin(); it != mat_no.end(); ++it)
2945  {
2946  char path[3];
2947  int code = it->first;
2948  path[0] = code >> 3; // ref_type (see SetDerefMatrixCodes())
2949  path[1] = code & 7; // child
2950  path[2] = 0;
2951
2952  GetPointMatrix(geom, path, transforms.point_matrices(it->second-1));
2953  }
2954  }
2955  return transforms;
2956 }
2957
2959 {
2961  transforms.embeddings.DeleteAll();
2963 }
2964
2965
2966 //// Utility ///////////////////////////////////////////////////////////////////
2967
2968 void NCMesh::GetEdgeVertices(const MeshId &edge_id, int vert_index[2]) const
2969 {
2970  const Element &el = elements[edge_id.element];
2971  const GeomInfo& gi = GI[(int) el.geom];
2972  const int* ev = gi.edges[edge_id.local];
2973
2974  int n0 = el.node[ev[0]], n1 = el.node[ev[1]];
2975  if (n0 > n1) { std::swap(n0, n1); }
2976
2977  vert_index[0] = nodes[n0].vert_index;
2978  vert_index[1] = nodes[n1].vert_index;
2979
2980  if (vert_index[0] > vert_index[1])
2981  {
2982  std::swap(vert_index[0], vert_index[1]);
2983  }
2984 }
2985
2987 {
2988  const Element &el = elements[edge_id.element];
2989  const GeomInfo& gi = GI[(int) el.geom];
2990  const int* ev = gi.edges[edge_id.local];
2991
2992  int v0 = nodes[el.node[ev[0]]].vert_index;
2993  int v1 = nodes[el.node[ev[1]]].vert_index;
2994
2995  return ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1])) ? -1 : 1;
2996 }
2997
2999  int vert_index[4], int edge_index[4],
3000  int edge_orientation[4]) const
3001 {
3002  const Element &el = elements[face_id.element];
3003  const int* fv = GI[(int) el.geom].faces[face_id.local];
3004
3005  for (int i = 0; i < 4; i++)
3006  {
3007  vert_index[i] = nodes[el.node[fv[i]]].vert_index;
3008  }
3009
3010  for (int i = 0; i < 4; i++)
3011  {
3012  int j = (i+1) & 0x3;
3013  int n1 = el.node[fv[i]];
3014  int n2 = el.node[fv[j]];
3015
3016  const Node* en = nodes.Find(n1, n2);
3018
3019  edge_index[i] = en->edge_index;
3020  edge_orientation[i] = (vert_index[i] < vert_index[j]) ? 1 : -1;
3021  }
3022 }
3023
3024 int NCMesh::GetEdgeMaster(int node) const
3025 {
3027  const Node &nd = nodes[node];
3028
3029  int p1 = nd.p1, p2 = nd.p2;
3030  MFEM_ASSERT(p1 != p2, "invalid edge node.");
3031
3032  const Node &n1 = nodes[p1], &n2 = nodes[p2];
3033
3034  int n1p1 = n1.p1, n1p2 = n1.p2;
3035  int n2p1 = n2.p1, n2p2 = n2.p2;
3036
3037  if ((n2p1 != n2p2) && (p1 == n2p1 || p1 == n2p2))
3038  {
3039  // n1 is parent of n2:
3040  // (n1)--(nd)--(n2)------(*)
3041  if (n2.HasEdge()) { return p2; }
3042  else { return GetEdgeMaster(p2); }
3043  }
3044
3045  if ((n1p1 != n1p2) && (p2 == n1p1 || p2 == n1p2))
3046  {
3047  // n2 is parent of n1:
3048  // (n2)--(nd)--(n1)------(*)
3049  if (n1.HasEdge()) { return p1; }
3050  else { return GetEdgeMaster(p1); }
3051  }
3052
3053  return -1;
3054 }
3055
3056 int NCMesh::GetEdgeMaster(int v1, int v2) const
3057 {
3058  int node = nodes.FindId(vertex_nodeId[v1], vertex_nodeId[v2]);
3059  MFEM_ASSERT(node >= 0 && nodes[node].HasEdge(), "(v1, v2) is not an edge.");
3060
3061  int master = GetEdgeMaster(node);
3062  return (master >= 0) ? nodes[master].edge_index : -1;
3063 }
3064
3065 int NCMesh::GetElementDepth(int i) const
3066 {
3067  int elem = leaf_elements[i];
3068  int depth = 0, parent;
3069  while ((parent = elements[elem].parent) != -1)
3070  {
3071  elem = parent;
3072  depth++;
3073  }
3074  return depth;
3075 }
3076
3077 void NCMesh::FindFaceNodes(int face, int node[4])
3078 {
3079  // Obtain face nodes from one of its elements (note that face->p1, p2, p3
3080  // cannot be used directly since they are not in order and p4 is missing).
3081
3082  Face &fa = faces[face];
3083
3084  int elem = fa.elem[0];
3085  if (elem < 0) { elem = fa.elem[1]; }
3086  MFEM_ASSERT(elem >= 0, "Face has no elements?");
3087
3088  Element &el = elements[elem];
3089  int f = find_hex_face(find_node(el, fa.p1),
3090  find_node(el, fa.p2),
3091  find_node(el, fa.p3));
3092
3093  const int* fv = GI[Geometry::CUBE].faces[f];
3094  for (int i = 0; i < 4; i++)
3095  {
3096  node[i] = el.node[fv[i]];
3097  }
3098 }
3099
3100 void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
3101  Array<int> &bdr_vertices, Array<int> &bdr_edges)
3102 {
3103  bdr_vertices.SetSize(0);
3104  bdr_edges.SetSize(0);
3105
3106  if (Dim == 3)
3107  {
3108  GetFaceList(); // make sure 'boundary_faces' is up to date
3109
3110  for (int i = 0; i < boundary_faces.Size(); i++)
3111  {
3112  int face = boundary_faces[i];
3113  if (bdr_attr_is_ess[faces[face].attribute - 1])
3114  {
3115  int node[4];
3116  FindFaceNodes(face, node);
3117
3118  for (int j = 0; j < 4; j++)
3119  {
3120  bdr_vertices.Append(nodes[node[j]].vert_index);
3121
3122  int enode = nodes.FindId(node[j], node[(j+1) % 4]);
3124  bdr_edges.Append(nodes[enode].edge_index);
3125
3126  while ((enode = GetEdgeMaster(enode)) >= 0)
3127  {
3128  // append master edges that may not be accessible from any
3129  // boundary element, this happens in 3D in re-entrant corners
3130  bdr_edges.Append(nodes[enode].edge_index);
3131  }
3132  }
3133  }
3134  }
3135  }
3136  else if (Dim == 2)
3137  {
3138  GetEdgeList(); // make sure 'boundary_faces' is up to date
3139
3140  for (int i = 0; i < boundary_faces.Size(); i++)
3141  {
3142  int face = boundary_faces[i];
3143  Face &fc = faces[face];
3144  if (bdr_attr_is_ess[fc.attribute - 1])
3145  {
3146  bdr_vertices.Append(nodes[fc.p1].vert_index);
3147  bdr_vertices.Append(nodes[fc.p3].vert_index);
3148  }
3149  }
3150  }
3151
3152  bdr_vertices.Sort();
3153  bdr_vertices.Unique();
3154
3155  bdr_edges.Sort();
3156  bdr_edges.Unique();
3157 }
3158
3159 int NCMesh::EdgeSplitLevel(int vn1, int vn2) const
3160 {
3161  int mid = nodes.FindId(vn1, vn2);
3162  if (mid < 0 || !nodes[mid].HasVertex()) { return 0; }
3163  return 1 + std::max(EdgeSplitLevel(vn1, mid), EdgeSplitLevel(mid, vn2));
3164 }
3165
3166 void NCMesh::FaceSplitLevel(int vn1, int vn2, int vn3, int vn4,
3167  int& h_level, int& v_level) const
3168 {
3169  int hl1, hl2, vl1, vl2;
3170  int mid[4];
3171
3172  switch (FaceSplitType(vn1, vn2, vn3, vn4, mid))
3173  {
3174  case 0: // not split
3175  h_level = v_level = 0;
3176  break;
3177
3178  case 1: // vertical
3179  FaceSplitLevel(vn1, mid[0], mid[2], vn4, hl1, vl1);
3180  FaceSplitLevel(mid[0], vn2, vn3, mid[2], hl2, vl2);
3181  h_level = std::max(hl1, hl2);
3182  v_level = std::max(vl1, vl2) + 1;
3183  break;
3184
3185  default: // horizontal
3186  FaceSplitLevel(vn1, vn2, mid[1], mid[3], hl1, vl1);
3187  FaceSplitLevel(mid[3], mid[1], vn3, vn4, hl2, vl2);
3188  h_level = std::max(hl1, hl2) + 1;
3189  v_level = std::max(vl1, vl2);
3190  }
3191 }
3192
3193 static int max8(int a, int b, int c, int d, int e, int f, int g, int h)
3194 {
3195  return std::max(std::max(std::max(a, b), std::max(c, d)),
3196  std::max(std::max(e, f), std::max(g, h)));
3197 }
3198
3199 void NCMesh::CountSplits(int elem, int splits[3]) const
3200 {
3201  const Element &el = elements[elem];
3202  const int* node = el.node;
3203  GeomInfo& gi = GI[(int) el.geom];
3204
3205  int elevel[12];
3206  for (int i = 0; i < gi.ne; i++)
3207  {
3208  const int* ev = gi.edges[i];
3209  elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
3210  }
3211
3212  if (el.geom == Geometry::CUBE)
3213  {
3214  int flevel[6][2];
3215  for (int i = 0; i < gi.nf; i++)
3216  {
3217  const int* fv = gi.faces[i];
3218  FaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]],
3219  flevel[i][1], flevel[i][0]);
3220  }
3221
3222  splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
3223  elevel[0], elevel[2], elevel[4], elevel[6]);
3224
3225  splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
3226  elevel[1], elevel[3], elevel[5], elevel[7]);
3227
3228  splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
3229  elevel[8], elevel[9], elevel[10], elevel[11]);
3230  }
3231  else if (el.geom == Geometry::SQUARE)
3232  {
3233  splits[0] = std::max(elevel[0], elevel[2]);
3234  splits[1] = std::max(elevel[1], elevel[3]);
3235  }
3236  else if (el.geom == Geometry::TRIANGLE)
3237  {
3238  splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
3239  splits[1] = splits[0];
3240  }
3241  else
3242  {
3243  MFEM_ABORT("Unsupported element geometry.");
3244  }
3245 }
3246
3247 void NCMesh::GetLimitRefinements(Array<Refinement> &refinements, int max_level)
3248 {
3249  for (int i = 0; i < leaf_elements.Size(); i++)
3250  {
3251  if (IsGhost(elements[leaf_elements[i]])) { break; }
3252
3253  int splits[3];
3254  CountSplits(leaf_elements[i], splits);
3255
3256  char ref_type = 0;
3257  for (int k = 0; k < Dim; k++)
3258  {
3259  if (splits[k] > max_level)
3260  {
3261  ref_type |= (1 << k);
3262  }
3263  }
3264
3265  if (ref_type)
3266  {
3267  if (Iso)
3268  {
3269  // iso meshes should only be modified by iso refinements
3270  ref_type = 7;
3271  }
3272  refinements.Append(Refinement(i, ref_type));
3273  }
3274  }
3275 }
3276
3277 void NCMesh::LimitNCLevel(int max_nc_level)
3278 {
3279  MFEM_VERIFY(max_nc_level >= 1, "'max_nc_level' must be 1 or greater.");
3280
3281  while (1)
3282  {
3283  Array<Refinement> refinements;
3284  GetLimitRefinements(refinements, max_nc_level);
3285
3286  if (!refinements.Size()) { break; }
3287
3288  Refine(refinements);
3289  }
3290 }
3291
3292 void NCMesh::PrintVertexParents(std::ostream &out) const
3293 {
3294  // count vertices with parents
3295  int nv = 0;
3296  for (node_const_iterator node = nodes.cbegin(); node != nodes.cend(); ++node)
3297  {
3298  if (node->HasVertex() && node->p1 != node->p2) { nv++; }
3299  }
3300  out << nv << "\n";
3301
3302  // print the relations
3303  for (node_const_iterator node = nodes.cbegin(); node != nodes.cend(); ++node)
3304  {
3305  if (node->HasVertex() && node->p1 != node->p2)
3306  {
3307  const Node &p1 = nodes[node->p1];
3308  const Node &p2 = nodes[node->p2];
3309
3310  MFEM_ASSERT(p1.HasVertex(), "");
3311  MFEM_ASSERT(p2.HasVertex(), "");
3312
3313  out << node->vert_index << " "
3314  << p1.vert_index << " " << p2.vert_index << "\n";
3315  }
3316  }
3317 }
3318
3320 {
3321  int nv;
3322  input >> nv;
3323  while (nv--)
3324  {
3325  int id, p1, p2;
3326  input >> id >> p1 >> p2;
3327  MFEM_VERIFY(input, "problem reading vertex parents.");
3328
3332
3333  // assign new parents for the node
3334  nodes.Reparent(id, p1, p2);
3335
3336  // NOTE: when loading an AMR mesh, node indices are guaranteed to have
3337  // the same indices as vertices, see NCMesh::NCMesh.
3338  }
3339 }
3340
3342 {
3343  int num_top_level = 0;
3344  for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
3345  {
3346  if (node->p1 == node->p2) // see NCMesh::NCMesh
3347  {
3348  MFEM_VERIFY(node.index() == node->p1, "invalid top-level vertex.");
3350  MFEM_VERIFY(node->vert_index == node->p1, "bad top-level vertex index");
3351  num_top_level = std::max(num_top_level, node->p1 + 1);
3352  }
3353  }
3354
3355  top_vertex_pos.SetSize(3*num_top_level);
3356  for (int i = 0; i < num_top_level; i++)
3357  {
3358  memcpy(&top_vertex_pos[3*i], mvertices[i](), 3*sizeof(double));
3359  }
3360 }
3361
3362 static int ref_type_num_children[8] = { 0, 2, 2, 4, 2, 4, 4, 8 };
3363
3364 int NCMesh::PrintElements(std::ostream &out, int elem, int &coarse_id) const
3365 {
3366  const Element &el = elements[elem];
3367  if (el.ref_type)
3368  {
3369  int child_id[8], nch = 0;
3370  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
3371  {
3372  child_id[nch++] = PrintElements(out, el.child[i], coarse_id);
3373  }
3374  MFEM_ASSERT(nch == ref_type_num_children[(int) el.ref_type], "");
3375
3376  out << (int) el.ref_type;
3377  for (int i = 0; i < nch; i++)
3378  {
3379  out << " " << child_id[i];
3380  }
3381  out << "\n";
3382  return coarse_id++; // return new id for this coarse element
3383  }
3384  else
3385  {
3386  return el.index;
3387  }
3388 }
3389
3390 void NCMesh::PrintCoarseElements(std::ostream &out) const
3391 {
3392  // print the number of non-leaf elements
3393  out << (elements.Size() - free_element_ids.Size() - leaf_elements.Size())
3394  << "\n";
3395
3396  // print the hierarchy recursively
3397  int coarse_id = leaf_elements.Size();
3398  for (int i = 0; i < root_count; i++)
3399  {
3400  PrintElements(out, i, coarse_id);
3401  }
3402 }
3403
3404 void NCMesh::CopyElements(int elem,
3405  const BlockArray<Element> &tmp_elements,
3406  Array<int> &index_map)
3407 {
3408  Element &el = elements[elem];
3409  if (el.ref_type)
3410  {
3411  for (int i = 0; i < 8 && el.child[i] >= 0; i++)
3412  {
3413  int old_id = el.child[i];
3414  // here, we do not use the content of 'free_element_ids', if any
3415  int new_id = elements.Append(tmp_elements[old_id]);
3416  index_map[old_id] = new_id;
3417  el.child[i] = new_id;
3418  elements[new_id].parent = elem;
3419  CopyElements(new_id, tmp_elements, index_map);
3420  }
3421  }
3422 }
3423
3425 {
3426  int ne;
3427  input >> ne;
3428
3429  bool iso = true;
3430
3431  // load the coarse elements
3432  while (ne--)
3433  {
3434  int ref_type;
3435  input >> ref_type;
3436
3437  int elem = AddElement(Element(0, 0));
3438  Element &el = elements[elem];
3439  el.ref_type = ref_type;
3440
3441  if (Dim == 3 && ref_type != 7) { iso = false; }
3442
3444  int nch = ref_type_num_children[ref_type];
3445  for (int i = 0, id; i < nch; i++)
3446  {
3447  input >> id;
3448  MFEM_VERIFY(id >= 0, "");
3449  MFEM_VERIFY(id < leaf_elements.Size() ||
3450  id < elements.Size()-free_element_ids.Size(),
3451  "coarse element cannot be referenced before it is "
3452  "defined (id=" << id << ").");
3453
3454  Element &child = elements[id];
3455  MFEM_VERIFY(child.parent == -1,
3456  "element " << id << " cannot have two parents.");
3457
3458  el.child[i] = id;
3459  child.parent = elem;
3460
3461  if (!i) // copy geom and attribute from first child
3462  {
3463  el.geom = child.geom;
3464  el.attribute = child.attribute;
3465  }
3466  }
3467  }
3468
3469  // prepare for reordering the elements
3470  BlockArray<Element> tmp_elements;
3471  elements.Swap(tmp_elements);
3473
3474  Array<int> index_map(tmp_elements.Size());
3475  index_map = -1;
3476
3477  // copy roots, they need to be at the beginning of 'elements'
3478  root_count = 0;
3479  for (elem_iterator el = tmp_elements.begin(); el != tmp_elements.end(); ++el)
3480  {
3481  if (el->parent == -1)
3482  {
3483  int new_id = elements.Append(*el); // same as AddElement()
3484  index_map[el.index()] = new_id;
3485  root_count++;
3486  }
3487  }
3488
3489  // copy the rest of the hierarchy
3490  for (int i = 0; i < root_count; i++)
3491  {
3492  CopyElements(i, tmp_elements, index_map);
3493  }
3494
3495  // we also need to renumber element links in Face::elem[]
3496  for (face_iterator face = faces.begin(); face != faces.end(); ++face)
3497  {
3498  for (int i = 0; i < 2; i++)
3499  {
3500  if (face->elem[i] >= 0)
3501  {
3502  face->elem[i] = index_map[face->elem[i]];
3503  MFEM_ASSERT(face->elem[i] >= 0, "");
3504  }
3505  }
3506  }
3507
3508  // set the Iso flag (must be false if there are 3D aniso refinements)
3509  Iso = iso;
3510
3511  Update();
3512 }
3513
3515 {
3516  vertex_list.Clear(true);
3517  face_list.Clear(true);
3518  edge_list.Clear(true);
3519
3522
3523  ClearTransforms();
3524 }
3525
3527 {
3528  int pmsize = 0;
3529  if (slaves.size())
3530  {
3531  pmsize = slaves[0].point_matrix.MemoryUsage();
3532  }
3533
3534  return conforming.capacity() * sizeof(MeshId) +
3535  masters.capacity() * sizeof(Master) +
3536  slaves.capacity() * sizeof(Slave) +
3537  slaves.size() * pmsize;
3538 }
3539
3541 {
3542  return point_matrices.MemoryUsage() + embeddings.MemoryUsage();
3543 }
3544
3546 {
3547  return nodes.MemoryUsage() +
3548  faces.MemoryUsage() +
3549  elements.MemoryUsage() +
3559  ref_stack.MemoryUsage() +
3563  sizeof(*this);
3564 }
3565
3567 {
3568  nodes.PrintMemoryDetail(); mfem::out << " nodes\n";
3569  faces.PrintMemoryDetail(); mfem::out << " faces\n";
3570
3571  mfem::out << elements.MemoryUsage() << " elements\n"
3572  << free_element_ids.MemoryUsage() << " free_element_ids\n"
3573  << top_vertex_pos.MemoryUsage() << " top_vertex_pos\n"
3574  << leaf_elements.MemoryUsage() << " leaf_elements\n"
3575  << vertex_nodeId.MemoryUsage() << " vertex_nodeId\n"
3576  << face_list.MemoryUsage() << " face_list\n"
3577  << edge_list.MemoryUsage() << " edge_list\n"
3578  << vertex_list.MemoryUsage() << " vertex_list\n"
3579  << boundary_faces.MemoryUsage() << " boundary_faces\n"
3580  << element_vertex.MemoryUsage() << " element_vertex\n"
3581  << ref_stack.MemoryUsage() << " ref_stack\n"
3582  << derefinements.MemoryUsage() << " derefinements\n"
3583  << transforms.MemoryUsage() << " transforms\n"
3584  << coarse_elements.MemoryUsage() << " coarse_elements\n"
3585  << sizeof(*this) << " NCMesh"
3586  << std::endl;
3587
3588  return elements.Size() - free_element_ids.Size();
3589 }
3590
3591 void NCMesh::PrintStats(std::ostream &out) const
3592 {
3593  static const double MiB = 1024.*1024.;
3594  out <<
3595  "NCMesh statistics:\n"
3596  "------------------\n"
3597  " mesh and space dimensions : " << Dim << ", " << spaceDim << "\n"
3598  " isotropic only : " << (Iso ? "yes" : "no") << "\n"
3599  " number of Nodes : " << std::setw(9)
3600  << nodes.Size() << " + [ " << std::setw(9)
3601  << nodes.MemoryUsage()/MiB << " MiB ]\n"
3602  " free " << std::setw(9)
3603  << nodes.NumFreeIds() << "\n"
3604  " number of Faces : " << std::setw(9)
3605  << faces.Size() << " + [ " << std::setw(9)
3606  << faces.MemoryUsage()/MiB << " MiB ]\n"
3607  " free " << std::setw(9)
3608  << faces.NumFreeIds() << "\n"
3609  " number of Elements : " << std::setw(9)
3610  << elements.Size()-free_element_ids.Size() << " + [ " << std::setw(9)
3611  << (elements.MemoryUsage() +
3612  free_element_ids.MemoryUsage())/MiB << " MiB ]\n"
3613  " free " << std::setw(9)
3614  << free_element_ids.Size() << "\n"
3615  " number of root elements : " << std::setw(9) << root_count << "\n"
3616  " number of leaf elements : " << std::setw(9)
3617  << leaf_elements.Size() << "\n"
3618  " number of vertices : " << std::setw(9)
3619  << vertex_nodeId.Size() << "\n"
3620  " number of faces : " << std::setw(9)
3621  << face_list.TotalSize() << " = [ " << std::setw(9)
3622  << face_list.MemoryUsage()/MiB << " MiB ]\n"
3623  " conforming " << std::setw(9)
3624  << face_list.conforming.size() << " +\n"
3625  " master " << std::setw(9)
3626  << face_list.masters.size() << " +\n"
3627  " slave " << std::setw(9)
3628  << face_list.slaves.size() << "\n"
3629  " number of edges : " << std::setw(9)
3630  << edge_list.TotalSize() << " = [ " << std::setw(9)
3631  << edge_list.MemoryUsage()/MiB << " MiB ]\n"
3632  " conforming " << std::setw(9)
3633  << edge_list.conforming.size() << " +\n"
3634  " master " << std::setw(9)
3635  << edge_list.masters.size() << " +\n"
3636  " slave " << std::setw(9)
3637  << edge_list.slaves.size() << "\n"
3638  " total memory : " << std::setw(17)
3639  << "[ " << std::setw(9) << MemoryUsage()/MiB << " MiB ]\n"
3640  ;
3641 }
3642
3643 #ifdef MFEM_DEBUG
3644 void NCMesh::DebugDump(std::ostream &out) const
3645 {
3646  // dump nodes
3647  tmp_vertex = new TmpVertex[nodes.NumIds()];
3648  out << nodes.Size() << "\n";
3649  for (node_const_iterator node = nodes.cbegin(); node != nodes.cend(); ++node)
3650  {
3651  const double *pos = CalcVertexPos(node.index());
3652  out << node.index() << " "
3653  << pos[0] << " " << pos[1] << " " << pos[2] << " "
3654  << node->p1 << " " << node->p2 << " "
3655  << node->vert_index << " " << node->edge_index << " "
3656  << 0 << "\n";
3657  }
3658  delete [] tmp_vertex;
3659  out << "\n";
3660
3661  // dump elements
3662  int nleaves = 0;
3663  for (int i = 0; i < elements.Size(); i++)
3664  {
3665  const Element &el = elements[i];
3666  if (!el.ref_type && el.parent != -2 /*freed*/) { nleaves++; }
3667  }
3668  out << nleaves << "\n";
3669  for (int i = 0; i < elements.Size(); i++)
3670  {
3671  const Element &el = elements[i];
3672  if (el.ref_type || el.parent == -2) { continue; }
3673  const GeomInfo& gi = GI[(int) el.geom];
3674  out << gi.nv << " ";
3675  for (int j = 0; j < gi.nv; j++)
3676  {
3677  out << el.node[j] << " ";
3678  }
3679  out << el.attribute << " " << el.rank << " " << i << "\n";
3680  }
3681  out << "\n";
3682
3683  // dump faces
3684  out << faces.Size() << "\n";
3685  for (face_const_iterator face = faces.cbegin(); face != faces.cend(); ++face)
3686  {
3687  int elem = face->elem[0];
3688  if (elem < 0) { elem = face->elem[1]; }
3689  MFEM_ASSERT(elem >= 0, "");
3690  const Element &el = elements[elem];
3691
3692  int lf = find_hex_face(find_node(el, face->p1),
3693  find_node(el, face->p2),
3694  find_node(el, face->p3));
3695
3696  out << "4";
3697  const int* fv = GI[Geometry::CUBE].faces[lf];
3698  for (int i = 0; i < 4; i++)
3699  {
3700  out << " " << el.node[fv[i]];
3701  }
3702  //out << " # face " << face.index() << ", index " << face->index << "\n";
3703  out << "\n";
3704  }
3705 }
3706 #endif
3707
3708 } // namespace mfem
NCList face_list
lazy-initialized list of faces, see GetFaceList
Definition: ncmesh.hpp:430
void CollectLeafElements(int elem, int state)
Definition: ncmesh.cpp:1423
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
Definition: ncmesh.hpp:431
NCMesh(const Mesh *mesh, std::istream *vertex_parents=NULL)
Definition: ncmesh.cpp:67
int Size() const
Logical size of the array.
Definition: array.hpp:133
CoarseFineTransformations transforms
storage for data returned by Get[De]RefinementTransforms()
Definition: ncmesh.hpp:672
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:358
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:651
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:37
TmpVertex * tmp_vertex
Definition: ncmesh.hpp:690
Array< double > top_vertex_pos
Definition: ncmesh.hpp:406
void Unique()
Definition: array.hpp:245
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
Definition: mesh.cpp:3819
void SetDerefMatrixCodes(int parent, Array< int > &fine_coarse)
Definition: ncmesh.cpp:1353
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:2915
virtual int GetNEdges() const =0
bool HasEdge() const
Definition: ncmesh.hpp:343
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:3292
void GetPointMatrix(int geom, const char *ref_path, DenseMatrix &matrix)
Definition: ncmesh.cpp:2537
virtual void LimitNCLevel(int max_nc_level)
Definition: ncmesh.cpp:3277
void GetFaceVerticesEdges(const MeshId &face_id, int vert_index[4], int edge_index[4], int edge_orientation[4]) const
Return Mesh vertex and edge indices of a face identified by &#39;face_id&#39;.
Definition: ncmesh.cpp:2998
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:621
virtual void Trim()
Save memory by releasing all non-essential and cached data.
Definition: ncmesh.cpp:3514
int index
element number in the Mesh, -1 if refined
Definition: ncmesh.hpp:381
const Geometry::Type geom
Definition: ex1.cpp:40
void CopyElements(int elem, const BlockArray< Element > &tmp_elements, Array< int > &index_map)
Definition: ncmesh.cpp:3404
void ForceRefinement(int vn1, int vn2, int vn3, int vn4)
Definition: ncmesh.cpp:549
virtual int GetNumGhostElements() const
Definition: ncmesh.hpp:447
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:169
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:1528
int Size() const
Return the number of items actually stored.
Definition: array.hpp:437
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:680
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:190
void TraverseEdge(int vn0, int vn1, double t0, double t1, int flags, int level)
Definition: ncmesh.cpp:1874
T * GetData()
Returns the data.
Definition: array.hpp:115
int NVertices
Definition: ncmesh.hpp:424
void ClearTransforms()
Free all internal data created by the above three functions.
Definition: ncmesh.cpp:2958
Table derefinements
possible derefinements, see GetDerefinementTable
Definition: ncmesh.hpp:455
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:2157
NCList vertex_list
lazy-initialized list of vertices, see GetVertexList
Definition: ncmesh.hpp:432
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:189
char geom
Geometry::Type of the element.
Definition: ncmesh.hpp:378
long TotalSize() const
Definition: ncmesh.cpp:2087
void GetEdgeVertices(const MeshId &edge_id, int vert_index[2]) const
Return Mesh vertex indices of an edge identified by &#39;edge_id&#39;.
Definition: ncmesh.cpp:2968
int GetEdgeNCOrientation(const MeshId &edge_id) const
Definition: ncmesh.cpp:2986
void SetSize(int i, int j, int k)
Definition: densemat.hpp:682
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:618
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
Definition: ncmesh.cpp:2048
const Element * GetFace(int i) const
Definition: mesh.hpp:684
static GeomInfo & gi_tri
Definition: ncmesh.hpp:729
void CheckAnisoFace(int vn1, int vn2, int vn3, int vn4, int mid12, int mid34, int level=0)
Definition: ncmesh.cpp:583
virtual int GetNumGhostVertices() const
Definition: ncmesh.hpp:448
static PointMatrix pm_tri_identity
Definition: ncmesh.hpp:658
void FreeElement(int id)
Definition: ncmesh.hpp:471
void CountSplits(int elem, int splits[3]) const
Definition: ncmesh.cpp:3199
void CollectDerefinements(int elem, Array< Connection > &list)
Definition: ncmesh.cpp:1225
std::vector< MeshId > conforming
Definition: ncmesh.hpp:171
void GetLimitRefinements(Array< Refinement > &refinements, int max_level)
Definition: ncmesh.cpp:3247
int PrintMemoryDetail() const
Definition: ncmesh.cpp:3566
int spaceDim
dimensions of the elements and the vertex coordinates
Definition: ncmesh.hpp:323
Array< int > vertex_nodeId
Definition: ncmesh.hpp:428
int attribute
boundary element attribute, -1 if internal face
Definition: ncmesh.hpp:356
void DebugDump(std::ostream &out) const
Definition: ncmesh.cpp:3644
void FindFaceNodes(int face, int node[4])
Definition: ncmesh.cpp:3077
int GetGeometryType() const
Definition: element.hpp:47
void DeleteAll()
Delete whole array.
Definition: array.hpp:709
int index
Mesh element number.
Definition: ncmesh.hpp:36
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
Definition: ncmesh.hpp:188
const MeshId & LookUp(int index, int *type=NULL) const
Definition: ncmesh.cpp:2092
virtual void ElementSharesEdge(int elem, int enode)
Definition: ncmesh.hpp:545
virtual void OnMeshUpdated(Mesh *mesh)
Definition: ncmesh.cpp:1603
void CollectEdgeVertices(int v0, int v1, Array< int > &indices)
Definition: ncmesh.cpp:2145
Array< Refinement > ref_stack
stack of scheduled refinements (temporary)
Definition: ncmesh.hpp:453
iterator begin()
Definition: array.hpp:528
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3374
virtual bool IsGhost(const Element &el) const
Definition: ncmesh.hpp:446
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:403
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:2848
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:3364
int index
face number in the Mesh
Definition: ncmesh.hpp:357
const Table & GetDerefinementTable()
Definition: ncmesh.cpp:1258
void MakeFromList(int nrows, const Array< Connection > &list)
Definition: table.cpp:274
Table element_vertex
leaf-element to vertex table, see FindSetNeighbors
Definition: ncmesh.hpp:436
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:138
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:614
void GetMatrix(DenseMatrix &point_matrix) const
Definition: ncmesh.cpp:2501
Definition: ncmesh.hpp:460
void Clear()
Definition: table.cpp:373
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:2356
std::vector< Master > masters
Definition: ncmesh.hpp:172
Array< int > coarse_elements
state of leaf_elements before Refine(), set by MarkCoarseLevel()
Definition: ncmesh.hpp:675
virtual void UpdateVertices()
update Vertex::index and vertex_nodeId
Definition: ncmesh.cpp:1372
long MemoryUsage() const
Definition: table.cpp:406
void RegisterFaces(int elem, int *fattr=NULL)
Definition: ncmesh.cpp:305
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:146
Array< Embedding > embeddings
fine element positions in their parents
Definition: ncmesh.hpp:55
virtual void BuildEdgeList()
Definition: ncmesh.cpp:1913
const CoarseFineTransformations & GetRefinementTransforms()
Definition: ncmesh.cpp:2878
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:3424
Identifies a vertex/edge/face in both Mesh and NCMesh.
Definition: ncmesh.hpp:134
int parent
parent element, -1 if this is a root element, -2 if free
Definition: ncmesh.hpp:389
Data type triangle element.
Definition: triangle.hpp:23
void MarkCoarseLevel()
Definition: ncmesh.cpp:2834
const Element * GetElement(int i) const
Definition: mesh.hpp:676
long MemoryUsage() const
Definition: array.hpp:269
void Sort()
Sorts the array. This requires operator&lt; to be defined for T.
Definition: array.hpp:237
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
Definition: ncmesh.hpp:83
int Dimension() const
Definition: mesh.hpp:645
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:3100
int element
NCMesh::Element containing this vertex/edge/face.
Definition: ncmesh.hpp:137
static GeomInfo & gi_hex
Definition: ncmesh.hpp:729
static GeomInfo GI[Geometry::NumGeom]
Definition: ncmesh.hpp:727
void Clear(bool hard=false)
Definition: ncmesh.cpp:2069
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:544
int SpaceDimension() const
Definition: mesh.hpp:646
bool Iso
true if the mesh only contains isotropic refinements
Definition: ncmesh.hpp:324
bool HasVertex() const
Definition: ncmesh.hpp:342
std::map< std::string, int > RefPathMap
Definition: ncmesh.hpp:666
int FaceSplitType(int v1, int v2, int v3, int v4, int mid[4]=NULL) const
Definition: ncmesh.cpp:1652
std::vector< Slave > slaves
Definition: ncmesh.hpp:173
virtual void BuildFaceList()
Definition: ncmesh.cpp:1803
Nonconforming edge/face within a bigger edge/face.
Definition: ncmesh.hpp:155
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:1721
void SetVertexPositions(const Array< mfem::Vertex > &vertices)
I/O: Set positions of all vertices (used by mesh loader).
Definition: ncmesh.cpp:3341
int edge_flags
edge orientation flags
Definition: ncmesh.hpp:158
bool Boundary() const
Definition: ncmesh.hpp:362
void RegisterElement(int e)
Definition: ncmesh.cpp:330
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
static const PointMatrix & GetGeomIdentity(int geom)
Definition: ncmesh.cpp:2524
mfem::Element * NewMeshElement(int geom) const
Definition: ncmesh.cpp:1487
void DeleteLast()
Delete the last entry.
Definition: array.hpp:181
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
Definition: ncmesh.hpp:729
long MemoryUsage() const
Definition: ncmesh.cpp:3526
void BuildElementToVertexTable()
Definition: ncmesh.cpp:2181
HashTable< Node > nodes
Definition: ncmesh.hpp:396
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:615
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:434
int child[8]
2-8 children (if ref_type != 0)
Definition: ncmesh.hpp:387
Array< int > leaf_elements
Definition: ncmesh.hpp:427
Definition: ncmesh.hpp:659
void TraverseFace(int vn0, int vn1, int vn2, int vn3, const PointMatrix &pm, int level)
Definition: ncmesh.cpp:1754
int EdgeSplitLevel(int vn1, int vn2) const
Definition: ncmesh.cpp:3159
void SetIJ(int *newI, int *newJ, int newsize=-1)
Definition: table.cpp:206
static int find_element_edge(const Element &el, int vn0, int vn1)
Definition: ncmesh.cpp:1689
void Initialize(const mfem::Element *elem)
Definition: ncmesh.cpp:28
Array< int > free_element_ids
Definition: ncmesh.hpp:400
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
virtual void ElementSharesVertex(int elem, int vnode)
Definition: ncmesh.hpp:546
static PointMatrix pm_hex_identity
Definition: ncmesh.hpp:660
static int find_hex_face(int a, int b, int c)
Definition: ncmesh.cpp:1705
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:647
void NeighborExpand(const Array< int > &elems, Array< int > &expanded, const Array< int > *search_set=NULL)
Definition: ncmesh.cpp:2439
BlockArray< Element > elements
Definition: ncmesh.hpp:399
void FindSetNeighbors(const Array< char > &elem_set, Array< int > *neighbors, Array< char > *neighbor_set=NULL)
Definition: ncmesh.cpp:2253
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:624
virtual void Update()
Definition: ncmesh.cpp:188
iterator end()
Definition: array.hpp:529
long MemoryUsage() const
Definition: ncmesh.cpp:3545
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:3390
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:52
virtual void BuildVertexList()
Definition: ncmesh.cpp:2013
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
Definition: ncmesh.hpp:195
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
void UpdateElementToVertexTable()
Definition: ncmesh.hpp:583
int RowSize(int i) const
Definition: table.hpp:108
int GetElementDepth(int i) const
Return the distance of leaf &#39;i&#39; from the root.
Definition: ncmesh.cpp:3065
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:379
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:1679
int index
Mesh number.
Definition: ncmesh.hpp:136
Definition: ncmesh.cpp:3319
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:3591
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:382
bool NodeSetX1(int node, int *n)
Definition: ncmesh.cpp:530
int node[8]
element corners (if ref_type == 0)
Definition: ncmesh.hpp:386
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:680
virtual int GetType() const =0
Returns element&#39;s type.
DenseMatrix point_matrix
position within the master edge/face
Definition: ncmesh.hpp:159
HashTable< Face > faces
Definition: ncmesh.hpp:397
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:3166
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:3056