13 #include "../fem/fem.hpp"
14 #include "../general/sort_pairs.hpp"
28 std::map<Geometry::Type, DenseTensor>::const_iterator pm_it;
31 "cannot find point matrices for geometry type \"" << geom
47 :
geom(g), num_children(n), children(c) { }
49 bool operator<(
const RefType &other)
const
51 if (
geom < other.geom) {
return true; }
52 if (
geom > other.geom) {
return false; }
53 if (num_children < other.num_children) {
return true; }
54 if (num_children > other.num_children) {
return false; }
55 for (
int i = 0; i < num_children; i++)
57 if (children[i].one < other.children[i].one) {
return true; }
58 if (children[i].one > other.children[i].one) {
return false; }
73 for (
int i = 0; i < fine_ne; i++)
75 coarse_ne = std::max(coarse_ne,
embeddings[i].parent);
79 coarse_to_ref_type.
SetSize(coarse_ne);
80 coarse_to_fine.
SetDims(coarse_ne, fine_ne);
85 for (
int i = 0; i < fine_ne; i++)
90 MFEM_ASSERT(cf_i.Last() == cf_j.
Size(),
"internal error");
91 for (
int i = 0; i < fine_ne; i++)
95 cf_j[cf_i[e.
parent]].two = i;
98 std::copy_backward(cf_i.begin(), cf_i.end()-1, cf_i.end());
100 for (
int i = 0; i < coarse_ne; i++)
102 std::sort(&cf_j[cf_i[i]], cf_j.
GetData() + cf_i[i+1]);
104 for (
int i = 0; i < fine_ne; i++)
106 coarse_to_fine.
GetJ()[i] = cf_j[i].two;
109 using internal::RefType;
113 map<RefType,int> ref_type_map;
114 for (
int i = 0; i < coarse_ne; i++)
116 const int num_children = cf_i[i+1]-cf_i[i];
117 MFEM_ASSERT(num_children > 0,
"");
118 const int fine_el = cf_j[cf_i[i]].two;
121 const RefType ref_type(geom, num_children, &cf_j[cf_i[i]]);
122 pair<map<RefType,int>::iterator,
bool> res =
124 pair<const RefType,int>(ref_type, (
int)ref_type_map.size()));
125 coarse_to_ref_type[i] = res.first->second;
127 ref_type_to_matrix.
MakeI((
int)ref_type_map.size());
128 ref_type_to_geom.
SetSize((
int)ref_type_map.size());
129 for (map<RefType,int>::iterator it = ref_type_map.begin();
130 it != ref_type_map.end(); ++it)
132 ref_type_to_matrix.
AddColumnsInRow(it->second, it->first.num_children);
133 ref_type_to_geom[it->second] = it->first.geom;
135 ref_type_to_matrix.
MakeJ();
136 for (map<RefType,int>::iterator it = ref_type_map.begin();
137 it != ref_type_map.end(); ++it)
139 const RefType &rt = it->first;
140 for (
int j = 0; j < rt.num_children; j++)
142 ref_type_to_matrix.
AddConnection(it->second, rt.children[j].one);
162 for (
int i = 0; i <
ne; i++)
164 for (
int j = 0; j < 2; j++)
169 for (
int i = 0; i <
nf; i++)
171 for (
int j = 0; j <
nfv; j++)
180 for (
int i = 0; i <
ne; i++)
199 Iso = vertex_parents ?
false :
true;
204 for (
int i = 0; i < mesh->
GetNE(); i++)
208 for (
int j = 0; j < nv; j++)
210 max_id = std::max(max_id, v[j]);
213 for (
int id = 0;
id <= max_id;
id++)
216 int node =
nodes.GetId(
id,
id);
217 MFEM_CONTRACT_VAR(node);
218 MFEM_ASSERT(node ==
id,
"");
230 for (
int i = 0; i < mesh->
GetNV(); i++)
237 for (
int i = 0; i < mesh->
GetNE(); i++)
246 MFEM_ABORT(
"only triangles, quads and hexes are supported by NCMesh.");
254 MFEM_ASSERT(root_id == i,
"");
258 for (
int j = 0; j <
GI[
geom].
nv; j++)
260 root_elem.
node[j] = v[j];
272 for (
int i = 0; i < mesh->
GetNBE(); i++)
279 Face* face =
faces.Find(v[0], v[1], v[2], v[3]);
280 MFEM_VERIFY(face,
"boundary face not found.");
285 Face* face =
faces.Find(v[0], v[0], v[1], v[1]);
286 MFEM_VERIFY(face,
"boundary face not found.");
291 MFEM_ABORT(
"only segment and quadrilateral boundary "
292 "elements are supported by NCMesh.");
335 for (
int i = 0; i <
elements.Size(); i++)
357 "vert_refc: " << (
int)
vert_refc <<
", edge_refc: "
368 for (
int i = 0; i < gi.
nv; i++)
370 nodes[node[i]].vert_refc++;
374 for (
int i = 0; i < gi.
ne; i++)
376 const int* ev = gi.
edges[i];
377 nodes.Get(node[ev[0]], node[ev[1]])->edge_refc++;
381 for (
int i = 0; i < gi.
nf; i++)
383 const int* fv = gi.
faces[i];
384 faces.GetId(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
399 for (
int i = 0; i < gi.
nf; i++)
401 const int* fv = gi.
faces[i];
402 int face =
faces.FindId(node[fv[0]], node[fv[1]],
403 node[fv[2]], node[fv[3]]);
404 MFEM_ASSERT(face >= 0,
"face not found.");
405 faces[face].ForgetElement(elem);
413 for (
int i = 0; i < gi.
ne; i++)
415 const int* ev = gi.
edges[i];
417 MFEM_ASSERT(enode >= 0,
"edge not found.");
418 MFEM_ASSERT(
nodes.IdExists(enode),
"edge does not exist.");
419 if (!
nodes[enode].UnrefEdge())
426 for (
int i = 0; i < gi.
nv; i++)
428 if (!
nodes[node[i]].UnrefVertex())
430 nodes.Delete(node[i]);
440 for (
int i = 0; i < gi.
nf; i++)
443 MFEM_ASSERT(face,
"face not found.");
445 if (fattr) { face->
attribute = fattr[i]; }
451 for (
int i = 0; i < elemFaces.
Size(); i++)
453 if (
faces[elemFaces[i]].Unused())
455 faces.Delete(elemFaces[i]);
462 if (elem[0] < 0) { elem[0] = e; }
463 else if (elem[1] < 0) { elem[1] = e; }
464 else { MFEM_ABORT(
"can't have 3 elements in Face::elem[]."); }
469 if (elem[0] == e) { elem[0] = -1; }
470 else if (elem[1] == e) { elem[1] = -1; }
471 else { MFEM_ABORT(
"element " << e <<
" not found in Face::elem[]."); }
477 const int* fv = gi.
faces[face_no];
478 int* node = elem.
node;
479 return faces.Find(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]]);
486 MFEM_ASSERT(elem[1] < 0,
"not a single element face.");
491 MFEM_ASSERT(elem[1] >= 0,
"no elements in face.");
498 int mid =
nodes.FindId(node1, node2);
499 if (mid < 0 && Dim >= 3 && !
Iso)
526 int n1p1 = n1.
p1, n1p2 = n1.
p2;
527 int n2p1 = n2.p1, n2p2 = n2.p2;
529 if ((n1p1 != n1p2) && (n2p1 != n2p2))
534 if (a1 < 0 || a2 < 0)
541 if (a1 >= 0 && a2 >= 0)
543 mid =
nodes.FindId(a1, a2);
554 : geom(geom), ref_type(0), flag(0), index(-1), rank(0), attribute(attr)
557 for (
int i = 0; i < 8; i++) {
node[i] = -1; }
566 int n4,
int n5,
int n6,
int n7,
568 int fattr0,
int fattr1,
int fattr2,
569 int fattr3,
int fattr4,
int fattr5)
596 int eattr0,
int eattr1,
int eattr2,
int eattr3)
620 int attr,
int eattr0,
int eattr1,
int eattr2)
647 if (mid < 0) { mid =
nodes.GetId(vn1, vn2); }
654 int midf =
nodes.FindId(en1, en3);
655 if (midf >= 0) {
return midf; }
656 return nodes.GetId(en2, en4);
661 {
return node == n[0] || node == n[3] || node == n[4] || node == n[7]; }
664 {
return node == n[1] || node == n[2] || node == n[5] || node == n[6]; }
667 {
return node == n[0] || node == n[1] || node == n[4] || node == n[5]; }
670 {
return node == n[2] || node == n[3] || node == n[6] || node == n[7]; }
673 {
return node == n[0] || node == n[1] || node == n[2] || node == n[3]; }
676 {
return node == n[4] || node == n[5] || node == n[6] || node == n[7]; }
682 Face* face =
faces.Find(vn1, vn2, vn3, vn4);
683 if (!face) {
return; }
686 MFEM_ASSERT(!
elements[elem].ref_type,
"element already refined.");
708 MFEM_ABORT(
"inconsistent element/face structure.");
714 int mid12,
int mid34,
int level)
740 int mid23 =
nodes.FindId(vn2, vn3);
741 int mid41 =
nodes.FindId(vn4, vn1);
742 if (mid23 >= 0 && mid41 >= 0)
744 int midf =
nodes.FindId(mid23, mid41);
747 nodes.Reparent(midf, mid12, mid34);
780 int en1,
int en2,
int en3,
int en4,
int midf)
798 if (!ref_type) {
return; }
804 char remaining = ref_type & ~el.
ref_type;
807 for (
int i = 0; i < 8; i++)
818 for (
int i = 0; i < 8; i++) { child[i] = -1; }
823 for (
int i = 0; i < gi.
nf; i++)
825 const int* fv = gi.
faces[i];
826 Face* face =
faces.Find(no[fv[0]], no[fv[1]], no[fv[2]], no[fv[3]]);
855 no[4], mid45, mid67, no[7], attr,
856 fa[0], fa[1], -1, fa[3], fa[4], fa[5]);
859 mid45, no[5], no[6], mid67, attr,
860 fa[0], fa[1], fa[2], fa[3], -1, fa[5]);
867 else if (ref_type == 2)
875 no[4], no[5], mid56, mid74, attr,
876 fa[0], fa[1], fa[2], -1, fa[4], fa[5]);
879 mid74, mid56, no[6], no[7], attr,
880 fa[0], -1, fa[2], fa[3], fa[4], fa[5]);
887 else if (ref_type == 4)
895 mid04, mid15, mid26, mid37, attr,
896 fa[0], fa[1], fa[2], fa[3], fa[4], -1);
899 no[4], no[5], no[6], no[7], attr,
900 -1, fa[1], fa[2], fa[3], fa[4], fa[5]);
907 else if (ref_type == 3)
923 no[4], mid45, midf5, mid74, attr,
924 fa[0], fa[1], -1, -1, fa[4], fa[5]);
927 mid45, no[5], mid56, midf5, attr,
928 fa[0], fa[1], fa[2], -1, -1, fa[5]);
931 midf5, mid56, no[6], mid67, attr,
932 fa[0], -1, fa[2], fa[3], -1, fa[5]);
935 mid74, midf5, mid67, no[7], attr,
936 fa[0], -1, -1, fa[3], fa[4], fa[5]);
943 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
944 CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
946 else if (ref_type == 5)
962 mid04, midf1, midf3, mid37, attr,
963 fa[0], fa[1], -1, fa[3], fa[4], -1);
966 midf1, mid15, mid26, midf3, attr,
967 fa[0], fa[1], fa[2], fa[3], -1, -1);
970 mid45, no[5], no[6], mid67, attr,
971 -1, fa[1], fa[2], fa[3], -1, fa[5]);
974 no[4], mid45, mid67, no[7], attr,
975 -1, fa[1], -1, fa[3], fa[4], fa[5]);
982 CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
983 CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
985 else if (ref_type == 6)
1001 mid04, mid15, midf2, midf4, attr,
1002 fa[0], fa[1], fa[2], -1, fa[4], -1);
1005 midf4, midf2, mid26, mid37, attr,
1006 fa[0], -1, fa[2], fa[3], fa[4], -1);
1009 no[4], no[5], mid56, mid74, attr,
1010 -1, fa[1], fa[2], -1, fa[4], fa[5]);
1013 mid74, mid56, no[6], no[7], attr,
1014 -1, -1, fa[2], fa[3], fa[4], fa[5]);
1021 CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1022 CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1024 else if (ref_type == 7)
1051 mid04, midf1, midel, midf4, attr,
1052 fa[0], fa[1], -1, -1, fa[4], -1);
1055 midf1, mid15, midf2, midel, attr,
1056 fa[0], fa[1], fa[2], -1, -1, -1);
1059 midel, midf2, mid26, midf3, attr,
1060 fa[0], -1, fa[2], fa[3], -1, -1);
1063 midf4, midel, midf3, mid37, attr,
1064 fa[0], -1, -1, fa[3], fa[4], -1);
1067 no[4], mid45, midf5, mid74, attr,
1068 -1, fa[1], -1, -1, fa[4], fa[5]);
1071 mid45, no[5], mid56, midf5, attr,
1072 -1, fa[1], fa[2], -1, -1, fa[5]);
1075 midf5, mid56, no[6], mid67, attr,
1076 -1, -1, fa[2], fa[3], -1, fa[5]);
1079 mid74, midf5, mid67, no[7], attr,
1080 -1, -1, -1, fa[3], fa[4], fa[5]);
1082 CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid30, midf0);
1083 CheckIsoFace(no[0], no[1], no[5], no[4], mid01, mid15, mid45, mid04, midf1);
1084 CheckIsoFace(no[1], no[2], no[6], no[5], mid12, mid26, mid56, mid15, midf2);
1085 CheckIsoFace(no[2], no[3], no[7], no[6], mid23, mid37, mid67, mid26, midf3);
1086 CheckIsoFace(no[3], no[0], no[4], no[7], mid30, mid04, mid74, mid37, midf4);
1087 CheckIsoFace(no[4], no[5], no[6], no[7], mid45, mid56, mid67, mid74, midf5);
1091 MFEM_ABORT(
"invalid refinement type.");
1094 if (ref_type != 7) {
Iso =
false; }
1102 int mid01 =
nodes.GetId(no[0], no[1]);
1103 int mid23 =
nodes.GetId(no[2], no[3]);
1106 attr, fa[0], -1, fa[2], fa[3]);
1109 attr, fa[0], fa[1], fa[2], -1);
1111 else if (ref_type == 2)
1113 int mid12 =
nodes.GetId(no[1], no[2]);
1114 int mid30 =
nodes.GetId(no[3], no[0]);
1117 attr, fa[0], fa[1], -1, fa[3]);
1120 attr, -1, fa[1], fa[2], fa[3]);
1122 else if (ref_type == 3)
1124 int mid01 =
nodes.GetId(no[0], no[1]);
1125 int mid12 =
nodes.GetId(no[1], no[2]);
1126 int mid23 =
nodes.GetId(no[2], no[3]);
1127 int mid30 =
nodes.GetId(no[3], no[0]);
1129 int midel =
nodes.GetId(mid01, mid23);
1132 attr, fa[0], -1, -1, fa[3]);
1135 attr, fa[0], fa[1], -1, -1);
1138 attr, -1, fa[1], fa[2], -1);
1141 attr, -1, -1, fa[2], fa[3]);
1145 MFEM_ABORT(
"Invalid refinement type.");
1148 if (ref_type != 3) {
Iso =
false; }
1155 int mid01 =
nodes.GetId(no[0], no[1]);
1156 int mid12 =
nodes.GetId(no[1], no[2]);
1157 int mid20 =
nodes.GetId(no[2], no[0]);
1159 child[0] =
NewTriangle(no[0], mid01, mid20, attr, fa[0], -1, fa[2]);
1160 child[1] =
NewTriangle(mid01, no[1], mid12, attr, fa[0], fa[1], -1);
1161 child[2] =
NewTriangle(mid20, mid12, no[2], attr, -1, fa[1], fa[2]);
1162 child[3] =
NewTriangle(mid01, mid12, mid20, attr, -1, -1, -1);
1166 MFEM_ABORT(
"Unsupported element geometry.");
1170 for (
int i = 0; i < 8 && child[i] >= 0; i++)
1183 for (
int i = 0; i < 8 && child[i] >= 0; i++)
1192 for (
int i = 0; i < 8 && child[i] >= 0; i++)
1209 for (
int i = refinements.
Size()-1; i >= 0; i--)
1237 #if defined(MFEM_DEBUG) && !defined(MFEM_USE_MPI)
1238 mfem::out <<
"Refined " << refinements.
Size() <<
" + " << nforced
1239 <<
" elements" << std::endl;
1249 static int quad_deref_table[3][4 + 4] =
1251 { 0, 1, 1, 0, 1, 1, 0, 0 },
1252 { 0, 0, 1, 1, 0, 0, 1, 1 },
1253 { 0, 1, 2, 3, 1, 1, 3, 3 }
1255 static int hex_deref_table[7][8 + 6] =
1257 { 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0 },
1258 { 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1 },
1259 { 0, 1, 2, 3, 0, 1, 2, 3, 1, 1, 1, 3, 3, 3 },
1260 { 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1 },
1261 { 0, 1, 1, 0, 3, 2, 2, 3, 1, 1, 1, 3, 3, 3 },
1262 { 0, 0, 1, 1, 2, 2, 3, 3, 0, 0, 0, 3, 3, 3 },
1263 { 0, 1, 2, 3, 4, 5, 6, 7, 1, 1, 1, 7, 7, 7 }
1285 ch = el.
child[index];
1290 MFEM_ABORT(
"Unsupported element geometry.");
1302 memcpy(child, el.
child,
sizeof(child));
1305 for (
int i = 0; i < 8 && child[i] >= 0; i++)
1317 for (
int i = 0; i < 8; i++)
1322 for (
int i = 0; i < 6; i++)
1327 ch.
node[fv[2]], ch.
node[fv[3]])->attribute;
1332 for (
int i = 0; i < 4; i++)
1337 for (
int i = 0; i < 4; i++)
1342 ch.
node[fv[2]], ch.
node[fv[3]])->attribute;
1347 for (
int i = 0; i < 3; i++)
1353 ch.
node[fv[2]], ch.
node[fv[3]])->attribute;
1358 MFEM_ABORT(
"Unsupported element geometry.");
1370 for (
int i = 0; i < 8 && child[i] >= 0; i++)
1393 int total = 0, ref = 0, ghost = 0;
1394 for (
int i = 0; i < 8 && el.
child[i] >= 0; i++)
1402 if (!ref && ghost < total)
1405 int next_row = list.
Size() ? (list.
Last().from + 1) : 0;
1406 for (
int i = 0; i < 8 && el.
child[i] >= 0; i++)
1414 for (
int i = 0; i < 8 && el.
child[i] >= 0; i++)
1431 int size = list.
Size() ? (list.
Last().from + 1) : 0;
1440 for (
int i = 0; i < deref_table.
Size(); i++)
1442 const int* fine = deref_table.
GetRow(i), size = deref_table.
RowSize(i);
1446 for (
int j = 0; j < size; j++)
1451 for (
int k = 0; k <
Dim; k++)
1453 if ((parent.
ref_type & (1 << k)) &&
1454 splits[k] >= max_nc_level)
1467 MFEM_VERIFY(
Dim < 3 ||
Iso,
1468 "derefinement of 3D anisotropic meshes not implemented yet.");
1476 for (
int i = 0; i < derefs.
Size(); i++)
1478 int row = derefs[i];
1480 "invalid derefinement number.");
1495 for (
int i = 0; i < fine_coarse.
Size(); i++)
1506 for (
int i = 0; i < nfine; i++)
1513 std::map<Geometry::Type, DenseTensor>::iterator it;
1517 it->second.SetSize(0, 0, 0);
1525 for (
int i = 0; i < 8 && prn.
child[i] >= 0; i++)
1530 int code = (prn.
ref_type << 3) + i;
1532 fine_coarse[ch.
index] = parent;
1546 if (node->HasVertex()) { node->vert_index =
NVertices++; }
1558 static char quad_hilbert_child_order[8][4] =
1560 {0,1,2,3}, {0,3,2,1}, {1,2,3,0}, {1,0,3,2},
1561 {2,3,0,1}, {2,1,0,3}, {3,0,1,2}, {3,2,1,0}
1563 static char quad_hilbert_child_state[8][4] =
1565 {1,0,0,5}, {0,1,1,4}, {3,2,2,7}, {2,3,3,6},
1566 {5,4,4,1}, {4,5,5,0}, {7,6,6,3}, {6,7,7,2}
1568 static char hex_hilbert_child_order[24][8] =
1570 {0,1,2,3,7,6,5,4}, {0,3,7,4,5,6,2,1}, {0,4,5,1,2,6,7,3},
1571 {1,0,3,2,6,7,4,5}, {1,2,6,5,4,7,3,0}, {1,5,4,0,3,7,6,2},
1572 {2,1,5,6,7,4,0,3}, {2,3,0,1,5,4,7,6}, {2,6,7,3,0,4,5,1},
1573 {3,0,4,7,6,5,1,2}, {3,2,1,0,4,5,6,7}, {3,7,6,2,1,5,4,0},
1574 {4,0,1,5,6,2,3,7}, {4,5,6,7,3,2,1,0}, {4,7,3,0,1,2,6,5},
1575 {5,1,0,4,7,3,2,6}, {5,4,7,6,2,3,0,1}, {5,6,2,1,0,3,7,4},
1576 {6,2,3,7,4,0,1,5}, {6,5,1,2,3,0,4,7}, {6,7,4,5,1,0,3,2},
1577 {7,3,2,6,5,1,0,4}, {7,4,0,3,2,1,5,6}, {7,6,5,4,0,1,2,3}
1579 static char hex_hilbert_child_state[24][8] =
1581 {1,2,2,7,7,21,21,17}, {2,0,0,22,22,16,16,8}, {0,1,1,15,15,6,6,23},
1582 {4,5,5,10,10,18,18,14}, {5,3,3,19,19,13,13,11}, {3,4,4,12,12,9,9,20},
1583 {8,7,7,17,17,23,23,2}, {6,8,8,0,0,15,15,22}, {7,6,6,21,21,1,1,16},
1584 {11,10,10,14,14,20,20,5}, {9,11,11,3,3,12,12,19}, {10,9,9,18,18,4,4,13},
1585 {13,14,14,5,5,19,19,10}, {14,12,12,20,20,11,11,4}, {12,13,13,9,9,3,3,18},
1586 {16,17,17,2,2,22,22,7}, {17,15,15,23,23,8,8,1}, {15,16,16,6,6,0,0,21},
1587 {20,19,19,11,11,14,14,3}, {18,20,20,4,4,10,10,12}, {19,18,18,13,13,5,5,9},
1588 {23,22,22,8,8,17,17,0}, {21,23,23,1,1,7,7,15}, {22,21,21,16,16,2,2,6}
1605 for (
int i = 0; i < 4; i++)
1607 int ch = quad_hilbert_child_order[state][i];
1608 int st = quad_hilbert_child_state[state][i];
1614 for (
int i = 0; i < 8; i++)
1616 int ch = hex_hilbert_child_order[state][i];
1617 int st = hex_hilbert_child_state[state][i];
1623 for (
int i = 0; i < 8; i++)
1667 node_order = (
char*) quad_hilbert_child_order;
1672 node_order = (
char*) hex_hilbert_child_order;
1679 int entry_node = -2;
1682 for (
int i = 0; i < root_count; i++)
1687 if (v_in < 0) { v_in = 0; }
1690 bool shared[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
1691 if (i+1 < root_count)
1694 for (
int j = 0; j < nch; j++)
1697 if (node >= 0) { shared[node] =
true; }
1702 int state =
Dim*v_in;
1703 for (
int j = 0; j <
Dim; j++)
1705 if (shared[(
int) node_order[nch*(state + j) + nch-1]])
1714 entry_node =
RetrieveNode(el, node_order[nch*state + nch-1]);
1726 MFEM_ABORT(
"invalid geometry");
1745 MFEM_VERIFY(tv.
visited ==
false,
"cyclic vertex dependencies.");
1751 for (
int i = 0; i < 3; i++)
1753 tv.
pos[i] = (pos1[i] + pos2[i]) * 0.5;
1768 for (
int i = 0; i < mvertices.
Size(); i++)
1787 if (
IsGhost(nc_elem)) {
continue; }
1789 const int* node = nc_elem.
node;
1796 for (
int j = 0; j < gi.
nv; j++)
1798 elem->GetVertices()[j] =
nodes[node[j]].vert_index;
1802 for (
int k = 0; k < gi.
nf; k++)
1804 const int* fv = gi.
faces[k];
1805 const Face* face =
faces.Find(node[fv[0]], node[fv[1]],
1806 node[fv[2]], node[fv[3]]);
1813 for (
int j = 0; j < 4; j++)
1823 for (
int j = 0; j < 2; j++)
1827 mboundary.
Append(segment);
1839 for (
int i = 0; i < edge_vertex->
Size(); i++)
1841 const int *ev = edge_vertex->
GetRow(i);
1844 MFEM_ASSERT(node && node->
HasEdge(),
"edge not found.");
1851 const int* fv = mesh->
GetFace(i)->GetVertices();
1855 MFEM_ASSERT(mesh->
GetFace(i)->GetNVertices() == 4,
"");
1861 MFEM_ASSERT(mesh->
GetFace(i)->GetNVertices() == 2,
"");
1863 face =
faces.Find(n0, n0, n1, n1);
1867 const int *ev = edge_vertex->
GetRow(i);
1868 MFEM_ASSERT((ev[0] == fv[0] && ev[1] == fv[1]) ||
1869 (ev[1] == fv[0] && ev[0] == fv[1]),
"");
1872 MFEM_ASSERT(face,
"face not found.");
1886 MFEM_ASSERT(
Dim >= 3,
"");
1889 int e1 =
nodes.FindId(v1, v2);
1890 int e2 =
nodes.FindId(v2, v3);
1891 int e3 = (e1 >= 0 &&
nodes[e1].HasVertex()) ?
nodes.FindId(v3, v4) : -1;
1892 int e4 = (e2 >= 0 &&
nodes[e2].HasVertex()) ?
nodes.FindId(v4, v1) : -1;
1895 if (mid) { mid[0] = e1, mid[1] = e2, mid[2] = e3, mid[3] = e4; }
1898 int midf1 = -1, midf2 = -1;
1899 if (e1 >= 0 && e3 >= 0) { midf1 =
nodes.FindId(e1, e3); }
1900 if (e2 >= 0 && e4 >= 0) { midf2 =
nodes.FindId(e2, e4); }
1903 MFEM_ASSERT(!(midf1 >= 0 && midf2 >= 0),
"incorrectly split face!");
1905 if (midf1 < 0 && midf2 < 0) {
return 0; }
1906 else if (midf1 >= 0) {
return 1; }
1912 for (
int i = 0; i < 8; i++)
1914 if (el.
node[i] == node) {
return i; }
1916 MFEM_ABORT(
"Node not found.");
1922 for (
int i = 0; i <
GI[(int) el.
geom].
nv; i++)
1926 if (abort) { MFEM_ABORT(
"Node not found."); }
1934 for (
int i = 0; i < gi.
ne; i++)
1936 const int* ev = gi.
edges[i];
1937 int n0 = el.
node[ev[0]];
1938 int n1 = el.
node[ev[1]];
1939 if ((n0 == vn0 && n1 == vn1) ||
1940 (n0 == vn1 && n1 == vn0)) {
return i; }
1942 MFEM_ABORT(
"Edge not found");
1948 for (
int i = 0; i < 6; i++)
1951 if ((a == fv[0] || a == fv[1] || a == fv[2] || a == fv[3]) &&
1952 (b == fv[0] || b == fv[1] || b == fv[2] || b == fv[3]) &&
1953 (c == fv[0] || c == fv[1] || c == fv[2] || c == fv[3]))
1958 MFEM_ABORT(
"Face not found.");
1976 for (
int i = 0, j; i < 4; i++)
1978 for (j = 0; j < 4; j++)
1980 if (fv[i] == master[j])
1983 for (
int k = 0; k < mat.
Height(); k++)
1985 mat(k,i) = tmp(k,j);
1990 MFEM_ASSERT(j != 4,
"node not found.");
2001 Face* fa =
faces.Find(vn0, vn1, vn2, vn3);
2024 Point mid0(pm(0), pm(1)), mid2(pm(2), pm(3));
2032 else if (split == 2)
2034 Point mid1(pm(1), pm(2)), mid3(pm(3), pm(0));
2049 if (
Dim < 3) {
return; }
2052 processed_faces = 0;
2059 MFEM_ASSERT(!el.
ref_type,
"not a leaf element.");
2062 for (
int j = 0; j < gi.
nf; j++)
2066 for (
int k = 0; k < 4; k++)
2071 int face =
faces.FindId(node[0], node[1], node[2], node[3]);
2072 MFEM_ASSERT(face >= 0,
"face not found!");
2078 if (processed_faces[face]) {
continue; }
2079 processed_faces[face] = 1;
2082 if (fa.
elem[0] >= 0 && fa.
elem[1] >= 0)
2094 TraverseFace(node[0], node[1], node[2], node[3], pm, 0);
2103 for (
int i = sb; i < se; i++)
2118 int mid =
nodes.FindId(vn0, vn1);
2119 if (mid < 0) {
return; }
2122 if (nd.
HasEdge() && level > 0)
2134 int v0index =
nodes[vn0].vert_index;
2135 int v1index =
nodes[vn1].vert_index;
2136 if (v0index > v1index) { sl.
edge_flags |= 2; }
2141 Face* fa =
faces.Find(vn0, vn0, vn1, vn1);
2142 MFEM_ASSERT(fa != NULL,
"");
2149 double tmid = (t0 + t1) / 2;
2163 processed_edges = 0;
2174 MFEM_ASSERT(!el.
ref_type,
"not a leaf element.");
2177 for (
int j = 0; j < gi.
ne; j++)
2180 const int* ev = gi.
edges[j];
2181 int node[2] = { el.
node[ev[0]], el.
node[ev[1]] };
2183 int enode =
nodes.FindId(node[0], node[1]);
2184 MFEM_ASSERT(enode >= 0,
"edge node not found!");
2187 MFEM_ASSERT(nd.
HasEdge(),
"edge not found!");
2195 int face =
faces.FindId(node[0], node[0], node[1], node[1]);
2196 MFEM_ASSERT(face >= 0,
"face not found!");
2208 if (processed_edges[enode]) {
continue; }
2209 processed_edges[enode] = 1;
2212 double t0 = 0.0, t1 = 1.0;
2213 int v0index =
nodes[node[0]].vert_index;
2214 int v1index =
nodes[node[1]].vert_index;
2215 int flags = (v0index > v1index) ? 1 : 0;
2228 for (
int i = sb; i < se; i++)
2245 int local = edge_local[sl.
index];
2262 processed_vertices = 0;
2270 for (
int j = 0; j <
GI[(int) el.
geom].
nv; j++)
2272 int node = el.
node[j];
2280 if (processed_vertices[index]) {
continue; }
2281 processed_vertices[index] = 1;
2291 oriented_matrix = point_matrix;
2295 MFEM_ASSERT(oriented_matrix.
Height() == 1 &&
2296 oriented_matrix.
Width() == 2,
"not an edge point matrix");
2300 oriented_matrix(0,0) = 1.0 - oriented_matrix(0,0);
2301 oriented_matrix(0,1) = 1.0 - oriented_matrix(0,1);
2305 std::swap(oriented_matrix(0,0), oriented_matrix(0,1));
2323 slaves.swap(empty.
slaves);
2325 inv_index.DeleteAll();
2330 return conforming.size() + masters.size() + slaves.size();
2335 if (!inv_index.Size())
2338 for (
unsigned i = 0; i < conforming.size(); i++)
2340 max_index = std::max(conforming[i].index, max_index);
2342 for (
unsigned i = 0; i < masters.size(); i++)
2344 max_index = std::max(masters[i].index, max_index);
2346 for (
unsigned i = 0; i < slaves.size(); i++)
2348 max_index = std::max(slaves[i].index, max_index);
2351 inv_index.SetSize(max_index + 1);
2354 for (
unsigned i = 0; i < conforming.size(); i++)
2356 inv_index[conforming[i].index] = (i << 2);
2358 for (
unsigned i = 0; i < masters.size(); i++)
2360 inv_index[masters[i].index] = (i << 2) + 1;
2362 for (
unsigned i = 0; i < slaves.size(); i++)
2364 inv_index[slaves[i].index] = (i << 2) + 2;
2368 MFEM_ASSERT(index >= 0 && index < inv_index.Size(),
"");
2369 int key = inv_index[index];
2370 MFEM_VERIFY(key >= 0,
"entity not found.");
2372 if (type) { *type = key & 0x3; }
2376 case 0:
return conforming[key >> 2];
2377 case 1:
return masters[key >> 2];
2378 case 2:
return slaves[key >> 2];
2379 default: MFEM_ABORT(
"internal error");
return conforming[0];
2388 int mid =
nodes.FindId(v0, v1);
2389 if (mid >= 0 &&
nodes[mid].HasVertex())
2425 int* I =
new int[nrows + 1];
2426 int** JJ =
new int*[nrows];
2436 MFEM_ASSERT(!el.
ref_type,
"not a leaf element.");
2439 int* node = el.
node;
2442 for (
int j = 0; j < gi.
ne; j++)
2444 const int* ev = gi.
edges[j];
2450 for (
int j = 0; j < gi.
nf; j++)
2452 const int* fv = gi.
faces[j];
2454 node[fv[2]], node[fv[3]], indices);
2461 int size = indices.
Size();
2463 JJ[i] =
new int[size];
2464 std::memcpy(JJ[i], indices.
GetData(), size *
sizeof(int));
2469 for (
int i = 0; i < nrows; i++)
2478 int *J =
new int[nnz];
2480 for (
int i = 0; i < nrows; i++)
2482 int cnt = I[i+1] - I[i];
2483 std::memcpy(J+nnz, JJ[i], cnt *
sizeof(
int));
2510 MFEM_VERIFY(elem_set.
Size() == nleaves,
"");
2519 for (
int i = 0; i < nleaves; i++)
2525 for (
int j = 0; j < nv; j++)
2532 for (
int j = 0; j < nv; j++)
2534 vmark[el.
node[j]] = 1;
2545 neighbor_set->
SetSize(nleaves);
2549 for (
int i = 0; i < nleaves; i++)
2557 for (
int j = 0; j < nv; j++)
2559 if (vmark[v[j]]) { hit =
true;
break; }
2566 for (
int j = 0; j < nv; j++)
2568 if (vmark[el.
node[j]]) { hit =
true;
break; }
2575 if (neighbor_set) { (*neighbor_set)[i] = 1; }
2581 static bool sorted_lists_intersect(
const int* a,
const int* b,
int na,
int nb)
2583 if (!na || !nb) {
return false; }
2584 int a_last = a[na-1], b_last = b[nb-1];
2585 if (*b < *a) {
goto l2; }
2587 if (a_last < *b) {
return false; }
2588 while (*a < *b) { a++; }
2589 if (*a == *b) {
return true; }
2591 if (b_last < *a) {
return false; }
2592 while (*b < *a) { b++; }
2593 if (*a == *b) {
return true; }
2616 while (stack.
Size())
2625 for (
int i = 0; i < nv; i++)
2631 for (
int i = 0; i < nv; i++)
2638 for (
int i = 0; i < 8 && el.
child[i] >= 0; i++)
2649 int nv1 = vert.
Size();
2654 for (
int i = 0; i < search_set->
Size(); i++)
2656 int testme = (*search_set)[i];
2663 bool hit = sorted_lists_intersect(v1, v2, nv1, nv2);
2668 for (
int j = 0; j < nv; j++)
2670 hit = sorted_lists_intersect(&el.
node[j], v1, 1, nv1);
2675 if (hit) { neighbors.
Append(testme); }
2689 for (
int i = 0; i < elems.
Size(); i++)
2695 for (
int j = 0; j < nv; j++)
2701 for (
int j = 0; j < nv; j++)
2703 vmark[el.
node[j]] = 1;
2713 for (
int i = 0; i < search_set->
Size(); i++)
2715 int testme = (*search_set)[i];
2721 for (
int j = 0; j < nv; j++)
2723 if (vmark[v[j]]) { hit =
true;
break; }
2729 for (
int j = 0; j < nv; j++)
2731 if (vmark[el.
node[j]]) { hit =
true;
break; }
2735 if (hit) { expanded.
Append(testme); }
2745 for (
int i = 0; i < np; i++)
2747 for (
int j = 0; j < points[0].dim; j++)
2749 point_matrix(j, i) = points[i].coord[j];
2761 Point(0, 0, 0),
Point(1, 0, 0),
Point(1, 1, 0),
Point(0, 1, 0),
2762 Point(0, 0, 1),
Point(1, 0, 1),
Point(1, 1, 1),
Point(0, 1, 1)
2773 MFEM_ABORT(
"unsupported geometry.");
2784 int ref_type = *ref_path++;
2785 int child = *ref_path++;
2791 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2792 Point mid67(pm(6), pm(7)), mid45(pm(4), pm(5));
2797 pm(4), mid45, mid67, pm(7));
2799 else if (child == 1)
2802 mid45, pm(5), pm(6), mid67);
2805 else if (ref_type == 2)
2807 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2808 Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2813 pm(4), pm(5), mid56, mid74);
2815 else if (child == 1)
2818 mid74, mid56, pm(6), pm(7));
2821 else if (ref_type == 4)
2823 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2824 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2829 mid04, mid15, mid26, mid37);
2831 else if (child == 1)
2834 pm(4), pm(5), pm(6), pm(7));
2837 else if (ref_type == 3)
2839 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2840 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2841 Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2842 Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2844 Point midf0(mid23, mid12, mid01, mid30);
2845 Point midf5(mid45, mid56, mid67, mid74);
2850 pm(4), mid45, midf5, mid74);
2852 else if (child == 1)
2855 mid45, pm(5), mid56, midf5);
2857 else if (child == 2)
2860 midf5, mid56, pm(6), mid67);
2862 else if (child == 3)
2865 mid74, midf5, mid67, pm(7));
2868 else if (ref_type == 5)
2870 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
2871 Point mid45(pm(4), pm(5)), mid67(pm(6), pm(7));
2872 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2873 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2875 Point midf1(mid01, mid15, mid45, mid04);
2876 Point midf3(mid23, mid37, mid67, mid26);
2881 mid04, midf1, midf3, mid37);
2883 else if (child == 1)
2886 midf1, mid15, mid26, midf3);
2888 else if (child == 2)
2891 mid45, pm(5), pm(6), mid67);
2893 else if (child == 3)
2896 pm(4), mid45, mid67, pm(7));
2899 else if (ref_type == 6)
2901 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
2902 Point mid56(pm(5), pm(6)), mid74(pm(7), pm(4));
2903 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2904 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2906 Point midf2(mid12, mid26, mid56, mid15);
2907 Point midf4(mid30, mid04, mid74, mid37);
2912 mid04, mid15, midf2, midf4);
2914 else if (child == 1)
2917 midf4, midf2, mid26, mid37);
2919 else if (child == 2)
2922 pm(4), pm(5), mid56, mid74);
2924 else if (child == 3)
2927 mid74, mid56, pm(6), pm(7));
2930 else if (ref_type == 7)
2932 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
2933 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
2934 Point mid45(pm(4), pm(5)), mid56(pm(5), pm(6));
2935 Point mid67(pm(6), pm(7)), mid74(pm(7), pm(4));
2936 Point mid04(pm(0), pm(4)), mid15(pm(1), pm(5));
2937 Point mid26(pm(2), pm(6)), mid37(pm(3), pm(7));
2939 Point midf0(mid23, mid12, mid01, mid30);
2940 Point midf1(mid01, mid15, mid45, mid04);
2941 Point midf2(mid12, mid26, mid56, mid15);
2942 Point midf3(mid23, mid37, mid67, mid26);
2943 Point midf4(mid30, mid04, mid74, mid37);
2944 Point midf5(mid45, mid56, mid67, mid74);
2946 Point midel(midf1, midf3);
2951 mid04, midf1, midel, midf4);
2953 else if (child == 1)
2956 midf1, mid15, midf2, midel);
2958 else if (child == 2)
2961 midel, midf2, mid26, midf3);
2963 else if (child == 3)
2966 midf4, midel, midf3, mid37);
2968 else if (child == 4)
2971 pm(4), mid45, midf5, mid74);
2973 else if (child == 5)
2976 mid45, pm(5), mid56, midf5);
2978 else if (child == 6)
2981 midf5, mid56, pm(6), mid67);
2983 else if (child == 7)
2986 mid74, midf5, mid67, pm(7));
2994 Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
3000 else if (child == 1)
3005 else if (ref_type == 2)
3007 Point mid12(pm(1), pm(2)), mid30(pm(3), pm(0));
3013 else if (child == 1)
3018 else if (ref_type == 3)
3020 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2));
3021 Point mid23(pm(2), pm(3)), mid30(pm(3), pm(0));
3022 Point midel(mid01, mid23);
3028 else if (child == 1)
3032 else if (child == 2)
3036 else if (child == 3)
3044 Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid20(pm(2), pm(0));
3050 else if (child == 1)
3054 else if (child == 2)
3058 else if (child == 3)
3066 for (
int i = 0; i < pm.
np; i++)
3068 for (
int j = 0; j < pm(i).dim; j++)
3070 matrix(j, i) = pm(i).coord[j];
3095 int &matrix = map[ref_path];
3096 if (!matrix) { matrix = map.size(); }
3099 emb.
parent = coarse_index;
3105 ref_path.push_back(0);
3107 for (
int i = 0; i < 8; i++)
3109 if (el.
child[i] >= 0)
3111 ref_path[ref_path.length()-1] = i;
3115 ref_path.resize(ref_path.length()-2);
3122 "GetRefinementTransforms() must be preceded by MarkCoarseLevel()"
3129 std::string ref_path;
3130 ref_path.reserve(100);
3147 for (RefPathMap::iterator it = map.begin(); it != map.end(); ++it)
3159 "GetDerefinementTransforms() must be preceded by Derefine().");
3165 std::map<int, int> mat_no;
3174 int &matrix = mat_no[code];
3175 if (!matrix) { matrix = mat_no.size(); }
3185 std::map<int, int>::iterator it;
3186 for (it = mat_no.begin(); it != mat_no.end(); ++it)
3189 int code = it->first;
3190 path[0] = code >> 3;
3205 std::map<Geometry::Type, DenseTensor>::iterator it;
3209 it->second.SetSize(0, 0, 0);
3216 static int sgn(
int x)
3218 return (x < 0) ? -1 : (x > 0) ? 1 : 0;
3221 static void HilbertSfc2D(
int x,
int y,
int ax,
int ay,
int bx,
int by,
3224 int w = std::abs(ax + ay);
3225 int h = std::abs(bx + by);
3227 int dax = sgn(ax), day = sgn(ay);
3228 int dbx = sgn(bx), dby = sgn(by);
3232 for (
int i = 0; i < w; i++, x += dax, y += day)
3241 for (
int i = 0; i < h; i++, x += dbx, y += dby)
3249 int ax2 = ax/2, ay2 = ay/2;
3250 int bx2 = bx/2, by2 = by/2;
3252 int w2 = std::abs(ax2 + ay2);
3253 int h2 = std::abs(bx2 + by2);
3257 if ((w2 & 0x1) && (w > 2))
3259 ax2 += dax, ay2 += day;
3262 HilbertSfc2D(x, y, ax2, ay2, bx, by, coords);
3263 HilbertSfc2D(x+ax2, y+ay2, ax-ax2, ay-ay2, bx, by, coords);
3267 if ((h2 & 0x1) && (h > 2))
3269 bx2 += dbx, by2 += dby;
3272 HilbertSfc2D(x, y, bx2, by2, ax2, ay2, coords);
3273 HilbertSfc2D(x+bx2, y+by2, ax, ay, bx-bx2, by-by2, coords);
3274 HilbertSfc2D(x+(ax-dax)+(bx2-dbx), y+(ay-day)+(by2-dby),
3275 -bx2, -by2, -(ax-ax2), -(ay-ay2), coords);
3279 static void HilbertSfc3D(
int x,
int y,
int z,
3280 int ax,
int ay,
int az,
3281 int bx,
int by,
int bz,
3282 int cx,
int cy,
int cz,
3285 int w = std::abs(ax + ay + az);
3286 int h = std::abs(bx + by + bz);
3287 int d = std::abs(cx + cy + cz);
3289 int dax = sgn(ax), day = sgn(ay), daz = sgn(az);
3290 int dbx = sgn(bx), dby = sgn(by), dbz = sgn(bz);
3291 int dcx = sgn(cx), dcy = sgn(cy), dcz = sgn(cz);
3294 if (h == 1 && d == 1)
3296 for (
int i = 0; i < w; i++, x += dax, y += day, z += daz)
3304 if (w == 1 && d == 1)
3306 for (
int i = 0; i < h; i++, x += dbx, y += dby, z += dbz)
3314 if (w == 1 && h == 1)
3316 for (
int i = 0; i < d; i++, x += dcx, y += dcy, z += dcz)
3325 int ax2 = ax/2, ay2 = ay/2, az2 = az/2;
3326 int bx2 = bx/2, by2 = by/2, bz2 = bz/2;
3327 int cx2 = cx/2, cy2 = cy/2, cz2 = cz/2;
3329 int w2 = std::abs(ax2 + ay2 + az2);
3330 int h2 = std::abs(bx2 + by2 + bz2);
3331 int d2 = std::abs(cx2 + cy2 + cz2);
3334 if ((w2 & 0x1) && (w > 2))
3336 ax2 += dax, ay2 += day, az2 += daz;
3338 if ((h2 & 0x1) && (h > 2))
3340 bx2 += dbx, by2 += dby, bz2 += dbz;
3342 if ((d2 & 0x1) && (d > 2))
3344 cx2 += dcx, cy2 += dcy, cz2 += dcz;
3348 if ((2*w > 3*h) && (2*w > 3*d))
3350 HilbertSfc3D(x, y, z,
3353 cx, cy, cz, coords);
3355 HilbertSfc3D(x+ax2, y+ay2, z+az2,
3356 ax-ax2, ay-ay2, az-az2,
3358 cx, cy, cz, coords);
3363 HilbertSfc3D(x, y, z,
3366 ax2, ay2, az2, coords);
3368 HilbertSfc3D(x+bx2, y+by2, z+bz2,
3370 bx-bx2, by-by2, bz-bz2,
3371 cx, cy, cz, coords);
3373 HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
3374 y+(ay-day)+(by2-dby),
3375 z+(az-daz)+(bz2-dbz),
3378 -(ax-ax2), -(ay-ay2), -(az-az2), coords);
3383 HilbertSfc3D(x, y, z,
3386 bx, by, bz, coords);
3388 HilbertSfc3D(x+cx2, y+cy2, z+cz2,
3391 cx-cx2, cy-cy2, cz-cz2, coords);
3393 HilbertSfc3D(x+(ax-dax)+(cx2-dcx),
3394 y+(ay-day)+(cy2-dcy),
3395 z+(az-daz)+(cz2-dcz),
3397 -(ax-ax2), -(ay-ay2), -(az-az2),
3398 bx, by, bz, coords);
3403 HilbertSfc3D(x, y, z,
3406 ax2, ay2, az2, coords);
3408 HilbertSfc3D(x+bx2, y+by2, z+bz2,
3411 bx-bx2, by-by2, bz-bz2, coords);
3413 HilbertSfc3D(x+(bx2-dbx)+(cx-dcx),
3414 y+(by2-dby)+(cy-dcy),
3415 z+(bz2-dbz)+(cz-dcz),
3418 -(cx-cx2), -(cy-cy2), -(cz-cz2), coords);
3420 HilbertSfc3D(x+(ax-dax)+bx2+(cx-dcx),
3421 y+(ay-day)+by2+(cy-dcy),
3422 z+(az-daz)+bz2+(cz-dcz),
3424 -(ax-ax2), -(ay-ay2), -(az-az2),
3425 bx-bx2, by-by2, bz-bz2, coords);
3427 HilbertSfc3D(x+(ax-dax)+(bx2-dbx),
3428 y+(ay-day)+(by2-dby),
3429 z+(az-daz)+(bz2-dbz),
3432 -(ax-ax2), -(ay-ay2), -(az-az2), coords);
3439 coords.
Reserve(2*width*height);
3441 if (width >= height)
3443 HilbertSfc2D(0, 0, width, 0, 0, height, coords);
3447 HilbertSfc2D(0, 0, 0, height, width, 0, coords);
3455 coords.
Reserve(3*width*height*depth);
3457 if (width >= height && width >= depth)
3459 HilbertSfc3D(0, 0, 0,
3462 0, 0, depth, coords);
3464 else if (height >= width && height >= depth)
3466 HilbertSfc3D(0, 0, 0,
3469 0, 0, depth, coords);
3473 HilbertSfc3D(0, 0, 0,
3476 0, height, 0, coords);
3484 bool oriented)
const
3490 int n0 = el.
node[ev[0]], n1 = el.
node[ev[1]];
3491 if (n0 > n1) { std::swap(n0, n1); }
3493 vert_index[0] =
nodes[n0].vert_index;
3494 vert_index[1] =
nodes[n1].vert_index;
3496 if (oriented && vert_index[0] > vert_index[1])
3498 std::swap(vert_index[0], vert_index[1]);
3508 int v0 =
nodes[el.
node[ev[0]]].vert_index;
3509 int v1 =
nodes[el.
node[ev[1]]].vert_index;
3511 return ((v0 < v1 && ev[0] > ev[1]) || (v0 > v1 && ev[0] < ev[1])) ? -1 : 1;
3515 int vert_index[4],
int edge_index[4],
3516 int edge_orientation[4])
const
3521 for (
int i = 0; i < 4; i++)
3523 vert_index[i] =
nodes[el.
node[fv[i]]].vert_index;
3526 for (
int i = 0; i < 4; i++)
3528 int j = (i+1) & 0x3;
3529 int n1 = el.
node[fv[i]];
3530 int n2 = el.
node[fv[j]];
3533 MFEM_ASSERT(en != NULL,
"edge not found.");
3536 edge_orientation[i] = (vert_index[i] < vert_index[j]) ? 1 : -1;
3542 MFEM_ASSERT(node >= 0,
"edge node not found.");
3545 int p1 = nd.
p1, p2 = nd.
p2;
3546 MFEM_ASSERT(p1 != p2,
"invalid edge node.");
3550 int n1p1 = n1.
p1, n1p2 = n1.
p2;
3551 int n2p1 = n2.p1, n2p2 = n2.p2;
3553 if ((n2p1 != n2p2) && (p1 == n2p1 || p1 == n2p2))
3557 if (n2.HasEdge()) {
return p2; }
3561 if ((n1p1 != n1p2) && (p2 == n1p1 || p2 == n1p2))
3565 if (n1.
HasEdge()) {
return p1; }
3575 MFEM_ASSERT(node >= 0 &&
nodes[node].HasEdge(),
"(v1, v2) is not an edge.");
3578 return (master >= 0) ?
nodes[master].edge_index : -1;
3584 int depth = 0, parent;
3585 while ((parent =
elements[elem].parent) != -1)
3600 int elem = fa.
elem[0];
3601 if (elem < 0) { elem = fa.
elem[1]; }
3602 MFEM_ASSERT(elem >= 0,
"Face has no elements?");
3610 for (
int i = 0; i < 4; i++)
3612 node[i] = el.
node[fv[i]];
3629 if (bdr_attr_is_ess[
faces[face].attribute - 1])
3634 for (
int j = 0; j < 4; j++)
3638 int enode =
nodes.FindId(node[j], node[(j+1) % 4]);
3639 MFEM_ASSERT(enode >= 0 &&
nodes[enode].HasEdge(),
"Edge not found.");
3668 bdr_vertices.
Sort();
3677 int mid =
nodes.FindId(vn1, vn2);
3678 if (mid < 0 || !
nodes[mid].HasVertex()) {
return 0; }
3683 int& h_level,
int& v_level)
const
3685 int hl1, hl2, vl1, vl2;
3691 h_level = v_level = 0;
3697 h_level = std::max(hl1, hl2);
3698 v_level = std::max(vl1, vl2) + 1;
3704 h_level = std::max(hl1, hl2) + 1;
3705 v_level = std::max(vl1, vl2);
3709 static int max8(
int a,
int b,
int c,
int d,
int e,
int f,
int g,
int h)
3711 return std::max(std::max(std::max(a, b), std::max(c, d)),
3712 std::max(std::max(e, f), std::max(g, h)));
3718 const int* node = el.
node;
3722 for (
int i = 0; i < gi.
ne; i++)
3724 const int* ev = gi.
edges[i];
3731 for (
int i = 0; i < gi.
nf; i++)
3733 const int* fv = gi.
faces[i];
3734 FaceSplitLevel(node[fv[0]], node[fv[1]], node[fv[2]], node[fv[3]],
3735 flevel[i][1], flevel[i][0]);
3738 splits[0] = max8(flevel[0][0], flevel[1][0], flevel[3][0], flevel[5][0],
3739 elevel[0], elevel[2], elevel[4], elevel[6]);
3741 splits[1] = max8(flevel[0][1], flevel[2][0], flevel[4][0], flevel[5][1],
3742 elevel[1], elevel[3], elevel[5], elevel[7]);
3744 splits[2] = max8(flevel[1][1], flevel[2][1], flevel[3][1], flevel[4][1],
3745 elevel[8], elevel[9], elevel[10], elevel[11]);
3749 splits[0] = std::max(elevel[0], elevel[2]);
3750 splits[1] = std::max(elevel[1], elevel[3]);
3754 splits[0] = std::max(elevel[0], std::max(elevel[1], elevel[2]));
3755 splits[1] = splits[0];
3759 MFEM_ABORT(
"Unsupported element geometry.");
3773 for (
int k = 0; k <
Dim; k++)
3775 if (splits[k] > max_level)
3777 ref_type |= (1 << k);
3795 MFEM_VERIFY(max_nc_level >= 1,
"'max_nc_level' must be 1 or greater.");
3802 if (!refinements.
Size()) {
break; }
3814 if (node->HasVertex() && node->p1 != node->p2) { nv++; }
3821 if (node->HasVertex() && node->p1 != node->p2)
3829 out << node->vert_index <<
" "
3842 input >>
id >> p1 >> p2;
3843 MFEM_VERIFY(input,
"problem reading vertex parents.");
3845 MFEM_VERIFY(
nodes.IdExists(
id),
"vertex " <<
id <<
" not found.");
3846 MFEM_VERIFY(
nodes.IdExists(p1),
"parent " << p1 <<
" not found.");
3847 MFEM_VERIFY(
nodes.IdExists(p2),
"parent " << p2 <<
" not found.");
3850 nodes.Reparent(
id, p1, p2);
3859 int num_top_level = 0;
3862 if (node->p1 == node->p2)
3864 MFEM_VERIFY(node.index() == node->p1,
"invalid top-level vertex.");
3865 MFEM_VERIFY(node->HasVertex(),
"top-level vertex not found.");
3866 MFEM_VERIFY(node->vert_index == node->p1,
"bad top-level vertex index");
3867 num_top_level = std::max(num_top_level, node->p1 + 1);
3872 for (
int i = 0; i < num_top_level; i++)
3878 static int ref_type_num_children[8] = { 0, 2, 2, 4, 2, 4, 4, 8 };
3885 int child_id[8], nch = 0;
3886 for (
int i = 0; i < 8 && el.
child[i] >= 0; i++)
3890 MFEM_ASSERT(nch == ref_type_num_children[(
int) el.
ref_type],
"");
3893 for (
int i = 0; i < nch; i++)
3895 out <<
" " << child_id[i];
3927 for (
int i = 0; i < 8 && el.
child[i] >= 0; i++)
3929 int old_id = el.
child[i];
3931 int new_id =
elements.Append(tmp_elements[old_id]);
3932 index_map[old_id] = new_id;
3933 el.
child[i] = new_id;
3957 if (
Dim == 3 && ref_type != 7) { iso =
false; }
3960 int nch = ref_type_num_children[ref_type];
3961 for (
int i = 0,
id; i < nch; i++)
3964 MFEM_VERIFY(
id >= 0,
"");
3967 "coarse element cannot be referenced before it is "
3968 "defined (id=" <<
id <<
").");
3971 MFEM_VERIFY(child.parent == -1,
3972 "element " <<
id <<
" cannot have two parents.");
3975 child.parent = elem;
3979 el.
geom = child.geom;
3997 if (el->parent == -1)
4000 index_map[el.index()] = new_id;
4006 for (
int i = 0; i < root_count; i++)
4014 for (
int i = 0; i < 2; i++)
4016 if (face->elem[i] >= 0)
4018 face->elem[i] = index_map[face->elem[i]];
4019 MFEM_ASSERT(face->elem[i] >= 0,
"");
4049 pmsize = slaves[0].point_matrix.MemoryUsage();
4052 return conforming.capacity() *
sizeof(
MeshId) +
4053 masters.capacity() *
sizeof(
Master) +
4054 slaves.capacity() *
sizeof(
Slave) +
4055 slaves.size() * pmsize;
4060 long mem = embeddings.MemoryUsage();
4062 std::map<Geometry::Type, DenseTensor>::const_iterator it;
4063 for (it=point_matrices.begin();
4064 it!=point_matrices.end(); it++)
4066 mem += it->second.MemoryUsage();
4074 return nodes.MemoryUsage() +
4075 faces.MemoryUsage() +
4110 <<
ref_stack.MemoryUsage() <<
" ref_stack\n"
4114 <<
sizeof(*this) <<
" NCMesh"
4122 static const double MiB = 1024.*1024.;
4124 "NCMesh statistics:\n"
4125 "------------------\n"
4126 " mesh and space dimensions : " <<
Dim <<
", " <<
spaceDim <<
"\n"
4127 " isotropic only : " << (
Iso ?
"yes" :
"no") <<
"\n"
4128 " number of Nodes : " << std::setw(9)
4129 <<
nodes.Size() <<
" + [ " << std::setw(9)
4130 <<
nodes.MemoryUsage()/MiB <<
" MiB ]\n"
4131 " free " << std::setw(9)
4132 <<
nodes.NumFreeIds() <<
"\n"
4133 " number of Faces : " << std::setw(9)
4134 <<
faces.Size() <<
" + [ " << std::setw(9)
4135 <<
faces.MemoryUsage()/MiB <<
" MiB ]\n"
4136 " free " << std::setw(9)
4137 <<
faces.NumFreeIds() <<
"\n"
4138 " number of Elements : " << std::setw(9)
4142 " free " << std::setw(9)
4144 " number of root elements : " << std::setw(9)
4146 " number of leaf elements : " << std::setw(9)
4148 " number of vertices : " << std::setw(9)
4150 " number of faces : " << std::setw(9)
4153 " conforming " << std::setw(9)
4155 " master " << std::setw(9)
4157 " slave " << std::setw(9)
4159 " number of edges : " << std::setw(9)
4162 " conforming " << std::setw(9)
4164 " master " << std::setw(9)
4166 " slave " << std::setw(9)
4168 " total memory : " << std::setw(17)
4169 <<
"[ " << std::setw(9) <<
MemoryUsage()/MiB <<
" MiB ]\n"
4180 for (
int j = 0; j <
Dim; j++)
4184 for (
int k = 0; k < 8; k++)
4186 if (elem->
node[k] >= 0)
4192 out << sum / count <<
" ";
4203 out <<
nodes.Size() <<
"\n";
4207 out << node.index() <<
" "
4208 << pos[0] <<
" " << pos[1] <<
" " << pos[2] <<
" "
4209 << node->p1 <<
" " << node->p2 <<
" "
4210 << node->vert_index <<
" " << node->edge_index <<
" "
4218 for (
int i = 0; i <
elements.Size(); i++)
4223 out << nleaves <<
"\n";
4224 for (
int i = 0; i <
elements.Size(); i++)
4229 out << gi.
nv <<
" ";
4230 for (
int j = 0; j < gi.
nv; j++)
4232 out << el.
node[j] <<
" ";
4239 out <<
faces.Size() <<
"\n";
4242 int elem = face->elem[0];
4243 if (elem < 0) { elem = face->elem[1]; }
4244 MFEM_ASSERT(elem >= 0,
"");
4253 for (
int i = 0; i < 4; i++)
4255 out <<
" " << el.
node[fv[i]];
NCList face_list
lazy-initialized list of faces, see GetFaceList
void CollectLeafElements(int elem, int state)
NCList edge_list
lazy-initialized list of edges, see GetEdgeList
NCMesh(const Mesh *mesh, std::istream *vertex_parents=NULL)
Geometry::Type GetGeometryType() const
int Size() const
Logical size of the array.
static void GridSfcOrdering3D(int width, int height, int depth, Array< int > &coords)
CoarseFineTransformations transforms
storage for data returned by Get[De]RefinementTransforms()
static void GridSfcOrdering2D(int width, int height, Array< int > &coords)
int NewQuadrilateral(int n0, int n1, int n2, int n3, int attr, int eattr0, int eattr1, int eattr2, int eattr3)
int elem[2]
up to 2 elements sharing the face
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Array< double > top_vertex_pos
coordinates of top-level vertices (organized as triples)
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
void SetDerefMatrixCodes(int parent, Array< int > &fine_coarse)
const CoarseFineTransformations & GetDerefinementTransforms()
virtual int GetNEdges() const =0
void AddColumnsInRow(int r, int ncol)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
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)
virtual void LimitNCLevel(int max_nc_level)
void GetFaceVerticesEdges(const MeshId &face_id, int vert_index[4], int edge_index[4], int edge_orientation[4]) const
Return Mesh vertex and edge indices of a face identified by 'face_id'.
int GetNBE() const
Returns number of boundary elements.
virtual void Trim()
Save memory by releasing all non-essential and cached data.
Geometry::Type GetElementGeometry() const
Return the type of elements in the mesh.
int index
element number in the Mesh, -1 if refined
void CopyElements(int elem, const BlockArray< Element > &tmp_elements, Array< int > &index_map)
void ForceRefinement(int vn1, int vn2, int vn3, int vn4)
virtual int GetNumGhostElements() const
Lists all edges/faces in the nonconforming mesh.
void GetMeshComponents(Array< mfem::Vertex > &mvertices, Array< mfem::Element * > &melements, Array< mfem::Element * > &mboundary) const
Return the basic Mesh arrays for the current finest level.
int Size() const
Return the number of items actually stored.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void SetDims(int rows, int nnz)
void Copy(Array ©) const
Create a copy of the current array.
void TraverseEdge(int vn0, int vn1, double t0, double t1, int flags, int level)
T * GetData()
Returns the data.
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 CollectFaceVertices(int v0, int v1, int v2, int v3, Array< int > &indices)
NCList vertex_list
lazy-initialized list of vertices, see GetVertexList
void InitDerefTransforms()
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
void InitRootState(int root_count)
int GetEdgeNCOrientation(const MeshId &edge_id) const
int GetNE() const
Returns number of elements.
void OrientedPointMatrix(DenseMatrix &oriented_matrix) const
Return the point matrix oriented according to the master and slave edges.
const Element * GetFace(int i) const
void CheckAnisoFace(int vn1, int vn2, int vn3, int vn4, int mid12, int mid34, int level=0)
virtual int GetNumGhostVertices() const
static PointMatrix pm_tri_identity
bool operator<(const Pair< A, B > &p, const Pair< A, B > &q)
Comparison operator for class Pair, based on the first element only.
void CountSplits(int elem, int splits[3]) const
void CollectDerefinements(int elem, Array< Connection > &list)
std::vector< MeshId > conforming
Data type quadrilateral element.
void GetLimitRefinements(Array< Refinement > &refinements, int max_level)
int PrintMemoryDetail() const
Geometry::Type GetElementBaseGeometry(int i) const
int spaceDim
dimensions of the elements and the vertex coordinates
Array< int > vertex_nodeId
int attribute
boundary element attribute, -1 if internal face
void DebugDump(std::ostream &out) const
void FindFaceNodes(int face, int node[4])
void DeleteAll()
Delete whole array.
int index
Mesh element number.
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
const MeshId & LookUp(int index, int *type=NULL) const
virtual void OnMeshUpdated(Mesh *mesh)
void DebugLeafOrder(std::ostream &out) const
void CollectEdgeVertices(int v0, int v1, Array< int > &indices)
Array< Refinement > ref_stack
stack of scheduled refinements (temporary)
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
virtual bool IsGhost(const Element &el) const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void CheckIsoFace(int vn1, int vn2, int vn3, int vn4, int en1, int en2, int en3, int en4, int midf)
Data type hexahedron element.
Face * GetFace(Element &elem, int face_no)
void TraverseRefinements(int elem, int coarse_index, std::string &ref_path, RefPathMap &map)
virtual const int * GetEdgeVertices(int) const =0
Element(Geometry::Type geom, int attr)
int PrintElements(std::ostream &out, int elem, int &coarse_id) const
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 FindAltParents(int node1, int node2)
int Append(const T &el)
Append element to array, resize if necessary.
void GetMatrix(DenseMatrix &point_matrix) const
int AddElement(const Element &el)
bool NodeSetY1(int node, int *n)
virtual void Derefine(const Array< int > &derefs)
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
std::vector< Master > masters
Array< int > coarse_elements
state of leaf_elements before Refine(), set by MarkCoarseLevel()
virtual void UpdateVertices()
update Vertex::index and vertex_nodeId
void AddConnection(int r, int c)
void RegisterFaces(int elem, int *fattr=NULL)
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.
Geometry::Type geom
Geometry::Type of the element.
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.
int parent
parent element, -1 if this is a root element, -2 if free
Data type triangle element.
const Element * GetElement(int i) const
virtual void ElementSharesFace(int elem, int local, int face)
void Sort()
Sorts the array. This requires operator< to be defined for T.
virtual void ElementSharesVertex(int elem, int local, int vnode)
int Size() const
Returns the number of TYPE I elements.
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
int element
NCMesh::Element containing this vertex/edge/face.
static GeomInfo GI[Geometry::NumGeom]
void Clear(bool hard=false)
bool NodeSetY2(int node, int *n)
int NewTriangle(int n0, int n1, int n2, int attr, int eattr0, int eattr1, int eattr2)
void DeleteUnusedFaces(const Array< int > &elemFaces)
int SpaceDimension() const
bool Iso
true if the mesh only contains isotropic refinements
std::map< std::string, int > RefPathMap
int FaceSplitType(int v1, int v2, int v3, int v4, int mid[4]=NULL) const
std::vector< Slave > slaves
virtual void BuildFaceList()
Nonconforming edge/face within a bigger edge/face.
void DerefineElement(int elem)
int ReorderFacePointMat(int v0, int v1, int v2, int v3, int elem, DenseMatrix &mat) const
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 RegisterElement(int e)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
static const PointMatrix & GetGeomIdentity(int geom)
mfem::Element * NewMeshElement(int geom) const
void DeleteLast()
Delete the last entry.
Helper struct for defining a connectivity table, see Table::MakeFromList.
void UpdateLeafElements()
int RetrieveNode(const Element &el, int index)
Return el.node[index] correctly, even if the element is refined.
bool NodeSetZ2(int node, int *n)
static GeomInfo & gi_quad
void BuildElementToVertexTable()
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
bool NodeSetZ1(int node, int *n)
void ForgetElement(int e)
int GetMidFaceNode(int en1, int en2, int en3, int en4)
void RefineElement(int elem, char ref_type)
virtual int GetNFaces(int &nFaceVertices) const =0
int NewHexahedron(int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int attr, int fattr0, int fattr1, int fattr2, int fattr3, int fattr4, int fattr5)
Array< int > boundary_faces
subset of all faces, set by BuildFaceList
int child[8]
2-8 children (if ref_type != 0)
Array< int > leaf_elements
static PointMatrix pm_quad_identity
void TraverseFace(int vn0, int vn1, int vn2, int vn3, const PointMatrix &pm, int level)
int EdgeSplitLevel(int vn1, int vn2) const
void SetIJ(int *newI, int *newJ, int newsize=-1)
Replace the I and J arrays with the given newI and newJ arrays.
static int find_element_edge(const Element &el, int vn0, int vn1)
void Initialize(const mfem::Element *elem)
Array< int > free_element_ids
bool NodeSetX2(int node, int *n)
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)
T & Last()
Return the last element in the array.
void NeighborExpand(const Array< int > &elems, Array< int > &expanded, const Array< int > *search_set=NULL)
BlockArray< Element > elements
void FindSetNeighbors(const Array< char > &elem_set, Array< int > *neighbors, Array< char > *neighbor_set=NULL)
int GetNEdges() const
Return the number of edges.
void GetEdgeVertices(const MeshId &edge_id, int vert_index[2], bool oriented=true) const
Return Mesh vertex indices of an edge identified by 'edge_id'.
long MemoryUsage() const
Return total number of bytes allocated.
void PrintCoarseElements(std::ostream &out) const
I/O: Print the "coarse_elements" section of the mesh file (ver. >= 1.1).
virtual void BuildVertexList()
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void UpdateElementToVertexTable()
virtual void ElementSharesEdge(int elem, int local, int enode)
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.
char ref_type
bit mask of X,Y,Z refinements (bits 0,1,2 respectively)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
static int find_node(const Element &el, int node)
void LoadVertexParents(std::istream &input)
int GetMidEdgeNode(int vn1, int vn2)
Rank 3 tensor (array of matrices)
Abstract data type element.
void UnrefElement(int elem, Array< int > &elemFaces)
Data type line segment element.
void PrintStats(std::ostream &out=mfem::out) const
int GetSingleElement() const
Return one of elem[0] or elem[1] and make sure the other is -1.
int FindNodeExt(const Element &el, int node, bool abort=false)
Extended version of find_node: works if 'el' is refined; optional abort.
virtual Type GetType() const =0
Returns element's type.
const double * CalcVertexPos(int node) const
int rank
processor number (ParNCMesh), -1 if undefined/unknown
bool NodeSetX1(int node, int *n)
int node[8]
element corners (if ref_type == 0)
const Element * GetBdrElement(int i) const
DenseMatrix point_matrix
position within the master edge/face
Defines the position of a fine element within a coarse element.
void RefElement(int elem)
void FaceSplitLevel(int vn1, int vn2, int vn3, int vn4, int &h_level, int &v_level) const
virtual void Refine(const Array< Refinement > &refinements)
int matrix
Index into the DenseTensor corresponding to the parent Geometry::Type stored in CoarseFineTransformat...
int GetEdgeMaster(int v1, int v2) const