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