MFEM  v3.0
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
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 "ncmesh.hpp"
15
16 #include <algorithm>
17 #include <cmath>
18
19 namespace mfem
20 {
21
24 struct GeomInfo
25 {
26  int nv, ne, nf, nfv; // number of: vertices, edge, faces, face vertices
27  int edges[12][2]; // edge vertices (up to 12 edges)
28  int faces[6][4]; // face vertices (up to 6 faces)
29
30  bool initialized;
31  GeomInfo() : initialized(false) {}
32  void Initialize(const Element* elem);
33 };
34
35 static GeomInfo GI[Geometry::NumGeom];
36
37 static GeomInfo& gi_hex = GI[Geometry::CUBE];
38 static GeomInfo& gi_quad = GI[Geometry::SQUARE];
39 static GeomInfo& gi_tri = GI[Geometry::TRIANGLE];
40
41 void GeomInfo::Initialize(const Element* elem)
42 {
43  if (initialized) return;
44
45  nv = elem->GetNVertices();
46  ne = elem->GetNEdges();
47  nf = elem->GetNFaces(nfv);
48
49  for (int i = 0; i < ne; i++)
50  for (int j = 0; j < 2; j++)
51  edges[i][j] = elem->GetEdgeVertices(i)[j];
52
53  for (int i = 0; i < nf; i++)
54  for (int j = 0; j < nfv; j++)
55  faces[i][j] = elem->GetFaceVertices(i)[j];
56
57  initialized = true;
58 }
59
60
61 NCMesh::NCMesh(const Mesh *mesh)
62 {
63  Dim = mesh->Dimension();
64
65  vertex_nodeId.SetSize(mesh->GetNV());
66  vertex_nodeId = -1;
67
68  // create the NCMesh::Element struct for each Mesh element
69  for (int i = 0; i < mesh->GetNE(); i++)
70  {
71  const mfem::Element *elem = mesh->GetElement(i);
72  const int *v = elem->GetVertices();
73
74  int geom = elem->GetGeometryType();
75  if (geom != Geometry::TRIANGLE &&
76  geom != Geometry::SQUARE &&
77  geom != Geometry::CUBE)
78  {
79  mfem_error("NCMesh: only triangles, quads and hexes are supported.");
80  }
81
82  // initialize edge/face tables for this type of element
83  GI[geom].Initialize(elem);
84
85  // create our Element struct for this element
86  Element* nc_elem = new Element(geom, elem->GetAttribute());
87  root_elements.Append(nc_elem);
88
89  for (int j = 0; j < GI[geom].nv; j++)
90  {
91  // root nodes are special: they have p1 == p2 == orig. mesh vertex id
92  Node* node = nodes.Get(v[j], v[j]);
93
94  if (!node->vertex)
95  {
96  // create a vertex in the node and initialize its position
97  const double* pos = mesh->GetVertex(v[j]);
98  node->vertex = new Vertex(pos[0], pos[1], pos[2]);
99  vertex_nodeId[v[j]] = node->id;
100  }
101
102  nc_elem->node[j] = node;
103  }
104
105  // increase reference count of all nodes the element is using
106  // (NOTE: this will also create and reference all edge and face nodes)
107  RefElementNodes(nc_elem);
108
109  // make links from faces back to the element
110  RegisterFaces(nc_elem);
111  }
112
113  // store boundary element attributes
114  for (int i = 0; i < mesh->GetNBE(); i++)
115  {
116  const mfem::Element *be = mesh->GetBdrElement(i);
117  const int *v = be->GetVertices();
118
119  Node* node[4];
120  for (int i = 0; i < be->GetNVertices(); i++)
121  {
122  node[i] = nodes.Peek(v[i], v[i]);
123  if (!node[i])
124  MFEM_ABORT("Boundary elements inconsistent.");
125  }
126
128  {
129  Face* face = faces.Peek(node[0], node[1], node[2], node[3]);
130  if (!face)
132
133  face->attribute = be->GetAttribute();
134  }
135  else if (be->GetType() == mfem::Element::SEGMENT)
136  {
137  Edge* edge = nodes.Peek(node[0], node[1])->edge;
138  if (!edge)
140
141  edge->attribute = be->GetAttribute();
142  }
143  else
144  mfem_error("NCMesh: only segment and quadrilateral boundary "
145  "elements are supported.");
146  }
147
149 }
150
152 {
153  for (int i = 0; i < root_elements.Size(); i++)
155 }
156
158 {
159  if (elem->ref_type)
160  {
161  for (int i = 0; i < 8; i++)
162  if (elem->child[i])
163  DeleteHierarchy(elem->child[i]);
164  }
165  else
166  {
167  UnrefElementNodes(elem);
168  }
169  delete elem;
170 }
171
172
174
176 {
177  MFEM_ASSERT(vertex, "NCMesh::Node::RefVertex: can't create vertex here.");
178  vertex->Ref();
179 }
180
182 {
183  if (!edge) edge = new Edge;
184  edge->Ref();
185 }
186
188 {
189  MFEM_ASSERT(vertex, "Cannot unref a nonexistent vertex.");
190  if (!vertex->Unref()) vertex = NULL;
191  if (!vertex && !edge) nodes.Delete(this);
192 }
193
195 {
197  MFEM_ASSERT(edge, "Cannot unref a nonexistent edge.");
198  if (!edge->Unref()) edge = NULL;
199  if (!vertex && !edge) nodes.Delete(this);
200 }
201
203 {
204  MFEM_ASSERT(!vertex && !edge, "Node was not unreffed properly.");
205  if (vertex) delete vertex;
206  if (edge) delete edge;
207 }
208
210 {
211  Node** node = elem->node;
212  GeomInfo& gi = GI[elem->geom];
213
214  // ref all vertices
215  for (int i = 0; i < gi.nv; i++)
216  node[i]->RefVertex();
217
218  // ref all edges (possibly creating them)
219  for (int i = 0; i < gi.ne; i++)
220  {
221  const int* ev = gi.edges[i];
222  nodes.Get(node[ev[0]], node[ev[1]])->RefEdge();
223  }
224
225  // ref all faces (possibly creating them)
226  for (int i = 0; i < gi.nf; i++)
227  {
228  const int* fv = gi.faces[i];
229  faces.Get(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]])->Ref();
230  // NOTE: face->RegisterElement called elsewhere to avoid having
231  // to store 3 element pointers temporarily in the face when refining
232  }
233 }
234
236 {
237  Node** node = elem->node;
238  GeomInfo& gi = GI[elem->geom];
239
240  // unref all faces (possibly destroying them)
241  for (int i = 0; i < gi.nf; i++)
242  {
243  const int* fv = gi.faces[i];
244  Face* face = faces.Peek(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
245  face->ForgetElement(elem);
246  if (!face->Unref()) faces.Delete(face);
247  }
248
249  // unref all edges (possibly destroying them)
250  for (int i = 0; i < gi.ne; i++)
251  {
252  const int* ev = gi.edges[i];
253  //nodes.Peek(node[ev[0]], node[ev[1]])->UnrefEdge(nodes); -- pre-aniso
254  PeekAltParents(node[ev[0]], node[ev[1]])->UnrefEdge(nodes);
255  }
256
257  // unref all vertices (possibly destroying them)
258  for (int i = 0; i < gi.nv; i++)
259  elem->node[i]->UnrefVertex(nodes);
260 }
261
263 {
264  if (elem[0] == NULL)
265  elem[0] = e;
266  else if (elem[1] == NULL)
267  elem[1] = e;
268  else
269  MFEM_ASSERT(0, "Can't have 3 elements in Face::elem[2].");
270 }
271
273 {
274  if (elem[0] == e)
275  elem[0] = NULL;
276  else if (elem[1] == e)
277  elem[1] = NULL;
278  else
280 }
281
283 {
284  Node** node = elem->node;
285  GeomInfo& gi = GI[elem->geom];
286
287  for (int i = 0; i < gi.nf; i++)
288  {
289  const int* fv = gi.faces[i];
290  Face* face = faces.Peek(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
291  face->RegisterElement(elem);
292  }
293 }
294
296 {
297  if (elem[0])
298  {
299  MFEM_ASSERT(!elem[1], "Not a single element face.");
300  return elem[0];
301  }
302  else
303  {
304  MFEM_ASSERT(elem[1], "No elements in face.");
305  return elem[1];
306  }
307 }
308
310 {
311  Node* mid = nodes.Peek(v1, v2);
312  if (!mid)
313  {
314  // In rare cases, a mid-face node exists under alternate parents w1, w2
315  // (see picture) instead of the requested parents v1, v2. This is an
316  // inconsistent situation that may exist temporarily as a result of
317  // "nodes.Reparent" while doing anisotropic splits, before forced
318  // refinements are all processed. This function attempts to retrieve such
319  // a node. An extra twist is that w1 and w2 may themselves need to be
320  // obtained using this very function.
321  //
322  // v1->p1 v1 v1->p2
323  // *------*------*
324  // | | |
325  // | |mid |
326  // w1 *------*------* w2
327  // | | |
328  // | | |
329  // *------*------*
330  // v2->p1 v2 v2->p2
331  //
332  // NOTE: this function would not be needed if the elements remembered
333  // pointers to their edge nodes. We have however opted to save memory
334  // at the cost of this computation, which is only necessary when forced
335  // refinements are being done.
336
337  if ((v1->p1 != v1->p2) && (v2->p1 != v2->p2)) // non-top-level nodes?
338  {
339  Node *v1p1 = nodes.Peek(v1->p1), *v1p2 = nodes.Peek(v1->p2);
340  Node *v2p1 = nodes.Peek(v2->p1), *v2p2 = nodes.Peek(v2->p2);
341
342  Node* w1 = PeekAltParents(v1p1, v2p1);
343  Node* w2 = w1 ? PeekAltParents(v1p2, v2p2) : NULL /* optimization */;
344
345  if (!w1 || !w2) // one more try may be needed as p1, p2 are unordered
346  w1 = PeekAltParents(v1p1, v2p2),
347  w2 = w1 ? PeekAltParents(v1p2, v2p1) : NULL /* optimization */;
348
349  if (w1 && w2) // got both alternate parents?
350  mid = nodes.Peek(w1, w2);
351  }
352  }
353  return mid;
354 }
355
356
358
359 NCMesh::Element::Element(int geom, int attr)
360  : geom(geom), attribute(attr), ref_type(0), index(-1)
361 {
362  memset(node, 0, sizeof(node));
363
364  // NOTE: in 2D the 8-element node/child arrays are not optimal, however,
365  // testing shows we would only save 17% of the total NCMesh memory if
366  // 4-element arrays were used (e.g. through templates); we thus prefer to
367  // keep the code as simple as possible.
368 }
369
372  Node* n4, Node* n5, Node* n6, Node* n7,
373  int attr,
374  int fattr0, int fattr1, int fattr2,
375  int fattr3, int fattr4, int fattr5)
376 {
377  // create new unrefined element, initialize nodes
378  Element* e = new Element(Geometry::CUBE, attr);
379  e->node[0] = n0, e->node[1] = n1, e->node[2] = n2, e->node[3] = n3;
380  e->node[4] = n4, e->node[5] = n5, e->node[6] = n6, e->node[7] = n7;
381
382  // get face nodes and assign face attributes
383  Face* f[6];
384  for (int i = 0; i < gi_hex.nf; i++)
385  {
386  const int* fv = gi_hex.faces[i];
387  f[i] = faces.Get(e->node[fv[0]], e->node[fv[1]],
388  e->node[fv[2]], e->node[fv[3]]);
389  }
390
391  f[0]->attribute = fattr0, f[1]->attribute = fattr1;
392  f[2]->attribute = fattr2, f[3]->attribute = fattr3;
393  f[4]->attribute = fattr4, f[5]->attribute = fattr5;
394
395  return e;
396 }
397
400  int attr,
401  int eattr0, int eattr1, int eattr2, int eattr3)
402 {
403  // create new unrefined element, initialize nodes
404  Element* e = new Element(Geometry::SQUARE, attr);
405  e->node[0] = n0, e->node[1] = n1, e->node[2] = n2, e->node[3] = n3;
406
407  // get edge nodes and assign edge attributes
408  Edge* edge[4];
409  for (int i = 0; i < gi_quad.ne; i++)
410  {
411  const int* ev = gi_quad.edges[i];
412  Node* node = nodes.Get(e->node[ev[0]], e->node[ev[1]]);
413  if (!node->edge) node->edge = new Edge;
414  edge[i] = node->edge;
415  }
416
417  edge[0]->attribute = eattr0;
418  edge[1]->attribute = eattr1;
419  edge[2]->attribute = eattr2;
420  edge[3]->attribute = eattr3;
421
422  return e;
423 }
424
427  int attr, int eattr0, int eattr1, int eattr2)
428 {
429  // create new unrefined element, initialize nodes
430  Element* e = new Element(Geometry::TRIANGLE, attr);
431  e->node[0] = n0, e->node[1] = n1, e->node[2] = n2;
432
433  // get edge nodes and assign edge attributes
434  Edge* edge[3];
435  for (int i = 0; i < gi_tri.ne; i++)
436  {
437  const int* ev = gi_tri.edges[i];
438  Node* node = nodes.Get(e->node[ev[0]], e->node[ev[1]]);
439  if (!node->edge) node->edge = new Edge;
440  edge[i] = node->edge;
441  }
442
443  edge[0]->attribute = eattr0;
444  edge[1]->attribute = eattr1;
445  edge[2]->attribute = eattr2;
446
447  return e;
448 }
449
451 {
452  MFEM_ASSERT(v1->vertex && v2->vertex,
453  "NCMesh::NewVertex: missing parent vertices.");
454
455  // get the midpoint between v1 and v2
456  Vertex* v = new Vertex;
457  for (int i = 0; i < 3; i++)
458  v->pos[i] = (v1->vertex->pos[i] + v2->vertex->pos[i]) * 0.5;
459
460  return v;
461 }
462
464 {
465  // in 3D we must be careful about getting the mid-edge node
466  Node* mid = PeekAltParents(v1, v2);
467  if (!mid) mid = nodes.Get(v1, v2);
468  if (!mid->vertex) mid->vertex = NewVertex(v1, v2);
469  return mid;
470 }
471
473 {
474  // simple version for 2D cases
475  Node* mid = nodes.Get(v1, v2);
476  if (!mid->vertex) mid->vertex = NewVertex(v1, v2);
477  return mid;
478 }
479
482 {
483  // mid-face node can be created either from (e1, e3) or from (e2, e4)
484  Node* midf = nodes.Peek(e1, e3);
485  if (midf)
486  {
487  if (!midf->vertex) midf->vertex = NewVertex(e1, e3);
488  return midf;
489  }
490  else
491  {
492  midf = nodes.Get(e2, e4);
493  if (!midf->vertex) midf->vertex = NewVertex(e2, e4);
494  return midf;
495  }
496 }
497
498 //
499 inline bool NCMesh::NodeSetX1(Node* node, Node** n)
500 { return node == n[0] || node == n[3] || node == n[4] || node == n[7]; }
501
502 inline bool NCMesh::NodeSetX2(Node* node, Node** n)
503 { return node == n[1] || node == n[2] || node == n[5] || node == n[6]; }
504
505 inline bool NCMesh::NodeSetY1(Node* node, Node** n)
506 { return node == n[0] || node == n[1] || node == n[4] || node == n[5]; }
507
508 inline bool NCMesh::NodeSetY2(Node* node, Node** n)
509 { return node == n[2] || node == n[3] || node == n[6] || node == n[7]; }
510
511 inline bool NCMesh::NodeSetZ1(Node* node, Node** n)
512 { return node == n[0] || node == n[1] || node == n[2] || node == n[3]; }
513
514 inline bool NCMesh::NodeSetZ2(Node* node, Node** n)
515 { return node == n[4] || node == n[5] || node == n[6] || node == n[7]; }
516
517
518 void NCMesh::ForceRefinement(Node* v1, Node* v2, Node* v3, Node* v4)
519 {
520  // get the element this face belongs to
521  Face* face = faces.Peek(v1, v2, v3, v4);
522  if (!face) return;
523
524  Element* elem = face->GetSingleElement();
526
527  Node** nodes = elem->node;
528
529  // schedule the right split depending on face orientation
530  if ((NodeSetX1(v1, nodes) && NodeSetX2(v2, nodes)) ||
531  (NodeSetX1(v2, nodes) && NodeSetX2(v1, nodes)))
532  {
533  ref_stack.Append(RefStackItem(elem, 1)); // X split
534  }
535  else if ((NodeSetY1(v1, nodes) && NodeSetY2(v2, nodes)) ||
536  (NodeSetY1(v2, nodes) && NodeSetY2(v1, nodes)))
537  {
538  ref_stack.Append(RefStackItem(elem, 2)); // Y split
539  }
540  else if ((NodeSetZ1(v1, nodes) && NodeSetZ2(v2, nodes)) ||
541  (NodeSetZ1(v2, nodes) && NodeSetZ2(v1, nodes)))
542  {
543  ref_stack.Append(RefStackItem(elem, 4)); // Z split
544  }
545  else
546  MFEM_ASSERT(0, "Inconsistent element/face structure.");
547 }
548
549
550 void NCMesh::CheckAnisoFace(Node* v1, Node* v2, Node* v3, Node* v4,
551  Node* mid12, Node* mid34, int level)
552 {
553  // When a face is getting split anisotropically (without loss of generality
554  // we assume a "vertical" split here, see picture), it is important to make
555  // sure that the mid-face vertex (midf) has mid34 and mid12 as parents.
556  // This is necessary for the face traversal algorithm and at places like
557  // Refine() that assume the mid-edge nodes to be accessible through the right
558  // parents. However, midf may already exist under the parents mid41 and
559  // mid23. In that case we need to "reparent" midf, i.e., reinsert it to the
560  // hash-table under the correct parents. This doesn't affect other nodes as
561  // all IDs stay the same, only the face refinement "tree" is affected.
562  //
563  // v4 mid34 v3
564  // *------*------*
565  // | | |
566  // | |midf |
567  // mid41 *- - - *- - - * mid23
568  // | | |
569  // | | |
570  // *------*------*
571  // v1 mid12 v2
572  //
573  // This function is recusive, because the above applies to any node along the
574  // middle vertical edge. The function calls itself again for the bottom and
575  // upper half of the above picture.
576
577  Node* mid23 = nodes.Peek(v2, v3);
578  Node* mid41 = nodes.Peek(v4, v1);
579  if (mid23 && mid41)
580  {
581  Node* midf = nodes.Peek(mid23, mid41);
582  if (midf)
583  {
584  nodes.Reparent(midf, mid12->id, mid34->id);
585
586  CheckAnisoFace(v1, v2, mid23, mid41, mid12, midf, level+1);
587  CheckAnisoFace(mid41, mid23, v3, v4, midf, mid34, level+1);
588  return;
589  }
590  }
591
592  /* Also, this is the place where forced refinements begin. In the picture,
593  the edges mid12-midf and midf-mid34 should actually exist in the
594  neighboring elements, otherwise the mesh is inconsistent and needs to be
595  fixed. */
596
597  if (level > 0)
598  ForceRefinement(v1, v2, v3, v4);
599 }
600
601 void NCMesh::CheckIsoFace(Node* v1, Node* v2, Node* v3, Node* v4,
602  Node* e1, Node* e2, Node* e3, Node* e4, Node* midf)
603 {
604  /* If anisotropic refinements are present in the mesh, we need to check
605  isotropically split faces as well. The iso face can be thought to contain
606  four anisotropic cases as in the function CheckAnisoFace, that still need
607  to be checked for the correct parents. */
608
609  CheckAnisoFace(v1, v2, e2, e4, e1, midf);
610  CheckAnisoFace(e4, e2, v3, v4, midf, e3);
611  CheckAnisoFace(v4, v1, e1, e3, e4, midf);
612  CheckAnisoFace(e3, e1, v2, v3, midf, e2);
613 }
614
615
616 void NCMesh::Refine(Element* elem, int ref_type)
617 {
618  if (!ref_type) return;
619
620  // handle elements that may have been (force-) refined already
621  if (elem->ref_type)
622  {
623  int remaining = ref_type & ~elem->ref_type;
624
625  // do the remaining splits on the children
626  for (int i = 0; i < 8; i++)
627  if (elem->child[i])
628  Refine(elem->child[i], remaining);
629
630  return;
631  }
632
633  Node** no = elem->node;
634  int attr = elem->attribute;
635
636  Element* child[8];
637  memset(child, 0, sizeof(child));
638
639  // create child elements
640  if (elem->geom == Geometry::CUBE)
641  {
642  // get parent's face attributes
643  int fa[6];
644  for (int i = 0; i < gi_hex.nf; i++)
645  {
646  const int* fv = gi_hex.faces[i];
647  Face* face = faces.Peek(no[fv[0]], no[fv[1]], no[fv[2]], no[fv[3]]);
648  fa[i] = face->attribute;
649  }
650
651  // Vertex numbering is assumed to be as follows:
652  //
653  // 7 6
654  // +------------+ Faces: 0 bottom
655  // /| /| 1 front
656  // 4 / | 5 / | 2 right
657  // +------------+ | 3 back
658  // | | | | 4 left
659  // | +---------|--+ 5 top
660  // | / 3 | / 2 Z Y
661  // |/ |/ |/
662  // +------------+ *--X
663  // 0 1
664
665  if (ref_type == 1) // split along X axis
666  {
667  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
668  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
669  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
670  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
671
672  child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
673  no[4], mid45, mid67, no[7], attr,
674  fa[0], fa[1], -1, fa[3], fa[4], fa[5]);
675
676  child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
677  mid45, no[5], no[6], mid67, attr,
678  fa[0], fa[1], fa[2], fa[3], -1, fa[5]);
679
680  CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
681  CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
682  CheckAnisoFace(no[4], no[5], no[6], no[7], mid45, mid67);
683  CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
684  }
685  else if (ref_type == 2) // split along Y axis
686  {
687  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
688  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
689  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
690  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
691
692  child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
693  no[4], no[5], mid56, mid74, attr,
694  fa[0], fa[1], fa[2], -1, fa[4], fa[5]);
695
696  child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
697  mid74, mid56, no[6], no[7], attr,
698  fa[0], -1, fa[2], fa[3], fa[4], fa[5]);
699
700  CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
701  CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
702  CheckAnisoFace(no[5], no[6], no[7], no[4], mid56, mid74);
703  CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
704  }
705  else if (ref_type == 4) // split along Z axis
706  {
707  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
708  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
709  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
710  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
711
712  child[0] = NewHexahedron(no[0], no[1], no[2], no[3],
713  mid04, mid15, mid26, mid37, attr,
714  fa[0], fa[1], fa[2], fa[3], fa[4], -1);
715
716  child[1] = NewHexahedron(mid04, mid15, mid26, mid37,
717  no[4], no[5], no[6], no[7], attr,
718  -1, fa[1], fa[2], fa[3], fa[4], fa[5]);
719
720  CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
721  CheckAnisoFace(no[5], no[1], no[2], no[6], mid15, mid26);
722  CheckAnisoFace(no[6], no[2], no[3], no[7], mid26, mid37);
723  CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
724  }
725  else if (ref_type == 3) // XY split
726  {
727  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
728  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
729  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
730  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
731
732  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
733  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
734  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
735  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
736
737  Node* midf0 = GetMidFaceVertex(mid23, mid12, mid01, mid30);
738  Node* midf5 = GetMidFaceVertex(mid45, mid56, mid67, mid74);
739
740  child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
741  no[4], mid45, midf5, mid74, attr,
742  fa[0], fa[1], -1, -1, fa[4], fa[5]);
743
744  child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
745  mid45, no[5], mid56, midf5, attr,
746  fa[0], fa[1], fa[2], -1, -1, fa[5]);
747
748  child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
749  midf5, mid56, no[6], mid67, attr,
750  fa[0], -1, fa[2], fa[3], -1, fa[5]);
751
752  child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
753  mid74, midf5, mid67, no[7], attr,
754  fa[0], -1, -1, fa[3], fa[4], fa[5]);
755
756  CheckAnisoFace(no[0], no[1], no[5], no[4], mid01, mid45);
757  CheckAnisoFace(no[1], no[2], no[6], no[5], mid12, mid56);
758  CheckAnisoFace(no[2], no[3], no[7], no[6], mid23, mid67);
759  CheckAnisoFace(no[3], no[0], no[4], no[7], mid30, mid74);
760
761  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
762  CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
763  }
764  else if (ref_type == 5) // XZ split
765  {
766  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
767  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
768  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
769  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
770
771  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
772  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
773  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
774  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
775
776  Node* midf1 = GetMidFaceVertex(mid01, mid15, mid45, mid04);
777  Node* midf3 = GetMidFaceVertex(mid23, mid37, mid67, mid26);
778
779  child[0] = NewHexahedron(no[0], mid01, mid23, no[3],
780  mid04, midf1, midf3, mid37, attr,
781  fa[0], fa[1], -1, fa[3], fa[4], -1);
782
783  child[1] = NewHexahedron(mid01, no[1], no[2], mid23,
784  midf1, mid15, mid26, midf3, attr,
785  fa[0], fa[1], fa[2], fa[3], -1, -1);
786
787  child[2] = NewHexahedron(midf1, mid15, mid26, midf3,
788  mid45, no[5], no[6], mid67, attr,
789  -1, fa[1], fa[2], fa[3], -1, fa[5]);
790
791  child[3] = NewHexahedron(mid04, midf1, midf3, mid37,
792  no[4], mid45, mid67, no[7], attr,
793  -1, fa[1], -1, fa[3], fa[4], fa[5]);
794
795  CheckAnisoFace(no[3], no[2], no[1], no[0], mid23, mid01);
796  CheckAnisoFace(no[2], no[6], no[5], no[1], mid26, mid15);
797  CheckAnisoFace(no[6], no[7], no[4], no[5], mid67, mid45);
798  CheckAnisoFace(no[7], no[3], no[0], no[4], mid37, mid04);
799
800  CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
801  CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
802  }
803  else if (ref_type == 6) // YZ split
804  {
805  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
806  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
807  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
808  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
809
810  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
811  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
812  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
813  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
814
815  Node* midf2 = GetMidFaceVertex(mid12, mid26, mid56, mid15);
816  Node* midf4 = GetMidFaceVertex(mid30, mid04, mid74, mid37);
817
818  child[0] = NewHexahedron(no[0], no[1], mid12, mid30,
819  mid04, mid15, midf2, midf4, attr,
820  fa[0], fa[1], fa[2], -1, fa[4], -1);
821
822  child[1] = NewHexahedron(mid30, mid12, no[2], no[3],
823  midf4, midf2, mid26, mid37, attr,
824  fa[0], -1, fa[2], fa[3], fa[4], -1);
825
826  child[2] = NewHexahedron(mid04, mid15, midf2, midf4,
827  no[4], no[5], mid56, mid74, attr,
828  -1, fa[1], fa[2], -1, fa[4], fa[5]);
829
830  child[3] = NewHexahedron(midf4, midf2, mid26, mid37,
831  mid74, mid56, no[6], no[7], attr,
832  -1, -1, fa[2], fa[3], fa[4], fa[5]);
833
834  CheckAnisoFace(no[4], no[0], no[1], no[5], mid04, mid15);
835  CheckAnisoFace(no[0], no[3], no[2], no[1], mid30, mid12);
836  CheckAnisoFace(no[3], no[7], no[6], no[2], mid37, mid26);
837  CheckAnisoFace(no[7], no[4], no[5], no[6], mid74, mid56);
838
839  CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
840  CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
841  }
842  else if (ref_type == 7) // full isotropic refinement
843  {
844  Node* mid01 = GetMidEdgeVertex(no[0], no[1]);
845  Node* mid12 = GetMidEdgeVertex(no[1], no[2]);
846  Node* mid23 = GetMidEdgeVertex(no[2], no[3]);
847  Node* mid30 = GetMidEdgeVertex(no[3], no[0]);
848
849  Node* mid45 = GetMidEdgeVertex(no[4], no[5]);
850  Node* mid56 = GetMidEdgeVertex(no[5], no[6]);
851  Node* mid67 = GetMidEdgeVertex(no[6], no[7]);
852  Node* mid74 = GetMidEdgeVertex(no[7], no[4]);
853
854  Node* mid04 = GetMidEdgeVertex(no[0], no[4]);
855  Node* mid15 = GetMidEdgeVertex(no[1], no[5]);
856  Node* mid26 = GetMidEdgeVertex(no[2], no[6]);
857  Node* mid37 = GetMidEdgeVertex(no[3], no[7]);
858
859  Node* midf0 = GetMidFaceVertex(mid23, mid12, mid01, mid30);
860  Node* midf1 = GetMidFaceVertex(mid01, mid15, mid45, mid04);
861  Node* midf2 = GetMidFaceVertex(mid12, mid26, mid56, mid15);
862  Node* midf3 = GetMidFaceVertex(mid23, mid37, mid67, mid26);
863  Node* midf4 = GetMidFaceVertex(mid30, mid04, mid74, mid37);
864  Node* midf5 = GetMidFaceVertex(mid45, mid56, mid67, mid74);
865
866  Node* midel = GetMidEdgeVertex(midf1, midf3);
867
868  child[0] = NewHexahedron(no[0], mid01, midf0, mid30,
869  mid04, midf1, midel, midf4, attr,
870  fa[0], fa[1], -1, -1, fa[4], -1);
871
872  child[1] = NewHexahedron(mid01, no[1], mid12, midf0,
873  midf1, mid15, midf2, midel, attr,
874  fa[0], fa[1], fa[2], -1, -1, -1);
875
876  child[2] = NewHexahedron(midf0, mid12, no[2], mid23,
877  midel, midf2, mid26, midf3, attr,
878  fa[0], -1, fa[2], fa[3], -1, -1);
879
880  child[3] = NewHexahedron(mid30, midf0, mid23, no[3],
881  midf4, midel, midf3, mid37, attr,
882  fa[0], -1, -1, fa[3], fa[4], -1);
883
884  child[4] = NewHexahedron(mid04, midf1, midel, midf4,
885  no[4], mid45, midf5, mid74, attr,
886  -1, fa[1], -1, -1, fa[4], fa[5]);
887
888  child[5] = NewHexahedron(midf1, mid15, midf2, midel,
889  mid45, no[5], mid56, midf5, attr,
890  -1, fa[1], fa[2], -1, -1, fa[5]);
891
892  child[6] = NewHexahedron(midel, midf2, mid26, midf3,
893  midf5, mid56, no[6], mid67, attr,
894  -1, -1, fa[2], fa[3], -1, fa[5]);
895
896  child[7] = NewHexahedron(midf4, midel, midf3, mid37,
897  mid74, midf5, mid67, no[7], attr,
898  -1, -1, -1, fa[3], fa[4], fa[5]);
899
900  CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
901  CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
902  CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
903  CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
904  CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
905  CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
906  }
907  else
908  MFEM_ABORT("Invalid refinement type.");
909  }
910  else if (elem->geom == Geometry::SQUARE)
911  {
912  // get parent's edge attributes
913  int ea0 = nodes.Peek(no[0], no[1])->edge->attribute;
914  int ea1 = nodes.Peek(no[1], no[2])->edge->attribute;
915  int ea2 = nodes.Peek(no[2], no[3])->edge->attribute;
916  int ea3 = nodes.Peek(no[3], no[0])->edge->attribute;
917
918  ref_type &= ~4; // ignore Z bit
919
920  if (ref_type == 1) // X split
921  {
922  Node* mid01 = GetMidEdgeVertexSimple(no[0], no[1]);
923  Node* mid23 = GetMidEdgeVertexSimple(no[2], no[3]);
924
925  child[0] = NewQuadrilateral(no[0], mid01, mid23, no[3],
926  attr, ea0, -1, ea2, ea3);
927
928  child[1] = NewQuadrilateral(mid01, no[1], no[2], mid23,
929  attr, ea0, ea1, ea2, -1);
930  }
931  else if (ref_type == 2) // Y split
932  {
933  Node* mid12 = GetMidEdgeVertexSimple(no[1], no[2]);
934  Node* mid30 = GetMidEdgeVertexSimple(no[3], no[0]);
935
936  child[0] = NewQuadrilateral(no[0], no[1], mid12, mid30,
937  attr, ea0, ea1, -1, ea3);
938
939  child[1] = NewQuadrilateral(mid30, mid12, no[2], no[3],
940  attr, -1, ea1, ea2, ea3);
941  }
942  else if (ref_type == 3) // iso split
943  {
944  Node* mid01 = GetMidEdgeVertexSimple(no[0], no[1]);
945  Node* mid12 = GetMidEdgeVertexSimple(no[1], no[2]);
946  Node* mid23 = GetMidEdgeVertexSimple(no[2], no[3]);
947  Node* mid30 = GetMidEdgeVertexSimple(no[3], no[0]);
948
949  Node* midel = GetMidEdgeVertexSimple(mid01, mid23);
950
951  child[0] = NewQuadrilateral(no[0], mid01, midel, mid30,
952  attr, ea0, -1, -1, ea3);
953
954  child[1] = NewQuadrilateral(mid01, no[1], mid12, midel,
955  attr, ea0, ea1, -1, -1);
956
957  child[2] = NewQuadrilateral(midel, mid12, no[2], mid23,
958  attr, -1, ea1, ea2, -1);
959
960  child[3] = NewQuadrilateral(mid30, midel, mid23, no[3],
961  attr, -1, -1, ea2, ea3);
962  }
963  else
964  MFEM_ABORT("Invalid refinement type.");
965  }
966  else if (elem->geom == Geometry::TRIANGLE)
967  {
968  // get parent's edge attributes
969  int ea0 = nodes.Peek(no[0], no[1])->edge->attribute;
970  int ea1 = nodes.Peek(no[1], no[2])->edge->attribute;
971  int ea2 = nodes.Peek(no[2], no[0])->edge->attribute;
972
973  // isotropic split - the only ref_type available for triangles
974  Node* mid01 = GetMidEdgeVertexSimple(no[0], no[1]);
975  Node* mid12 = GetMidEdgeVertexSimple(no[1], no[2]);
976  Node* mid20 = GetMidEdgeVertexSimple(no[2], no[0]);
977
978  child[0] = NewTriangle(no[0], mid01, mid20, attr, ea0, -1, ea2);
979  child[1] = NewTriangle(mid01, no[1], mid12, attr, ea0, ea1, -1);
980  child[2] = NewTriangle(mid20, mid12, no[2], attr, -1, ea1, ea2);
981  child[3] = NewTriangle(mid01, mid12, mid20, attr, -1, -1, -1);
982  }
983  else
984  MFEM_ABORT("Unsupported element geometry.");
985
986  // start using the nodes of the children, create edges & faces
987  for (int i = 0; i < 8; i++)
988  if (child[i])
989  RefElementNodes(child[i]);
990
991  // sign off of all nodes of the parent, clean up unused nodes
992  UnrefElementNodes(elem);
993
994  // register the children in their faces once the parent is out of the way
995  for (int i = 0; i < 8; i++)
996  if (child[i])
997  RegisterFaces(child[i]);
998
999  // finish the refinement
1000  elem->ref_type = ref_type;
1001  memcpy(elem->child, child, sizeof(elem->child));
1002 }
1003
1004
1005 void NCMesh::Refine(const Array<Refinement>& refinements)
1006 {
1007  // push all refinements on the stack in reverse order
1008  for (int i = refinements.Size()-1; i >= 0; i--)
1009  {
1010  const Refinement& ref = refinements[i];
1011  ref_stack.Append(RefStackItem(leaf_elements[ref.index], ref.ref_type));
1012  }
1013
1014  // keep refining as long as the stack contains something
1015  int nforced = 0;
1016  while (ref_stack.Size())
1017  {
1018  RefStackItem ref = ref_stack.Last();
1019  ref_stack.DeleteLast();
1020
1021  int size = ref_stack.Size();
1022  Refine(ref.elem, ref.ref_type);
1023  nforced += ref_stack.Size() - size;
1024  }
1025
1026  /* TODO: the current algorithm of forced refinements is not optimal. As
1027  forced refinements spread through the mesh, some may not be necessary
1028  in the end, since the affected elements may still be scheduled for
1029  refinement that could stop the propagation. We should introduce the
1030  member Element::ref_pending that would show the intended refinement in
1031  the batch. A forced refinement would be combined with ref_pending to
1032  (possibly) stop the propagation earlier. */
1033
1034 #ifdef MFEM_DEBUG
1035  std::cout << "Refined " << refinements.Size() << " + " << nforced
1036  << " elements" << std::endl;
1037 #endif
1038
1040  UpdateVertices();
1041 }
1042
1043
1044 /*void NCMesh::Derefine(Element* elem)
1045  {
1046
1047  }*/
1048
1049
1051 {
1052  int num_vert = 0;
1053  for (HashTable<Node>::Iterator it(nodes); it; ++it)
1054  if (it->vertex)
1055  it->vertex->index = num_vert++;
1056
1057  vertex_nodeId.SetSize(num_vert);
1058
1059  num_vert = 0;
1060  for (HashTable<Node>::Iterator it(nodes); it; ++it)
1061  if (it->vertex)
1062  vertex_nodeId[num_vert++] = it->id;
1063 }
1064
1065
1067
1069 {
1070  if (!e->ref_type)
1071  {
1072  e->index = leaf_elements.Size();
1073  leaf_elements.Append(e);
1074  }
1075  else
1076  {
1077  e->index = -1;
1078  for (int i = 0; i < 8; i++)
1079  if (e->child[i])
1080  GetLeafElements(e->child[i]);
1081  }
1082 }
1083
1085 {
1086  // collect leaf elements
1087  leaf_elements.SetSize(0);
1088  for (int i = 0; i < root_elements.Size(); i++)
1090 }
1091
1093  Array<mfem::Element*>& elements,
1094  Array<mfem::Element*>& boundary)
1095 {
1096  // copy vertex coordinates
1097  vertices.SetSize(vertex_nodeId.Size());
1098  int i = 0;
1099  for (HashTable<Node>::Iterator it(nodes); it; ++it)
1100  if (it->vertex)
1101  vertices[i++].SetCoords(it->vertex->pos);
1102
1103  elements.SetSize(leaf_elements.Size());
1104  boundary.SetSize(0);
1105
1106  for (int i = 0; i < leaf_elements.Size(); i++)
1107  {
1108  Element* nc_elem = leaf_elements[i];
1109  Node** node = nc_elem->node;
1110  GeomInfo& gi = GI[nc_elem->geom];
1111
1112  // create an mfem::Element for each leaf Element
1113  mfem::Element* elem = NULL;
1114  switch (nc_elem->geom)
1115  {
1116  case Geometry::CUBE: elem = new Hexahedron; break;
1117  case Geometry::SQUARE: elem = new Quadrilateral; break;
1118  case Geometry::TRIANGLE: elem = new Triangle; break;
1119  }
1120
1121  elements[i] = elem;
1122  elem->SetAttribute(nc_elem->attribute);
1123  for (int j = 0; j < gi.nv; j++)
1124  elem->GetVertices()[j] = node[j]->vertex->index;
1125
1126  // create boundary elements
1127  if (nc_elem->geom == Geometry::CUBE)
1128  {
1129  for (int k = 0; k < gi.nf; k++)
1130  {
1131  const int* fv = gi.faces[k];
1132  Face* face = faces.Peek(node[fv[0]], node[fv[1]],
1133  node[fv[2]], node[fv[3]]);
1134  if (face->Boundary())
1135  {
1138  for (int j = 0; j < 4; j++)
1140
1142  }
1143  }
1144  }
1145  else // quad & triangle boundary elements
1146  {
1147  for (int k = 0; k < gi.ne; k++)
1148  {
1149  const int* ev = gi.edges[k];
1150  Edge* edge = nodes.Peek(node[ev[0]], node[ev[1]])->edge;
1151  if (edge->Boundary())
1152  {
1153  Segment* segment = new Segment;
1154  segment->SetAttribute(edge->attribute);
1155  for (int j = 0; j < 2; j++)
1156  segment->GetVertices()[j] = node[ev[j]]->vertex->index;
1157
1158  boundary.Append(segment);
1159  }
1160  }
1161  }
1162  }
1163 }
1164
1165
1167 {
1168  Table *edge_vertex = mesh->GetEdgeVertexTable();
1169
1170  for (int i = 0; i < edge_vertex->Size(); i++)
1171  {
1172  const int *ev = edge_vertex->GetRow(i);
1173  Node* node = nodes.Peek(vertex_nodeId[ev[0]], vertex_nodeId[ev[1]]);
1174
1176  node->edge->index = i;
1177  }
1178 }
1179
1181 {
1182  for (int i = 0; i < mesh->GetNFaces(); i++)
1183  {
1184  const int* fv = mesh->GetFace(i)->GetVertices();
1185  Face* face = faces.Peek(vertex_nodeId[fv[0]], vertex_nodeId[fv[1]],
1186  vertex_nodeId[fv[2]], vertex_nodeId[fv[3]]);
1187
1189  face->index = i;
1190  }
1191 }
1192
1193
1195
1196 int NCMesh::FaceSplitType(Node* v1, Node* v2, Node* v3, Node* v4,
1197  Node* mid[4])
1198 {
1199  // find edge nodes
1200  Node* e1 = nodes.Peek(v1, v2);
1201  Node* e2 = nodes.Peek(v2, v3);
1202  Node* e3 = e1 ? nodes.Peek(v3, v4) : NULL;
1203  Node* e4 = e2 ? nodes.Peek(v4, v1) : NULL;
1204
1205  // optional: return the mid-edge nodes if requested
1206  if (mid) mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4;
1207
1208  // try to get a mid-face node, either by (e1, e3) or by (e2, e4)
1209  Node *midf1 = NULL, *midf2 = NULL;
1210  if (e1 && e3) midf1 = nodes.Peek(e1, e3);
1211  if (e2 && e4) midf2 = nodes.Peek(e2, e4);
1212
1213  // only one way to access the mid-face node must always exist
1214  MFEM_ASSERT(!(midf1 && midf2), "Incorrectly split face!");
1215
1216  if (!midf1 && !midf2)
1217  return 0; // face not split
1218
1219  if (midf1)
1220  return 1; // face split "vertically"
1221  else
1222  return 2; // face split "horizontally"
1223 }
1224
1225 int NCMesh::find_node(Element* elem, Node* node)
1226 {
1227  for (int i = 0; i < 8; i++)
1228  if (elem->node[i] == node)
1229  return i;
1230
1232  return -1;
1233 }
1234
1235 static int find_hex_face(int a, int b, int c)
1236 {
1237  for (int i = 0; i < 6; i++)
1238  {
1239  const int* fv = gi_hex.faces[i];
1240  if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
1241  (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
1242  (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
1243  {
1244  return i;
1245  }
1246  }
1248  return -1;
1249 }
1250
1252  Element* elem, DenseMatrix& pm)
1253 {
1254  int master[4] = {
1255  find_node(elem, v0), find_node(elem, v1),
1256  find_node(elem, v2), find_node(elem, v3)
1257  };
1258
1259  int fi = find_hex_face(master[0], master[1], master[2]);
1260  const int* fv = gi_hex.faces[fi];
1261
1262  DenseMatrix tmp(pm);
1263  for (int i = 0, j; i < 4; i++)
1264  {
1265  for (j = 0; j < 4; j++)
1266  if (fv[i] == master[j])
1267  {
1268  // pm.column(i) = tmp.column(j)
1269  for (int k = 0; k < pm.Height(); k++)
1270  pm(k,i) = tmp(k,j);
1271  break;
1272  }
1273
1275  }
1276 }
1277
1278 inline int decode_dof(int dof, double& sign)
1279 {
1280  if (dof >= 0)
1281  return (sign = 1.0, dof);
1282  else
1283  return (sign = -1.0, -1 - dof);
1284 }
1285
1286 void NCMesh::AddDependencies(Array<int>& master_dofs, Array<int>& slave_dofs,
1287  DenseMatrix& I)
1288 {
1289  // make each slave DOF dependent on all master DOFs
1290  for (int i = 0; i < slave_dofs.Size(); i++)
1291  {
1292  double ssign;
1293  int sdof = decode_dof(slave_dofs[i], ssign);
1294  DofData& sdata = dof_data[sdof];
1295
1296  if (!sdata.dep_list.Size()) // not processed yet?
1297  {
1298  for (int j = 0; j < master_dofs.Size(); j++)
1299  {
1300  double coef = I(i, j);
1301  if (std::abs(coef) > 1e-12)
1302  {
1303  double msign;
1304  int mdof = decode_dof(master_dofs[j], msign);
1305  if (mdof != sdof)
1306  sdata.dep_list.Append(Dependency(mdof, coef * ssign * msign));
1307  }
1308  }
1309  }
1310  }
1311 }
1312
1313 void NCMesh::ConstrainEdge(Node* v0, Node* v1, double t0, double t1,
1314  Array<int>& master_dofs, int level)
1315 {
1316  Node* mid = nodes.Peek(v0, v1);
1317  if (!mid) return;
1318
1319  if (mid->edge && level > 0)
1320  {
1321  // we need to make this edge constrained; get its DOFs
1322  Array<int> slave_dofs;
1323  space->GetEdgeDofs(mid->edge->index, slave_dofs);
1324
1325  if (slave_dofs.Size() > 0)
1326  {
1327  // prepare edge transformation
1329  edge_T.SetFE(&SegmentFE);
1330
1331  DenseMatrix& pm = edge_T.GetPointMat();
1332  pm.SetSize(1, 2);
1333  pm(0,0) = t0, pm(0,1) = t1;
1334
1335  // handle slave edge orientation
1336  if (v0->vertex->index > v1->vertex->index)
1337  std::swap(pm(0,0), pm(0,1));
1338
1339  // obtain the local interpolation matrix
1340  const FiniteElement* edge_fe =
1342
1343  DenseMatrix I(edge_fe->GetDof());
1344  edge_fe->GetLocalInterpolation(edge_T, I);
1345
1346  // make each slave DOF dependent on all master edge DOFs
1348  }
1349  }
1350
1351  // recurse deeper
1352  double tmid = (t0 + t1) / 2;
1353  ConstrainEdge(v0, mid, t0, tmid, master_dofs, level+1);
1354  ConstrainEdge(mid, v1, tmid, t1, master_dofs, level+1);
1355 }
1356
1357 void NCMesh::ConstrainFace(Node* v0, Node* v1, Node* v2, Node* v3,
1358  const PointMatrix& pm,
1359  Array<int>& master_dofs, int level)
1360 {
1361  if (level > 0)
1362  {
1363  // check if we made it to a face that is not split further
1364  Face* face = faces.Peek(v0, v1, v2, v3);
1365  if (face)
1366  {
1367  // yes, we need to make this face constrained; get its DOFs
1368  Array<int> slave_dofs;
1369  space->GetFaceDofs(face->index, slave_dofs);
1370
1371  if (slave_dofs.Size() > 0)
1372  {
1373  // prepare face transformation
1376  pm.GetMatrix(face_T.GetPointMat());
1377
1378  // reorder face_T point matrix according to slave face orientation
1379  ReorderFacePointMat(v0, v1, v2, v3, face->GetSingleElement(),
1380  face_T.GetPointMat());
1381
1382  // obtain the local interpolation matrix
1383  const FiniteElement* face_fe =
1385
1386  DenseMatrix I(face_fe->GetDof());
1387  face_fe->GetLocalInterpolation(face_T, I);
1388
1389  // make each slave DOF dependent on all master face DOFs
1391  }
1392  return;
1393  }
1394  }
1395
1396  // we need to recurse deeper
1397  Node* mid[4];
1398  int split = FaceSplitType(v0, v1, v2, v3, mid);
1399
1400  if (split == 1) // "X" split face
1401  {
1402  Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
1403
1404  ConstrainFace(v0, mid[0], mid[2], v3,
1405  PointMatrix(pm(0), mid0, mid2, pm(3)),
1406  master_dofs, level+1);
1407
1408  ConstrainFace(mid[0], v1, v2, mid[2],
1409  PointMatrix(mid0, pm(1), pm(2), mid2),
1410  master_dofs, level+1);
1411  }
1412  else if (split == 2) // "Y" split face
1413  {
1414  Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
1415
1416  ConstrainFace(v0, v1, mid[1], mid[3],
1417  PointMatrix(pm(0), pm(1), mid1, mid3),
1418  master_dofs, level+1);
1419
1420  ConstrainFace(mid[3], mid[1], v2, v3,
1421  PointMatrix(mid3, mid1, pm(2), pm(3)),
1422  master_dofs, level+1);
1423  }
1424 }
1425
1426 void NCMesh::ProcessMasterEdge(Node* node[2], Node* edge)
1427 {
1428  // get a list of DOFs on the master edge
1429  Array<int> master_dofs;
1430  space->GetEdgeDofs(edge->edge->index, master_dofs);
1431
1432  if (master_dofs.Size() > 0)
1433  {
1434  // we'll keep track of our position within the master edge;
1435  // the initial transformation is identity (interval 0..1)
1436  double t0 = 0.0, t1 = 1.0;
1437  if (node[0]->vertex->index > node[1]->vertex->index)
1438  std::swap(t0, t1);
1439
1440  ConstrainEdge(node[0], node[1], t0, t1, master_dofs, 0);
1441  }
1442 }
1443
1444 void NCMesh::ProcessMasterFace(Node* node[4], Face* face)
1445 {
1446  // get a list of DOFs on the master face
1447  Array<int> master_dofs;
1448  space->GetFaceDofs(face->index, master_dofs);
1449
1450  if (master_dofs.Size() > 0)
1451  {
1452  // we'll keep track of our position within the master face;
1453  // the initial transformation is identity (vertices of the unit square)
1454  PointMatrix pm(Point(0,0), Point(1,0), Point(1,1), Point(0,1));
1455
1456  ConstrainFace(node[0], node[1], node[2], node[3], pm, master_dofs, 0);
1457  }
1458 }
1459
1461 {
1462  // are all constraining DOFs finalized?
1463  for (int i = 0; i < dd.dep_list.Size(); i++)
1464  if (!dof_data[dd.dep_list[i].dof].finalized)
1465  return false;
1466
1467  return true;
1468 }
1469
1471  SparseMatrix **cR_ptr)
1472 {
1473  // allocate temporary data for each DOF
1474  int n_dofs = space->GetNDofs();
1475  dof_data = new DofData[n_dofs];
1476
1477  this->space = space;
1478
1479  // visit edges and faces of leaf elements
1480  for (int i = 0; i < leaf_elements.Size(); i++)
1481  {
1482  Element* elem = leaf_elements[i];
1483  MFEM_ASSERT(!elem->ref_type, "Not a leaf element.");
1484  GeomInfo& gi = GI[elem->geom];
1485
1486  // visit edges of 'elem'
1487  for (int j = 0; j < gi.ne; j++)
1488  {
1489  const int* ev = gi.edges[j];
1490  Node* node[2] = { elem->node[ev[0]], elem->node[ev[1]] };
1491
1492  Node* edge = nodes.Peek(node[0], node[1]);
1494
1495  // this edge could contain slave edges that need constraining; traverse
1496  // them recursively and make them dependent on this master edge
1497  ProcessMasterEdge(node, edge);
1498  }
1499
1500  // visit faces of 'elem'
1501  for (int j = 0; j < gi.nf; j++)
1502  {
1503  Node* node[4];
1504  const int* fv = gi.faces[j];
1505  for (int k = 0; k < 4; k++)
1506  node[k] = elem->node[fv[k]];
1507
1508  Face* face = faces.Peek(node[0], node[1], node[2], node[3]);
1510
1511  if (face->ref_count == 1 && !face->Boundary())
1512  {
1513  // this is a potential master face that could be constraining
1514  // smaller faces adjacent to it; traverse them recursively and
1515  // make them dependent on this master face
1516  ProcessMasterFace(node, face);
1517  }
1518  }
1519  }
1520
1521  // DOFs that stayed independent are true DOFs
1522  int n_true_dofs = 0;
1523  for (int i = 0; i < n_dofs; i++)
1524  {
1525  DofData& dd = dof_data[i];
1526  if (dd.Independent())
1527  n_true_dofs++;
1528  }
1529
1530  // if all dofs are true dofs return empty (NULL) matrices
1531  if (n_true_dofs == n_dofs)
1532  {
1533  delete [] dof_data;
1534  if (cR_ptr)
1535  *cR_ptr = NULL;
1536  return NULL;
1537  }
1538
1539  // create the conforming restiction matrix
1540  int *cR_J = NULL;
1541  SparseMatrix *cR = NULL;
1542  if (cR_ptr)
1543  {
1544  int *cR_I = new int[n_true_dofs+1];
1545  double *cR_A = new double[n_true_dofs];
1546  cR_J = new int[n_true_dofs];
1547  for (int i = 0; i < n_true_dofs; i++)
1548  {
1549  cR_I[i] = i;
1550  cR_A[i] = 1.0;
1551  }
1552  cR_I[n_true_dofs] = n_true_dofs;
1553  cR = new SparseMatrix(cR_I, cR_J, cR_A, n_true_dofs, n_dofs);
1554  }
1555
1556  // create the conforming prolongation matrix
1557  SparseMatrix* cP = new SparseMatrix(n_dofs, n_true_dofs);
1558
1559  // put identity in the restiction and prolongation matrices for true DOFs
1560  for (int i = 0, true_dof = 0; i < n_dofs; i++)
1561  {
1562  DofData& dd = dof_data[i];
1563  if (dd.Independent())
1564  {
1565  if (cR_ptr)
1566  cR_J[true_dof] = i;
1568  dd.finalized = true;
1569  }
1570  }
1571
1572  // resolve dependencies of slave DOFs
1573  bool finished;
1574  int n_finalized = n_true_dofs;
1575  Array<int> cols;
1576  Vector srow;
1577  do
1578  {
1579  finished = true;
1580  for (int i = 0; i < n_dofs; i++)
1581  {
1582  DofData& dd = dof_data[i];
1583  if (!dd.finalized && DofFinalizable(dd))
1584  {
1585  for (int j = 0; j < dd.dep_list.Size(); j++)
1586  {
1587  Dependency& dep = dd.dep_list[j];
1588
1589  cP->GetRow(dep.dof, cols, srow);
1590  srow *= dep.coef;
1592  }
1593
1594  dd.finalized = true;
1595  n_finalized++;
1596  finished = false;
1597  }
1598  }
1599  }
1600  while (!finished);
1601
1602  // if everything is consistent (mesh, face orientations, etc.), we should
1603  // be able to finalize all slave DOFs, otherwise it's a serious error
1604  if (n_finalized != n_dofs)
1605  MFEM_ABORT("Error creating cP matrix.");
1606
1607  delete [] dof_data;
1608
1609  if (cR_ptr)
1610  *cR_ptr = cR;
1611
1612  cP->Finalize();
1613  return cP;
1614 }
1615
1616
1618
1620 {
1621  point_matrix.SetSize(points[0].dim, np);
1622  for (int i = 0; i < np; i++)
1623  for (int j = 0; j < points[0].dim; j++)
1624  point_matrix(j, i) = points[i].coord[j];
1625 }
1626
1627 void NCMesh::GetFineTransforms(Element* elem, int coarse_index,
1628  FineTransform* transforms,
1629  const PointMatrix& pm)
1630 {
1631  if (!elem->ref_type)
1632  {
1633  // we got to a leaf, store the fine element transformation
1634  FineTransform& ft = transforms[elem->index];
1635  ft.coarse_index = coarse_index;
1636  pm.GetMatrix(ft.point_matrix);
1637  return;
1638  }
1639
1640  // recurse into the finer children, adjusting the point matrix
1641  if (elem->geom == Geometry::CUBE)
1642  {
1643  if (elem->ref_type == 1) // split along X axis
1644  {
1645  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
1646  Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
1647
1648  GetFineTransforms(elem->child[0], coarse_index, transforms,
1649  PointMatrix(pm(0), mid01, mid23, pm(3),
1650  pm(4), mid45, mid67, pm(7)));
1651
1652  GetFineTransforms(elem->child[1], coarse_index, transforms,
1653  PointMatrix(mid01, pm(1), pm(2), mid23,
1654  mid45, pm(5), pm(6), mid67));
1655  }
1656  else if (elem->ref_type == 2) // split along Y axis
1657  {
1658  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
1659  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
1660
1661  GetFineTransforms(elem->child[0], coarse_index, transforms,
1662  PointMatrix(pm(0), pm(1), mid12, mid30,
1663  pm(4), pm(5), mid56, mid74));
1664
1665  GetFineTransforms(elem->child[1], coarse_index, transforms,
1666  PointMatrix(mid30, mid12, pm(2), pm(3),
1667  mid74, mid56, pm(6), pm(7)));
1668  }
1669  else if (elem->ref_type == 4) // split along Z axis
1670  {
1671  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
1672  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
1673
1674  GetFineTransforms(elem->child[0], coarse_index, transforms,
1675  PointMatrix(pm(0), pm(1), pm(2), pm(3),
1676  mid04, mid15, mid26, mid37));
1677
1678  GetFineTransforms(elem->child[1], coarse_index, transforms,
1679  PointMatrix(mid04, mid15, mid26, mid37,
1680  pm(4), pm(5), pm(6), pm(7)));
1681  }
1682  else if (elem->ref_type == 3) // XY split
1683  {
1684  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
1685  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
1686  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
1687  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
1688
1689  Point midf0(mid23, mid12, mid01, mid30);
1690  Point midf5(mid45, mid56, mid67, mid74);
1691
1692  GetFineTransforms(elem->child[0], coarse_index, transforms,
1693  PointMatrix(pm(0), mid01, midf0, mid30,
1694  pm(4), mid45, midf5, mid74));
1695
1696  GetFineTransforms(elem->child[1], coarse_index, transforms,
1697  PointMatrix(mid01, pm(1), mid12, midf0,
1698  mid45, pm(5), mid56, midf5));
1699
1700  GetFineTransforms(elem->child[2], coarse_index, transforms,
1701  PointMatrix(midf0, mid12, pm(2), mid23,
1702  midf5, mid56, pm(6), mid67));
1703
1704  GetFineTransforms(elem->child[3], coarse_index, transforms,
1705  PointMatrix(mid30, midf0, mid23, pm(3),
1706  mid74, midf5, mid67, pm(7)));
1707  }
1708  else if (elem->ref_type == 5) // XZ split
1709  {
1710  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
1711  Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
1712  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
1713  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
1714
1715  Point midf1(mid01, mid15, mid45, mid04);
1716  Point midf3(mid23, mid37, mid67, mid26);
1717
1718  GetFineTransforms(elem->child[0], coarse_index, transforms,
1719  PointMatrix(pm(0), mid01, mid23, pm(3),
1720  mid04, midf1, midf3, mid37));
1721
1722  GetFineTransforms(elem->child[1], coarse_index, transforms,
1723  PointMatrix(mid01, pm(1), pm(2), mid23,
1724  midf1, mid15, mid26, midf3));
1725
1726  GetFineTransforms(elem->child[2], coarse_index, transforms,
1727  PointMatrix(midf1, mid15, mid26, midf3,
1728  mid45, pm(5), pm(6), mid67));
1729
1730  GetFineTransforms(elem->child[3], coarse_index, transforms,
1731  PointMatrix(mid04, midf1, midf3, mid37,
1732  pm(4), mid45, mid67, pm(7)));
1733  }
1734  else if (elem->ref_type == 6) // YZ split
1735  {
1736  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
1737  Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
1738  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
1739  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
1740
1741  Point midf2(mid12, mid26, mid56, mid15);
1742  Point midf4(mid30, mid04, mid74, mid37);
1743
1744  GetFineTransforms(elem->child[0], coarse_index, transforms,
1745  PointMatrix(pm(0), pm(1), mid12, mid30,
1746  mid04, mid15, midf2, midf4));
1747
1748  GetFineTransforms(elem->child[1], coarse_index, transforms,
1749  PointMatrix(mid30, mid12, pm(2), pm(3),
1750  midf4, midf2, mid26, mid37));
1751
1752  GetFineTransforms(elem->child[2], coarse_index, transforms,
1753  PointMatrix(mid04, mid15, midf2, midf4,
1754  pm(4), pm(5), mid56, mid74));
1755
1756  GetFineTransforms(elem->child[3], coarse_index, transforms,
1757  PointMatrix(midf4, midf2, mid26, mid37,
1758  mid74, mid56, pm(6), pm(7)));
1759  }
1760  else if (elem->ref_type == 7) // full isotropic refinement
1761  {
1762  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
1763  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
1764  Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
1765  Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
1766  Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
1767  Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
1768
1769  Point midf0(mid23, mid12, mid01, mid30);
1770  Point midf1(mid01, mid15, mid45, mid04);
1771  Point midf2(mid12, mid26, mid56, mid15);
1772  Point midf3(mid23, mid37, mid67, mid26);
1773  Point midf4(mid30, mid04, mid74, mid37);
1774  Point midf5(mid45, mid56, mid67, mid74);
1775
1776  Point midel(midf1, midf3);
1777
1778  GetFineTransforms(elem->child[0], coarse_index, transforms,
1779  PointMatrix(pm(0), mid01, midf0, mid30,
1780  mid04, midf1, midel, midf4));
1781
1782  GetFineTransforms(elem->child[1], coarse_index, transforms,
1783  PointMatrix(mid01, pm(1), mid12, midf0,
1784  midf1, mid15, midf2, midel));
1785
1786  GetFineTransforms(elem->child[2], coarse_index, transforms,
1787  PointMatrix(midf0, mid12, pm(2), mid23,
1788  midel, midf2, mid26, midf3));
1789
1790  GetFineTransforms(elem->child[3], coarse_index, transforms,
1791  PointMatrix(mid30, midf0, mid23, pm(3),
1792  midf4, midel, midf3, mid37));
1793
1794  GetFineTransforms(elem->child[4], coarse_index, transforms,
1795  PointMatrix(mid04, midf1, midel, midf4,
1796  pm(4), mid45, midf5, mid74));
1797
1798  GetFineTransforms(elem->child[5], coarse_index, transforms,
1799  PointMatrix(midf1, mid15, midf2, midel,
1800  mid45, pm(5), mid56, midf5));
1801
1802  GetFineTransforms(elem->child[6], coarse_index, transforms,
1803  PointMatrix(midel, midf2, mid26, midf3,
1804  midf5, mid56, pm(6), mid67));
1805
1806  GetFineTransforms(elem->child[7], coarse_index, transforms,
1807  PointMatrix(midf4, midel, midf3, mid37,
1808  mid74, midf5, mid67, pm(7)));
1809  }
1810  }
1811  else if (elem->geom == Geometry::SQUARE)
1812  {
1813  if (elem->ref_type == 1) // X split
1814  {
1815  Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
1816
1817  GetFineTransforms(elem->child[0], coarse_index, transforms,
1818  PointMatrix(pm(0), mid01, mid23, pm(3)));
1819
1820  GetFineTransforms(elem->child[1], coarse_index, transforms,
1821  PointMatrix(mid01, pm(1), pm(2), mid23));
1822  }
1823  else if (elem->ref_type == 2) // Y split
1824  {
1825  Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
1826
1827  GetFineTransforms(elem->child[0], coarse_index, transforms,
1828  PointMatrix(pm(0), pm(1), mid12, mid30));
1829
1830  GetFineTransforms(elem->child[1], coarse_index, transforms,
1831  PointMatrix(mid30, mid12, pm(2), pm(3)));
1832  }
1833  else if (elem->ref_type == 3) // iso split
1834  {
1835  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
1836  Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
1837  Point midel(mid01, mid23);
1838
1839  GetFineTransforms(elem->child[0], coarse_index, transforms,
1840  PointMatrix(pm(0), mid01, midel, mid30));
1841
1842  GetFineTransforms(elem->child[1], coarse_index, transforms,
1843  PointMatrix(mid01, pm(1), mid12, midel));
1844
1845  GetFineTransforms(elem->child[2], coarse_index, transforms,
1846  PointMatrix(midel, mid12, pm(2), mid23));
1847
1848  GetFineTransforms(elem->child[3], coarse_index, transforms,
1849  PointMatrix(mid30, midel, mid23, pm(3)));
1850  }
1851  }
1852  else if (elem->geom == Geometry::TRIANGLE)
1853  {
1854  Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
1855
1856  GetFineTransforms(elem->child[0], coarse_index, transforms,
1857  PointMatrix(pm(0), mid01, mid20));
1858
1859  GetFineTransforms(elem->child[1], coarse_index, transforms,
1860  PointMatrix(mid01, pm(1), mid12));
1861
1862  GetFineTransforms(elem->child[2], coarse_index, transforms,
1863  PointMatrix(mid20, mid12, pm(2)));
1864
1865  GetFineTransforms(elem->child[3], coarse_index, transforms,
1866  PointMatrix(mid01, mid12, mid20));
1867  }
1868 }
1869
1871 {
1872  if (!coarse_elements.Size())
1873  MFEM_ABORT("You need to call MarkCoarseLevel before calling Refine and "
1874  "GetFineTransformations.");
1875
1876  FineTransform* transforms = new FineTransform[leaf_elements.Size()];
1877
1878  // get transformations for fine elements, starting from coarse elements
1879  for (int i = 0; i < coarse_elements.Size(); i++)
1880  {
1881  Element* c_elem = coarse_elements[i];
1882  if (c_elem->ref_type)
1883  {
1884  if (c_elem->geom == Geometry::CUBE)
1885  {
1886  PointMatrix pm
1887  (Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0),
1888  Point(0,0,1), Point(1,0,1), Point(1,1,1), Point(0,1,1));
1889  GetFineTransforms(c_elem, i, transforms, pm);
1890  }
1891  else if (c_elem->geom == Geometry::SQUARE)
1892  {
1893  PointMatrix pm(Point(0,0), Point(1,0), Point(1,1), Point(0,1));
1894  GetFineTransforms(c_elem, i, transforms, pm);
1895  }
1896  else if (c_elem->geom == Geometry::TRIANGLE)
1897  {
1898  PointMatrix pm(Point(0,0), Point(1,0), Point(0,1));
1899  GetFineTransforms(c_elem, i, transforms, pm);
1900  }
1901  else
1903  }
1904  else
1905  {
1906  // element not refined, return identity transform
1907  transforms[c_elem->index].coarse_index = i;
1908  // leave point_matrix empty...
1909  }
1910  }
1911
1912  return transforms;
1913 }
1914
1915
1917
1919 {
1920  MFEM_ASSERT(n != NULL && n->p1 != n->p2, "Invalid node.");
1921
1922  Node *n1 = nodes.Peek(n->p1);
1923  Node *n2 = nodes.Peek(n->p2);
1924
1925  if ((n2->p1 != n2->p2) && (n1->id == n2->p1 || n1->id == n2->p2))
1926  {
1927  // (n1) is parent of (n2):
1928  // (n1)--(n)--(n2)----(*) or (*)----(n2)--(n)--(n1)
1929  if (n2->edge)
1930  return n2->edge->index;
1931  return GetEdgeMaster(n2);
1932  }
1933
1934  if ((n1->p1 != n1->p2) && (n2->id == n1->p1 || n2->id == n1->p2))
1935  {
1936  // (n2) is parent of (n1):
1937  // (n2)--(n)--(n1)----(*) or (*)----(n1)--(n)--(n2)
1938  if (n1->edge)
1939  return n1->edge->index;
1940  return GetEdgeMaster(n1);
1941  }
1942
1943  return n->edge ? n->edge->index : -1;
1944 }
1945
1946 int NCMesh::GetEdgeMaster(int v1, int v2) const
1947 {
1948  Node *n = nodes.Peek(vertex_nodeId[v1], vertex_nodeId[v2]);
1949  int master_edge = GetEdgeMaster(n);
1950  MFEM_ASSERT(n->edge != NULL, "Invalid edge.");
1951  return (n->edge->index != master_edge) ? master_edge : -1;
1952 }
1953
1954 void NCMesh::FaceSplitLevel(Node* v1, Node* v2, Node* v3, Node* v4,
1955  int& h_level, int& v_level)
1956 {
1957  int hl1, hl2, vl1, vl2;
1958  Node* mid[4];
1959
1960  switch (FaceSplitType(v1, v2, v3, v4, mid))
1961  {
1962  case 0: // not split
1963  h_level = v_level = 0;
1964  break;
1965
1966  case 1: // vertical
1967  FaceSplitLevel(v1, mid[0], mid[2], v4, hl1, vl1);
1968  FaceSplitLevel(mid[0], v2, v3, mid[2], hl2, vl2);
1969  h_level = std::max(hl1, hl2);
1970  v_level = std::max(vl1, vl2) + 1;
1971  break;
1972
1973  default: // horizontal
1974  FaceSplitLevel(v1, v2, mid[1], mid[3], hl1, vl1);
1975  FaceSplitLevel(mid[3], mid[1], v3, v4, hl2, vl2);
1976  h_level = std::max(hl1, hl2) + 1;
1977  v_level = std::max(vl1, vl2);
1978  }
1979 }
1980
1981 static int max4(int a, int b, int c, int d)
1982 {
1983  return std::max(std::max(a, b), std::max(c, d));
1984 }
1985
1986 void NCMesh::CountSplits(Element* elem, int splits[3])
1987 {
1988  Node** node = elem->node;
1989  GeomInfo& gi = GI[elem->geom];
1990
1991  MFEM_ASSERT(elem->geom == Geometry::CUBE, "TODO");
1992  // TODO: triangles and quads
1993
1994  int level[6][2];
1995  for (int i = 0; i < gi.nf; i++)
1996  {
1997  const int* fv = gi.faces[i];
1998  FaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]],
1999  level[i][1], level[i][0]);
2000  }
2001
2002  splits[0] = max4(level[0][0], level[1][0], level[3][0], level[5][0]);
2003  splits[1] = max4(level[0][1], level[2][0], level[4][0], level[5][1]);
2004  splits[2] = max4(level[1][1], level[2][1], level[3][1], level[4][1]);
2005 }
2006
2007 void NCMesh::LimitNCLevel(int max_level)
2008 {
2009  if (max_level < 1)
2010  MFEM_ABORT("'max_level' must be 1 or greater.");
2011
2012  while (1)
2013  {
2014  Array<Refinement> refinements;
2015  for (int i = 0; i < leaf_elements.Size(); i++)
2016  {
2017  int splits[3];
2018  CountSplits(leaf_elements[i], splits);
2019
2020  int ref_type = 0;
2021  for (int k = 0; k < 3; k++)
2022  if (splits[k] > max_level)
2023  ref_type |= (1 << k);
2024
2025  // TODO: isotropic meshes should only be modified with iso refinements
2026
2027  if (ref_type)
2028  refinements.Append(Refinement(i, ref_type));
2029  }
2030
2031  if (!refinements.Size()) break;
2032
2033  Refine(refinements);
2034  }
2035 }
2036
2038 {
2039  int n = 1;
2040  if (elem->ref_type)
2041  {
2042  for (int i = 0; i < 8; i++)
2043  if (elem->child[i])
2044  n += CountElements(elem->child[i]);
2045  }
2046  return n;
2047 }
2048
2050 {
2051  int num_elem = 0;
2052  for (int i = 0; i < root_elements.Size(); i++)
2053  num_elem += CountElements(root_elements[i]);
2054
2055  int num_vert = 0, num_edges = 0;
2056  for (HashTable<Node>::Iterator it(nodes); it; ++it)
2057  {
2058  if (it->vertex) num_vert++;
2059  if (it->edge) num_edges++;
2060  }
2061
2062  return num_elem * sizeof(Element) +
2063  num_vert * sizeof(Vertex) +
2064  num_edges * sizeof(Edge) +
2065  nodes.MemoryUsage() +
2066  faces.MemoryUsage() +
2067  root_elements.Capacity() * sizeof(Element*) +
2068  leaf_elements.Capacity() * sizeof(Element*) +
2069  vertex_nodeId.Capacity() * sizeof(int) +
2070  sizeof(*this);
2071 }
2072
2073 } // namespace mfem
Abstract class for Finite Elements.
Definition: fe.hpp:42
void ReorderFacePointMat(Node *v0, Node *v1, Node *v2, Node *v3, Element *elem, DenseMatrix &pm)
Definition: ncmesh.cpp:1251
void LimitNCLevel(int max_level)
Definition: ncmesh.cpp:2007
int Size() const
Logical size of the array.
Definition: array.hpp:108
Element * NewHexahedron(Node *n0, Node *n1, Node *n2, Node *n3, Node *n4, Node *n5, Node *n6, Node *n7, int attr, int fattr0, int fattr1, int fattr2, int fattr3, int fattr4, int fattr5)
Definition: ncmesh.cpp:371
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:1525
DepList dep_list
list of other DOFs this DOF depends on
Definition: ncmesh.hpp:357
void GetVerticesElementsBoundary(Array< mfem::Vertex > &vertices, Array< mfem::Element * > &elements, Array< mfem::Element * > &boundary)
Definition: ncmesh.cpp:1092
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:149
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:421
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
Definition: mesh.cpp:3451
int index
edge number in the Mesh
Definition: ncmesh.hpp:183
Array< Element * > coarse_elements
Definition: ncmesh.hpp:268
Iterator over items contained in the HashTable.
Definition: hash.hpp:154
bool finalized
true if cP matrix row is known for this DOF
Definition: ncmesh.hpp:356
FiniteElementSpace * space
Definition: ncmesh.hpp:365
FineTransform * GetFineTransforms()
Definition: ncmesh.cpp:1870
SparseMatrix * GetInterpolation(FiniteElementSpace *space, SparseMatrix **cR_ptr=NULL)
Definition: ncmesh.cpp:1470
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
static const int NumGeom
Definition: geom.hpp:34
Node * PeekAltParents(Node *v1, Node *v2)
Definition: ncmesh.cpp:309
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:399
void SetEdgeIndicesFromMesh(Mesh *mesh)
Definition: ncmesh.cpp:1166
Array< RefStackItem > ref_stack
stack of scheduled refinements
Definition: ncmesh.hpp:284
int attribute
boundary element attribute, -1 if internal edge
Definition: ncmesh.hpp:182
bool Boundary() const
Definition: ncmesh.hpp:186
Data type dense matrix.
Definition: densemat.hpp:22
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:148
void RefElementNodes(Element *elem)
Definition: ncmesh.cpp:209
static int find_node(Element *elem, Node *node)
Definition: ncmesh.cpp:1225
void GetLeafElements(Element *e)
Definition: ncmesh.cpp:1068
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
const Element * GetFace(int i) const
Definition: mesh.hpp:432
int index
vertex number in the Mesh
Definition: ncmesh.hpp:172
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe.cpp:69
void ConstrainEdge(Node *v0, Node *v1, double t0, double t1, Array< int > &master_dofs, int level)
Definition: ncmesh.cpp:1313
bool DofFinalizable(DofData &vd)
Definition: ncmesh.cpp:1460
void ProcessMasterEdge(Node *node[2], Node *edge)
Definition: ncmesh.cpp:1426
void Delete(ItemT *item)
Remove an item from the hash table and also delete the item itself.
Definition: hash.hpp:403
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:1728
Vertex * NewVertex(Node *v1, Node *v2)
Definition: ncmesh.cpp:450
bool NodeSetY1(Node *node, Node **n)
Definition: ncmesh.cpp:505
Array< int > vertex_nodeId
Definition: ncmesh.hpp:270
int attribute
boundary element attribute, -1 if internal face
Definition: ncmesh.hpp:228
int GetGeometryType() const
Definition: element.hpp:49
int index
Mesh element number.
Definition: ncmesh.hpp:37
Element * NewQuadrilateral(Node *n0, Node *n1, Node *n2, Node *n3, int attr, int eattr0, int eattr1, int eattr2, int eattr3)
Definition: ncmesh.cpp:399
int coarse_index
coarse Mesh element index
Definition: ncmesh.hpp:100
virtual void GetFaceDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:888
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
Definition: segment.cpp:35
Data type hexahedron element.
Definition: hexahedron.hpp:22
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Definition: sparsemat.cpp:1661
void ConstrainFace(Node *v0, Node *v1, Node *v2, Node *v3, const PointMatrix &pm, Array< int > &master_dofs, int level)
Definition: ncmesh.cpp:1357
int index
face number in the Mesh
Definition: ncmesh.hpp:229
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
Definition: operator.hpp:35
void ProcessMasterFace(Node *node[4], Face *face)
Definition: ncmesh.cpp:1444
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:330
bool Independent() const
Definition: ncmesh.hpp:360
void GetMatrix(DenseMatrix &point_matrix) const
Definition: ncmesh.cpp:1619
void CheckAnisoFace(Node *v1, Node *v2, Node *v3, Node *v4, Node *mid12, Node *mid34, int level=0)
Definition: ncmesh.cpp:550
DofData * dof_data
DOF temporary data.
Definition: ncmesh.hpp:363
void UpdateVertices()
Definition: ncmesh.cpp:1050
void RegisterElement(Element *e)
Definition: ncmesh.cpp:262
DenseMatrix point_matrix
for use in IsoparametricTransformation
Definition: ncmesh.hpp:101
void AddDependencies(Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I)
Definition: ncmesh.cpp:1286
Element(int geom, int attr)
Definition: ncmesh.cpp:359
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:79
Data type triangle element.
Definition: triangle.hpp:22
const Element * GetElement(int i) const
Definition: mesh.hpp:424
int decode_dof(int dof, double &sign)
Definition: ncmesh.cpp:1278
int CountElements(Element *elem)
Definition: ncmesh.cpp:2037
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:66
virtual void Finalize(int skip_zeros=1)
Definition: sparsemat.cpp:529
bool NodeSetX1(Node *node, Node **n)
Definition: ncmesh.cpp:499
int Dimension() const
Definition: mesh.hpp:417
void RegisterFaces(Element *elem)
Definition: ncmesh.cpp:282
void GetEdgeDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:925
void DeleteHierarchy(Element *elem)
Definition: ncmesh.cpp:157
Vertex * vertex
Definition: ncmesh.hpp:201
void UnrefEdge(HashTable< Node > &nodes)
Definition: ncmesh.cpp:194
void ForceRefinement(Node *v1, Node *v2, Node *v3, Node *v4)
Definition: ncmesh.cpp:518
Abstract finite element space.
Definition: fespace.hpp:61
int GetDof() const
Returns the degrees of freedom in the FE space.
Definition: fe.hpp:86
bool Boundary() const
Definition: ncmesh.hpp:235
double pos[3]
3D position
Definition: ncmesh.hpp:171
void swap(Vector *v1, Vector *v2)
Definition: vector.cpp:213
void mfem_error(const char *msg)
Definition: error.cpp:23
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:293
Element * child[8]
Definition: ncmesh.hpp:260
Array< Element * > leaf_elements
Definition: ncmesh.hpp:267
void CountSplits(Element *elem, int splits[3])
Definition: ncmesh.cpp:1986
void UpdateLeafElements()
Definition: ncmesh.cpp:1084
bool NodeSetZ1(Node *node, Node **n)
Definition: ncmesh.cpp:511
Node * GetMidEdgeVertexSimple(Node *v1, Node *v2)
Definition: ncmesh.cpp:472
HashTable< Node > nodes
Definition: ncmesh.hpp:272
int GetNV() const
Definition: mesh.hpp:393
void UnrefElementNodes(Element *elem)
Definition: ncmesh.cpp:235
void CheckIsoFace(Node *v1, Node *v2, Node *v3, Node *v4, Node *e1, Node *e2, Node *e3, Node *e4, Node *midf)
Definition: ncmesh.cpp:601
int FaceSplitType(Node *v1, Node *v2, Node *v3, Node *v4, Node *mid[4]=NULL)
Definition: ncmesh.cpp:1196
Node * GetMidFaceVertex(Node *e1, Node *e2, Node *e3, Node *e4)
Definition: ncmesh.cpp:481
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:67
bool NodeSetZ2(Node *node, Node **n)
Definition: ncmesh.cpp:514
int Capacity() const
Definition: array.hpp:112
bool NodeSetX2(Node *node, Node **n)
Definition: ncmesh.cpp:502
void SetFaceIndicesFromMesh(Mesh *mesh)
Definition: ncmesh.cpp:1180
long MemoryUsage()
Definition: ncmesh.cpp:2049
Element * NewTriangle(Node *n0, Node *n1, Node *n2, int attr, int eattr0, int eattr1, int eattr2)
Definition: ncmesh.cpp:426
NCMesh(const Mesh *mesh)
Definition: ncmesh.cpp:61
virtual int GetNVertices() const =0
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:82
Array< Element * > root_elements
Definition: ncmesh.hpp:266
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:405
Vector data type.
Definition: vector.hpp:29
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:162
Node * GetMidEdgeVertex(Node *v1, Node *v2)
Definition: ncmesh.cpp:463
Element * GetSingleElement() const
Definition: ncmesh.cpp:295
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element&#39;s vertices.
void SetSize(int s)
If the matrix is not a square matrix of size s then recreate it.
Definition: densemat.cpp:73
void ForgetElement(Element *e)
Definition: ncmesh.cpp:272
bool NodeSetY2(Node *node, Node **n)
Definition: ncmesh.cpp:508
int ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:38
Abstract data type element.
Definition: element.hpp:27
Data type line segment element.
Definition: segment.hpp:22
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:41
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:428
virtual int GetType() const =0
Returns element&#39;s type.
void FaceSplitLevel(Node *v1, Node *v2, Node *v3, Node *v4, int &h_level, int &v_level)
Definition: ncmesh.cpp:1954
HashTable< Face > faces
Definition: ncmesh.hpp:273
void Refine(const Array< Refinement > &refinements)
Definition: ncmesh.cpp:1005
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const =0
int GetEdgeMaster(int v1, int v2) const
Definition: ncmesh.cpp:1946
void UnrefVertex(HashTable< Node > &nodes)
Definition: ncmesh.cpp:187