12 #include "../fem/fem.hpp"
37 for (
int i = 0; i <
ne; i++)
39 for (
int j = 0; j < 2; j++)
44 for (
int i = 0; i <
nf; i++)
46 for (
int j = 0; j <
nfv; j++)
55 for (
int i = 0; i <
ne; i++)
74 Iso = vertex_parents ?
false :
true;
79 for (
int i = 0; i < mesh->
GetNE(); i++)
83 for (
int j = 0; j < nv; j++)
85 max_id = std::max(max_id, v[j]);
88 for (
int i = 0; i <= max_id; i++)
92 MFEM_CONTRACT_VAR(node);
93 MFEM_ASSERT(node->
id == i,
"");
104 for (
int i = 0; i < mesh->
GetNE(); i++)
114 MFEM_ABORT(
"only triangles, quads and hexes are supported by NCMesh.");
124 for (
int j = 0; j <
GI[
geom].
nv; j++)
129 if (v[j] < mesh->
GetNV())
132 const double* pos = mesh->
GetVertex(v[j]);
141 nc_elem->
node[j] = node;
153 for (
int i = 0; i < mesh->
GetNBE(); i++)
161 node[i] =
nodes.Peek(v[i]);
162 MFEM_VERIFY(node[i],
"boundary elements inconsistent.");
167 Face* face =
faces.Peek(node[0], node[1], node[2], node[3]);
168 MFEM_VERIFY(face,
"boundary face not found.");
173 Edge* edge =
nodes.Peek(node[0], node[1])->edge;
174 MFEM_VERIFY(edge,
"boundary edge not found.");
179 MFEM_ABORT(
"only segment and quadrilateral boundary "
180 "elements are supported by NCMesh.");
203 for (
int i = 0; i < 8; i++)
215 for (
int i = 0; i < gi.
nv; i++)
228 for (
int i = 0; i < 8; i++)
267 MFEM_ASSERT(
vertex,
"can't create vertex here.");
273 if (!edge) { edge =
new Edge; }
279 MFEM_ASSERT(vertex,
"cannot unref a nonexistent vertex.");
280 if (!vertex->Unref()) { vertex = NULL; }
281 if (!vertex && !edge) { nodes.
Delete(
this); }
286 MFEM_ASSERT(node,
"node not found.");
287 MFEM_ASSERT(node->
edge,
"cannot unref a nonexistent edge.");
294 std::memcpy(
this, &other,
sizeof(*
this));
295 if (vertex) { vertex =
new Vertex(*vertex); }
296 if (edge) { edge =
new Edge(*edge); }
301 MFEM_ASSERT(!vertex && !edge,
"node was not unreffed properly.");
302 if (vertex) {
delete vertex; }
303 if (edge) {
delete edge; }
312 for (
int i = 0; i < gi.
nv; i++)
318 for (
int i = 0; i < gi.
ne; i++)
320 const int* ev = gi.
edges[i];
321 nodes.Get(node[ev[0]], node[ev[1]])->RefEdge();
325 for (
int i = 0; i < gi.
nf; i++)
327 const int* fv = gi.
faces[i];
328 faces.Get(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]])->Ref();
341 for (
int i = 0; i < gi.
nf; i++)
343 const int* fv = gi.
faces[i];
344 Face* face =
faces.Peek(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
350 for (
int i = 0; i < gi.
ne; i++)
352 const int* ev = gi.
edges[i];
358 for (
int i = 0; i < gi.
nv; i++)
373 std::memcpy(
this, &other,
sizeof(*
this));
374 elem[0] = elem[1] = NULL;
379 if (elem[0] == NULL) { elem[0] = e; }
380 else if (elem[1] == NULL) { elem[1] = e; }
381 else { MFEM_ABORT(
"can't have 3 elements in Face::elem[]."); }
386 if (elem[0] == e) { elem[0] = NULL; }
387 else if (elem[1] == e) { elem[1] = NULL; }
388 else { MFEM_ABORT(
"element not found in Face::elem[]."); }
394 const int* fv = gi.
faces[face_no];
396 return faces.Peek(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
402 for (
int i = 0; i < gi.
nf; i++)
406 if (fattr) { face->
attribute = fattr[i]; }
414 MFEM_ASSERT(!elem[1],
"not a single element face.");
419 MFEM_ASSERT(elem[1],
"no elements in face.");
452 if ((v1->
p1 != v1->
p2) && (v2->
p1 != v2->
p2))
468 mid =
nodes.Peek(w1, w2);
479 : geom(geom), ref_type(0), flag(0), index(-1), rank(0), attribute(attr)
494 int fattr0,
int fattr1,
int fattr2,
495 int fattr3,
int fattr4,
int fattr5)
504 for (
int i = 0; i < gi_hex.
nf; i++)
506 const int* fv = gi_hex.
faces[i];
521 int eattr0,
int eattr1,
int eattr2,
int eattr3)
529 for (
int i = 0; i < gi_quad.
ne; i++)
531 const int* ev = gi_quad.
edges[i];
534 edge[i] = node->
edge;
547 int attr,
int eattr0,
int eattr1,
int eattr2)
555 for (
int i = 0; i < gi_tri.
ne; i++)
557 const int* ev = gi_tri.
edges[i];
560 edge[i] = node->
edge;
572 MFEM_ASSERT(v1->
vertex && v2->
vertex,
"missing parent vertices.");
576 for (
int i = 0; i < 3; i++)
588 if (!mid) { mid =
nodes.Get(v1, v2); }
613 midf =
nodes.Get(e2, e4);
621 {
return node == n[0] || node == n[3] || node == n[4] || node == n[7]; }
624 {
return node == n[1] || node == n[2] || node == n[5] || node == n[6]; }
627 {
return node == n[0] || node == n[1] || node == n[4] || node == n[5]; }
630 {
return node == n[2] || node == n[3] || node == n[6] || node == n[7]; }
633 {
return node == n[0] || node == n[1] || node == n[2] || node == n[3]; }
636 {
return node == n[4] || node == n[5] || node == n[6] || node == n[7]; }
643 if (!face) {
return; }
646 MFEM_ASSERT(!elem->
ref_type,
"element already refined.");
668 MFEM_ABORT(
"inconsistent element/face structure.");
674 Node* mid12,
Node* mid34,
int level)
707 nodes.Reparent(midf, mid12->
id, mid34->
id);
758 if (!ref_type) {
return; }
763 char remaining = ref_type & ~elem->
ref_type;
766 for (
int i = 0; i < 8; i++)
777 memset(child, 0,
sizeof(child));
784 for (
int i = 0; i < gi_hex.
nf; i++)
786 const int* fv = gi_hex.
faces[i];
787 Face* face =
faces.Peek(no[fv[0]], no[fv[1]], no[fv[2]], no[fv[3]]);
813 no[4], mid45, mid67, no[7], attr,
814 fa[0], fa[1], -1, fa[3], fa[4], fa[5]);
817 mid45, no[5], no[6], mid67, attr,
818 fa[0], fa[1], fa[2], fa[3], -1, fa[5]);
825 else if (ref_type == 2)
833 no[4], no[5], mid56, mid74, attr,
834 fa[0], fa[1], fa[2], -1, fa[4], fa[5]);
837 mid74, mid56, no[6], no[7], attr,
838 fa[0], -1, fa[2], fa[3], fa[4], fa[5]);
845 else if (ref_type == 4)
853 mid04, mid15, mid26, mid37, attr,
854 fa[0], fa[1], fa[2], fa[3], fa[4], -1);
857 no[4], no[5], no[6], no[7], attr,
858 -1, fa[1], fa[2], fa[3], fa[4], fa[5]);
865 else if (ref_type == 3)
881 no[4], mid45, midf5, mid74, attr,
882 fa[0], fa[1], -1, -1, fa[4], fa[5]);
885 mid45, no[5], mid56, midf5, attr,
886 fa[0], fa[1], fa[2], -1, -1, fa[5]);
889 midf5, mid56, no[6], mid67, attr,
890 fa[0], -1, fa[2], fa[3], -1, fa[5]);
893 mid74, midf5, mid67, no[7], attr,
894 fa[0], -1, -1, fa[3], fa[4], fa[5]);
901 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
902 CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
904 else if (ref_type == 5)
920 mid04, midf1, midf3, mid37, attr,
921 fa[0], fa[1], -1, fa[3], fa[4], -1);
924 midf1, mid15, mid26, midf3, attr,
925 fa[0], fa[1], fa[2], fa[3], -1, -1);
928 mid45, no[5], no[6], mid67, attr,
929 -1, fa[1], fa[2], fa[3], -1, fa[5]);
932 no[4], mid45, mid67, no[7], attr,
933 -1, fa[1], -1, fa[3], fa[4], fa[5]);
940 CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
941 CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
943 else if (ref_type == 6)
959 mid04, mid15, midf2, midf4, attr,
960 fa[0], fa[1], fa[2], -1, fa[4], -1);
963 midf4, midf2, mid26, mid37, attr,
964 fa[0], -1, fa[2], fa[3], fa[4], -1);
967 no[4], no[5], mid56, mid74, attr,
968 -1, fa[1], fa[2], -1, fa[4], fa[5]);
971 mid74, mid56, no[6], no[7], attr,
972 -1, -1, fa[2], fa[3], fa[4], fa[5]);
979 CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
980 CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
982 else if (ref_type == 7)
1009 mid04, midf1, midel, midf4, attr,
1010 fa[0], fa[1], -1, -1, fa[4], -1);
1013 midf1, mid15, midf2, midel, attr,
1014 fa[0], fa[1], fa[2], -1, -1, -1);
1017 midel, midf2, mid26, midf3, attr,
1018 fa[0], -1, fa[2], fa[3], -1, -1);
1021 midf4, midel, midf3, mid37, attr,
1022 fa[0], -1, -1, fa[3], fa[4], -1);
1025 no[4], mid45, midf5, mid74, attr,
1026 -1, fa[1], -1, -1, fa[4], fa[5]);
1029 mid45, no[5], mid56, midf5, attr,
1030 -1, fa[1], fa[2], -1, -1, fa[5]);
1033 midf5, mid56, no[6], mid67, attr,
1034 -1, -1, fa[2], fa[3], -1, fa[5]);
1037 mid74, midf5, mid67, no[7], attr,
1038 -1, -1, -1, fa[3], fa[4], fa[5]);
1040 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1041 CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1042 CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1043 CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1044 CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1045 CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1049 MFEM_ABORT(
"invalid refinement type.");
1052 if (ref_type != 7) {
Iso =
false; }
1057 int ea0 =
nodes.Peek(no[0], no[1])->edge->attribute;
1058 int ea1 =
nodes.Peek(no[1], no[2])->edge->attribute;
1059 int ea2 =
nodes.Peek(no[2], no[3])->edge->attribute;
1060 int ea3 =
nodes.Peek(no[3], no[0])->edge->attribute;
1070 attr, ea0, -1, ea2, ea3);
1073 attr, ea0, ea1, ea2, -1);
1075 else if (ref_type == 2)
1081 attr, ea0, ea1, -1, ea3);
1084 attr, -1, ea1, ea2, ea3);
1086 else if (ref_type == 3)
1096 attr, ea0, -1, -1, ea3);
1099 attr, ea0, ea1, -1, -1);
1102 attr, -1, ea1, ea2, -1);
1105 attr, -1, -1, ea2, ea3);
1109 MFEM_ABORT(
"Invalid refinement type.");
1112 if (ref_type != 3) {
Iso =
false; }
1117 int ea0 =
nodes.Peek(no[0], no[1])->edge->attribute;
1118 int ea1 =
nodes.Peek(no[1], no[2])->edge->attribute;
1119 int ea2 =
nodes.Peek(no[2], no[0])->edge->attribute;
1128 child[0] =
NewTriangle(no[0], mid01, mid20, attr, ea0, -1, ea2);
1129 child[1] =
NewTriangle(mid01, no[1], mid12, attr, ea0, ea1, -1);
1130 child[2] =
NewTriangle(mid20, mid12, no[2], attr, -1, ea1, ea2);
1131 child[3] =
NewTriangle(mid01, mid12, mid20, attr, -1, -1, -1);
1135 MFEM_ABORT(
"Unsupported element geometry.");
1139 for (
int i = 0; i < 8; i++)
1148 for (
int i = 0; i < 8; i++)
1154 for (
int i = 0; i < 8; i++)
1165 memcpy(elem->
child, child,
sizeof(elem->
child));
1172 for (
int i = refinements.
Size()-1; i >= 0; i--)
1201 std::cout <<
"Refined " << refinements.
Size() <<
" + " << nforced
1202 <<
" elements" << std::endl;
1216 std::memcpy(child, elem->
child,
sizeof(child));
1219 for (
int i = 0; i < 8; i++)
1221 if (child[i] && child[i]->ref_type)
1231 const int table[7][8 + 6] =
1233 { 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0 },
1234 { 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1 },
1235 { 0, 1, 2, 3, 0, 1, 2, 3, 1, 1, 1, 3, 3, 3 },
1236 { 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1 },
1237 { 0, 1, 1, 0, 3, 2, 2, 3, 1, 1, 1, 3, 3, 3 },
1238 { 0, 0, 1, 1, 2, 2, 3, 3, 0, 0, 0, 3, 3, 3 },
1239 { 0, 1, 2, 3, 4, 5, 6, 7, 1, 1, 1, 7, 7, 7 }
1241 for (
int i = 0; i < 8; i++)
1245 for (
int i = 0; i < 6; i++)
1248 const int* fv = gi_hex.
faces[i];
1250 ch->
node[fv[2]], ch->
node[fv[3]])->attribute;
1255 const int table[3][4 + 4] =
1257 { 0, 1, 1, 0, 1, 1, 0, 0 },
1258 { 0, 0, 1, 1, 0, 0, 1, 1 },
1259 { 0, 1, 2, 3, 1, 1, 3, 3 }
1261 for (
int i = 0; i < 4; i++)
1265 for (
int i = 0; i < 4; i++)
1268 const int* ev = gi_quad.
edges[i];
1269 fa[i] =
nodes.Peek(ch->
node[ev[0]], ch->
node[ev[1]])->edge->attribute;
1274 for (
int i = 0; i < 3; i++)
1278 const int* ev = gi_tri.
edges[i];
1279 fa[i] =
nodes.Peek(ch->
node[ev[0]], ch->
node[ev[1]])->edge->attribute;
1284 MFEM_ABORT(
"Unsupported element geometry.");
1291 elem->
rank = INT_MAX;
1292 for (
int i = 0; i < 8; i++)
1309 for (
int i = 0; i < gi.
ne; i++)
1311 const int* ev = gi.
edges[i];
1312 nodes.Peek(node[ev[0]], node[ev[1]])->edge->attribute = fa[i];
1324 int total = 0, ref = 0, ghost = 0;
1325 for (
int i = 0; i < 8; i++)
1331 if (ch->
ref_type) { ref++;
break; }
1336 if (!ref && ghost < total)
1339 int next_row = list.
Size() ? (list.
Last().from + 1) : 0;
1340 for (
int i = 0; i < 8; i++)
1348 for (
int i = 0; i < 8; i++)
1366 int size = list.
Size() ? (list.
Last().from + 1) : 0;
1375 for (
int i = 0; i < deref_table.
Size(); i++)
1377 const int* fine = deref_table.
GetRow(i), size = deref_table.
RowSize(i);
1381 for (
int j = 0; j < size; j++)
1386 for (
int k = 0; k <
Dim; k++)
1388 if ((parent->
ref_type & (1 << k)) &&
1389 splits[k] >= max_nc_level)
1402 MFEM_VERIFY(
Dim < 3 ||
Iso,
1403 "derefinement of 3D anisotropic meshes not implemented yet.");
1411 for (
int i = 0; i < derefs.
Size(); i++)
1413 int row = derefs[i];
1415 "invalid derefinement number.");
1430 for (
int i = 0; i < coarse.
Size(); i++)
1441 for (
int i = 0; i < nfine; i++)
1454 for (
int i = 0; i < 8; i++)
1457 if (ch && ch->
index >= 0)
1459 int code = (parent->
ref_type << 3) + i;
1461 coarse[ch->
index] = parent;
1475 if (it->vertex) { it->vertex->index = num_vert++; }
1487 static char quad_hilbert_child_order[8][4] =
1489 {0,1,2,3}, {0,3,2,1}, {1,2,3,0}, {1,0,3,2},
1490 {2,3,0,1}, {2,1,0,3}, {3,0,1,2}, {3,2,1,0}
1492 static char quad_hilbert_child_state[8][4] =
1494 {1,0,0,5}, {0,1,1,4}, {3,2,2,7}, {2,3,3,6},
1495 {5,4,4,1}, {4,5,5,0}, {7,6,6,3}, {6,7,7,2}
1497 static char hex_hilbert_child_order[24][8] =
1499 {0,1,2,3,7,6,5,4}, {0,3,7,4,5,6,2,1}, {0,4,5,1,2,6,7,3},
1500 {1,0,3,2,6,7,4,5}, {1,2,6,5,4,7,3,0}, {1,5,4,0,3,7,6,2},
1501 {2,1,5,6,7,4,0,3}, {2,3,0,1,5,4,7,6}, {2,6,7,3,0,4,5,1},
1502 {3,0,4,7,6,5,1,2}, {3,2,1,0,4,5,6,7}, {3,7,6,2,1,5,4,0},
1503 {4,0,1,5,6,2,3,7}, {4,5,6,7,3,2,1,0}, {4,7,3,0,1,2,6,5},
1504 {5,1,0,4,7,3,2,6}, {5,4,7,6,2,3,0,1}, {5,6,2,1,0,3,7,4},
1505 {6,2,3,7,4,0,1,5}, {6,5,1,2,3,0,4,7}, {6,7,4,5,1,0,3,2},
1506 {7,3,2,6,5,1,0,4}, {7,4,0,3,2,1,5,6}, {7,6,5,4,0,1,2,3}
1508 static char hex_hilbert_child_state[24][8] =
1510 {1,2,2,7,7,21,21,17}, {2,0,0,22,22,16,16,8}, {0,1,1,15,15,6,6,23},
1511 {4,5,5,10,10,18,18,14}, {5,3,3,19,19,13,13,11}, {3,4,4,12,12,9,9,20},
1512 {8,7,7,17,17,23,23,2}, {6,8,8,0,0,15,15,22}, {7,6,6,21,21,1,1,16},
1513 {11,10,10,14,14,20,20,5}, {9,11,11,3,3,12,12,19}, {10,9,9,18,18,4,4,13},
1514 {13,14,14,5,5,19,19,10}, {14,12,12,20,20,11,11,4}, {12,13,13,9,9,3,3,18},
1515 {16,17,17,2,2,22,22,7}, {17,15,15,23,23,8,8,1}, {15,16,16,6,6,0,0,21},
1516 {20,19,19,11,11,14,14,3}, {18,20,20,4,4,10,10,12}, {19,18,18,13,13,5,5,9},
1517 {23,22,22,8,8,17,17,0}, {21,23,23,1,1,7,7,15}, {22,21,21,16,16,2,2,6}
1524 if (elem->
rank >= 0)
1533 for (
int i = 0; i < 4; i++)
1535 int ch = quad_hilbert_child_order[state][i];
1536 int st = quad_hilbert_child_state[state][i];
1542 for (
int i = 0; i < 8; i++)
1544 int ch = hex_hilbert_child_order[state][i];
1545 int st = hex_hilbert_child_state[state][i];
1551 for (
int i = 0; i < 8; i++)
1591 MFEM_ABORT(
"invalid geometry");
1601 for (
int i = 0; i < vertices.
Size(); i++)
1604 vertices[i].SetCoords(node->
vertex->
pos);
1616 if (
IsGhost(nc_elem)) {
continue; }
1625 for (
int j = 0; j < gi.
nv; j++)
1633 for (
int k = 0; k < gi.
nf; k++)
1635 const int* fv = gi.
faces[k];
1636 Face* face =
faces.Peek(node[fv[0]], node[fv[1]],
1637 node[fv[2]], node[fv[3]]);
1642 for (
int j = 0; j < 4; j++)
1652 for (
int k = 0; k < gi.
ne; k++)
1654 const int* ev = gi.
edges[k];
1655 Edge* edge =
nodes.Peek(node[ev[0]], node[ev[1]])->edge;
1660 for (
int j = 0; j < 2; j++)
1664 boundary.
Append(segment);
1676 for (
int i = 0; i < edge_vertex->
Size(); i++)
1678 const int *ev = edge_vertex->
GetRow(i);
1681 MFEM_ASSERT(node && node->
edge,
"edge not found.");
1688 const int* fv = mesh->
GetFace(i)->GetVertices();
1692 MFEM_ASSERT(mesh->
GetFace(i)->GetNVertices() == 4,
"");
1698 MFEM_ASSERT(mesh->
GetFace(i)->GetNVertices() == 2,
"");
1700 face =
faces.Peek(n0, n0, n1, n1);
1704 const int *ev = edge_vertex->
GetRow(i);
1705 MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
1706 (ev[1] == fv[0] && ev[0] == fv[1]),
"");
1709 MFEM_ASSERT(face,
"face not found.");
1720 MFEM_ASSERT(
Dim >= 3,
"");
1725 Node* e3 = e1 ?
nodes.Peek(v3, v4) : NULL;
1726 Node* e4 = e2 ?
nodes.Peek(v4, v1) : NULL;
1729 if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
1732 Node *midf1 = NULL, *midf2 = NULL;
1733 if (e1 && e3) { midf1 =
nodes.Peek(e1, e3); }
1734 if (e2 && e4) { midf2 =
nodes.Peek(e2, e4); }
1737 MFEM_ASSERT(!(midf1 && midf2),
"incorrectly split face!");
1739 if (!midf1 && !midf2) {
return 0; }
1741 if (midf1) {
return 1; }
1747 for (
int i = 0; i < 8; i++)
1748 if (elem->
node[i] == node) {
return i; }
1750 MFEM_ABORT(
"Node not found.");
1756 for (
int i = 0; i < 8; i++)
1757 if (elem->
node[i]->
id == node_id) {
return i; }
1759 MFEM_ABORT(
"Node not found.");
1767 for (
int i = 0; i < gi.
ne; i++)
1769 const int* ev = gi.
edges[i];
1770 int n0 = elem->
node[ev[0]]->
id;
1771 int n1 = elem->
node[ev[1]]->
id;
1772 if ((n0 == v0 && n1 == v1) ||
1773 (n0 == v1 && n1 == v0)) {
return i; }
1775 MFEM_ABORT(
"Edge not found");
1781 for (
int i = 0; i < 6; i++)
1783 const int* fv = gi_hex.
faces[i];
1784 if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
1785 (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
1786 (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
1791 MFEM_ABORT(
"Face not found.");
1805 const int* fv = gi_hex.
faces[local];
1808 for (
int i = 0, j; i < 4; i++)
1810 for (j = 0; j < 4; j++)
1812 if (fv[i] == master[j])
1815 for (
int k = 0; k < mat.
Height(); k++)
1817 mat(k,i) = tmp(k,j);
1822 MFEM_ASSERT(j != 4,
"node not found.");
1833 Face* face =
faces.Peek(v0, v1, v2, v3);
1856 Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
1864 else if (split == 2)
1866 Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
1881 if (
Dim < 3) {
return; }
1888 MFEM_ASSERT(!elem->
ref_type,
"not a leaf element.");
1891 for (
int j = 0; j < gi.
nf; j++)
1895 for (
int k = 0; k < 4; k++)
1900 Face* face =
faces.Peek(node[0], node[1], node[2], node[3]);
1901 MFEM_ASSERT(face,
"face not found!");
1907 int index = face->
index;
1908 if (index >= processed_faces.
Size())
1910 processed_faces.
SetSize(index + 1, 0);
1912 else if (processed_faces[index])
1916 processed_faces[index] = 1;
1930 TraverseFace(node[0], node[1], node[2], node[3], pm, 0);
1939 for (
int i = sb; i < se; i++)
1955 if (!mid) {
return; }
1957 if (mid->
edge && level > 0)
1974 Face* face =
faces.Peek(v0, v0, v1, v1);
1975 MFEM_ASSERT(face != NULL,
"");
1982 double tmid = (t0 + t1) / 2;
1997 MFEM_ASSERT(!elem->
ref_type,
"not a leaf element.");
2000 for (
int j = 0; j < gi.
ne; j++)
2003 const int* ev = gi.
edges[j];
2004 Node* node[2] = { elem->
node[ev[0]], elem->
node[ev[1]] };
2006 Node* edge =
nodes.Peek(node[0], node[1]);
2007 MFEM_ASSERT(edge && edge->
edge,
"edge not found!");
2020 if (index >= processed_edges.
Size())
2022 processed_edges.
SetSize(index + 1, 0);
2024 else if (processed_edges[index])
2028 processed_edges[index] = 1;
2031 double t0 = 0.0, t1 = 1.0;
2045 for (
int i = sb; i < se; i++)
2061 oriented_matrix = point_matrix;
2065 MFEM_ASSERT(oriented_matrix.
Height() == 1 &&
2066 oriented_matrix.
Width() == 2,
"not an edge point matrix");
2070 oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
2071 oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
2075 std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
2102 indices.
Append(mid[0]->vertex->index);
2103 indices.
Append(mid[2]->vertex->index);
2110 indices.
Append(mid[1]->vertex->index);
2111 indices.
Append(mid[3]->vertex->index);
2122 int* I =
new int[nrows + 1];
2123 int** JJ =
new int*[nrows];
2132 MFEM_ASSERT(!elem->
ref_type,
"not a leaf element.");
2138 for (
int j = 0; j < gi.
nv; j++)
2140 indices.
Append(node[j]->vertex->index);
2143 for (
int j = 0; j < gi.
ne; j++)
2145 const int* ev = gi.
edges[j];
2151 for (
int j = 0; j < gi.
nf; j++)
2153 const int* fv = gi.
faces[j];
2155 node[fv[2]], node[fv[3]], indices);
2162 int size = indices.
Size();
2164 JJ[i] =
new int[size];
2165 memcpy(JJ[i], indices.
GetData(), size *
sizeof(int));
2170 for (
int i = 0; i < nrows; i++)
2179 int *J =
new int[nnz];
2181 for (
int i = 0; i < nrows; i++)
2183 int cnt = I[i+1] - I[i];
2184 memcpy(J+nnz, JJ[i], cnt *
sizeof(
int));
2213 MFEM_VERIFY(elem_set.
Size() == nleaves,
"");
2222 for (
int i = 0; i < nleaves; i++)
2228 for (
int j = 0; j < nv; j++)
2241 neighbor_set->
SetSize(nleaves);
2245 for (
int i = 0; i < nleaves; i++)
2251 for (
int j = 0; j < nv; j++)
2256 if (neighbor_set) { (*neighbor_set)[i] = 1; }
2264 static bool sorted_lists_intersect(
const int* a,
const int* b,
int na,
int nb)
2266 if (!na || !nb) {
return false; }
2267 int a_last = a[na-1], b_last = b[nb-1];
2268 if (*b < *a) {
goto l2; }
2270 if (a_last < *b) {
return false; }
2271 while (*a < *b) { a++; }
2272 if (*a == *b) {
return true; }
2274 if (b_last < *a) {
return false; }
2275 while (*b < *a) { b++; }
2276 if (*a == *b) {
return true; }
2300 while (stack.
Size())
2311 for (
int i = 0; i < 8; i++)
2326 for (
int i = 0; i < search_set->
Size(); i++)
2328 Element* testme = (*search_set)[i];
2334 if (sorted_lists_intersect(v1, v2, nv1, nv2))
2336 neighbors.
Append(testme);
2351 for (
int i = 0; i < elements.
Size(); i++)
2353 int index = elements[i]->index;
2356 for (
int j = 0; j < nv; j++)
2368 for (
int i = 0; i < search_set->
Size(); i++)
2370 Element* testme = (*search_set)[i];
2373 for (
int j = 0; j < nv; j++)
2390 for (
int i = 0; i < neighbors.
Size(); i++)
2392 elem_set[neighbors[i]->index] = 2;
2403 for (
int i = 0; i < np; i++)
2405 for (
int j = 0; j < points[0].dim; j++)
2407 point_matrix(j, i) = points[i].coord[j];
2419 Point(0, 0, 0),
Point(1, 0, 0),
Point(1, 1, 0),
Point(0, 1, 0),
2420 Point(0, 0, 1),
Point(1, 0, 1),
Point(1, 1, 1),
Point(0, 1, 1)
2431 MFEM_ABORT(
"unsupported geometry.");
throw;
2441 int ref_type = *ref_path++;
2442 int child = *ref_path++;
2448 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2449 Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
2454 pm(4), mid45, mid67, pm(7));
2456 else if (child == 1)
2459 mid45, pm(5), pm(6), mid67);
2462 else if (ref_type == 2)
2464 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2465 Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2470 pm(4), pm(5), mid56, mid74);
2472 else if (child == 1)
2475 mid74, mid56, pm(6), pm(7));
2478 else if (ref_type == 4)
2480 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2481 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2486 mid04, mid15, mid26, mid37);
2488 else if (child == 1)
2491 pm(4), pm(5), pm(6), pm(7));
2494 else if (ref_type == 3)
2496 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2497 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2498 Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2499 Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2501 Point midf0(mid23, mid12, mid01, mid30);
2502 Point midf5(mid45, mid56, mid67, mid74);
2507 pm(4), mid45, midf5, mid74);
2509 else if (child == 1)
2512 mid45, pm(5), mid56, midf5);
2514 else if (child == 2)
2517 midf5, mid56, pm(6), mid67);
2519 else if (child == 3)
2522 mid74, midf5, mid67, pm(7));
2525 else if (ref_type == 5)
2527 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2528 Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
2529 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2530 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2532 Point midf1(mid01, mid15, mid45, mid04);
2533 Point midf3(mid23, mid37, mid67, mid26);
2538 mid04, midf1, midf3, mid37);
2540 else if (child == 1)
2543 midf1, mid15, mid26, midf3);
2545 else if (child == 2)
2548 mid45, pm(5), pm(6), mid67);
2550 else if (child == 3)
2553 pm(4), mid45, mid67, pm(7));
2556 else if (ref_type == 6)
2558 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2559 Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2560 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2561 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2563 Point midf2(mid12, mid26, mid56, mid15);
2564 Point midf4(mid30, mid04, mid74, mid37);
2569 mid04, mid15, midf2, midf4);
2571 else if (child == 1)
2574 midf4, midf2, mid26, mid37);
2576 else if (child == 2)
2579 pm(4), pm(5), mid56, mid74);
2581 else if (child == 3)
2584 mid74, mid56, pm(6), pm(7));
2587 else if (ref_type == 7)
2589 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2590 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2591 Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2592 Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2593 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2594 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2596 Point midf0(mid23, mid12, mid01, mid30);
2597 Point midf1(mid01, mid15, mid45, mid04);
2598 Point midf2(mid12, mid26, mid56, mid15);
2599 Point midf3(mid23, mid37, mid67, mid26);
2600 Point midf4(mid30, mid04, mid74, mid37);
2601 Point midf5(mid45, mid56, mid67, mid74);
2603 Point midel(midf1, midf3);
2608 mid04, midf1, midel, midf4);
2610 else if (child == 1)
2613 midf1, mid15, midf2, midel);
2615 else if (child == 2)
2618 midel, midf2, mid26, midf3);
2620 else if (child == 3)
2623 midf4, midel, midf3, mid37);
2625 else if (child == 4)
2628 pm(4), mid45, midf5, mid74);
2630 else if (child == 5)
2633 mid45, pm(5), mid56, midf5);
2635 else if (child == 6)
2638 midf5, mid56, pm(6), mid67);
2640 else if (child == 7)
2643 mid74, midf5, mid67, pm(7));
2651 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2657 else if (child == 1)
2662 else if (ref_type == 2)
2664 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2670 else if (child == 1)
2675 else if (ref_type == 3)
2677 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2678 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2679 Point midel(mid01, mid23);
2685 else if (child == 1)
2689 else if (child == 2)
2693 else if (child == 3)
2701 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
2707 else if (child == 1)
2711 else if (child == 2)
2715 else if (child == 3)
2723 for (
int i = 0; i < pm.
np; i++)
2725 for (
int j = 0; j < pm(i).dim; j++)
2727 matrix(j, i) = pm(i).coord[j];
2751 int &matrix = map[ref_path];
2752 if (!matrix) { matrix = map.size(); }
2755 emb.
parent = coarse_index;
2760 ref_path.push_back(elem->
ref_type);
2761 ref_path.push_back(0);
2763 for (
int i = 0; i < 8; i++)
2767 ref_path[ref_path.length()-1] = i;
2771 ref_path.resize(ref_path.length()-2);
2778 "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
2785 std::string ref_path;
2786 ref_path.reserve(100);
2803 for (RefPathMap::iterator it = map.begin(); it != map.end(); ++it)
2815 "GetDerefinementTransforms() must be preceded by Derefine().");
2819 std::map<int, int> mat_no;
2828 int &matrix = mat_no[code];
2829 if (!matrix) { matrix = mat_no.size(); }
2840 std::map<int, int>::iterator it;
2841 for (it = mat_no.begin(); it != mat_no.end(); ++it)
2844 int code = it->first;
2845 path[0] = code >> 3;
2867 MFEM_ASSERT(node != NULL && node->
p1 != node->
p2,
"Invalid edge node.");
2872 if ((n2->
p1 != n2->
p2) && (n1->
id == n2->
p1 || n1->
id == n2->
p2))
2876 if (n2->
edge) {
return n2; }
2880 if ((n1->
p1 != n1->
p2) && (n2->
id == n1->
p1 || n2->
id == n1->
p2))
2884 if (n1->
edge) {
return n1; }
2894 MFEM_ASSERT(node->
edge != NULL,
"(v1, v2) is not an edge.");
2897 return master ? master->
edge->
index : -1;
2918 if (!elem) { elem = face->
elem[1]; }
2919 MFEM_ASSERT(elem,
"Face has no elements?");
2926 for (
int i = 0; i < 4; i++)
2928 node[i] = elem->
node[fv[i]];
2944 if (bdr_attr_is_ess[face->
attribute - 1])
2949 for (
int j = 0; j < 4; j++)
2951 bdr_vertices.
Append(node[j]->vertex->index);
2953 Node* edge =
nodes.Peek(node[j], node[(j+1) % 4]);
2954 MFEM_ASSERT(edge && edge->
edge,
"Edge not found.");
2981 bdr_vertices.
Sort();
2991 if (!mid || !mid->
vertex) {
return 0; }
2996 int& h_level,
int& v_level)
const
2998 int hl1, hl2, vl1, vl2;
3004 h_level = v_level = 0;
3010 h_level = std::max(hl1, hl2);
3011 v_level = std::max(vl1, vl2) + 1;
3017 h_level = std::max(hl1, hl2) + 1;
3018 v_level = std::max(vl1, vl2);
3022 static int max8(
int a,
int b,
int c,
int d,
int e,
int f,
int g,
int h)
3024 return std::max(std::max(std::max(a, b), std::max(c, d)),
3025 std::max(std::max(e, f), std::max(g, h)));
3034 for (
int i = 0; i < gi.
ne; i++)
3036 const int* ev = gi.
edges[i];
3043 for (
int i = 0; i < gi.
nf; i++)
3045 const int* fv = gi.
faces[i];
3046 FaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]],
3047 flevel[i][1], flevel[i][0]);
3050 splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
3051 elevel[0], elevel[2], elevel[4], elevel[6]);
3053 splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
3054 elevel[1], elevel[3], elevel[5], elevel[7]);
3056 splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
3057 elevel[8], elevel[9], elevel[10], elevel[11]);
3061 splits[0] = std::max(elevel[0], elevel[2]);
3062 splits[1] = std::max(elevel[1], elevel[3]);
3066 splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
3067 splits[1] = splits[0];
3071 MFEM_ABORT(
"Unsupported element geometry.");
3085 for (
int k = 0; k <
Dim; k++)
3087 if (splits[k] > max_level)
3089 ref_type |= (1 << k);
3107 MFEM_VERIFY(max_nc_level >= 1,
"'max_nc_level' must be 1 or greater.");
3114 if (!refinements.
Size()) {
break; }
3126 if (it->vertex && it->p1 != it->p2) { nv++; }
3133 if (it->vertex && it->p1 != it->p2)
3138 MFEM_ASSERT(p1 && p1->
vertex,
"");
3139 MFEM_ASSERT(p2 && p2->
vertex,
"");
3141 out << it->vertex->index <<
" "
3155 input >>
id >> p1 >> p2;
3156 MFEM_VERIFY(input,
"problem reading vertex parents.");
3159 MFEM_VERIFY(node,
"vertex " <<
id <<
" not found.");
3160 MFEM_VERIFY(
nodes.Peek(p1),
"parent " << p1 <<
" not found.");
3161 MFEM_VERIFY(
nodes.Peek(p2),
"parent " << p2 <<
" not found.");
3164 nodes.Reparent(node, p1, p2);
3170 for (
int i = 0; i < vertices.
Size(); i++)
3173 MFEM_ASSERT(node && node->
vertex,
"");
3175 const double* pos = vertices[i]();
3180 static int ref_type_num_children[8] = {0, 2, 2, 4, 2, 4, 4, 8 };
3183 int &coarse_id)
const
3187 int child_id[8], nch = 0;
3188 for (
int i = 0; i < 8; i++)
3195 MFEM_ASSERT(nch == ref_type_num_children[(
int) elem->
ref_type],
"");
3198 for (
int i = 0; i < nch; i++)
3200 out <<
" " << child_id[i];
3237 int nleaf = leaves.
Size();
3249 if (
Dim == 3 && ref_type != 7) { iso =
false; }
3252 int nch = ref_type_num_children[ref_type];
3253 for (
int i = 0,
id; i < nch; i++)
3256 MFEM_VERIFY(
id >= 0,
"");
3257 MFEM_VERIFY(
id < nleaf ||
id - nleaf < coarse.
Size(),
3258 "coarse element cannot be referenced before it is "
3259 "defined (id=" <<
id <<
").");
3261 Element* &child = (
id < nleaf) ? leaves[
id] : coarse[
id - nleaf];
3263 MFEM_VERIFY(child,
"element " <<
id <<
" cannot have two parents.");
3264 elem->
child[i] = child;
3281 for (
int i = 0; i < coarse.
Size(); i++)
3285 for (
int i = 0; i < leaves.
Size(); i++)
3299 for (
int i = 0; i < 8; i++)
3313 pmsize = pm.
Width() * pm.
Height() *
sizeof(double);
3316 return conforming.capacity() *
sizeof(
MeshId) +
3317 masters.capacity() *
sizeof(
Master) +
3318 slaves.capacity() *
sizeof(
Slave) +
3319 slaves.size() * pmsize;
3324 nelem = nvert = nedges = 0;
3331 if (it->vertex) { nvert++; }
3332 if (it->edge) { nedges++; }
3338 int nelem, nvert, nedges;
3341 return nelem *
sizeof(
Element) +
3343 nedges *
sizeof(
Edge) +
3344 nodes.MemoryUsage() +
3345 faces.MemoryUsage() +
3361 int nelem, nvert, nedges;
3364 std::cout << nelem *
sizeof(
Element) <<
" elements\n"
3365 << nvert *
sizeof(
Vertex) <<
" vertices\n"
3366 << nedges *
sizeof(
Edge) <<
" edges\n";
3368 nodes.PrintMemoryDetail(); std::cout <<
" nodes\n";
3369 faces.PrintMemoryDetail(); std::cout <<
" faces\n";
3371 std::cout <<
root_elements.MemoryUsage() <<
" root_elements\n"
3379 <<
ref_stack.MemoryUsage() <<
" ref_stack\n"
3381 <<
sizeof(*this) <<
" NCMesh" << std::endl;
3390 for (
int j = 0; j <
Dim; j++)
3394 for (
int k = 0; k < 8; k++)
3402 std::cout << sum / count <<
" ";
NCList face_list
lazy-initialized list of faces, see GetFaceList
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
NCMesh(const Mesh *mesh, std::istream *vertex_parents=NULL)
int Size() const
Logical size of the array.
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)
CoarseFineTransformations transforms
storage for data returned by Get[De]RefinementTransforms()
void FindNeighbors(const Element *elem, Array< Element * > &neighbors, const Array< Element * > *search_set=NULL)
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
const CoarseFineTransformations & GetDerefinementTransforms()
int index
edge number in the Mesh
Array< Element * > coarse_elements
state of leaf_elements before Refine(), set by MarkCoarseLevel()
virtual int GetNEdges() const =0
Iterator over items contained in the HashTable.
int PrintElements(std::ostream &out, Element *elem, int &coarse_id) const
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
void PrintVertexParents(std::ostream &out) const
I/O: Print the "vertex_parents" section of the mesh file (ver. >= 1.1).
void GetPointMatrix(int geom, const char *ref_path, DenseMatrix &matrix)
Node * PeekAltParents(Node *v1, Node *v2)
virtual void LimitNCLevel(int max_nc_level)
void PrintMemoryDetail() const
int GetNBE() const
Returns number of boundary elements.
int index
element number in the Mesh, -1 if refined
virtual int GetNumGhosts() const
void DebugLeafOrder() const
Print the space-filling curve formed by the sequence of leaf elements.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
T * GetData()
Returns the data.
int attribute
boundary element attribute, -1 if internal edge (2D)
void ClearTransforms()
Free all internal data created by the above three functions.
Table derefinements
possible derefinements, see GetDerefinementTable
Data type dense matrix using column-major storage.
void InitDerefTransforms()
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
void RefElementNodes(Element *elem)
char geom
Geometry::Type of the element.
static int find_node(Element *elem, Node *node)
void SetDerefMatrixCodes(Element *parent, Array< Element * > &coarse)
int FaceSplitType(Node *v1, Node *v2, Node *v3, Node *v4, Node *mid[4]=NULL) const
void SetSize(int i, int j, int k)
int GetNE() const
Returns number of elements.
virtual void ElementSharesFace(Element *elem, Face *face)
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
const Element * GetFace(int i) const
static PointMatrix pm_tri_identity
int index
vertex number in the Mesh
void NeighborExpand(const Array< Element * > &elements, Array< Element * > &expanded, const Array< Element * > *search_set=NULL)
void CountObjects(int &nelem, int &nvert, int &nedges) const
std::vector< MeshId > conforming
Data type quadrilateral element.
void GetLimitRefinements(Array< Refinement > &refinements, int max_level)
void Delete(ItemT *item)
Remove an item from the hash table and also delete the item itself.
Element * parent
parent element, NULL if this is a root element
Vertex * NewVertex(Node *v1, Node *v2)
bool NodeSetY1(Node *node, Node **n)
int spaceDim
dimensions of the elements and the vertex coordinates
Array< int > vertex_nodeId
int attribute
boundary element attribute, -1 if internal face
int GetGeometryType() const
int index
Mesh element number.
Element * element
NCMesh::Element containing this vertex/edge/face.
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
virtual void OnMeshUpdated(Mesh *mesh)
int Width() const
Returns the number of TYPE II elements (after Finalize() is called).
Element * NewQuadrilateral(Node *n0, Node *n1, Node *n2, Node *n3, int attr, int eattr0, int eattr1, int eattr2, int eattr3)
void DebugNeighbors(Array< char > &elem_set)
void CollectEdgeVertices(Node *v0, Node *v1, Array< int > &indices)
virtual bool IsGhost(const Element *elem) const
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
Data type hexahedron element.
virtual const int * GetEdgeVertices(int) const =0
int num_vertices
width of the table
int index
face number in the Mesh
const Table & GetDerefinementTable()
void MakeFromList(int nrows, const Array< Connection > &list)
Table element_vertex
leaf-element to vertex table, see FindSetNeighbors
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
int local
local number within 'element'
int Append(const T &el)
Append element to array, resize if necessary.
Element * elem[2]
up to 2 elements sharing the face
static void UnrefEdge(Node *node, HashTable< Node > &nodes)
void GetMatrix(DenseMatrix &point_matrix) const
void CheckAnisoFace(Node *v1, Node *v2, Node *v3, Node *v4, Node *mid12, Node *mid34, int level=0)
virtual void Derefine(const Array< int > &derefs)
void RegisterFaces(Element *elem, int *fattr=NULL)
std::vector< Master > masters
virtual void UpdateVertices()
update Vertex::index and vertex_nodeId
void RegisterElement(Element *e)
int ReorderFacePointMat(Node *v0, Node *v1, Node *v2, Node *v3, Element *elem, DenseMatrix &mat) const
Element * CopyHierarchy(Element *elem)
void RefineElement(Element *elem, char ref_type)
void FaceSplitLevel(Node *v1, Node *v2, Node *v3, Node *v4, int &h_level, int &v_level) const
Element(int geom, int attr)
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
virtual void BuildEdgeList()
const CoarseFineTransformations & GetRefinementTransforms()
int GetAttribute() const
Return element's attribute.
void LoadCoarseElements(std::istream &input)
I/O: Load the element refinement hierarchy from a mesh file.
Identifies a vertex/edge/face in both Mesh and NCMesh.
Data type triangle element.
const Element * GetElement(int i) const
virtual void ElementSharesEdge(Element *elem, Edge *edge)
void Sort()
Sorts the array. This requires operator< to be defined for T.
int Size() const
Returns the number of TYPE I elements.
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
bool NodeSetX1(Node *node, Node **n)
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
static GeomInfo GI[Geometry::NumGeom]
void DeleteHierarchy(Element *elem)
int SpaceDimension() const
bool Iso
true if the mesh only contains isotropic refinements
void TraverseEdge(Node *v0, Node *v1, double t0, double t1, int flags, int level)
void DerefineElement(Element *elem)
std::map< std::string, int > RefPathMap
Array< Node * > boundary_edges
subset of all edges, set by BuildEdgeList
void FindSetNeighbors(const Array< char > &elem_set, Array< Element * > *neighbors, Array< char > *neighbor_set=NULL)
std::vector< Slave > slaves
virtual void BuildFaceList()
void ForceRefinement(Node *v1, Node *v2, Node *v3, Node *v4)
void SetVertexPositions(const Array< mfem::Vertex > &vertices)
I/O: Set positions of all vertices (used by mesh loader).
int edge_flags
edge orientation flags
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
static const PointMatrix & GetGeomIdentity(int geom)
Element * child[8]
2-8 children (if ref_type != 0)
Array< Element * > leaf_elements
Face * GetFace(Element *elem, int face_no)
mfem::Element * NewMeshElement(int geom) const
void DeleteLast()
Delete the last entry.
Helper struct for defining a connectivity table, see Table::MakeFromList.
void UpdateLeafElements()
bool NodeSetZ1(Node *node, Node **n)
Node * GetMidEdgeVertexSimple(Node *v1, Node *v2)
void BuildElementToVertexTable()
void UnrefElementNodes(Element *elem)
virtual int GetNFaces(int &nFaceVertices) const =0
static PointMatrix pm_quad_identity
void CheckIsoFace(Node *v1, Node *v2, Node *v3, Node *v4, Node *e1, Node *e2, Node *e3, Node *e4, Node *midf)
void SetIJ(int *newI, int *newJ, int newsize=-1)
Node * GetMidFaceVertex(Node *e1, Node *e2, Node *e3, Node *e4)
void CollectFaceVertices(Node *v0, Node *v1, Node *v2, Node *v3, Array< int > &indices)
bool NodeSetZ2(Node *node, Node **n)
void Initialize(const mfem::Element *elem)
bool NodeSetX2(Node *node, Node **n)
void GetMeshComponents(Array< mfem::Vertex > &vertices, Array< mfem::Element * > &elements, Array< mfem::Element * > &boundary) const
Return the basic Mesh arrays for the current finest level.
Node * node[8]
element corners (if ref_type == 0)
void CollectLeafElements(Element *elem, int state)
Element * NewTriangle(Node *n0, Node *n1, Node *n2, int attr, int eattr0, int eattr1, int eattr2)
virtual int GetNVertices() const =0
int parent
element index in the coarse mesh
void SetAttribute(const int attr)
Set element's attribute.
virtual void AssignLeafIndices()
static PointMatrix pm_hex_identity
static int find_hex_face(int a, int b, int c)
virtual const int * GetFaceVertices(int fi) const =0
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
Array< ElemRefType > ref_stack
stack of scheduled refinements (temporary)
T & Last()
Return the last element in the array.
static void find_face_nodes(const Face *face, Node *node[4])
Array< Element * > root_elements
long MemoryUsage() const
Return total number of bytes allocated.
Node * GetMidEdgeVertex(Node *v1, Node *v2)
int EdgeSplitLevel(Node *v1, Node *v2) const
Element * GetSingleElement() const
Return one of elem[0] or elem[1] and make sure the other is NULL.
void PrintCoarseElements(std::ostream &out) const
I/O: Print the "coarse_elements" section of the mesh file (ver. >= 1.1).
static int find_element_edge(Element *elem, int v0, int v1)
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
void TraverseFace(Node *v0, Node *v1, Node *v2, Node *v3, const PointMatrix &pm, int level)
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void UpdateElementToVertexTable()
int GetElementDepth(int i) const
Return the distance of leaf 'i' from the root.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void CountSplits(Element *elem, int splits[3]) const
void ForgetElement(Element *e)
char ref_type
bit mask of X,Y,Z refinements (bits 0,1,2 respectively)
bool NodeSetY2(Node *node, Node **n)
void LoadVertexParents(std::istream &input)
int CountElements(Element *elem) const
Abstract data type element.
Data type line segment element.
Array< Face * > boundary_faces
subset of all faces, set by BuildFaceList
void TraverseRefinements(Element *elem, int coarse_index, std::string &ref_path, RefPathMap &map)
int rank
processor number (ParNCMesh)
void CollectDerefinements(Element *elem, Array< Connection > &list)
const Element * GetBdrElement(int i) const
virtual int GetType() const =0
Returns element's type.
DenseMatrix point_matrix
position within the master edge/face
Defines the position of a fine element within a coarse element.
virtual void Refine(const Array< Refinement > &refinements)
int matrix
index into CoarseFineTransformations::point_matrices
int GetEdgeMaster(int v1, int v2) const
void UnrefVertex(HashTable< Node > &nodes)