34#include <unordered_map>
35#include <unordered_set>
39#if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
44#if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
49 int*,
int*,
int*,
int*,
int*,
idxtype*);
51 int*,
int*,
int*,
int*,
int*,
idxtype*);
53 int*,
int*,
int*,
int*,
int*,
idxtype*);
118 return sqrt((d_hat * d_hat) / (dir * dir));
157 if (coord[d] < min(d)) { min(d) = coord[d]; }
158 if (coord[d] > max(d)) { max(d) = coord[d]; }
164 const bool use_boundary =
false;
173 for (
int i = 0; i < ne; i++)
190 for (
int j = 0; j < pointmat.
Width(); j++)
192 for (
int d = 0; d < pointmat.
Height(); d++)
194 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
195 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
217 h_max = kappa_max = -h_min;
218 if (
dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
226 if (Vh) { (*Vh)(i) = h; }
227 if (Vk) { (*Vk)(i) =
kappa; }
229 if (h < h_min) { h_min = h; }
230 if (h > h_max) { h_max = h; }
244 if (!num_elems_by_geom[g]) {
continue; }
245 if (!first) { os <<
" + "; }
253 real_t h_min, h_max, kappa_min, kappa_max;
255 os <<
"Mesh Characteristics:";
260 num_elems_by_geom = 0;
261 for (
int i = 0; i <
GetNE(); i++)
272 <<
"Number of vertices : " <<
GetNV() <<
'\n'
273 <<
"Number of elements : " <<
GetNE() <<
'\n'
274 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n';
279 <<
"Number of vertices : " <<
GetNV() <<
'\n'
280 <<
"Number of elements : " <<
GetNE() <<
'\n'
281 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
282 <<
"h_min : " << h_min <<
'\n'
283 <<
"h_max : " << h_max <<
'\n';
288 <<
"Number of vertices : " <<
GetNV() <<
'\n'
289 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
290 <<
"Number of elements : " <<
GetNE() <<
" -- ";
293 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
295 <<
"h_min : " << h_min <<
'\n'
296 <<
"h_max : " << h_max <<
'\n'
297 <<
"kappa_min : " << kappa_min <<
'\n'
298 <<
"kappa_max : " << kappa_max <<
'\n';
303 num_bdr_elems_by_geom = 0;
304 for (
int i = 0; i <
GetNBE(); i++)
309 num_faces_by_geom = 0;
316 <<
"Number of vertices : " <<
GetNV() <<
'\n'
317 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
318 <<
"Number of faces : " <<
GetNFaces() <<
" -- ";
321 <<
"Number of elements : " <<
GetNE() <<
" -- ";
324 <<
"Number of bdr elem : " <<
GetNBE() <<
" -- ";
328 <<
"h_min : " << h_min <<
'\n'
329 <<
"h_max : " << h_max <<
'\n'
330 <<
"kappa_min : " << kappa_min <<
'\n'
331 <<
"kappa_max : " << kappa_max <<
'\n';
333 os <<
'\n' << std::flush;
349 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
352 MFEM_ABORT(
"Unknown element type");
381 for (
int j = 0; j < n; j++)
383 pm(k,j) =
nodes(vdofs[n*k+j]);
432 int nv =
elements[i]->GetNVertices();
433 const int *v =
elements[i]->GetVertices();
438 for (
int j = 0; j < nv; j++)
440 pm(k, j) =
nodes(k*n+v[j]);
454 for (
int j = 0; j < n; j++)
456 pm(k,j) =
nodes(vdofs[n*k+j]);
490 for (
int j = 0; j < n; j++)
492 int idx = vdofs[n*k+j];
493 pm(k,j) =
nodes((idx<0)? -1-idx:idx);
500 int elem_id, face_info;
516 "Mesh requires nodal Finite Element.");
525 ElTr->
SetFE(face_el);
547 const int *v = (
Dim == 1) ? &FaceNo :
faces[FaceNo]->GetVertices();
548 const int nv = (
Dim == 1) ? 1 :
faces[FaceNo]->GetNVertices();
552 for (
int j = 0; j < nv; j++)
572 for (
int j = 0; j < n; j++)
574 pm(i, j) =
nodes(vdofs[n*i+j]);
593 "Mesh requires nodal Finite Element.");
623 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
640 for (
int j = 0; j < nv; j++)
660 for (
int j = 0; j < n; j++)
662 pm(i, j) =
nodes(vdofs[n*i+j]);
665 EdTr->
SetFE(edge_el);
669 MFEM_ABORT(
"Not implemented.");
709 for (
int j = 0; j < 2; j++)
711 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
712 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
729 for (
int j = 0; j < 2; j++)
731 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
732 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
751 for (
int j = 0; j < 3; j++)
754 locpm(0, j) = vert.
x;
755 locpm(1, j) = vert.
y;
756 locpm(2, j) = vert.
z;
768 MFEM_VERIFY(i < 128,
"Local face index " << i/64
769 <<
" is not a triangular face of a wedge.");
777 for (
int j = 0; j < 3; j++)
780 locpm(0, j) = vert.
x;
781 locpm(1, j) = vert.
y;
782 locpm(2, j) = vert.
z;
793 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
794 <<
" is not a triangular face of a pyramid.");
802 for (
int j = 0; j < 3; j++)
805 locpm(0, j) = vert.
x;
806 locpm(1, j) = vert.
y;
807 locpm(2, j) = vert.
z;
824 for (
int j = 0; j < 4; j++)
827 locpm(0, j) = vert.
x;
828 locpm(1, j) = vert.
y;
829 locpm(2, j) = vert.
z;
841 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
842 <<
" is not a quadrilateral face of a wedge.");
848 for (
int j = 0; j < 4; j++)
851 locpm(0, j) = vert.
x;
852 locpm(1, j) = vert.
y;
853 locpm(2, j) = vert.
z;
864 MFEM_VERIFY(i < 64,
"Local face index " << i/64
865 <<
" is not a quadrilateral face of a pyramid.");
871 for (
int j = 0; j < 4; j++)
874 locpm(0, j) = vert.
x;
875 locpm(1, j) = vert.
y;
876 locpm(2, j) = vert.
z;
975 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
976 "face type " << face_type
977 <<
" and element type " << elem_type <<
"\n");
996 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
997 "face type " << face_type
998 <<
" and element type " << elem_type <<
"\n");
1030 FElTr.
Elem1 = &ElTr1;
1043 { MFEM_ABORT(
"NURBS mesh not supported!"); }
1046 FElTr.
Elem2 = &ElTr2;
1094 mfem::out <<
"\nInternal error: face id = " << FaceNo
1095 <<
", dist = " << dist <<
'\n';
1097 MFEM_ABORT(
"internal error");
1155 const FaceInfo &fi,
bool is_ghost)
const
1157#ifdef MFEM_THREAD_SAFE
1162 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1174 std::swap(composition(0,0), composition(0,1));
1175 std::swap(composition(1,0), composition(1,1));
1232 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1321 res.
Elem1No = element[0].index;
1322 res.Elem2No = element[1].index;
1323 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1324 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1325 res.NCFace = ncface;
1328 res.Elem1No = element[0].index;
1329 res.Elem2No = element[1].index;
1330 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1331 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1332 res.NCFace = ncface;
1335 res.Elem1No = element[0].index;
1336 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1339 res.Elem1No = element[0].index;
1340 res.Elem2No = -1 - element[1].index;
1341 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1342 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1345 res.Elem1No = element[0].index;
1346 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1349 res.Elem1No = element[0].index;
1350 res.Elem2No = -1 - element[1].index;
1351 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1352 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1357 res.Elem1No = element[0].index;
1358 res.Elem2No = -1 - element[1].index;
1359 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1360 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1368 os <<
"face topology=";
1378 os <<
"Non-conforming";
1385 os <<
"element[0].location=";
1399 os <<
"element[1].location=";
1413 os <<
"element[0].conformity=";
1430 os <<
"element[1].conformity=";
1447 os <<
"element[0].index=" << info.
element[0].
index <<
'\n'
1448 <<
"element[1].index=" << info.
element[1].
index <<
'\n'
1453 <<
"ncface=" << info.
ncface << std::endl;
1485 return faces[Face]->GetGeometryType();
1488 const int nc_face_id =
faces_info[Face].NCFace;
1490 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1551 "Could not determine a typical element Geometry!");
1559 face_marker.
SetSize(num_faces);
1561 for (
int f = 0;
f < num_faces;
f++)
1578 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1579 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1581 Array<bool> interior_bdr(max_bdr_attr); interior_bdr =
false;
1582 Array<bool> exterior_bdr(max_bdr_attr); exterior_bdr =
false;
1586 for (
int be = 0; be <
boundary.Size(); be++)
1588 const int bea =
boundary[be]->GetAttribute();
1590 if (bdr_marker[bea-1] != 0)
1596 interior_bdr[bea-1] =
true;
1600 exterior_bdr[bea-1] =
true;
1607 for (
int b = 0;
b < max_bdr_attr;
b++)
1609 if (bdr_marker[
b] != 0 && interior_bdr[
b])
1611 if (!excl || !exterior_bdr[
b])
1625 "Named set is not defined in this mesh!");
1627 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1631 for (
int b = 0;
b < max_bdr_attr;
b++)
1644 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1645 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1647 Array<bool> interior_bdr(max_bdr_attr); interior_bdr =
false;
1648 Array<bool> exterior_bdr(max_bdr_attr); exterior_bdr =
false;
1652 for (
int be = 0; be <
boundary.Size(); be++)
1654 const int bea =
boundary[be]->GetAttribute();
1660 interior_bdr[bea-1] =
true;
1664 exterior_bdr[bea-1] =
true;
1670 for (
int b = 0;
b < max_bdr_attr;
b++)
1672 if (bdr_marker[
b] == 0 && exterior_bdr[
b])
1674 if (!excl || !interior_bdr[
b])
1688 "Named set is not defined in this mesh!");
1689 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1690 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1694 for (
int b = 0;
b < max_bdr_attr;
b++)
1772 for (
int i = 0; i <
faces.Size(); i++)
1800#ifdef MFEM_USE_MEMALLOC
1824 for (
int i = 0; i < attribs.
Size(); i++)
1833 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1837 for (
int i = 0; i < attribs.
Size(); i++)
1846 MFEM_WARNING(
"Non-positive attributes in the domain!");
1868static void CheckEnlarge(
Array<T> &array,
int size)
1870 if (size >= array.Size()) { array.SetSize(size + 1); }
1893 "invalid 'coords' size: " << coords.
Size());
1906 for (
int j = 0; j < 3; j++)
1908 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1917 for (
int i = 0; i < nverts; i++)
1920 for (
int j = 0; j <
dim; j++)
1974 int vi[4] = {v1, v2, v3, v4};
1981#ifdef MFEM_USE_MEMALLOC
2021int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
2026 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
2039 static const int hex_to_tet[6][4] =
2041 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
2042 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
2046 for (
int i = 0; i < 6; i++)
2048 for (
int j = 0; j < 4; j++)
2050 ti[j] = vi[hex_to_tet[i][j]];
2058 static const int hex_to_wdg[2][6] =
2060 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
2064 for (
int i = 0; i < 2; i++)
2066 for (
int j = 0; j < 6; j++)
2068 ti[j] = vi[hex_to_wdg[i][j]];
2076 static const int hex_to_pyr[6][5] =
2078 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
2079 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
2083 for (
int i = 0; i < 6; i++)
2085 for (
int j = 0; j < 5; j++)
2087 ti[j] = vi[hex_to_pyr[i][j]];
2096 static const int quad_to_tri[4][2] =
2098 {0, 1}, {1, 2}, {2, 3}, {3, 0}
2104 ti[2] = elem_center_index;
2105 for (
int i = 0; i < num_faces; i++)
2107 for (
int j = 0; j < 2; j++)
2109 ti[j] = vi[quad_to_tri[i][j]];
2118 static const int quad_faces[4][2] =
2120 {0, 1}, {1, 2}, {2, 3}, {3, 0}
2124 for (
int i = 0; i < 4; i++)
2133 real_t r = 0.25, s = 0.25;
2134 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2135 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2140 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2141 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2146 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2147 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2152 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2153 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2157 static const int quad_faces_new[4][2] =
2159 { 1, 0}, { 2, 1}, { 3, 2}, { 0, 3}
2163 for (
int i = 0; i < num_faces; i++)
2165 for (
int j = 0; j < 2; j++)
2167 ti[j] = vi[quad_faces[i][j]];
2168 ti[j+2] = vnew_index[quad_faces_new[i][j]];
2176 std::map<std::array<int, 4>,
int> &hex_face_verts,
2182 return std::array<int, 4> {v[0], v[1], v[2], v[3]};
2186 static const int hex_to_tet[6][4] =
2188 { 0, 1, 2, 3 }, { 1, 2, 6, 5 }, { 5, 4, 7, 6},
2189 { 0, 1, 5, 4 }, { 2, 3, 7, 6 }, { 0,3, 7, 4}
2197 static const int tet_face[4][2] =
2199 {0, 1}, {1, 2}, {3, 2}, {3, 0}
2202 for (
int i = 0; i < num_faces; i++)
2204 for (
int j = 0; j < 4; j++)
2206 flist[j] = vi[hex_to_tet[i][j]];
2208 int face_center_index;
2210 auto t = get4arraysorted(flist);
2211 auto it = hex_face_verts.find(t);
2212 if (it == hex_face_verts.end())
2215 flist.
Size(), 3) - 1;
2216 hex_face_verts.insert({t, face_center_index});
2220 face_center_index = it->second;
2223 fti[2] = face_center_index;
2224 fti[3] = elem_center_index;
2225 for (
int j = 0; j < 4; j++)
2227 for (
int k = 0; k < 2; k++)
2229 fti[k] = flist[tet_face[j][k]];
2254 MFEM_ASSERT(bdr_elems.
Size() == new_be_to_face.
Size(),
"wrong size");
2255 for (
int i = 0; i < bdr_elems.
Size(); i++)
2306 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
2309 for (
int i = 0; i < 2; i++)
2311 for (
int j = 0; j < 3; j++)
2313 ti[j] = vi[quad_to_tri[i][j]];
2349 for (
int i = 0, j = 0; i <
faces_info.Size(); i++)
2365 "incorrect number of vertices: preallocated: " <<
vertices.Size()
2368 "incorrect number of elements: preallocated: " <<
elements.Size()
2371 "incorrect number of boundary elements: preallocated: "
2405 bool fix_orientation)
2408 if (fix_orientation)
2438 GeckoProgress(
real_t limit) : limit(limit) { sw.
Start(); }
2439 bool quit()
const override {
return limit > 0 && sw.
UserTime() > limit; }
2442class GeckoVerboseProgress :
public GeckoProgress
2448 GeckoVerboseProgress(
real_t limit) : GeckoProgress(limit) {}
2450 void beginorder(
const Graph* graph, Float cost)
const override
2451 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2452 void endorder(
const Graph* graph, Float cost)
const override
2453 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2455 void beginiter(
const Graph* graph,
2456 uint iter, uint maxiter, uint window)
const override
2458 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
2459 << window << std::flush;
2461 void enditer(
const Graph* graph, Float mincost, Float cost)
const override
2462 {
mfem::out <<
", cost = " << cost << endl; }
2467 int iterations,
int window,
2468 int period,
int seed,
bool verbose,
2474 GeckoProgress progress(time_limit);
2475 GeckoVerboseProgress vprogress(time_limit);
2478 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2486 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2488 const int *neighid = my_el_to_el.
GetRow(elemid);
2489 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2491 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2496 graph.
order(&functional, iterations, window, period, seed,
2497 verbose ? &vprogress : &progress);
2503 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2506 return graph.
cost();
2518 : coord(coord), dir(dir), points(points), mid(mid) {}
2520 bool operator()(
int i)
const
2522 return (points[3*i + coord] < mid) != dir;
2526static void HilbertSort2D(
int coord1,
2529 const Array<real_t> &points,
int *beg,
int *end,
2532 if (end - beg <= 1) {
return; }
2534 real_t xmid = (xmin + xmax)*0.5;
2535 real_t ymid = (ymin + ymax)*0.5;
2537 int coord2 = (coord1 + 1) % 2;
2540 int *p0 = beg, *p4 = end;
2541 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2542 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2543 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2547 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2548 ymin, xmin, ymid, xmid);
2550 if (p1 != p0 || p2 != p4)
2552 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2553 xmin, ymid, xmid, ymax);
2555 if (p2 != p0 || p3 != p4)
2557 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2558 xmid, ymid, xmax, ymax);
2562 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2563 ymid, xmax, ymin, xmid);
2567static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2568 const Array<real_t> &points,
int *beg,
int *end,
2572 if (end - beg <= 1) {
return; }
2574 real_t xmid = (xmin + xmax)*0.5;
2575 real_t ymid = (ymin + ymax)*0.5;
2576 real_t zmid = (zmin + zmax)*0.5;
2578 int coord2 = (coord1 + 1) % 3;
2579 int coord3 = (coord1 + 2) % 3;
2582 int *p0 = beg, *p8 = end;
2583 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2584 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2585 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2586 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2587 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2588 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2589 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2593 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2594 zmin, xmin, ymin, zmid, xmid, ymid);
2596 if (p1 != p0 || p2 != p8)
2598 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2599 ymin, zmid, xmin, ymid, zmax, xmid);
2601 if (p2 != p0 || p3 != p8)
2603 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2604 ymid, zmid, xmin, ymax, zmax, xmid);
2606 if (p3 != p0 || p4 != p8)
2608 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2609 xmin, ymax, zmid, xmid, ymid, zmin);
2611 if (p4 != p0 || p5 != p8)
2613 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2614 xmid, ymax, zmid, xmax, ymid, zmin);
2616 if (p5 != p0 || p6 != p8)
2618 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2619 ymax, zmid, xmax, ymid, zmax, xmid);
2621 if (p6 != p0 || p7 != p8)
2623 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2624 ymid, zmid, xmax, ymin, zmax, xmid);
2628 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2629 zmid, xmax, ymin, zmin, xmid, ymid);
2643 if (
spaceDim < 3) { points = 0.0; }
2646 for (
int i = 0; i <
GetNE(); i++)
2651 points[3*i + j] = center(j);
2658 indices.
Sort([&](
int a,
int b)
2659 {
return points[3*
a] < points[3*
b]; });
2664 HilbertSort2D(0,
false,
false,
2665 points, indices.
begin(), indices.
end(),
2666 min(0), min(1), max(0), max(1));
2671 HilbertSort3D(0,
false,
false,
false,
2672 points, indices.
begin(), indices.
end(),
2673 min(0), min(1), min(2), max(0), max(1), max(2));
2678 for (
int i = 0; i <
GetNE(); i++)
2680 ordering[indices[i]] = i;
2689 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2694 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
2698 MFEM_VERIFY(ordering.
Size() ==
GetNE(),
"invalid reordering array.")
2731 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2735 old_elem_node_vals[old_elid] =
new Vector(vals);
2741 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2743 int new_elid = ordering[old_elid];
2744 new_elements[new_elid] =
elements[old_elid];
2749 if (reorder_vertices)
2754 vertex_ordering = -1;
2756 int new_vertex_ind = 0;
2757 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2759 int *elem_vert =
elements[new_elid]->GetVertices();
2760 int nv =
elements[new_elid]->GetNVertices();
2761 for (
int vi = 0; vi < nv; ++vi)
2763 int old_vertex_ind = elem_vert[vi];
2764 if (vertex_ordering[old_vertex_ind] == -1)
2766 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2767 new_vertices[new_vertex_ind] =
vertices[old_vertex_ind];
2777 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2779 int *elem_vert =
elements[new_elid]->GetVertices();
2780 int nv =
elements[new_elid]->GetNVertices();
2781 for (
int vi = 0; vi < nv; ++vi)
2783 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2788 for (
int belid = 0; belid <
GetNBE(); ++belid)
2790 int *be_vert =
boundary[belid]->GetVertices();
2791 int nv =
boundary[belid]->GetNVertices();
2792 for (
int vi = 0; vi < nv; ++vi)
2794 be_vert[vi] = vertex_ordering[be_vert[vi]];
2823 nodes_fes->
Update(
false);
2826 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2828 int new_elid = ordering[old_elid];
2831 delete old_elem_node_vals[old_elid];
2880 length_idx[j].one =
GetLength(i, it.Column());
2881 length_idx[j].two = j;
2890 order[length_idx[i].two] = i;
2905 elements[i]->MarkEdge(v_to_v, order);
2912 boundary[i]->MarkEdge(v_to_v, order);
2919 if (*old_v_to_v && *old_elem_vert)
2926 if (*old_v_to_v == NULL)
2928 bool need_v_to_v =
false;
2936 if (dofs.
Size() > 0)
2948 if (*old_elem_vert == NULL)
2950 bool need_elem_vert =
false;
2952 for (
int i = 0; i <
GetNE(); i++)
2958 if (dofs.
Size() > 1)
2960 need_elem_vert =
true;
2966 *old_elem_vert =
new Table;
2968 for (
int i = 0; i <
GetNE(); i++)
2970 (*old_elem_vert)->AddColumnsInRow(i,
elements[i]->GetNVertices());
2972 (*old_elem_vert)->MakeJ();
2973 for (
int i = 0; i <
GetNE(); i++)
2978 (*old_elem_vert)->ShiftUpI();
2991 const int num_edge_dofs = old_dofs.
Size();
3002 if (num_edge_dofs > 0)
3011 const int old_i = (*old_v_to_v)(i, it.Column());
3012 const int new_i = it.Index();
3013 if (new_i == old_i) {
continue; }
3015 old_dofs.
SetSize(num_edge_dofs);
3016 new_dofs.
SetSize(num_edge_dofs);
3017 for (
int j = 0; j < num_edge_dofs; j++)
3019 old_dofs[j] = offset + old_i * num_edge_dofs + j;
3020 new_dofs[j] = offset + new_i * num_edge_dofs + j;
3024 for (
int j = 0; j < old_dofs.
Size(); j++)
3026 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3038 Table old_face_vertex;
3044 old_face_vertex.
MakeJ();
3047 faces[i]->GetNVertices());
3059 const int *old_v = old_face_vertex.
GetRow(i);
3061 switch (old_face_vertex.
RowSize(i))
3064 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
3068 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
3072 new_fdofs[new_i+1] = old_dofs.
Size();
3079 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
3082 switch (old_face_vertex.
RowSize(i))
3085 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
3086 new_v =
faces[new_i]->GetVertices();
3092 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
3093 new_v =
faces[new_i]->GetVertices();
3101 for (
int j = 0; j < old_dofs.
Size(); j++)
3104 const int old_j = dof_ord[j];
3105 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
3109 for (
int j = 0; j < old_dofs.
Size(); j++)
3111 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3133 for (
int i = 0; i <
GetNE(); i++)
3135 const int *old_v = old_elem_vert->
GetRow(i);
3136 const int *new_v =
elements[i]->GetVertices();
3143 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
3157 <<
" FE collection) are not supported yet!");
3161 MFEM_VERIFY(dof_ord != NULL,
3162 "FE collection '" << fec->
Name()
3163 <<
"' does not define reordering for "
3167 for (
int j = 0; j < new_dofs.
Size(); j++)
3170 const int old_j = dof_ord[j];
3171 new_dofs[old_j] = offset + j;
3173 offset += new_dofs.
Size();
3176 for (
int j = 0; j < old_dofs.
Size(); j++)
3178 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3216 MFEM_ASSERT(
NURBSext,
"SetPatchAttribute is only for NURBS meshes");
3219 for (
auto e : elems)
3227 MFEM_ASSERT(
NURBSext,
"GetPatchAttribute is only for NURBS meshes");
3233 MFEM_ASSERT(
NURBSext,
"SetPatchBdrAttribute is only for NURBS meshes");
3237 for (
auto be : bdryelems)
3245 MFEM_ASSERT(
NURBSext,
"GetBdrPatchBdrAttribute is only for NURBS meshes");
3273 if (generate_edges == 1)
3291 bool fix_orientation)
3308 if (generate_edges == 1)
3374 bool generate_edges =
true;
3383 MFEM_VERIFY(
ncmesh == NULL,
"");
3418 if (
Dim > 1 && generate_edges)
3482 const bool check_orientation =
true;
3483 const bool curved = (
Nodes != NULL);
3484 const bool may_change_topology =
3486 ( check_orientation && fix_orientation &&
3490 Table *old_elem_vert = NULL;
3492 if (curved && may_change_topology)
3497 if (check_orientation)
3507 if (may_change_topology)
3512 delete old_elem_vert;
3533 for (
int i = 0; i < num_faces; i++)
3536 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
3537 " Interior face with incompatible orientations.");
3548 int NVert, NElem, NBdrElem;
3550 NVert = (nx+1) * (ny+1) * (nz+1);
3551 NElem = nx * ny * nz;
3552 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3561 NBdrElem += 2*nx*ny;
3566 NVert += nx * ny * nz;
3569 InitMesh(3, 3, NVert, NElem, NBdrElem);
3575 for (z = 0; z <= nz; z++)
3577 coord[2] = ((
real_t) z / nz) * sz;
3578 for (y = 0; y <= ny; y++)
3580 coord[1] = ((
real_t) y / ny) * sy;
3581 for (x = 0; x <= nx; x++)
3583 coord[0] = ((
real_t) x / nx) * sx;
3590 for (z = 0; z < nz; z++)
3592 coord[2] = (((
real_t) z + 0.5) / nz) * sz;
3593 for (y = 0; y < ny; y++)
3595 coord[1] = (((
real_t) y + 0.5) / ny) * sy;
3596 for (x = 0; x < nx; x++)
3598 coord[0] = (((
real_t) x + 0.5) / nx) * sx;
3605#define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3606#define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3613 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3615 for (
int k = 0; k < nx*ny*nz; k++)
3622 ind[0] = VTX(x , y , z );
3623 ind[1] = VTX(x+1, y , z );
3624 ind[2] = VTX(x+1, y+1, z );
3625 ind[3] = VTX(x , y+1, z );
3626 ind[4] = VTX(x , y , z+1);
3627 ind[5] = VTX(x+1, y , z+1);
3628 ind[6] = VTX(x+1, y+1, z+1);
3629 ind[7] = VTX(x , y+1, z+1);
3637 for (z = 0; z < nz; z++)
3639 for (y = 0; y < ny; y++)
3641 for (x = 0; x < nx; x++)
3644 ind[0] = VTX(x , y , z );
3645 ind[1] = VTX(x+1, y , z );
3646 ind[2] = VTX(x+1, y+1, z );
3647 ind[3] = VTX(x , y+1, z );
3648 ind[4] = VTX(x , y , z+1);
3649 ind[5] = VTX(x+1, y , z+1);
3650 ind[6] = VTX(x+1, y+1, z+1);
3651 ind[7] = VTX( x, y+1, z+1);
3663 ind[8] = VTXP(x, y, z);
3677 for (y = 0; y < ny; y++)
3679 for (x = 0; x < nx; x++)
3682 ind[0] = VTX(x , y , 0);
3683 ind[1] = VTX(x , y+1, 0);
3684 ind[2] = VTX(x+1, y+1, 0);
3685 ind[3] = VTX(x+1, y , 0);
3702 for (y = 0; y < ny; y++)
3704 for (x = 0; x < nx; x++)
3707 ind[0] = VTX(x , y , nz);
3708 ind[1] = VTX(x+1, y , nz);
3709 ind[2] = VTX(x+1, y+1, nz);
3710 ind[3] = VTX(x , y+1, nz);
3727 for (z = 0; z < nz; z++)
3729 for (y = 0; y < ny; y++)
3732 ind[0] = VTX(0 , y , z );
3733 ind[1] = VTX(0 , y , z+1);
3734 ind[2] = VTX(0 , y+1, z+1);
3735 ind[3] = VTX(0 , y+1, z );
3748 for (z = 0; z < nz; z++)
3750 for (y = 0; y < ny; y++)
3753 ind[0] = VTX(nx, y , z );
3754 ind[1] = VTX(nx, y+1, z );
3755 ind[2] = VTX(nx, y+1, z+1);
3756 ind[3] = VTX(nx, y , z+1);
3769 for (x = 0; x < nx; x++)
3771 for (z = 0; z < nz; z++)
3774 ind[0] = VTX(x , 0, z );
3775 ind[1] = VTX(x+1, 0, z );
3776 ind[2] = VTX(x+1, 0, z+1);
3777 ind[3] = VTX(x , 0, z+1);
3790 for (x = 0; x < nx; x++)
3792 for (z = 0; z < nz; z++)
3795 ind[0] = VTX(x , ny, z );
3796 ind[1] = VTX(x , ny, z+1);
3797 ind[2] = VTX(x+1, ny, z+1);
3798 ind[3] = VTX(x+1, ny, z );
3815 ofstream test_stream(
"debug.mesh");
3817 test_stream.close();
3845 for (
real_t j = 0; j < ny+1; j++)
3847 real_t cy = (j / ny) * sy;
3848 for (
real_t i = 0; i < nx+1; i++)
3850 real_t cx = (i / nx) * sx;
3857 for (
int y = 0; y < ny; y++)
3859 for (
int x = 0; x < nx; x++)
3861 ind[0] = x + y*(nx+1);
3862 ind[1] = x + 1 +y*(nx+1);
3863 ind[2] = x + 1 + (y+1)*(nx+1);
3864 ind[3] = x + (y+1)*(nx+1);
3870 for (
int i = 0; i < nx; i++)
3876 for (
int j = 0; j < ny; j++)
3919 for (
real_t j = 0; j < ny+1; j++)
3921 real_t cy = (j / ny) * sy;
3922 for (
real_t i = 0; i < nx+1; i++)
3924 real_t cx = (i / nx) * sx;
3931 for (
int y = 0; y < ny; y++)
3933 for (
int x = 0; x < nx; x++)
3935 ind[0] = x + y*(nx+1);
3936 ind[1] = x + 1 +y*(nx+1);
3937 ind[2] = x + 1 + (y+1)*(nx+1);
3938 ind[3] = x + (y+1)*(nx+1);
3944 for (
int i = 0; i < nx; i++)
3950 for (
int j = 0; j < ny; j++)
3976 const int NVert = (nx+1) * (ny+1) * (nz+1);
3977 const int NElem = nx * ny * nz * 24;
3978 const int NBdrElem = 2*(nx*ny+nx*nz+ny*nz)*4;
3980 InitMesh(3, 3, NVert, NElem, NBdrElem);
3985 for (
real_t z = 0; z <= nz; z++)
3987 coord[2] = ( z / nz) * sz;
3988 for (
real_t y = 0; y <= ny; y++)
3990 coord[1] = (y / ny) * sy;
3991 for (
real_t x = 0; x <= nx; x++)
3993 coord[0] = (x / nx) * sx;
3999 std::map<std::array<int, 4>,
int> hex_face_verts;
4000 auto VertexIndex = [nx, ny](
int xc,
int yc,
int zc)
4002 return xc + (yc + zc*(ny+1))*(nx+1);
4006 for (
int z = 0; z < nz; z++)
4008 for (
int y = 0; y < ny; y++)
4010 for (
int x = 0; x < nx; x++)
4013 ind[0] = VertexIndex(x , y , z );
4014 ind[1] = VertexIndex(x+1, y , z );
4015 ind[2] = VertexIndex(x+1, y+1, z );
4016 ind[3] = VertexIndex(x , y+1, z );
4017 ind[4] = VertexIndex(x , y , z+1);
4018 ind[5] = VertexIndex(x+1, y , z+1);
4019 ind[6] = VertexIndex(x+1, y+1, z+1);
4020 ind[7] = VertexIndex( x, y+1, z+1);
4026 hex_face_verts.clear();
4035 std::map<std::array<int, 3>,
int> tet_face_count;
4037 std::map<std::array<int, 3>,
int> face_count_map;
4042 return std::array<int, 3>{v[0], v[1], v[2]};
4051 for (
int j = 0; j < el_faces.
Size(); j++)
4054 auto t = get3array(vertidxs);
4055 auto it = tet_face_count.find(t);
4056 if (it == tet_face_count.end())
4058 tet_face_count.insert({t, 1});
4059 face_count_map.insert({t, el_faces[j]});
4068 for (
const auto &edge : tet_face_count)
4070 if (edge.second == 1)
4072 int facenum = (face_count_map.find(edge.first))->second;
4079 ofstream test_stream(
"debug.mesh");
4081 test_stream.close();
4090 bool generate_edges,
bool sfc_ordering)
4114 for (j = 0; j < ny+1; j++)
4116 cy = ((
real_t) j / ny) * sy;
4117 for (i = 0; i < nx+1; i++)
4119 cx = ((
real_t) i / nx) * sx;
4131 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
4133 for (k = 0; k < nx*ny; k++)
4137 ind[0] = i + j*(nx+1);
4138 ind[1] = i + 1 +j*(nx+1);
4139 ind[2] = i + 1 + (j+1)*(nx+1);
4140 ind[3] = i + (j+1)*(nx+1);
4147 for (j = 0; j < ny; j++)
4149 for (i = 0; i < nx; i++)
4151 ind[0] = i + j*(nx+1);
4152 ind[1] = i + 1 +j*(nx+1);
4153 ind[2] = i + 1 + (j+1)*(nx+1);
4154 ind[3] = i + (j+1)*(nx+1);
4163 for (i = 0; i < nx; i++)
4169 for (j = 0; j < ny; j++)
4191 for (j = 0; j < ny+1; j++)
4193 cy = ((
real_t) j / ny) * sy;
4194 for (i = 0; i < nx+1; i++)
4196 cx = ((
real_t) i / nx) * sx;
4205 for (j = 0; j < ny; j++)
4207 for (i = 0; i < nx; i++)
4209 ind[0] = i + j*(nx+1);
4210 ind[1] = i + 1 + (j+1)*(nx+1);
4211 ind[2] = i + (j+1)*(nx+1);
4214 ind[1] = i + 1 + j*(nx+1);
4215 ind[2] = i + 1 + (j+1)*(nx+1);
4223 for (i = 0; i < nx; i++)
4229 for (j = 0; j < ny; j++)
4239 MFEM_ABORT(
"Unsupported element type.");
4245 if (generate_edges == 1)
4283 for (j = 0; j < n+1; j++)
4289 for (j = 0; j < n; j++)
4316 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4368 for (
int i = 0; i <
faces.Size(); i++)
4410 if (
dynamic_cast<const ParMesh*
>(&mesh))
4422 if (mesh.
Nodes && copy_nodes)
4454 int refine,
bool fix_orientation)
4458 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
4459 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
4476 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
4486 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
4522 ref_factors = ref_factor;
4535Mesh::Mesh(
const std::string &filename,
int generate_edges,
int refine,
4536 bool fix_orientation)
4537 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4546 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
4550 Load(imesh, generate_edges, refine, fix_orientation);
4555 bool fix_orientation)
4556 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4559 Load(input, generate_edges, refine, fix_orientation);
4568 "Not enough vertices in external array : "
4569 "len_vertex_data = "<< len_vertex_data <<
", "
4574 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
4579 memcpy(vertex_data,
vertices.GetData(),
4588 int *element_attributes,
int num_elements,
4590 int *boundary_attributes,
int num_boundary_elements,
4592 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4594 if (space_dimension == -1)
4600 num_boundary_elements);
4603 int boundary_index_stride = num_boundary_elements > 0 ?
4607 vertices.MakeRef(
reinterpret_cast<Vertex*
>(vertices_), num_vertices);
4610 for (
int i = 0; i < num_elements; i++)
4613 elements[i]->SetVertices(element_indices + i * element_index_stride);
4614 elements[i]->SetAttribute(element_attributes[i]);
4618 for (
int i = 0; i < num_boundary_elements; i++)
4621 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
4622 boundary[i]->SetAttribute(boundary_attributes[i]);
4630 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4667 MFEM_ABORT(
"NURBS mesh has no patches.");
4681#ifdef MFEM_USE_MEMALLOC
4690 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
4703 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
4706 for (
int i = 0; i < nv; i++)
4719 for (
int j = 0; j < nv; j++)
4791 MFEM_ABORT(
"invalid element type: " << type);
4798 std::string parse_tag)
4800 int curved = 0, read_gf = 1;
4801 bool finalize_topo =
true;
4805 MFEM_ABORT(
"Input stream is not open");
4812 getline(input, mesh_type);
4816 int mfem_version = 0;
4817 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
4818 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
4819 else if (mesh_type ==
"MFEM mesh v1.3") { mfem_version = 13; }
4823 int mfem_nc_version = 0;
4824 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
4825 else if (mesh_type ==
"MFEM NC mesh v1.1") { mfem_nc_version = 11; }
4826 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
4834 if (mfem_version >= 12 && parse_tag.empty())
4836 parse_tag =
"mfem_mesh_end";
4840 else if (mfem_nc_version)
4842 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4849 MFEM_VERIFY(mfem_nc_version >= 10,
4850 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
4851 "used to load a parallel nonconforming mesh, sorry.");
4854 input, mfem_nc_version, curved, is_nc);
4859 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
4872 else if (mesh_type ==
"linemesh")
4876 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
4878 if (mesh_type ==
"curved_areamesh2")
4884 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
4888 else if (mesh_type ==
"TrueGrid")
4892 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
4894 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
4896 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4897 "Unsupported VTK format");
4898 ReadVTKMesh(input, curved, read_gf, finalize_topo);
4900 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
4904 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
4908 else if (mesh_type ==
"MFEM NURBS mesh v1.1")
4912 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
4917 else if (mesh_type ==
"$MeshFormat")
4922 ((mesh_type.size() > 2 &&
4923 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
4924 (mesh_type.size() > 3 &&
4925 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
4930#ifdef MFEM_USE_NETCDF
4933 MFEM_ABORT(
"NetCDF support requires configuration with"
4934 " MFEM_USE_NETCDF=YES");
4940 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
4941 " Use mfem::named_ifgzstream for input.");
4947 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
4973 bool generate_bdr =
false;
4978 if (curved && read_gf)
4992 if (mfem_version >= 12)
4998 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
4999 getline(input, line);
5005 if (line ==
"mfem_mesh_end") {
break; }
5007 while (line != parse_tag);
5009 else if (mfem_nc_version >= 10)
5014 MFEM_VERIFY(ident ==
"mfem_mesh_end",
5015 "invalid mesh: end of file tag not found");
5022 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
5024 int i, j, ie, ib, iv, *v, nv;
5055 for (i = 0; i < num_pieces; i++)
5062 for (i = 0; i < num_pieces; i++)
5068 for (j = 0; j < m->
GetNE(); j++)
5073 for (j = 0; j < m->
GetNBE(); j++)
5078 for (
int k = 0; k < nv; k++)
5080 v[k] = lvert_vert[v[k]];
5085 for (j = 0; j < m->
GetNV(); j++)
5097 for (i = 0; i < num_pieces; i++)
5108 for (i = 0; i < num_pieces; i++)
5112 for (j = 0; j < m->
GetNE(); j++)
5117 for (
int k = 0; k < nv; k++)
5124 for (j = 0; j < m->
GetNBE(); j++)
5129 for (
int k = 0; k < nv; k++)
5136 for (j = 0; j < m->
GetNV(); j++)
5150 for (i = 0; i < num_pieces; i++)
5152 gf_array[i] = mesh_array[i]->
GetNodes();
5165 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
5168 ref_factors = ref_factor;
5179 int orig_ne = orig_mesh.
GetNE();
5180 MFEM_VERIFY(ref_factors.
Size() == orig_ne,
5181 "Number of refinement factors must equal number of elements")
5182 MFEM_VERIFY(orig_ne == 0 ||
5183 ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
5186 "Invalid refinement type. Must use closed basis type.");
5188 int min_ref = orig_ne > 0 ? ref_factors.
Min() : 1;
5189 int max_ref = orig_ne > 0 ? ref_factors.
Max() : 1;
5191 bool var_order = (min_ref != max_ref);
5204 for (
int i = 0; i < orig_ne; i++)
5220 for (
int el = 0; el < orig_ne; el++)
5232 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
5233 for (
int i = 0; i < phys_pts.
Width(); i++)
5242 for (
int k = 0; k < nvert; k++)
5245 v[k] = rdofs[c2h_map[cid]];
5258 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
5269 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
5275 for (
int k = 0; k < nvert; k++)
5278 v[k] = rdofs[c2h_map[cid]];
5300 for (
int iel = 0; iel < orig_ne; iel++)
5309 const int *node_map = NULL;
5312 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
5313 const int *vertex_map = vertex_fec.
GetDofMap(geom);
5314 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
5315 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
5318 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
5321 int iv = vertex_map[iv_lex];
5323 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
5325 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
5328 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
5339 using GeomRef = std::pair<Geometry::Type, int>;
5340 std::map<GeomRef, int> point_matrices_offsets;
5342 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5346 GeomRef id(geom, ref_factors[el_coarse]);
5347 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
5353 point_matrices_offsets[id] = n_point_matrices[geom];
5354 n_point_matrices[geom] += nref_el;
5361 int nmatrices = n_point_matrices[geom];
5368 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5371 int ref = ref_factors[el_coarse];
5372 int offset = point_matrices_offsets[GeomRef(geom, ref)];
5378 for (
int k = 0; k < nvert; k++)
5411 "Mesh::MakeSimplicial requires a properly oriented input mesh");
5413 "Mesh::MakeSimplicial does not support non-conforming meshes.")
5420 Mesh copy(orig_mesh);
5425 int nv = orig_mesh.
GetNV();
5426 int ne = orig_mesh.
GetNE();
5427 int nbe = orig_mesh.
GetNBE();
5440 int new_ne = 0, new_nbe = 0;
5441 for (
int i=0; i<ne; ++i)
5445 for (
int i=0; i<nbe; ++i)
5454 for (
int i=0; i<nv; ++i)
5464 if (vglobal == NULL)
5467 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
5468 vglobal = vglobal_id.
GetData();
5471 constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
5472 constexpr int quad_ntris = 2, prism_ntets = 3;
5473 static const int quad_trimap[2][nv_tri*quad_ntris] =
5485 static const int prism_rot[nv_prism*nv_prism] =
5494 static const int prism_f[nv_quad] = {1, 2, 5, 4};
5495 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
5509 static const int hex_rot[nv_hex*nv_hex] =
5511 0, 1, 2, 3, 4, 5, 6, 7,
5512 1, 0, 4, 5, 2, 3, 7, 6,
5513 2, 1, 5, 6, 3, 0, 4, 7,
5514 3, 0, 1, 2, 7, 4, 5, 6,
5515 4, 0, 3, 7, 5, 1, 2, 6,
5516 5, 1, 0, 4, 6, 2, 3, 7,
5517 6, 2, 1, 5, 7, 3, 0, 4,
5518 7, 3, 2, 6, 4, 0, 1, 5
5520 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
5521 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
5522 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
5523 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
5524 static const int hex_tetmap0[nv_tet*5] =
5531 static const int hex_tetmap1[nv_tet*6] =
5538 static const int hex_tetmap2[nv_tet*6] =
5545 static const int hex_tetmap3[nv_tet*6] =
5552 static const int *hex_tetmaps[4] =
5554 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
5557 auto find_min = [](
const int*
a,
int n) {
return std::min_element(
a,
a+n)-
a; };
5559 for (
int i=0; i<ne; ++i)
5561 const int *v = orig_mesh.
elements[i]->GetVertices();
5565 if (num_subdivisions[orig_geom] == 1)
5577 for (
int itri=0; itri<quad_ntris; ++itri)
5582 for (
int iv=0; iv<nv_tri; ++iv)
5584 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
5592 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
5595 int irot = find_min(vg, nv_prism);
5596 for (
int iv=0; iv<nv_prism; ++iv)
5598 int jv = prism_rot[iv + irot*nv_prism];
5603 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
5604 int j = find_min(q, nv_quad);
5605 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
5606 for (
int itet=0; itet<prism_ntets; ++itet)
5611 for (
int iv=0; iv<nv_tet; ++iv)
5613 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
5621 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
5625 int irot = find_min(vg, nv_hex);
5626 for (
int iv=0; iv<nv_hex; ++iv)
5628 int jv = hex_rot[iv + irot*nv_hex];
5638 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
5639 j = find_min(q, nv_quad);
5640 if (j == 0 || j == 2) { bitmask += 4; }
5642 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
5643 j = find_min(q, nv_quad);
5644 if (j == 1 || j == 3) { bitmask += 2; }
5646 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
5647 j = find_min(q, nv_quad);
5648 if (j == 0 || j == 2) { bitmask += 1; }
5651 int nrot = num_rot[bitmask];
5652 for (
int k=0; k<nrot; ++k)
5666 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
5667 int ntets = (ndiags == 0) ? 5 : 6;
5668 const int *tetmap = hex_tetmaps[ndiags];
5669 for (
int itet=0; itet<ntets; ++itet)
5674 for (
int iv=0; iv<nv_tet; ++iv)
5676 v2[iv] = vg[tetmap[itet + iv*ntets]];
5685 for (
int i=0; i<nbe; ++i)
5687 const int *v = orig_mesh.
boundary[i]->GetVertices();
5690 if (num_subdivisions[orig_geom] == 1)
5700 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
5702 int iv_min = find_min(vg, nv_quad);
5703 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
5704 for (
int itri=0; itri<quad_ntris; ++itri)
5709 for (
int iv=0; iv<nv_tri; ++iv)
5711 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
5718 MFEM_ABORT(
"Unreachable");
5732 Mesh periodic_mesh(orig_mesh,
true);
5738 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
5743 for (
int j = 0; j < nv; j++)
5749 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
5754 for (
int j = 0; j < nv; j++)
5761 return periodic_mesh;
5765 const std::vector<Vector> &translations,
real_t tol)
const
5769 Vector coord(sdim), at(sdim), dx(sdim);
5770 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
5771 xMax = xMin = xDiff = 0.0;
5774 unordered_set<int> bdr_v;
5775 for (
int be = 0; be <
GetNBE(); be++)
5780 for (
int i = 0; i < dofs.
Size(); i++)
5782 bdr_v.insert(dofs[i]);
5785 for (
int j = 0; j < sdim; j++)
5787 xMax[j] = max(xMax[j], coord[j]);
5788 xMin[j] = min(xMin[j], coord[j]);
5792 add(xMax, -1.0, xMin, xDiff);
5801 unordered_map<int, int> replica2primary;
5803 unordered_map<int, unordered_set<int>> primary2replicas;
5806 std::unique_ptr<KDTreeBase<int,real_t>> kdtree;
5807 if (sdim == 1) { kdtree.reset(
new KDTree1D); }
5808 else if (sdim == 2) { kdtree.reset(
new KDTree2D); }
5809 else if (sdim == 3) { kdtree.reset(
new KDTree3D); }
5810 else { MFEM_ABORT(
"Invalid space dimension."); }
5814 for (
const int v : bdr_v)
5816 primary2replicas[v];
5824 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
5826 if (r ==
p) {
return; }
5827 primary2replicas[
p].insert(r);
5828 replica2primary[r] =
p;
5829 for (
const int s : primary2replicas[r])
5831 primary2replicas[
p].insert(s);
5832 replica2primary[s] =
p;
5834 primary2replicas.erase(r);
5837 for (
unsigned int i = 0; i < translations.size(); i++)
5839 for (
int vi : bdr_v)
5842 add(coord, translations[i], at);
5844 const int vj = kdtree->FindClosestPoint(at.GetData());
5846 add(at, -1.0, coord, dx);
5848 if (dx.
Norml2() > dia*tol) {
continue; }
5853 const bool pi = primary2replicas.find(vi) != primary2replicas.end();
5854 const bool pj = primary2replicas.find(vj) != primary2replicas.end();
5860 make_replica(vj, vi);
5865 const int owner_of_vj = replica2primary[vj];
5867 make_replica(vi, owner_of_vj);
5873 const int owner_of_vi = replica2primary[vi];
5874 make_replica(vj, owner_of_vi);
5881 const int owner_of_vi = replica2primary[vi];
5882 const int owner_of_vj = replica2primary[vj];
5883 make_replica(owner_of_vj, owner_of_vi);
5888 std::vector<int> v2v(
GetNV());
5889 for (
size_t i = 0; i < v2v.size(); i++)
5891 v2v[i] =
static_cast<int>(i);
5893 for (
const auto &r2p : replica2primary)
5895 v2v[r2p.first] = r2p.second;
5902 MFEM_VERIFY(
NURBSext,
"Mesh::RefineNURBSFromFile: Not a NURBS mesh!");
5903 mfem::out<<
"Refining NURBS from refinement file: "<<ref_file<<endl;
5906 ifstream input(ref_file);
5913 mfem::out<<
"Knot vectors in ref_file: "<<nkv<<endl;
5915 MFEM_ABORT(
"Refine file does not have the correct number of knot vectors");
5920 for (
int kv = 0; kv < nkv; kv++)
5922 knotVec[kv] =
new Vector();
5923 knotVec[kv]->
Load(input);
5931 for (
int kv = 0; kv < nkv; kv++)
5941 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5946 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5963 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5968 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5985 mfem_error(
"Mesh::KnotRemove : Not a NURBS mesh!");
5990 mfem_error(
"Mesh::KnotRemove : KnotVector array size mismatch!");
6013 "Refinement factors must be defined for each dimension");
6015 MFEM_VERIFY(
NURBSext,
"NURBSUniformRefinement is only for NURBS meshes");
6025 cf1 = (cf1 &&
f == 1);
6042 for (
int i=0; i<cf.
Size(); ++i) { cf[i] *= rf[i]; }
6056 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
6080 for (
int i = 0; i <
elements.Size(); i++)
6090 for (
int i = 0; i <
boundary.Size(); i++)
6108 for (
int i = 0; i < vd; i++)
6175 input >> edge_to_knot[j] >> v[0] >> v[1];
6178 edge_to_knot[j] = -1 - edge_to_knot[j];
6198 if (edge_to_knot.
Size() == 0)
6202 constexpr int notset = -9999999;
6203 edge_to_knot = notset;
6215 edge0[0] = 0; edge1[0] = 2;
6216 edge0[1] = 1; edge1[1] = 3;
6224 edge0[0] = 0; edge1[0] = 2;
6225 edge0[1] = 0; edge1[1] = 4;
6226 edge0[2] = 0; edge1[2] = 6;
6228 edge0[3] = 1; edge1[3] = 3;
6229 edge0[4] = 1; edge1[4] = 5;
6230 edge0[5] = 1; edge1[5] = 7;
6232 edge0[6] = 8; edge1[6] = 9;
6233 edge0[7] = 8; edge1[7] = 10;
6234 edge0[8] = 8; edge1[8] = 11;
6242 int e0, e1, v0, v1, df;
6248 const int *v =
elements[
p]->GetVertices();
6249 for (j = 0; j < edges.
Size(); j++)
6252 const int *e =
elements[
p]->GetEdgeVertices(j);
6265 for (j = 0; j < edge1.
Size(); j++)
6267 e0 = edges[edge0[j]];
6268 e1 = edges[edge1[j]];
6269 v0 = edge_to_knot[e0];
6270 v1 = edge_to_knot[e1];
6271 df = flip*oedge[edge0[j]]*oedge[edge1[j]];
6274 if ((v0 == notset) && (v1 == notset))
6276 edge_to_knot[e0] = knot;
6277 edge_to_knot[e1] = knot;
6283 else if ((v0 != notset) && (v1 == notset))
6285 edge_to_knot[e1] = (df >= 0 ? -v0-1 : v0);
6287 else if ((v0 == notset) && (v1 != notset))
6289 edge_to_knot[e0] = (df >= 0 ? -v1-1 : v1);
6309 for (j = 0; j < edge1.
Size(); j++)
6311 e0 = edges[edge0[j]];
6312 e1 = edges[edge1[j]];
6313 v0 = edge_to_knot[e0];
6314 v1 = edge_to_knot[e1];
6315 v0 = ( v0 >= 0 ? v0 : -v0-1);
6316 v1 = ( v1 >= 0 ? v1 : -v1-1);
6322 edge_to_knot[e1] = (oedge[edge1[j]] >= 0 ? v0 : -v0-1);
6326 edge_to_knot[e0] = (oedge[edge0[j]] >= 0 ? v1 : -v1-1);
6334 while (corrections > 0 && passes <
GetNE() + 1);
6337 if (corrections > 0)
6339 mfem::err<<
"Edge_to_knot mapping potentially incorrect"<<endl;
6341 mfem::err<<
" corrections = "<<corrections<<endl;
6351 k = edge_to_knot[j];
6352 cnt[(k >= 0 ? k : -k-1)]++;
6356 for (j = 0; j < cnt.
Size(); j++)
6358 cnt[j] = (cnt[j] > 0 ? k++ : -1);
6363 k = edge_to_knot[j];
6364 edge_to_knot[j] = (k >= 0 ? cnt[k]:-cnt[-k-1]-1);
6368 mfem::out<<
"Generated edge to knot mapping:"<<endl;
6372 k = edge_to_knot[j];
6381 mfem::out<<(k >= 0 ? k:-k-1)<<
" "<< v[0] <<
" "<<v[1]<<endl;
6385 if (corrections > 0 ) {
mfem_error(
"Mesh::LoadPatchTopo");}
6391 if (
p.Size() >= v.
Size())
6393 for (
int d = 0; d < v.
Size(); d++)
6401 for (d = 0; d <
p.Size(); d++)
6405 for ( ; d < v.
Size(); d++)
6418 nodes.ProjectCoefficient(xyz);
6448 const bool warn =
true;
6451 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
6455 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
6456 " H1-continuous mesh!\n "
6457 "If this is the desired behavior, you can silence"
6458 " this warning by converting\n "
6459 "the NURBS mesh to high-order mesh in advance by"
6460 " calling the method\n "
6461 "Mesh::SetCurvature().");
6492 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
6506 const int old_space_dim =
spaceDim;
6519 MFEM_ASSERT(
nodes != NULL,
"");
6523 nodes->GetNodalValues(vert_val, i+1);
6535 case 1:
return GetNV();
6571#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
6572static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
6577 int i, j, k, wo = 0, fo = 0;
6586 int *vi =
elements[i]->GetVertices();
6589 for (j = 0; j < 3; j++)
6593 for (j = 0; j < 2; j++)
6594 for (k = 0; k < 2; k++)
6596 J(j, k) = v[j+1][k] - v[0][k];
6617 MFEM_ABORT(
"Invalid 2D element type \""
6634 int *vi =
elements[i]->GetVertices();
6640 for (j = 0; j < 4; j++)
6644 for (j = 0; j < 3; j++)
6645 for (k = 0; k < 3; k++)
6647 J(j, k) = v[j+1][k] - v[0][k];
6707 MFEM_ABORT(
"Invalid 3D element type \""
6713#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
6716 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
6717 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
6721 MFEM_CONTRACT_VAR(fo);
6734 if (test[0] == base[0])
6735 if (test[1] == base[1])
6743 else if (test[0] == base[1])
6744 if (test[1] == base[0])
6753 if (test[1] == base[0])
6764 for (
int j = 0; j < 3; j++)
6765 if (test[aor[j]] != base[j])
6767 mfem::err <<
"Mesh::GetTriOrientation(...)" << endl;
6769 for (
int k = 0; k < 3; k++)
6774 for (
int k = 0; k < 3; k++)
6795 const int oo[6][6] =
6805 int ori_a_c = oo[ori_a_b][ori_b_c];
6811 const int inv_ori[6] = {0, 1, 4, 3, 2, 5};
6812 return inv_ori[ori];
6819 for (i = 0; i < 4; i++)
6820 if (test[i] == base[0])
6827 if (test[(i+1)%4] == base[1])
6836 for (
int j = 0; j < 4; j++)
6837 if (test[aor[j]] != base[j])
6839 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
6841 for (
int k = 0; k < 4; k++)
6846 for (
int k = 0; k < 4; k++)
6855 if (test[(i+1)%4] == base[1])
6872 const int oo[8][8] =
6874 {0, 1, 2, 3, 4, 5, 6, 7},
6875 {1, 0, 3, 2, 5, 4, 7, 6},
6876 {2, 7, 4, 1, 6, 3, 0, 5},
6877 {3, 6, 5, 0, 7, 2, 1, 4},
6878 {4, 5, 6, 7, 0, 1, 2, 3},
6879 {5, 4, 7, 6, 1, 0, 3, 2},
6880 {6, 3, 0, 5, 2, 7, 4, 1},
6881 {7, 2, 1, 4, 3, 6, 5, 0}
6884 int ori_a_c = oo[ori_a_b][ori_b_c];
6890 const int inv_ori[8] = {0, 1, 6, 3, 4, 5, 2, 7};
6891 return inv_ori[ori];
6902 if (test[0] == base[0])
6903 if (test[1] == base[1])
6904 if (test[2] == base[2])
6912 else if (test[2] == base[1])
6913 if (test[3] == base[2])
6922 if (test[1] == base[2])
6930 else if (test[1] == base[0])
6931 if (test[2] == base[1])
6932 if (test[0] == base[2])
6940 else if (test[3] == base[1])
6941 if (test[2] == base[2])
6950 if (test[3] == base[2])
6958 else if (test[2] == base[0])
6959 if (test[3] == base[1])
6960 if (test[0] == base[2])
6968 else if (test[0] == base[1])
6969 if (test[1] == base[2])
6978 if (test[3] == base[2])
6987 if (test[0] == base[1])
6988 if (test[2] == base[2])
6996 else if (test[1] == base[1])
6997 if (test[0] == base[2])
7006 if (test[1] == base[2])
7017 for (
int j = 0; j < 4; j++)
7018 if (test[aor[j]] != base[j])
7043 int *bv =
boundary[i]->GetVertices();
7063 if (
faces_info[fi].Elem2No >= 0) {
continue; }
7066 int *bv =
boundary[i]->GetVertices();
7068 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
7069 const int *fv =
faces[fi]->GetVertices();
7086 MFEM_ABORT(
"Invalid 2D boundary element type \""
7087 << bdr_type <<
"\"");
7092 if (orientation % 2 == 0) {
continue; }
7094 if (!fix_it) {
continue; }
7130 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
7148 MFEM_ASSERT(o >= 0 && o < 2,
"Invalid orientation for Geometry::SEGMENT!");
7160 MFEM_ASSERT(o >= 0 && o < 6,
"Invalid orientation for Geometry::TRIANGLE!");
7173 fip.
x = 1.0 - ip.
x - ip.
y;
7178 fip.
x = 1.0 - ip.
x - ip.
y;
7184 fip.
y = 1.0 - ip.
x - ip.
y;
7189 fip.
y = 1.0 - ip.
x - ip.
y;
7194 MFEM_ASSERT(o >= 0 && o < 8,
"Invalid orientation for Geometry::SQUARE!");
7238 MFEM_ABORT(
"Unsupported face geometry for TransformBdrElementToFace!");
7245 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
7256 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
7275 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
7276 "is not generated.");
7279 const int *v =
elements[i]->GetVertices();
7280 const int ne =
elements[i]->GetNEdges();
7282 for (
int j = 0; j < ne; j++)
7284 const int *e =
elements[i]->GetEdgeVertices(j);
7285 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7296 const int *v =
boundary[i]->GetVertices();
7297 cor[0] = (v[0] < v[1]) ? (1) : (-1);
7310 const int *v =
boundary[i]->GetVertices();
7311 const int ne =
boundary[i]->GetNEdges();
7313 for (
int j = 0; j < ne; j++)
7315 const int *e =
boundary[i]->GetEdgeVertices(j);
7316 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7328 const int *v =
faces[i]->GetVertices();
7329 o[0] = (v[0] < v[1]) ? (1) : (-1);
7341 const int *v =
faces[i]->GetVertices();
7342 const int ne =
faces[i]->GetNEdges();
7344 for (
int j = 0; j < ne; j++)
7346 const int *e =
faces[i]->GetEdgeVertices(j);
7347 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7375 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
7422 const int nv =
elements[i]->GetNVertices();
7423 const int *v =
elements[i]->GetVertices();
7424 for (
int j = 0; j < nv; j++)
7434 const int nv =
elements[i]->GetNVertices();
7435 const int *v =
elements[i]->GetVertices();
7436 for (
int j = 0; j < nv; j++)
7455 const int nv =
boundary[i]->GetNVertices();
7456 const int *v =
boundary[i]->GetVertices();
7457 for (
int j = 0; j < nv; j++)
7463 vert_bdr_elem->
MakeJ();
7467 const int nv =
boundary[i]->GetNVertices();
7468 const int *v =
boundary[i]->GetVertices();
7469 for (
int j = 0; j < nv; j++)
7477 return vert_bdr_elem;
7516 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
7520 int n = el_faces.
Size();
7523 for (
int j = 0; j < n; j++)
7527 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
7531 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
7532 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
7549 for (
auto f : elem_faces)
7570 const int *bv =
boundary[i]->GetVertices();
7580 default: MFEM_ABORT(
"invalid geometry");
7589 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
7592 const int *bv =
boundary[bdr_el]->GetVertices();
7600 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
7607 int bdr_el,
int &el,
int &info)
const
7612 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
7615 const int *bv =
boundary[bdr_el]->GetVertices();
7623 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
7656 for (j = 0; j < nv; j++)
7658 pointmat(k, j) =
vertices[v[j]](k);
7673 for (j = 0; j < nv; j++)
7675 pointmat(k, j) =
vertices[v[j]](k);
7687 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
7690 return sqrt(length);
7698 for (
int i = 0; i < elem_array.
Size(); i++)
7703 for (
int i = 0; i < elem_array.
Size(); i++)
7705 const int *v = elem_array[i]->GetVertices();
7706 const int ne = elem_array[i]->GetNEdges();
7707 for (
int j = 0; j < ne; j++)
7709 const int *e = elem_array[i]->GetEdgeVertices(j);
7723 v_to_v.
Push(v[0], v[1]);
7730 const int *v =
elements[i]->GetVertices();
7731 const int ne =
elements[i]->GetNEdges();
7732 for (
int j = 0; j < ne; j++)
7734 const int *e =
elements[i]->GetEdgeVertices(j);
7735 v_to_v.
Push(v[e[0]], v[e[1]]);
7743 int i, NumberOfEdges;
7759 const int *v =
boundary[i]->GetVertices();
7773 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
7777 return NumberOfEdges;
7836 if (
faces[gf] == NULL)
7868 if (
faces[gf] == NULL)
7878 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
7879 "Interior edge found between 2D elements "
7881 <<
" and " << el <<
".");
7882 int *v =
faces[gf]->GetVertices();
7884 if (v[1] == v0 && v[0] == v1)
7888 else if (v[0] == v0 && v[1] == v1)
7898 MFEM_ABORT(
"internal error");
7904 int v0,
int v1,
int v2)
7906 if (
faces[gf] == NULL)
7916 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
7917 "Interior triangular face found connecting elements "
7919 <<
" and " << el <<
".");
7920 int orientation, vv[3] = { v0, v1, v2 };
7927 faces_info[gf].Elem2Inf = 64 * lf + orientation;
7932 int v0,
int v1,
int v2,
int v3)
7944 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
7945 "Interior quadrilateral face found connecting elements "
7947 <<
" and " << el <<
".");
7948 int vv[4] = { v0, v1, v2, v3 };
7967 faces.SetSize(nfaces);
7969 for (
int i = 0; i < nfaces; ++i)
7988 const int ne =
elements[i]->GetNEdges();
7989 for (
int j = 0; j < ne; j++)
7991 const int *e =
elements[i]->GetEdgeVertices(j);
8002 for (
int j = 0; j < 4; j++)
8006 v[fv[0]], v[fv[1]], v[fv[2]]);
8012 for (
int j = 0; j < 2; j++)
8016 v[fv[0]], v[fv[1]], v[fv[2]]);
8018 for (
int j = 2; j < 5; j++)
8022 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8028 for (
int j = 0; j < 1; j++)
8032 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8034 for (
int j = 1; j < 5; j++)
8038 v[fv[0]], v[fv[1]], v[fv[2]]);
8044 for (
int j = 0; j < 6; j++)
8048 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8053 MFEM_ABORT(
"Unexpected type of Element.");
8061 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
8072 nc_faces_info.Reserve(list.masters.Size() + list.slaves.Size());
8079 if (master.index >= nfaces) {
continue; }
8085 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
8086 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
8092 if (slave.index < 0 ||
8093 slave.index >= nfaces ||
8094 slave.master >= nfaces)
8113 list.point_matrices[slave.geom][slave.matrix]));
8122 const int *v =
elements[i]->GetVertices();
8127 for (
int j = 0; j < 4; j++)
8130 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8136 for (
int j = 0; j < 1; j++)
8139 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8141 for (
int j = 1; j < 5; j++)
8144 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8150 for (
int j = 0; j < 2; j++)
8153 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8155 for (
int j = 2; j < 5; j++)
8158 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8166 for (
int j = 0; j < 6; j++)
8169 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8198 for (
int j = 0; j < 4; j++)
8202 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
8208 for (
int j = 0; j < 2; j++)
8212 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
8214 for (
int j = 2; j < 5; j++)
8218 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
8224 for (
int j = 0; j < 1; j++)
8228 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
8230 for (
int j = 1; j < 5; j++)
8234 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
8242 for (
int j = 0; j < 6; j++)
8246 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
8251 MFEM_ABORT(
"Unexpected type of Element.");
8265 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
8270 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
8274 MFEM_ABORT(
"Unexpected type of boundary Element.");
8288void Rotate3(
int &
a,
int &
b,
int &c)
8320 Table *old_elem_vert = NULL;
8331 int *v =
elements[i]->GetVertices();
8333 Rotate3(v[0], v[1], v[2]);
8336 Rotate3(v[1], v[2], v[3]);
8349 int *v =
boundary[i]->GetVertices();
8351 Rotate3(v[0], v[1], v[2]);
8367 delete old_elem_vert;
8383 if (
p[i] < pmin[i]) { pmin[i] =
p[i]; }
8384 if (
p[i] > pmax[i]) { pmax[i] =
p[i]; }
8398 for (
int i =
spaceDim-1; i >= 0; i--)
8400 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
8401 if (idx < 0) { idx = 0; }
8402 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
8403 part = part * nxyz[i] + idx;
8405 partitioning[el] = part;
8408 return partitioning;
8418#ifdef MFEM_USE_METIS
8420 int print_messages = 1;
8423 int init_flag, fin_flag;
8424 MPI_Initialized(&init_flag);
8425 MPI_Finalized(&fin_flag);
8426 if (init_flag && !fin_flag)
8430 if (rank != 0) { print_messages = 0; }
8434 int i, *partitioning;
8444 partitioning[i] = 0;
8451 partitioning[i] = i;
8457#ifndef MFEM_USE_METIS_5
8469 bool freedata =
false;
8471 idx_t *mpartitioning;
8474 if (
sizeof(
idx_t) ==
sizeof(int))
8478 mpartitioning = (
idx_t*) partitioning;
8487 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
8488 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
8489 mpartitioning =
new idx_t[n];
8492#ifndef MFEM_USE_METIS_5
8495 METIS_SetDefaultOptions(options);
8496 options[METIS_OPTION_CONTIG] = 1;
8504 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
8509 if (part_method >= 0 && part_method <= 2)
8511 for (i = 0; i < n; i++)
8517 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
8523 if (part_method == 0 || part_method == 3)
8525#ifndef MFEM_USE_METIS_5
8554 " error in METIS_PartGraphRecursive!");
8561 if (part_method == 1 || part_method == 4)
8563#ifndef MFEM_USE_METIS_5
8592 " error in METIS_PartGraphKway!");
8599 if (part_method == 2 || part_method == 5)
8601#ifndef MFEM_USE_METIS_5
8614 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
8631 " error in METIS_PartGraphKway!");
8639 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
8643 nparts = (int) mparts;
8644 if (mpartitioning != (
idx_t*)partitioning)
8648 partitioning[k] = mpartitioning[k];
8655 delete[] mpartitioning;
8671 auto count_partition_elements = [&]()
8673 for (i = 0; i < nparts; i++)
8681 psize[partitioning[i]].one++;
8685 for (i = 0; i < nparts; i++)
8687 if (psize[i].one == 0) { empty_parts++; }
8691 count_partition_elements();
8699 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
8700 << empty_parts <<
" empty parts!"
8701 <<
" Applying a simple fix ..." << endl;
8706 for (i = nparts-1; i > nparts-1-empty_parts; i--)
8713 for (i = nparts-1; i > nparts-1-empty_parts; i--)
8715 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
8721 partitioning[j] = psize[nparts-1-i].two;
8728 count_partition_elements();
8732 return partitioning;
8736 mfem_error(
"Mesh::GeneratePartitioning(...): "
8737 "MFEM was compiled without Metis.");
8751 int num_elem, *i_elem_elem, *j_elem_elem;
8753 num_elem = elem_elem.
Size();
8754 i_elem_elem = elem_elem.
GetI();
8755 j_elem_elem = elem_elem.
GetJ();
8760 int stack_p, stack_top_p, elem;
8764 for (i = 0; i < num_elem; i++)
8766 if (partitioning[i] > num_part)
8768 num_part = partitioning[i];
8775 for (i = 0; i < num_part; i++)
8782 for (elem = 0; elem < num_elem; elem++)
8784 if (component[elem] >= 0)
8789 component[elem] = num_comp[partitioning[elem]]++;
8791 elem_stack[stack_top_p++] = elem;
8793 for ( ; stack_p < stack_top_p; stack_p++)
8795 i = elem_stack[stack_p];
8796 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
8799 if (partitioning[k] == partitioning[i])
8801 if (component[k] < 0)
8803 component[k] = component[i];
8804 elem_stack[stack_top_p++] = k;
8806 else if (component[k] != component[i])
8818 int i, n_empty, n_mcomp;
8826 n_empty = n_mcomp = 0;
8827 for (i = 0; i < num_comp.
Size(); i++)
8828 if (num_comp[i] == 0)
8832 else if (num_comp[i] > 1)
8839 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
8840 <<
"The following subdomains are empty :\n";
8841 for (i = 0; i < num_comp.
Size(); i++)
8842 if (num_comp[i] == 0)
8850 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
8851 <<
"The following subdomains are NOT connected :\n";
8852 for (i = 0; i < num_comp.
Size(); i++)
8853 if (num_comp[i] > 1)
8859 if (n_empty == 0 && n_mcomp == 0)
8860 mfem::out <<
"Mesh::CheckPartitioning(...) : "
8861 "All subdomains are connected." << endl;
8885 c(0) =
a[0]*
a[3]-
a[1]*
a[2];
8886 c(1) =
a[0]*
b[3]-
a[1]*
b[2]+
b[0]*
a[3]-
b[1]*
a[2];
8887 c(2) =
b[0]*
b[3]-
b[1]*
b[2];
8908 c(0) = (
a[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
8909 a[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
8910 a[2] * (
a[3] *
a[7] -
a[4] *
a[6]));
8912 c(1) = (
b[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
8913 b[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
8914 b[2] * (
a[3] *
a[7] -
a[4] *
a[6]) +
8916 a[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
8917 a[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
8918 a[2] * (
b[3] *
a[7] -
b[4] *
a[6]) +
8920 a[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
8921 a[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
8922 a[2] * (
a[3] *
b[7] -
a[4] *
b[6]));
8924 c(2) = (
a[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
8925 a[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
8926 a[2] * (
b[3] *
b[7] -
b[4] *
b[6]) +
8928 b[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
8929 b[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
8930 b[2] * (
a[3] *
b[7] -
a[4] *
b[6]) +
8932 b[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
8933 b[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
8934 b[2] * (
b[3] *
a[7] -
b[4] *
a[6]));
8936 c(3) = (
b[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
8937 b[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
8938 b[2] * (
b[3] *
b[7] -
b[4] *
b[6]));
8984 real_t a = z(2),
b = z(1), c = z(0);
8992 x(0) = x(1) = -0.5 *
b /
a;
8997 x(0) = -(x(1) = fabs(0.5 * sqrt(D) /
a));
9005 t = -0.5 * (
b + sqrt(D));
9009 t = -0.5 * (
b - sqrt(D));
9023 real_t a = z(2)/z(3),
b = z(1)/z(3), c = z(0)/z(3);
9027 real_t R = (2 *
a *
a *
a - 9 *
a *
b + 27 * c) / 54;
9035 x(0) = x(1) = x(2) = -
a / 3;
9043 x(0) = -2 * sqrtQ -
a / 3;
9044 x(1) = x(2) = sqrtQ -
a / 3;
9048 x(0) = x(1) = - sqrtQ -
a / 3;
9049 x(2) = 2 * sqrtQ -
a / 3;
9056 real_t theta = acos(R / sqrt(Q3));
9059 x0 = A * cos(theta / 3) -
a / 3;
9060 x1 = A * cos((theta + 2.0 * M_PI) / 3) -
a / 3;
9061 x2 = A * cos((theta - 2.0 * M_PI) / 3) -
a / 3;
9086 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
9090 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
9092 x(0) = A + Q / A -
a / 3;
9101 const real_t factor,
const int Dim)
9104 c(0) = c0 * (1.0 - pow(factor, -Dim));
9106 for (
int j = 0; j < nr; j++)
9118 c(0) = c0 * (1.0 - pow(factor, Dim));
9120 for (
int j = 0; j < nr; j++)
9139 const real_t factor = 2.0;
9154 for (
int k = 0; k < nv; k++)
9157 V(j, k) = displacements(v[k]+j*nvs);
9187 for (
int j = 0; j < nv; j++)
9213 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
9216 vertices[i](j) += displacements(j*nv+i);
9224 for (
int i = 0; i < nv; i++)
9227 vert_coord(j*nv+i) =
vertices[i](j);
9233 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
9236 vertices[i](j) = vert_coord(j*nv+i);
9266 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
9283 (*Nodes) += displacements;
9298 node_coord = (*Nodes);
9310 (*Nodes) = node_coord;
9373 for (j = 1; j < n; j++)
9410 int quad_counter = 0;
9432 const int attr =
elements[i]->GetAttribute();
9433 int *v =
elements[i]->GetVertices();
9439 for (
int ei = 0; ei < 3; ei++)
9441 for (
int k = 0; k < 2; k++)
9449 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
9451 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
9453 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
9455 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
9459 const int qe = quad_counter;
9463 for (
int ei = 0; ei < 4; ei++)
9465 for (
int k = 0; k < 2; k++)
9473 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
9475 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
9477 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
9479 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
9483 MFEM_ABORT(
"unknown element type: " << el_type);
9493 const int attr =
boundary[i]->GetAttribute();
9494 int *v =
boundary[i]->GetVertices();
9503 static const real_t A = 0.0, B = 0.5, C = 1.0;
9504 static real_t tri_children[2*3*4] =
9511 static real_t quad_children[2*4*4] =
9525 for (
int i = 0; i <
elements.Size(); i++)
9546 if (!
Nodes || update_nodes)
9576 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
9579 int NumOfQuadFaces = 0;
9585 for (
int i = 0; i <
faces.Size(); i++)
9589 f2qf[i] = NumOfQuadFaces;
9600 int hex_counter = 0;
9603 for (
int i = 0; i <
elements.Size(); i++)
9612 int pyr_counter = 0;
9615 for (
int i = 0; i <
elements.Size(); i++)
9633 DSTable *v_to_v_ptr = v_to_v_p;
9649 std::sort(row_start, J_v2v.
end());
9652 for (
int i = 0; i < J_v2v.
Size(); i++)
9654 e2v[J_v2v[i].two] = i;
9667 it.SetIndex(e2v[it.Index()]);
9677 const int oelem = oface + NumOfQuadFaces;
9690 const int attr =
elements[i]->GetAttribute();
9691 int *v =
elements[i]->GetVertices();
9698 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
9706 for (
int ei = 0; ei < 6; ei++)
9708 for (
int k = 0; k < 2; k++)
9718 const int rt_algo = 1;
9731 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
9732 sqr(J(1,0)-J(1,1)-J(1,2)) +
9733 sqr(J(2,0)-J(2,1)-J(2,2));
9736 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
9737 sqr(J(1,1)-J(1,0)-J(1,2)) +
9738 sqr(J(2,1)-J(2,0)-J(2,2));
9739 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
9741 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
9742 sqr(J(1,2)-J(1,0)-J(1,1)) +
9743 sqr(J(2,2)-J(2,0)-J(2,1));
9744 if (len_sqr < min_len) { rt = 2; }
9749 real_t Em_data[18], Js_data[9], Jp_data[9];
9754 for (
int s = 0; s < 3; s++)
9756 for (
int t = 0; t < 3; t++)
9758 Em(t,s) = 0.5*J(t,s);
9761 for (
int t = 0; t < 3; t++)
9763 Em(t,3) = 0.5*(J(t,0)+J(t,1));
9764 Em(t,4) = 0.5*(J(t,0)+J(t,2));
9765 Em(t,5) = 0.5*(J(t,1)+J(t,2));
9769 for (
int t = 0; t < 3; t++)
9771 Js(t,0) = Em(t,5)-Em(t,0);
9772 Js(t,1) = Em(t,1)-Em(t,0);
9773 Js(t,2) = Em(t,2)-Em(t,0);
9777 for (
int t = 0; t < 3; t++)
9779 Js(t,0) = Em(t,5)-Em(t,0);
9780 Js(t,1) = Em(t,2)-Em(t,0);
9781 Js(t,2) = Em(t,4)-Em(t,0);
9785 kappa_min = std::max(ar1, ar2);
9789 for (
int t = 0; t < 3; t++)
9791 Js(t,0) = Em(t,0)-Em(t,1);
9792 Js(t,1) = Em(t,4)-Em(t,1);
9793 Js(t,2) = Em(t,2)-Em(t,1);
9797 for (
int t = 0; t < 3; t++)
9799 Js(t,0) = Em(t,2)-Em(t,1);
9800 Js(t,1) = Em(t,4)-Em(t,1);
9801 Js(t,2) = Em(t,5)-Em(t,1);
9805 kappa = std::max(ar1, ar2);
9806 if (
kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
9809 for (
int t = 0; t < 3; t++)
9811 Js(t,0) = Em(t,0)-Em(t,2);
9812 Js(t,1) = Em(t,1)-Em(t,2);
9813 Js(t,2) = Em(t,3)-Em(t,2);
9817 for (
int t = 0; t < 3; t++)
9819 Js(t,0) = Em(t,1)-Em(t,2);
9820 Js(t,1) = Em(t,5)-Em(t,2);
9821 Js(t,2) = Em(t,3)-Em(t,2);
9825 kappa = std::max(ar1, ar2);
9826 if (
kappa < kappa_min) { rt = 2; }
9829 static const int mv_all[3][4][4] =
9831 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
9832 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
9833 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
9835 const int (&mv)[4][4] = mv_all[rt];
9837#ifndef MFEM_USE_MEMALLOC
9839 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
9841 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
9843 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
9845 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
9847 for (
int k = 0; k < 4; k++)
9849 new_elements[j+4+k] =
9850 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
9851 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
9855 new_elements[j+0] = tet =
TetMemory.Alloc();
9856 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
9858 new_elements[j+1] = tet =
TetMemory.Alloc();
9859 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
9861 new_elements[j+2] = tet =
TetMemory.Alloc();
9862 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
9864 new_elements[j+3] = tet =
TetMemory.Alloc();
9865 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
9867 for (
int k = 0; k < 4; k++)
9869 new_elements[j+4+k] = tet =
TetMemory.Alloc();
9870 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
9871 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
9874 for (
int k = 0; k < 4; k++)
9879 for (
int k = 0; k < 4; k++)
9893 for (
int fi = 2; fi < 5; fi++)
9895 for (
int k = 0; k < 4; k++)
9902 for (
int ei = 0; ei < 9; ei++)
9904 for (
int k = 0; k < 2; k++)
9911 const int qf2 = f2qf[
f[2]];
9912 const int qf3 = f2qf[
f[3]];
9913 const int qf4 = f2qf[
f[4]];
9916 new Wedge(v[0], oedge+e[0], oedge+e[2],
9917 oedge+e[6], oface+qf2, oface+qf4, attr);
9920 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
9921 oface+qf3, oface+qf4, oface+qf2, attr);
9924 new Wedge(oedge+e[0], v[1], oedge+e[1],
9925 oface+qf2, oedge+e[7], oface+qf3, attr);
9928 new Wedge(oedge+e[2], oedge+e[1], v[2],
9929 oface+qf4, oface+qf3, oedge+e[8], attr);
9932 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
9933 v[3], oedge+e[3], oedge+e[5], attr);
9936 new Wedge(oface+qf3, oface+qf4, oface+qf2,
9937 oedge+e[4], oedge+e[5], oedge+e[3], attr);
9940 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
9941 oedge+e[3], v[4], oedge+e[4], attr);
9944 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
9945 oedge+e[5], oedge+e[4], v[5], attr);
9954 for (
int fi = 0; fi < 1; fi++)
9956 for (
int k = 0; k < 4; k++)
9963 for (
int ei = 0; ei < 8; ei++)
9965 for (
int k = 0; k < 2; k++)
9972 const int qf0 = f2qf[
f[0]];
9975 new Pyramid(v[0], oedge+e[0], oface+qf0,
9976 oedge+e[3], oedge+e[4], attr);
9979 new Pyramid(oedge+e[0], v[1], oedge+e[1],
9980 oface+qf0, oedge+e[5], attr);
9983 new Pyramid(oface+qf0, oedge+e[1], v[2],
9984 oedge+e[2], oedge+e[6], attr);
9987 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
9988 v[3], oedge+e[7], attr);
9991 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
9992 oedge+e[7], v[4], attr);
9995 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
9996 oedge+e[4], oface+qf0, attr);
9998#ifndef MFEM_USE_MEMALLOC
10000 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
10003 new_elements[j++] =
10004 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
10007 new_elements[j++] =
10008 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
10011 new_elements[j++] =
10012 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
10016 new_elements[j++] = tet =
TetMemory.Alloc();
10017 tet->
Init(oedge+e[0], oedge+e[4], oedge+e[5],
10020 new_elements[j++] = tet =
TetMemory.Alloc();
10021 tet->
Init(oedge+e[1], oedge+e[5], oedge+e[6],
10024 new_elements[j++] = tet =
TetMemory.Alloc();
10025 tet->
Init(oedge+e[2], oedge+e[6], oedge+e[7],
10028 new_elements[j++] = tet =
TetMemory.Alloc();
10029 tet->
Init(oedge+e[3], oedge+e[7], oedge+e[4],
10042 const int he = hex_counter;
10047 if (f2qf.
Size() == 0)
10053 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[
f[k]]; }
10059 for (
int fi = 0; fi < 6; fi++)
10061 for (
int k = 0; k < 4; k++)
10068 for (
int ei = 0; ei < 12; ei++)
10070 for (
int k = 0; k < 2; k++)
10077 new_elements[j++] =
10078 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
10079 oedge+e[3], oedge+e[8], oface+qf[1],
10080 oelem+he, oface+qf[4], attr);
10081 new_elements[j++] =
10082 new Hexahedron(oedge+e[0], v[1], oedge+e[1],
10083 oface+qf[0], oface+qf[1], oedge+e[9],
10084 oface+qf[2], oelem+he, attr);
10085 new_elements[j++] =
10086 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
10087 oedge+e[2], oelem+he, oface+qf[2],
10088 oedge+e[10], oface+qf[3], attr);
10089 new_elements[j++] =
10090 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
10091 v[3], oface+qf[4], oelem+he,
10092 oface+qf[3], oedge+e[11], attr);
10093 new_elements[j++] =
10094 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
10095 oface+qf[4], v[4], oedge+e[4],
10096 oface+qf[5], oedge+e[7], attr);
10097 new_elements[j++] =
10098 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
10099 oelem+he, oedge+e[4], v[5],
10100 oedge+e[5], oface+qf[5], attr);
10101 new_elements[j++] =
10102 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
10103 oface+qf[3], oface+qf[5], oedge+e[5],
10104 v[6], oedge+e[6], attr);
10105 new_elements[j++] =
10106 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
10107 oedge+e[11], oedge+e[7], oface+qf[5],
10108 oedge+e[6], v[7], attr);
10113 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
10125 const int attr =
boundary[i]->GetAttribute();
10126 int *v =
boundary[i]->GetVertices();
10133 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
10139 new_boundary[j++] =
10140 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
10141 new_boundary[j++] =
10142 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
10143 new_boundary[j++] =
10144 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
10145 new_boundary[j++] =
10146 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
10153 new_boundary[j++] =
10154 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
10155 new_boundary[j++] =
10156 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
10157 new_boundary[j++] =
10158 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
10159 new_boundary[j++] =
10160 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
10164 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
10170 static const real_t A = 0.0, B = 0.5, C = 1.0, D = -1.0;
10171 static real_t tet_children[3*4*16] =
10173 A,A,A, B,A,A, A,B,A, A,A,B,
10174 B,A,A, C,A,A, B,B,A, B,A,B,
10175 A,B,A, B,B,A, A,C,A, A,B,B,
10176 A,A,B, B,A,B, A,B,B, A,A,C,
10181 B,A,A, A,B,B, A,B,A, A,A,B,
10182 B,A,A, A,B,B, A,A,B, B,A,B,
10183 B,A,A, A,B,B, B,A,B, B,B,A,
10184 B,A,A, A,B,B, B,B,A, A,B,A,
10186 A,B,A, B,A,A, B,A,B, A,A,B,
10187 A,B,A, A,A,B, B,A,B, A,B,B,
10188 A,B,A, A,B,B, B,A,B, B,B,A,
10189 A,B,A, B,B,A, B,A,B, B,A,A,
10191 A,A,B, B,A,A, A,B,A, B,B,A,
10192 A,A,B, A,B,A, A,B,B, B,B,A,
10193 A,A,B, A,B,B, B,A,B, B,B,A,
10194 A,A,B, B,A,B, B,A,A, B,B,A
10196 static real_t pyr_children[3*5*10] =
10198 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
10199 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
10200 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
10201 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
10202 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
10203 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
10204 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
10205 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
10206 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
10207 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
10209 static real_t pri_children[3*6*8] =
10211 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
10212 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
10213 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
10214 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
10215 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
10216 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
10217 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
10218 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
10220 static real_t hex_children[3*8*8] =
10222 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
10223 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
10224 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
10225 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
10226 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
10227 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
10228 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
10229 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
10241 for (
int i = 0; i <
elements.Size(); i++)
10272 int i, j, ind, nedges;
10279 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
10293 for (j = 0; j < marked_el.
Size(); j++)
10298 int new_v = cnv + j, new_e = cne + j;
10307 static real_t seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
10309 UseExternalData(seg_children, 1, 2, 3);
10322 int *edge1 =
new int[nedges];
10323 int *edge2 =
new int[nedges];
10324 int *middle =
new int[nedges];
10326 for (i = 0; i < nedges; i++)
10328 edge1[i] = edge2[i] = middle[i] = -1;
10334 for (j = 1; j < v.
Size(); j++)
10336 ind = v_to_v(v[j-1], v[j]);
10337 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
10339 ind = v_to_v(v[0], v[v.
Size()-1]);
10340 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
10344 for (i = 0; i < marked_el.
Size(); i++)
10350 int need_refinement;
10353 need_refinement = 0;
10354 for (i = 0; i < nedges; i++)
10356 if (middle[i] != -1 && edge1[i] != -1)
10358 need_refinement = 1;
10363 while (need_refinement == 1);
10366 int v1[2], v2[2],
bisect, temp;
10368 for (i = 0; i < temp; i++)
10371 bisect = v_to_v(v[0], v[1]);
10372 if (middle[
bisect] != -1)
10376 v1[0] = v[0]; v1[1] = middle[
bisect];
10377 v2[0] = middle[
bisect]; v2[1] = v[1];
10383 mfem_error(
"Only bisection of segment is implemented"
10407 MFEM_VERIFY(
GetNE() == 0 ||
10409 "tetrahedral mesh is not marked for refinement:"
10410 " call Finalize(true)");
10417 for (i = 0; i < marked_el.
Size(); i++)
10423 for (i = 0; i < marked_el.
Size(); i++)
10432 for (i = 0; i < marked_el.
Size(); i++)
10449 int need_refinement;
10454 need_refinement = 0;
10462 if (
elements[i]->NeedRefinement(v_to_v))
10464 need_refinement = 1;
10469 while (need_refinement == 1);
10476 need_refinement = 0;
10478 if (
boundary[i]->NeedRefinement(v_to_v))
10480 need_refinement = 1;
10484 while (need_refinement == 1);
10515 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
10516 "not supported. Project the NURBS to Nodes first.");
10526 if (!refinements.
Size())
10547 Swap(*mesh2,
false);
10559 const int *fine,
int nfine,
int op)
10561 real_t error = (op == 3) ? std::pow(elem_error[fine[0]],
10562 2.0) : elem_error[fine[0]];
10564 for (
int i = 1; i < nfine; i++)
10566 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
10568 real_t err_fine = elem_error[fine[i]];
10571 case 0: error = std::min(error, err_fine);
break;
10572 case 1: error += err_fine;
break;
10573 case 2: error = std::max(error, err_fine);
break;
10574 case 3: error += std::pow(err_fine, 2.0);
break;
10575 default: MFEM_ABORT(
"Invalid operation.");
10578 return (op == 3) ? std::sqrt(error) : error;
10582 real_t threshold,
int nc_limit,
int op)
10584 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
10585 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
10586 "Project the NURBS to Nodes first.");
10599 for (
int i = 0; i < dt.
Size(); i++)
10601 if (nc_limit > 0 && !level_ok[i]) {
continue; }
10606 if (error < threshold) { derefs.
Append(i); }
10609 if (!derefs.
Size()) {
return false; }
10616 Swap(*mesh2,
false);
10630 int nc_limit,
int op)
10640 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
10647 int nc_limit,
int op)
10650 for (
int i = 0; i < tmp.
Size(); i++)
10652 tmp[i] = elem_error(i);
10695 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
10741#ifdef MFEM_USE_MEMALLOC
10769 for (
int i = 0; i < elem_array.
Size(); i++)
10771 if (elem_array[i]->GetGeometryType() == geom)
10776 elem_vtx.
SetSize(nv*num_elems);
10780 for (
int i = 0; i < elem_array.
Size(); i++)
10786 elem_vtx.
Append(loc_vtx);
10794 for (
int i = 0; i < nelem; i++) { list[i] = i; }
10810 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
10822 default: MFEM_ABORT(
"internal error");
10836 bool noInitialCoarsening =
true;
10837 for (
auto f : initialCoarsening)
10839 noInitialCoarsening = (noInitialCoarsening &&
f == 1);
10842 if (noInitialCoarsening)
10862 bool divisible =
true;
10863 for (
int i=0; i<rf.
Size(); ++i)
10866 divisible = divisible && cf * rf[i] == initialCoarsening[i];
10869 MFEM_VERIFY(divisible,
"Invalid coarsening");
10884 int nonconforming,
int nc_limit)
10894 else if (nonconforming < 0)
10915 for (
int i = 0; i < refinements.
Size(); i++)
10917 el_to_refine[i] = refinements[i].index;
10921 int type, rt = (refinements.
Size() ? refinements[0].GetType() : 7);
10922 if (rt == 1 || rt == 2 || rt == 4)
10926 else if (rt == 3 || rt == 5 || rt == 6)
10944 for (
int i = 0; i < el_to_refine.
Size(); i++)
10946 refinements[i] =
Refinement(el_to_refine[i]);
10953 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
10954 "Please project the NURBS to Nodes first, with SetCurvature().");
10957 MFEM_VERIFY(
ncmesh != NULL ||
dynamic_cast<const ParMesh*
>(
this) == NULL,
10958 "Sorry, converting a conforming ParMesh to an NC mesh is "
10966 (simplices_nonconforming && (
meshgen & 0x1)) )
10979 for (
int i = 0; i <
GetNE(); i++)
10986 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
10998 for (
int i = 0; i <
GetNE(); i++)
11001 bool refine =
false;
11002 for (
int j = 0; j < v.
Size(); j++)
11005 for (
int l = 0; l <
spaceDim; l++)
11010 if (dist <= eps*eps) { refine =
true;
break; }
11021 int nonconforming,
int nc_limit)
11023 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
11025 for (
int i = 0; i <
GetNE(); i++)
11027 if (elem_error[i] > threshold)
11041 int nonconforming,
int nc_limit)
11044 elem_error.
Size());
11045 return RefineByError(tmp, threshold, nonconforming, nc_limit);
11050 int *edge1,
int *edge2,
int *middle)
11053 int v[2][4], v_new,
bisect, t;
11065 bisect = v_to_v(vert[0], vert[1]);
11066 MFEM_ASSERT(
bisect >= 0,
"");
11068 if (middle[
bisect] == -1)
11071 for (
int d = 0; d <
spaceDim; d++)
11096 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
11097 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
11118 bisect = v_to_v(v[1][0], v[1][1]);
11119 MFEM_ASSERT(
bisect >= 0,
"");
11125 else if (edge2[
bisect] == i)
11134 MFEM_ABORT(
"Bisection for now works only for triangles.");
11141 int v[2][4], v_new,
bisect, t;
11151 "TETRAHEDRON element is not marked for refinement.");
11160 for (
int j = 0; j < 3; j++)
11173 int type, old_redges[2], flag;
11176 int new_type, new_redges[2][2];
11179 new_redges[0][0] = 2;
11180 new_redges[0][1] = 1;
11181 new_redges[1][0] = 2;
11182 new_redges[1][1] = 1;
11183 int tr1 = -1, tr2 = -1;
11184 switch (old_redges[0])
11187 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
11192 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
11196 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
11199 switch (old_redges[1])
11202 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
11207 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
11211 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
11218#ifdef MFEM_USE_MEMALLOC
11254 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
11261 int v[2][3], v_new,
bisect, t;
11273 MFEM_ASSERT(
bisect >= 0,
"");
11275 MFEM_ASSERT(v_new != -1,
"");
11279 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
11280 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
11290 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
11296 int *edge1,
int *edge2,
int *middle)
11299 int j, v1[3], v2[3], v3[3], v4[3], v_new[3],
bisect[3];
11308 bisect[0] = v_to_v(v[0],v[1]);
11309 bisect[1] = v_to_v(v[1],v[2]);
11310 bisect[2] = v_to_v(v[0],v[2]);
11313 for (j = 0; j < 3; j++)
11315 if (middle[
bisect[j]] == -1)
11318 for (
int d = 0; d <
spaceDim; d++)
11326 if (edge1[
bisect[j]] == i)
11331 middle[
bisect[j]] = v_new[j];
11335 v_new[j] = middle[
bisect[j]];
11344 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
11345 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
11346 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
11347 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
11381 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
11417 for (
int i = 0; i < elem_geoms.
Size(); i++)
11425 std::map<unsigned, int> mat_no;
11429 for (
int j = 0; j <
elements.Size(); j++)
11432 unsigned code =
elements[j]->GetTransform();
11435 int &matrix = mat_no[code];
11436 if (!matrix) { matrix =
static_cast<int>(mat_no.size()); }
11443 pmats.
SetSize(
Dim,
Dim+1,
static_cast<int>((mat_no.size())));
11446 std::map<unsigned, int>::iterator it;
11447 for (it = mat_no.begin(); it != mat_no.end(); ++it)
11461 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
11462 " geom = " << geom);
11472 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
11481 os <<
"areamesh2\n\n";
11485 os <<
"curved_areamesh2\n\n";
11494 os <<
boundary[i]->GetAttribute();
11495 for (j = 0; j < v.
Size(); j++)
11497 os <<
' ' << v[j] + 1;
11509 for (j = 0; j < v.
Size(); j++)
11511 os <<
' ' << v[j] + 1;
11523 for (j = 1; j <
Dim; j++)
11540 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
11548 os <<
"NETGEN_Neutral_Format\n";
11553 for (j = 0; j <
Dim; j++)
11566 os <<
elements[i]->GetAttribute();
11567 for (j = 0; j < nv; j++)
11569 os <<
' ' << ind[j]+1;
11580 os <<
boundary[i]->GetAttribute();
11581 for (j = 0; j < nv; j++)
11583 os <<
' ' << ind[j]+1;
11595 <<
" 0 0 0 0 0 0 0\n"
11596 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
11598 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
11599 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11603 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
11609 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
11610 for (j = 0; j < nv; j++)
11612 os <<
' ' << ind[j]+1;
11621 os <<
boundary[i]->GetAttribute();
11622 for (j = 0; j < nv; j++)
11624 os <<
' ' << ind[j]+1;
11626 os <<
" 1.0 1.0 1.0 1.0\n";
11635 const std::string &comments)
const
11669 os <<
"\n# mesh curvature GridFunction";
11674 os <<
"\nmfem_mesh_end" << endl;
11681 os << (!set_names && section_delimiter.empty()
11682 ?
"MFEM mesh v1.0\n" :
11683 (!set_names ?
"MFEM mesh v1.2\n" :
"MFEM mesh v1.3\n"));
11685 if (set_names && section_delimiter.empty())
11687 section_delimiter =
"mfem_mesh_end";
11691 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
11694 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
11699 "# TETRAHEDRON = 4\n"
11705 os <<
"\ndimension\n" <<
Dim;
11715 os <<
"\nattribute_sets\n";
11727 os <<
"\nbdr_attribute_sets\n";
11752 if (!section_delimiter.empty())
11755 << section_delimiter << endl;
11760 const int version,
const std::string &comments)
const
11762 MFEM_VERIFY(version == 10 || version == 11,
"Invalid NURBS mesh version");
11767 os <<
"MFEM NURBS mesh v" << int(version / 10) <<
"." << version % 10 <<
"\n";
11770 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
11773 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
11779 os <<
"\ndimension\n" <<
Dim
11796 int ki = e_to_k[i];
11801 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
11808 ofstream ofs(fname);
11809 ofs.precision(precision);
11813#ifdef MFEM_USE_ADIOS2
11823 "# vtk DataFile Version 3.0\n"
11824 "Generated by MFEM\n"
11826 "DATASET UNSTRUCTURED_GRID\n";
11839 for ( ; j < 3; j++)
11850 for (
int i = 0; i <
Nodes->
FESpace()->GetNDofs(); i++)
11855 os << (*Nodes)(vdofs[0]);
11859 os <<
' ' << (*Nodes)(vdofs[j]);
11861 for ( ; j < 3; j++)
11875 size +=
elements[i]->GetNVertices() + 1;
11880 const int *v =
elements[i]->GetVertices();
11881 const int nv =
elements[i]->GetNVertices();
11885 for (
int j = 0; j < nv; j++)
11887 os <<
' ' << v[perm ? perm[j] : j];
11900 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
11901 "Point meshes should have a single dof per element");
11902 size += dofs.
Size() + 1;
11907 if (!strcmp(fec_name,
"Linear") ||
11908 !strcmp(fec_name,
"H1_0D_P1") ||
11909 !strcmp(fec_name,
"H1_1D_P1") ||
11910 !strcmp(fec_name,
"H1_2D_P1") ||
11911 !strcmp(fec_name,
"H1_3D_P1"))
11915 else if (!strcmp(fec_name,
"Quadratic") ||
11916 !strcmp(fec_name,
"H1_1D_P2") ||
11917 !strcmp(fec_name,
"H1_2D_P2") ||
11918 !strcmp(fec_name,
"H1_3D_P2"))
11924 mfem::err <<
"Mesh::PrintVTK : can not save '"
11925 << fec_name <<
"' elements!" << endl;
11934 for (
int j = 0; j < dofs.
Size(); j++)
11936 os <<
' ' << dofs[j];
11939 else if (order == 2)
11941 const int *vtk_mfem;
11942 switch (
elements[i]->GetGeometryType())
11956 for (
int j = 0; j < dofs.
Size(); j++)
11958 os <<
' ' << dofs[vtk_mfem[j]];
11968 int vtk_cell_type = 5;
11972 os << vtk_cell_type <<
'\n';
11977 <<
"SCALARS material int\n"
11978 <<
"LOOKUP_TABLE default\n";
11981 os <<
elements[i]->GetAttribute() <<
'\n';
11988 bool high_order_output,
11989 int compression_level,
11992 int ref = (high_order_output &&
Nodes)
11995 fname = fname +
".vtu";
11996 std::fstream os(fname.c_str(),std::ios::out);
11997 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
11998 if (compression_level != 0)
12000 os <<
" compressor=\"vtkZLibDataCompressor\"";
12003 os <<
"<UnstructuredGrid>\n";
12004 PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
12005 os <<
"</Piece>\n";
12006 os <<
"</UnstructuredGrid>\n";
12007 os <<
"</VTKFile>" << std::endl;
12014 bool high_order_output,
12015 int compression_level)
12017 PrintVTU(fname, format, high_order_output, compression_level,
true);
12021 bool high_order_output,
int compression_level,
12029 std::vector<char> buf;
12031 auto get_geom = [&](
int i)
12039 int np = 0, nc_ref = 0;
12040 for (
int i = 0; i < ne; i++)
12049 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
12050 << (high_order_output ? ne : nc_ref) <<
"\">\n";
12053 os <<
"<Points>\n";
12054 os <<
"<DataArray type=\"" << type_str
12055 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
12056 for (
int i = 0; i < ne; i++)
12069 for (
int j = 0; j < pmat.
Width(); j++)
12095 os <<
"</DataArray>" << std::endl;
12096 os <<
"</Points>" << std::endl;
12098 os <<
"<Cells>" << std::endl;
12099 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
12100 << fmt_str <<
"\">" << std::endl;
12102 std::vector<int> offset;
12105 if (high_order_output)
12108 for (
int iel = 0; iel < ne; iel++)
12112 int nnodes = local_connectivity.
Size();
12113 for (
int i=0; i<nnodes; ++i)
12120 offset.push_back(np);
12126 for (
int i = 0; i < ne; i++)
12132 for (
int j = 0; j < RG.
Size(); )
12135 offset.push_back(coff);
12137 for (
int k = 0; k < nv; k++, j++)
12151 os <<
"</DataArray>" << std::endl;
12153 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
12154 << fmt_str <<
"\">" << std::endl;
12156 for (
size_t ii=0; ii<offset.size(); ii++)
12164 os <<
"</DataArray>" << std::endl;
12165 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
12166 << fmt_str <<
"\">" << std::endl;
12168 const int *vtk_geom_map =
12170 for (
int i = 0; i < ne; i++)
12173 uint8_t vtk_cell_type = 5;
12175 vtk_cell_type = vtk_geom_map[geom];
12177 if (high_order_output)
12186 for (
int j = 0; j < RG.
Size(); j += nv)
12196 os <<
"</DataArray>" << std::endl;
12197 os <<
"</Cells>" << std::endl;
12199 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
12200 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
12201 << fmt_str <<
"\">" << std::endl;
12202 for (
int i = 0; i < ne; i++)
12205 if (high_order_output)
12224 os <<
"</DataArray>" << std::endl;
12225 os <<
"</CellData>" << std::endl;
12236 "# vtk DataFile Version 3.0\n"
12237 "Generated by MFEM\n"
12239 "DATASET UNSTRUCTURED_GRID\n";
12244 os <<
"FIELD FieldData 1\n"
12254 np = nc = size = 0;
12255 for (
int i = 0; i <
GetNE(); i++)
12264 os <<
"POINTS " << np <<
" double\n";
12266 for (
int i = 0; i <
GetNE(); i++)
12273 for (
int j = 0; j < pmat.
Width(); j++)
12275 os << pmat(0, j) <<
' ';
12278 os << pmat(1, j) <<
' ';
12290 os << 0.0 <<
' ' << 0.0;
12297 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
12299 for (
int i = 0; i <
GetNE(); i++)
12306 for (
int j = 0; j < RG.
Size(); )
12309 for (
int k = 0; k < nv; k++, j++)
12311 os <<
' ' << np + RG[j];
12317 os <<
"CELL_TYPES " << nc <<
'\n';
12318 for (
int i = 0; i <
GetNE(); i++)
12326 for (
int j = 0; j < RG.
Size(); j += nv)
12328 os << vtk_cell_type <<
'\n';
12332 os <<
"CELL_DATA " << nc <<
'\n'
12333 <<
"SCALARS material int\n"
12334 <<
"LOOKUP_TABLE default\n";
12335 for (
int i = 0; i <
GetNE(); i++)
12343 os << attr <<
'\n';
12350 srand((
unsigned)time(0));
12352 int el0 = (int)floor(
a *
GetNE());
12354 os <<
"SCALARS element_coloring int\n"
12355 <<
"LOOKUP_TABLE default\n";
12356 for (
int i = 0; i <
GetNE(); i++)
12363 os << coloring[i] + 1 <<
'\n';
12369 os <<
"POINT_DATA " << np <<
'\n' << flush;
12374 int delete_el_to_el = (
el_to_el) ? (0) : (1);
12376 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
12379 const int *i_el_el = el_el.
GetI();
12380 const int *j_el_el = el_el.
GetJ();
12385 stack_p = stack_top_p = 0;
12386 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
12388 if (colors[el] != -2)
12394 el_stack[stack_top_p++] = el;
12396 for ( ; stack_p < stack_top_p; stack_p++)
12398 int i = el_stack[stack_p];
12399 int num_nb = i_el_el[i+1] - i_el_el[i];
12400 if (max_num_col < num_nb + 1)
12402 max_num_col = num_nb + 1;
12404 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
12406 int k = j_el_el[j];
12407 if (colors[k] == -2)
12410 el_stack[stack_top_p++] = k;
12418 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
12420 int i = el_stack[stack_p], col;
12422 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
12424 col = colors[j_el_el[j]];
12427 col_marker[col] = 1;
12431 for (col = 0; col < max_num_col; col++)
12432 if (col_marker[col] == 0)
12440 if (delete_el_to_el)
12448 int elem_attr)
const
12450 if (
Dim != 3 &&
Dim != 2) {
return; }
12452 int i, j, k, l, nv, nbe, *v;
12454 os <<
"MFEM mesh v1.0\n";
12458 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12463 "# TETRAHEDRON = 4\n"
12468 os <<
"\ndimension\n" <<
Dim
12472 os << int((elem_attr) ? partitioning[i]+1 :
elements[i]->GetAttribute())
12473 <<
' ' <<
elements[i]->GetGeometryType();
12476 for (j = 0; j < nv; j++)
12488 l = partitioning[l];
12503 os <<
"\nboundary\n" << nbe <<
'\n';
12509 l = partitioning[l];
12512 nv =
faces[i]->GetNVertices();
12513 v =
faces[i]->GetVertices();
12514 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
12515 for (j = 0; j < nv; j++)
12522 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
12523 for (j = nv-1; j >= 0; j--)
12534 nv =
faces[i]->GetNVertices();
12535 v =
faces[i]->GetVertices();
12536 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
12537 for (j = 0; j < nv; j++)
12568 int interior_faces)
12570 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
12571 if (
Dim != 3 &&
Dim != 2) {
return; }
12580 int nv =
elements[i]->GetNVertices();
12581 const int *ind =
elements[i]->GetVertices();
12582 for (
int j = 0; j < nv; j++)
12592 voff[i] = vcount[i-1] + voff[i-1];
12598 vown[i] =
new int[vcount[i]];
12610 int nv =
elements[i]->GetNVertices();
12611 const int *ind =
elements[i]->GetVertices();
12612 for (
int j = 0; j < nv; j++)
12615 vown[ind[j]][vcount[ind[j]]] = i;
12621 vcount[i] = voff[i+1] - voff[i];
12625 for (
int i = 0; i < edge_el.
Size(); i++)
12627 const int *el = edge_el.
GetRow(i);
12630 int k = partitioning[el[0]];
12631 int l = partitioning[el[1]];
12632 if (interior_faces || k != l)
12644 os <<
"areamesh2\n\n" << nbe <<
'\n';
12646 for (
int i = 0; i < edge_el.
Size(); i++)
12648 const int *el = edge_el.
GetRow(i);
12651 int k = partitioning[el[0]];
12652 int l = partitioning[el[1]];
12653 if (interior_faces || k != l)
12658 for (
int j = 0; j < 2; j++)
12659 for (
int s = 0; s < vcount[ev[j]]; s++)
12660 if (vown[ev[j]][s] == el[0])
12662 os <<
' ' << voff[ev[j]]+s+1;
12666 for (
int j = 1; j >= 0; j--)
12667 for (
int s = 0; s < vcount[ev[j]]; s++)
12668 if (vown[ev[j]][s] == el[1])
12670 os <<
' ' << voff[ev[j]]+s+1;
12677 int k = partitioning[el[0]];
12681 for (
int j = 0; j < 2; j++)
12682 for (
int s = 0; s < vcount[ev[j]]; s++)
12683 if (vown[ev[j]][s] == el[0])
12685 os <<
' ' << voff[ev[j]]+s+1;
12695 int nv =
elements[i]->GetNVertices();
12696 const int *ind =
elements[i]->GetVertices();
12697 os << partitioning[i]+1 <<
' ';
12699 for (
int j = 0; j < nv; j++)
12701 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
12702 vown[ind[j]][vcount[ind[j]]] = i;
12709 vcount[i] = voff[i+1] - voff[i];
12715 for (
int k = 0; k < vcount[i]; k++)
12717 for (
int j = 0; j <
Dim; j++)
12727 os <<
"NETGEN_Neutral_Format\n";
12731 for (
int k = 0; k < vcount[i]; k++)
12733 for (
int j = 0; j <
Dim; j++)
12744 int nv =
elements[i]->GetNVertices();
12745 const int *ind =
elements[i]->GetVertices();
12746 os << partitioning[i]+1;
12747 for (
int j = 0; j < nv; j++)
12749 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
12750 vown[ind[j]][vcount[ind[j]]] = i;
12757 vcount[i] = voff[i+1] - voff[i];
12767 int k = partitioning[
faces_info[i].Elem1No];
12768 l = partitioning[l];
12769 if (interior_faces || k != l)
12786 int k = partitioning[
faces_info[i].Elem1No];
12787 l = partitioning[l];
12788 if (interior_faces || k != l)
12790 int nv =
faces[i]->GetNVertices();
12791 const int *ind =
faces[i]->GetVertices();
12793 for (
int j = 0; j < nv; j++)
12794 for (
int s = 0; s < vcount[ind[j]]; s++)
12795 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
12797 os <<
' ' << voff[ind[j]]+s+1;
12801 for (
int j = nv-1; j >= 0; j--)
12802 for (
int s = 0; s < vcount[ind[j]]; s++)
12803 if (vown[ind[j]][s] ==
faces_info[i].Elem2No)
12805 os <<
' ' << voff[ind[j]]+s+1;
12812 int k = partitioning[
faces_info[i].Elem1No];
12813 int nv =
faces[i]->GetNVertices();
12814 const int *ind =
faces[i]->GetVertices();
12816 for (
int j = 0; j < nv; j++)
12817 for (
int s = 0; s < vcount[ind[j]]; s++)
12818 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
12820 os <<
' ' << voff[ind[j]]+s+1;
12836 int k = partitioning[
faces_info[i].Elem1No];
12837 l = partitioning[l];
12838 if (interior_faces || k != l)
12851 <<
" 0 0 0 0 0 0 0\n"
12852 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
12853 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
12854 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
12855 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
12858 for (
int k = 0; k < vcount[i]; k++)
12859 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
12864 int nv =
elements[i]->GetNVertices();
12865 const int *ind =
elements[i]->GetVertices();
12866 os << i+1 <<
' ' << partitioning[i]+1;
12867 for (
int j = 0; j < nv; j++)
12869 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
12870 vown[ind[j]][vcount[ind[j]]] = i;
12877 vcount[i] = voff[i+1] - voff[i];
12886 int k = partitioning[
faces_info[i].Elem1No];
12887 l = partitioning[l];
12888 if (interior_faces || k != l)
12890 int nv =
faces[i]->GetNVertices();
12891 const int *ind =
faces[i]->GetVertices();
12893 for (
int j = 0; j < nv; j++)
12894 for (
int s = 0; s < vcount[ind[j]]; s++)
12895 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
12897 os <<
' ' << voff[ind[j]]+s+1;
12899 os <<
" 1.0 1.0 1.0 1.0\n";
12901 for (
int j = nv-1; j >= 0; j--)
12902 for (
int s = 0; s < vcount[ind[j]]; s++)
12903 if (vown[ind[j]][s] ==
faces_info[i].Elem2No)
12905 os <<
' ' << voff[ind[j]]+s+1;
12907 os <<
" 1.0 1.0 1.0 1.0\n";
12912 int k = partitioning[
faces_info[i].Elem1No];
12913 int nv =
faces[i]->GetNVertices();
12914 const int *ind =
faces[i]->GetVertices();
12916 for (
int j = 0; j < nv; j++)
12917 for (
int s = 0; s < vcount[ind[j]]; s++)
12918 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
12920 os <<
' ' << voff[ind[j]]+s+1;
12922 os <<
" 1.0 1.0 1.0 1.0\n";
12946 " NURBS mesh is not supported!");
12950 os <<
"MFEM mesh v1.0\n";
12954 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12959 "# TETRAHEDRON = 4\n"
12964 os <<
"\ndimension\n" <<
Dim
12972 const int *
const i_AF_f = Aface_face.
GetI();
12973 const int *
const j_AF_f = Aface_face.
GetJ();
12975 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
12976 for (
const int * iface = j_AF_f + i_AF_f[iAF];
12977 iface < j_AF_f + i_AF_f[iAF+1];
12980 os << iAF+1 <<
' ';
13013 int *nbea =
new int[na];
13020 for (i = 0; i < na; i++)
13032 for (k = 0; k < vert.
Size(); k++)
13044 for (k = 0; k < vert.
Size(); k++)
13045 if (vn[vert[k]] == 1)
13050 cg[bea*
spaceDim+j] += pointmat(j,k);
13061 for (k = 0; k < vert.
Size(); k++)
13066 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
13083 int *nbea =
new int[na];
13090 for (i = 0; i < na; i++)
13102 for (k = 0; k < vert.
Size(); k++)
13114 for (k = 0; k < vert.
Size(); k++)
13115 if (vn[vert[k]] == 1)
13120 cg[bea*
spaceDim+j] += pointmat(j,k);
13131 for (k = 0; k < vert.
Size(); k++)
13136 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
13152 for (
int i = 0; i <
vertices.Size(); i++)
13154 for (
int j = 0; j <
spaceDim; j++)
13175 "incompatible vector dimensions");
13183 for (
int d = 0; d <
spaceDim; d++)
13203 for (
int i = 0; i <
GetNE(); i++)
13208 for (
int j = 0; j < nv; j++)
13213 for (
int i = 0; i <
GetNBE(); i++)
13218 for (
int j = 0; j < nv; j++)
13224 for (
int i = 0; i < v2v.
Size(); i++)
13229 v2v[i] = num_vert++;
13233 if (num_vert == v2v.
Size()) {
return; }
13235 Vector nodes_by_element;
13240 for (
int i = 0; i <
GetNE(); i++)
13247 for (
int i = 0; i <
GetNE(); i++)
13256 for (
int i = 0; i <
GetNE(); i++)
13261 for (
int j = 0; j < nv; j++)
13266 for (
int i = 0; i <
GetNBE(); i++)
13271 for (
int j = 0; j < nv; j++)
13295 for (
int i = 0; i <
GetNE(); i++)
13308 int num_bdr_elem = 0;
13309 int new_bel_to_edge_nnz = 0;
13310 for (
int i = 0; i <
GetNBE(); i++)
13326 if (num_bdr_elem ==
GetNBE()) {
return; }
13330 Table *new_bel_to_edge = NULL;
13332 new_be_to_face.
Reserve(num_bdr_elem);
13335 new_bel_to_edge =
new Table;
13336 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
13338 for (
int i = 0; i <
GetNBE(); i++)
13343 int row = new_be_to_face.
Size();
13349 int *new_e = new_bel_to_edge->
GetRow(row);
13350 for (
int j = 0; j < ne; j++)
13354 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
13371 for (
int i = 0; i < attribs.
Size(); i++)
13383#ifdef MFEM_USE_MEMALLOC
13410 const int npts = point_mat.
Width();
13411 if (!npts) {
return 0; }
13412 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
13416 if (!
GetNE()) {
return 0; }
13425 min_dist = std::numeric_limits<real_t>::max();
13429 for (
int i = 0; i <
GetNE(); i++)
13433 for (
int k = 0; k < npts; k++)
13436 if (dist < min_dist(k))
13438 min_dist(k) = dist;
13447 for (
int k = 0; k < npts; k++)
13451 int res = inv_tr->
Transform(pt, ips[k]);
13454 elem_ids[k] = e_idx[k];
13458 if (pts_found != npts)
13462 for (
int k = 0; k < npts; k++)
13464 if (elem_ids[k] != -1) {
continue; }
13468 for (
int v = 0; v < elvertices.
Size(); v++)
13470 int vv = elvertices[v];
13472 const int* els = vtoel->
GetRow(vv);
13473 for (
int e = 0; e < ne; e++)
13475 if (els[e] == e_idx[k]) {
continue; }
13477 int res = inv_tr->
Transform(pt, ips[k]);
13480 elem_ids[k] = els[e];
13492 for (
int e = 0; e < neigh.
Size(); e++)
13498 int res = inv_tr->
Transform(pt, ips[k]);
13511 if (inv_trans == NULL) {
delete inv_tr; }
13513 if (warn && pts_found != npts)
13515 MFEM_WARNING((npts-pts_found) <<
" points were not found");
13530 MFEM_VERIFY(
Dim == 2 ||
Dim == 3,
"Only 2D/3D meshes supported right now.");
13531 MFEM_VERIFY(
Dim ==
spaceDim,
"Surface meshes not currently supported.");
13548 skew(0) = std::atan2(J.
Det(), col1 * col2);
13551 ori(0) = std::atan2(J(1,0), J(0,0));
13558 Vector col1, col2, col3;
13569 col1unit *= 1.0/len1;
13570 col2unit *= 1.0/len2;
13571 col3unit *= 1.0/len3;
13577 aspr(0) = len1/std::sqrt(len2*len3),
13578 aspr(1) = len2/std::sqrt(len1*len3);
13581 aspr(2) = std::sqrt(len1/(len2*len3)),
13582 aspr(3) = std::sqrt(len2/(len1*len3));
13585 Vector crosscol12, crosscol13;
13586 col1.
cross3D(col2, crosscol12);
13587 col1.
cross3D(col3, crosscol13);
13588 skew(0) = std::acos(col1unit*col2unit);
13589 skew(1) = std::acos(col1unit*col3unit);
13590 skew(2) = std::atan(len1*volume/(crosscol12*crosscol13));
13596 for (
int d=0; d<
Dim; d++) { rot(d, 0) = col1unit(d); }
13600 rot1 *= col1unit*col2unit;
13602 col1unit.
cross3D(col2unit, rot1);
13604 for (
int d=0; d <
Dim; d++) { rot(d, 1) = rot2(d); }
13607 for (
int d=0; d <
Dim; d++) { rot(d, 2) = rot1(d); }
13608 real_t delta = sqrt(pow(rot(2,1)-rot(1,2), 2.0) +
13609 pow(rot(0,2)-rot(2,0), 2.0) +
13610 pow(rot(1,0)-rot(0,1), 2.0));
13615 for (
int d = 0; d <
Dim; d++) { Iden(d, d) = 1.0; };
13621 MFEM_ABORT(
"Invalid rotation matrix. Contact TMOP Developers.");
13626 ori(0) = (1./
delta)*(rot(2,1)-rot(1,2));
13627 ori(1) = (1./
delta)*(rot(0,2)-rot(2,0));
13628 ori(2) = (1./
delta)*(rot(1,0)-rot(0,1));
13629 ori(3) = std::acos(0.5*(rot.
Trace()-1.0));
13638 entity_to_vertex(entity_to_vertex_)
13640 int geom_offset = 0;
13654 while (geom_offsets[geom+1] <= bytype_entity_id) { geom++; }
13658 const int geom_elem_id = bytype_entity_id - geom_offsets[geom];
13660 return { geom, nv, v };
13665 os <<
"MFEM mesh v1.2\n";
13669 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
13674 "# TETRAHEDRON = 4\n"
13681 os <<
"\ndimension\n" <<
dim;
13687 "invalid MeshPart state");
13690 "invalid MeshPart state");
13691 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
13693 const int bytype_elem_id = have_element_map ?
13698 for (
int i = 0; i < ent.
num_verts; i++)
13700 os <<
' ' << ent.
verts[i];
13710 "invalid MeshPart state");
13713 "invalid MeshPart state");
13716 const int bytype_bdr_id = have_boundary_map ?
13721 for (
int i = 0; i < ent.
num_verts; i++)
13723 os <<
' ' << ent.
verts[i];
13733 os << sdim <<
'\n';
13737 for (
int d = 1; d < sdim; d++)
13750 os <<
"\nmfem_serial_mesh_end\n";
13754 os <<
"\ncommunication_groups\n";
13755 os <<
"number_of_groups " << num_groups <<
"\n\n";
13757 os <<
"# number of entities in each group, followed by ranks in group\n";
13758 for (
int group_id = 0; group_id < num_groups; ++group_id)
13763 for (
int group_member_index = 0; group_member_index < group_size;
13764 ++group_member_index)
13766 os <<
' ' << group_ptr[group_member_index];
13777 MFEM_VERIFY(g2v.
RowSize(0) == 0,
"internal erroor");
13781 MFEM_VERIFY(g2ev.
RowSize(0) == 0,
"internal erroor");
13786 MFEM_VERIFY(g2tv.
RowSize(0) == 0,
"internal erroor");
13787 MFEM_VERIFY(g2qv.
RowSize(0) == 0,
"internal erroor");
13788 const int total_shared_faces =
13790 os <<
"total_shared_faces " << total_shared_faces <<
'\n';
13792 os <<
"\n# group 0 has no shared entities\n";
13793 for (
int gr = 1; gr < num_groups; gr++)
13796 const int nv = g2v.
RowSize(gr);
13797 const int *sv = g2v.
GetRow(gr);
13798 os <<
"\n# group " << gr <<
"\nshared_vertices " << nv <<
'\n';
13799 for (
int i = 0; i < nv; i++)
13801 os << sv[i] <<
'\n';
13806 const int ne = g2ev.
RowSize(gr)/2;
13807 const int *se = g2ev.
GetRow(gr);
13808 os <<
"\nshared_edges " << ne <<
'\n';
13809 for (
int i = 0; i < ne; i++)
13811 const int *v = se + 2*i;
13812 os << v[0] <<
' ' << v[1] <<
'\n';
13817 const int nt = g2tv.
RowSize(gr)/3;
13818 const int *st = g2tv.
GetRow(gr);
13819 const int nq = g2qv.
RowSize(gr)/4;
13821 os <<
"\nshared_faces " << nt+nq <<
'\n';
13822 for (
int i = 0; i < nt; i++)
13825 const int *v = st + 3*i;
13826 for (
int j = 0; j < 3; j++) { os <<
' ' << v[j]; }
13829 for (
int i = 0; i < nq; i++)
13832 const int *v =
sq + 4*i;
13833 for (
int j = 0; j < 4; j++) { os <<
' ' << v[j]; }
13840 os <<
"\nmfem_mesh_end" << endl;
13857 "invalid MeshPart state");
13860 "invalid MeshPart state");
13862 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
13864 const int bytype_elem_id = have_element_map ?
13875 static_cast<Tetrahedron*
>(el)->SetRefinementFlag(ref_flag);
13877 mesh->AddElement(el);
13885 "invalid MeshPart state");
13888 "invalid MeshPart state");
13891 const int bytype_bdr_id = have_boundary_map ?
13897 mesh->AddBdrElement(bdr);
13904 MFEM_ASSERT(!
nodes,
"invalid MeshPart state");
13905 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
13913 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
13915 mesh->AddVertex(0., 0., 0.);
13920 mesh->FinalizeTopology(
false);
13928 const int *partitioning_,
13960 for (
int i = 0; i < boundary_to_part.
Size(); i++)
13962 int face, o, el1, el2;
13965 boundary_to_part[i] =
13971 for (
int i = 0; i < boundary_to_part.
Size(); i++)
13980 for (
int i = 0; i < boundary_to_part.
Size(); i++)
13995 delete vert_element;
14002 MFEM_VERIFY(0 <= part_id && part_id < num_parts,
14003 "invalid part_id = " << part_id
14004 <<
", num_parts = " << num_parts);
14037 mesh_part.
nodes.reset(
nullptr);
14039 mesh_part.
mesh.reset(
nullptr);
14047 int geom_marker = 0, num_geom = 0;
14048 for (
int i = 0; i < num_elems; i++)
14054 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
14056 "overflow in 'entity_to_vertex[geom]', geom: "
14081 if ((geom_marker & (1 << geom)) == 0)
14083 geom_marker |= (1 << geom);
14098 offsets[g] = offset;
14102 for (
int i = 0; i < num_elems; i++)
14113 geom_marker = 0; num_geom = 0;
14114 for (
int i = 0; i < num_bdr_elems; i++)
14120 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
14122 "overflow in 'entity_to_vertex[geom]', geom: "
14126 if ((geom_marker & (1 << geom)) == 0)
14128 geom_marker |= (1 << geom);
14139 offsets[g] = offset;
14143 for (
int i = 0; i < num_bdr_elems; i++)
14154 std::unordered_set<int> vertex_set;
14155 for (
int i = 0; i < num_elems; i++)
14161 vertex_set.insert(v, v + nv);
14163 vertex_loc_to_glob.
SetSize(
static_cast<int>(vertex_set.size()));
14164 std::copy(vertex_set.begin(), vertex_set.end(),
14165 vertex_loc_to_glob.
begin());
14167 vertex_loc_to_glob.
Sort();
14177 for (
int i = 0; i < vert_array.
Size(); i++)
14179 const int glob_id = vert_array[i];
14180 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
14181 MFEM_ASSERT(loc_id >= 0,
"internal error: global vertex id not found");
14182 vert_array[i] = loc_id;
14189 MFEM_VERIFY(numeric_limits<int>::max()/sdim >= vertex_loc_to_glob.
Size(),
14190 "overflow in 'vertex_coordinates', num_vertices = "
14191 << vertex_loc_to_glob.
Size() <<
", sdim = " << sdim);
14193 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
14196 for (
int d = 0; d < sdim; d++)
14213 mesh_part.
mesh->NewNodes(*mesh_part.
nodes,
false);
14235 std::unordered_set<int> face_set;
14238 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14240 const int glob_elem_id = elem_list[loc_elem_id];
14241 const int nfaces = elem_to_face.
RowSize(glob_elem_id);
14242 const int *faces = elem_to_face.
GetRow(glob_elem_id);
14243 face_set.insert(faces, faces + nfaces);
14247 for (
int glob_face_id : face_set)
14251 if (el[1] < 0) {
continue; }
14254 MFEM_ASSERT(el[0] == part_id || el[1] == part_id,
"internal error");
14255 if (el[0] != part_id || el[1] != part_id)
14258 const int group_id = groups.
Insert(group);
14262 shared_faces.
Sort();
14273 std::unordered_set<int> edge_set;
14276 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14278 const int glob_elem_id = elem_list[loc_elem_id];
14279 const int nedges = elem_to_edge.
RowSize(glob_elem_id);
14280 const int *edges = elem_to_edge.
GetRow(glob_elem_id);
14281 edge_set.insert(edges, edges + nedges);
14285 for (
int glob_edge_id : edge_set)
14291 for (
int j = 0; j < nelem; j++)
14297 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
14298 if (group.
Size() > 1)
14300 const int group_id = groups.
Insert(group);
14304 shared_edges.
Sort();
14315 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
14318 const int glob_vertex_id = vertex_loc_to_glob[i];
14323 for (
int j = 0; j < nelem; j++)
14329 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
14330 if (group.
Size() > 1)
14332 const int group_id = groups.
Insert(group);
14339 const int num_groups = groups.
Size();
14345 Table &group__shared_vertex_to_vertex =
14347 group__shared_vertex_to_vertex.
MakeI(num_groups);
14348 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
14350 const int group_id = shared_verts[sv].two;
14353 group__shared_vertex_to_vertex.
MakeJ();
14354 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
14356 const int glob_vertex_id = shared_verts[sv].one;
14357 const int group_id = shared_verts[sv].two;
14358 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(glob_vertex_id);
14359 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
14360 group__shared_vertex_to_vertex.
AddConnection(group_id, loc_vertex_id);
14362 group__shared_vertex_to_vertex.
ShiftUpI();
14367 Table &group__shared_edge_to_vertex =
14369 group__shared_edge_to_vertex.
MakeI(num_groups);
14370 for (
int se = 0; se < shared_edges.
Size(); se++)
14372 const int group_id = shared_edges[se].two;
14375 group__shared_edge_to_vertex.
MakeJ();
14377 for (
int se = 0; se < shared_edges.
Size(); se++)
14379 const int glob_edge_id = shared_edges[se].one;
14380 const int group_id = shared_edges[se].two;
14381 const int *v = edge_to_vertex.
GetRow(glob_edge_id);
14382 for (
int i = 0; i < 2; i++)
14384 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(v[i]);
14385 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
14386 group__shared_edge_to_vertex.
AddConnection(group_id, loc_vertex_id);
14389 group__shared_edge_to_vertex.
ShiftUpI();
14396 Table &group__shared_tria_to_vertex =
14398 Table &group__shared_quad_to_vertex =
14401 group__shared_tria_to_vertex.
MakeI(num_groups);
14402 group__shared_quad_to_vertex.
MakeI(num_groups);
14403 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
14405 const int glob_face_id = shared_faces[sf].one;
14406 const int group_id = shared_faces[sf].two;
14411 group__shared_tria_to_vertex.
MakeJ();
14412 group__shared_quad_to_vertex.
MakeJ();
14413 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
14415 const int glob_face_id = shared_faces[sf].one;
14416 const int group_id = shared_faces[sf].two;
14452 for (
int i = 0; i < vertex_ids.
Size(); i++)
14454 const int glob_id = vertex_ids[i];
14455 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
14456 MFEM_ASSERT(loc_id >= 0,
"internal error");
14457 vertex_ids[i] = loc_id;
14460 AddConnections(group_id, vertex_ids, vertex_ids.
Size());
14462 group__shared_tria_to_vertex.
ShiftUpI();
14463 group__shared_quad_to_vertex.
ShiftUpI();
14467std::unique_ptr<FiniteElementSpace>
14475 return std::unique_ptr<FiniteElementSpace>(
14477 global_fespace.
FEColl(),
14482std::unique_ptr<GridFunction>
14487 std::unique_ptr<GridFunction> local_gf(
new GridFunction(&local_fespace));
14495 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14497 const int glob_elem_id = elem_list[loc_elem_id];
14500 if (glob_dt) { glob_dt->InvTransformPrimal(loc_vals); }
14503 local_gf->SetSubVector(lvdofs, loc_vals);
14517 "mixed meshes are not supported!");
14518 MFEM_ASSERT(
mesh->
GetNodes(),
"meshes without nodes are not supported!");
14531 Compute(
nodes, d_mt);
14541 const int vdim = fespace->
GetVDim();
14542 const int NE = fespace->
GetNE();
14543 const int ND = fe->
GetDof();
14546 unsigned eval_flags = 0;
14548 Device::GetDeviceMemoryType();
14571 qi->DisableTensorProducts(!use_tensor_products);
14579 Vector Enodes(vdim*ND*NE, my_d_mt);
14580 elem_restr->Mult(
nodes, Enodes);
14581 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
14601 const int vdim = fespace->
GetVDim();
14617 unsigned eval_flags = 0;
14664 V(1) = s * ((ip.
y + layer) / n);
14669 V(2) = s * ((ip.
z + layer) / n);
14678 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
14682 int nvy = (closed) ? (ny) : (ny + 1);
14683 int nvt = mesh->
GetNV() * nvy;
14692 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
14697 for (
int i = 0; i < mesh->
GetNV(); i++)
14700 for (
int j = 0; j < nvy; j++)
14702 vc[1] = sy * (
real_t(j) / ny);
14708 for (
int i = 0; i < mesh->
GetNE(); i++)
14713 for (
int j = 0; j < ny; j++)
14716 qv[0] = vert[0] * nvy + j;
14717 qv[1] = vert[1] * nvy + j;
14718 qv[2] = vert[1] * nvy + (j + 1) % nvy;
14719 qv[3] = vert[0] * nvy + (j + 1) % nvy;
14725 for (
int i = 0; i < mesh->
GetNBE(); i++)
14730 for (
int j = 0; j < ny; j++)
14733 sv[0] = vert[0] * nvy + j;
14734 sv[1] = vert[0] * nvy + (j + 1) % nvy;
14750 for (
int i = 0; i < mesh->
GetNE(); i++)
14756 sv[0] = vert[0] * nvy;
14757 sv[1] = vert[1] * nvy;
14761 sv[0] = vert[1] * nvy + ny;
14762 sv[1] = vert[0] * nvy + ny;
14778 string cname = name;
14779 if (cname ==
"Linear")
14783 else if (cname ==
"Quadratic")
14787 else if (cname ==
"Cubic")
14791 else if (!strncmp(name,
"H1_", 3))
14795 else if (!strncmp(name,
"L2_T", 4))
14799 else if (!strncmp(name,
"L2_", 3))
14806 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
14818 for (
int i = 0; i < mesh->
GetNE(); i++)
14821 for (
int j = ny-1; j >= 0; j--)
14838 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
14843 int nvt = mesh->
GetNV() * nvz;
14848 bool wdgMesh =
false;
14849 bool hexMesh =
false;
14853 for (
int i = 0; i < mesh->
GetNV(); i++)
14857 for (
int j = 0; j < nvz; j++)
14859 vc[2] = sz * (
real_t(j) / nz);
14865 for (
int i = 0; i < mesh->
GetNE(); i++)
14875 for (
int j = 0; j < nz; j++)
14878 pv[0] = vert[0] * nvz + j;
14879 pv[1] = vert[1] * nvz + j;
14880 pv[2] = vert[2] * nvz + j;
14881 pv[3] = vert[0] * nvz + (j + 1) % nvz;
14882 pv[4] = vert[1] * nvz + (j + 1) % nvz;
14883 pv[5] = vert[2] * nvz + (j + 1) % nvz;
14890 for (
int j = 0; j < nz; j++)
14893 hv[0] = vert[0] * nvz + j;
14894 hv[1] = vert[1] * nvz + j;
14895 hv[2] = vert[2] * nvz + j;
14896 hv[3] = vert[3] * nvz + j;
14897 hv[4] = vert[0] * nvz + (j + 1) % nvz;
14898 hv[5] = vert[1] * nvz + (j + 1) % nvz;
14899 hv[6] = vert[2] * nvz + (j + 1) % nvz;
14900 hv[7] = vert[3] * nvz + (j + 1) % nvz;
14902 mesh3d->
AddHex(hv, attr);
14906 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
14907 << geom <<
"\'" << endl;
14913 for (
int i = 0; i < mesh->
GetNBE(); i++)
14918 for (
int j = 0; j < nz; j++)
14921 qv[0] = vert[0] * nvz + j;
14922 qv[1] = vert[1] * nvz + j;
14923 qv[2] = vert[1] * nvz + (j + 1) % nvz;
14924 qv[3] = vert[0] * nvz + (j + 1) % nvz;
14933 for (
int i = 0; i < mesh->
GetNE(); i++)
14944 tv[0] = vert[0] * nvz;
14945 tv[1] = vert[2] * nvz;
14946 tv[2] = vert[1] * nvz;
14950 tv[0] = vert[0] * nvz + nz;
14951 tv[1] = vert[1] * nvz + nz;
14952 tv[2] = vert[2] * nvz + nz;
14960 qv[0] = vert[0] * nvz;
14961 qv[1] = vert[3] * nvz;
14962 qv[2] = vert[2] * nvz;
14963 qv[3] = vert[1] * nvz;
14967 qv[0] = vert[0] * nvz + nz;
14968 qv[1] = vert[1] * nvz + nz;
14969 qv[2] = vert[2] * nvz + nz;
14970 qv[3] = vert[3] * nvz + nz;
14976 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
14977 << geom <<
"\'" << endl;
14983 if ( hexMesh && wdgMesh )
14987 else if ( hexMesh )
14991 else if ( wdgMesh )
15004 string cname = name;
15005 if (cname ==
"Linear")
15009 else if (cname ==
"Quadratic")
15013 else if (cname ==
"Cubic")
15017 else if (!strncmp(name,
"H1_", 3))
15021 else if (!strncmp(name,
"L2_T", 4))
15025 else if (!strncmp(name,
"L2_", 3))
15032 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
15044 for (
int i = 0; i < mesh->
GetNE(); i++)
15047 for (
int j = nz-1; j >= 0; j--)
15068 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
15069 <<
" 0 0 " << i <<
" -1 0\n";
15076 real_t mid[3] = {0, 0, 0};
15077 for (
int j = 0; j < 2; j++)
15079 for (
int k = 0; k <
spaceDim; k++)
15085 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
15086 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
Node::Index insert_node(Float length=1)
Arc::Index insert_arc(Node::Index i, Node::Index j, Float w=1, Float b=1)
uint rank(Node::Index i) const
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int FindSorted(const T &el) const
Do bisection search for 'el' in a sorted array; return -1 if not found.
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
void DeleteAll()
Delete the whole array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
T * GetData()
Returns the data.
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
T * end()
STL-like end. Returns pointer after the last element of the array.
T * begin()
STL-like begin. Returns pointer to the first element of the array.
T & Last()
Return the last element in the array.
bool SetsExist() const
Return true if any named sets are currently defined.
bool AttributeSetExists(const std::string &name) const
Return true is the named attribute set is present.
Array< int > GetAttributeSetMarker(const std::string &set_name) const
Return a marker array corresponding to a named attribute set.
void Print(std::ostream &out=mfem::out, int width=-1) const
Print the contents of the container to an output stream.
void Copy(AttributeSets ©) const
Create a copy of the internal data to the provided copy.
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Piecewise-(bi)cubic continuous finite elements.
int NumberOfEntries() const
Data type dense matrix using column-major storage.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
const real_t * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
real_t * Data() const
Returns the matrix data array.
real_t * GetData() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
real_t Trace() const
Trace of a square matrix.
real_t CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
void Print(std::ostream &out=mfem::out, int width_=4) const override
Prints matrix to stream out.
void GetColumn(int c, Vector &col) const
real_t FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Rank 3 tensor (array of matrices)
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
void UseExternalData(real_t *ext_data, int i, int j, int k)
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Abstract data type element.
Geometry::Type GetGeometryType() const
virtual Element * Duplicate(Mesh *m) const =0
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
virtual Type GetType() const =0
Returns element's type.
Type
Constants for the classes derived from Element.
int GetAttribute() const
Return element's attribute.
static Type TypeFromGeometry(const Geometry::Type geom)
Return the Element::Type associated with the given Geometry::Type.
virtual int GetNVertices() const =0
virtual void SetVertices(const Array< int > &v)=0
Set the indices defining the vertices.
Structure for storing face geometric factors: coordinates, Jacobians, determinants of the Jacobians,...
Vector normal
Normals at all quadrature points.
Vector J
Jacobians of the element transformations at all quadrature points.
const IntegrationRule * IntRule
Vector X
Mapped (physical) coordinates of all quadrature points.
FaceGeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, FaceType type, MemoryType d_mt=MemoryType::DEFAULT)
Vector detJ
Determinants of the Jacobians at all quadrature points.
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
@ DERIVATIVES
Evaluate the derivatives at quadrature points.
@ DETERMINANTS
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
@ VALUES
Evaluate the values at quadrature points.
void DisableTensorProducts(bool disable=true) const
Disable the use of tensor product evaluations, for tensor-product elements, e.g. quads and hexes.
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
void Mult(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
Interpolate the E-vector e_vec to quadrature points.
Base class for operators that extracts Face degrees of freedom.
void Mult(const Vector &x, Vector &y) const override=0
Extract the face degrees of freedom from x into y.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i].
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name.
virtual const char * Name() const
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
void GetVertexDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the specified vertices.
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
virtual void UpdateMeshPointer(Mesh *new_mesh)
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Ordering::Type GetOrdering() const
Return the ordering method.
int GetNE() const
Returns number of elements in the mesh.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
void SetRelaxedHpConformity(bool relaxed=true)
int GetNFDofs() const
Number of all scalar face-interior dofs.
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
const FiniteElementCollection * FEColl() const
Mesh * GetMesh() const
Returns the mesh.
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
int GetVDim() const
Returns the vector dimension of the finite element space.
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified face.
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
virtual int GetMaxElementOrder() const
Return the maximum polynomial order over all elements.
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified element.
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering f_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Abstract class for all finite elements.
int GetDim() const
Returns the reference space dimension for the finite element.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Vector X
Mapped (physical) coordinates of all quadrature points.
const IntegrationRule * IntRule
Vector detJ
Determinants of the Jacobians at all quadrature points.
Vector J
Jacobians of the element transformations at all quadrature points.
GeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, MemoryType d_mt=MemoryType::DEFAULT)
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
static const int Dimension[NumGeom]
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
static const char * Name[NumGeom]
static const int NumVerts[NumGeom]
const IntegrationRule * GetVertices(int GeomType) const
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType.
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
static int GetInverseOrientation(Type geom_type, int orientation)
Return the inverse of the given orientation for the specified geometry type.
static const int DimStart[MaxDim+2]
Class for grid function - Vector with associated FE space.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec_owned and fes.
FiniteElementSpace * FESpace()
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int VectorDim() const
Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying fes.
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Arbitrary order H1-conforming (continuous) finite elements.
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
int GetId(int p1, int p2)
Get id of item whose parents are p1, p2... Create it if it doesn't exist.
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn't exist.
Data type hexahedron element.
void Recreate(const int n, const int *p)
Create an integer set from C-array 'p' of 'n' integers. Overwrites any existing set data.
Class for integration point with weight.
void Get(real_t *p, const int dim) const
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Arbitrary order "L2-conforming" discontinuous finite elements.
Piecewise-(bi/tri)linear continuous finite elements.
int Insert(const IntegerSet &s)
Check to see if set 's' is in the list. If not append it to the end of the list. Returns the index of...
void AsTable(Table &t) const
Write the list of sets into table 't'.
int Size() const
Return the number of integer sets in the list.
Class containing a minimal description of a part (a subset of the elements) of a Mesh and its connect...
Array< real_t > vertex_coordinates
int dimension
Reference space dimension of the elements.
int num_vertices
Number of vertices.
Table group_shared_entity_to_vertex[Geometry::NumGeom]
Array< int > entity_to_vertex[Geometry::NumGeom]
std::unique_ptr< Mesh > mesh
Array< int > boundary_map
Optional re-ordering for the boundary elements, similar to 'element_map'.
std::unique_ptr< FiniteElementSpace > nodal_fes
int num_parts
Total number of MeshParts.
Array< int > tet_refine_flags
Store the refinement flags for tetraheral elements. If all tets have zero refinement flags then this ...
int space_dimension
Dimension of the physical space into which the MeshPart is embedded.
int num_bdr_elements
Number of boundary elements with reference space dimension equal to 'dimension'-1.
int num_elements
Number of elements with reference space dimension equal to 'dimension'.
Mesh & GetMesh()
Construct a serial Mesh object from the MeshPart.
int my_part_id
Index of the part described by this MeshPart: 0 <= 'my_part_id' < 'num_parts'.
std::unique_ptr< GridFunction > nodes
void Print(std::ostream &os) const
Write the MeshPart to a stream using the parallel format "MFEM mesh v1.2".
Array< int > bdr_attributes
Array< int > partitioning
std::unique_ptr< FiniteElementSpace > ExtractFESpace(MeshPart &mesh_part, const FiniteElementSpace &global_fespace) const
Construct a local version of the given FiniteElementSpace global_fespace corresponding to the given m...
MeshPartitioner(Mesh &mesh_, int num_parts_, const int *partitioning_=nullptr, int part_method=1)
Construct a MeshPartitioner.
std::unique_ptr< GridFunction > ExtractGridFunction(const MeshPart &mesh_part, const GridFunction &global_gf, FiniteElementSpace &local_fespace) const
Construct a local version of the given GridFunction, global_gf, corresponding to the given mesh_part....
void ExtractPart(int part_id, MeshPart &mesh_part) const
Construct a MeshPart corresponding to the given part_id.
List of mesh geometries stored as Array<Geometry::Type>.
int CheckElementOrientation(bool fix_it=true)
Check (and optionally attempt to fix) the orientation of the elements.
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
void GetEdgeOrdering(const DSTable &v_to_v, Array< int > &order)
void GetLocalFaceTransformation(int face_type, int elem_type, IsoparametricTransformation &Transf, int info) const
A helper method that constructs a transformation from the reference space of a face to the reference ...
void NURBSCoarsening(int cf=2, real_t tol=1.0e-12)
void SetVerticesFromNodes(const GridFunction *nodes)
Helper to set vertex coordinates given a high-order curvature function.
void ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
int GetPatchBdrAttribute(int i) const
Return the attribute of patch boundary element i, for a NURBS mesh.
int GetElementToEdgeTable(Table &)
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
void SetVertices(const Vector &vert_coord)
Element * NewElement(int geom)
Operation GetLastOperation() const
Return type of last modification of the mesh.
IsoparametricTransformation Transformation2
int GetNEdges() const
Return the number of edges.
void GetBdrElementFace(int i, int *f, int *o) const
void InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
Begin construction of a mesh.
Table * GetVertexToBdrElementTable()
static void PrintElement(const Element *el, std::ostream &os)
Array< FaceInfo > faces_info
int EulerNumber() const
Equals 1 + num_holes - num_loops.
CoarseFineTransformations CoarseFineTr
void GetElementJacobian(int i, DenseMatrix &J, const IntegrationPoint *ip=NULL)
int AddSegment(int v1, int v2, int attr=1)
Adds a segment to the mesh given by 2 vertices v1 and v2.
int AddBdrElement(Element *elem)
void GetElementColoring(Array< int > &colors, int el0=0)
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
void FinalizeMesh(int refine=0, bool fix_orientation=true)
Finalize the construction of any type of Mesh.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
static void PrintElementWithoutAttr(const Element *el, std::ostream &os)
MemAlloc< Tetrahedron, 1024 > TetMemory
void RedRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
void ReadTrueGridMesh(std::istream &input)
static const int vtk_quadratic_tet[10]
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
virtual void GetExteriorFaceMarker(Array< int > &face_marker) const
Populate a marker array identifying exterior faces.
IsoparametricTransformation EdgeTransformation
static FiniteElement * GetTransformationFEforElementType(Element::Type)
Return FiniteElement for reference element of the specified type.
int AddBdrQuad(int v1, int v2, int v3, int v4, int attr=1)
int * CartesianPartitioning(int nxyz[])
Array< int > FindFaceNeighbors(const int elem) const
Returns the sorted, unique indices of elements sharing a face with element elem, including elem.
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
Element::Type GetElementType(int i) const
Returns the type of element i.
void GetLocalSegToQuadTransformation(IsoparametricTransformation &loc, int i) const
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
void BdrBisection(int i, const HashTable< Hashed2 > &)
Bisect a boundary triangle: boundary element with index i is bisected.
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
const Table & ElementToEdgeTable() const
bool Conforming() const
Return a bool indicating whether this mesh is conforming.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
virtual void UnmarkNamedBoundaries(const std::string &set_name, Array< int > &bdr_marker) const
Unmark boundary attributes in the named set.
void ReadNetgen3DMesh(std::istream &input)
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Geometry::Type GetElementGeometry(int i) const
Geometry::Type GetBdrElementGeometry(int i) const
int AddTri(const int *vi, int attr=1)
Adds a triangle to the mesh given by 3 vertices vi.
static Mesh MakeCartesian1D(int n, real_t sx=1.0)
Creates 1D mesh, divided into n equal intervals.
int GetAttribute(int i) const
Return the attribute of element i.
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally,...
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
void UniformRefinement3D_base(Array< int > *f2qf=NULL, DSTable *v_to_v_p=NULL, bool update_nodes=true)
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
long nodes_sequence
Counter for geometric factor invalidation.
virtual void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
IsoparametricTransformation FaceTransformation
Array< NCFaceInfo > nc_faces_info
void MakeRefined_(Mesh &orig_mesh, const Array< int > &ref_factors, int ref_type)
Internal function used in Mesh::MakeRefined.
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
Adds a wedge to the mesh given by 6 vertices v1 through v6.
Array< int > GetFaceToBdrElMap() const
static Mesh MakeCartesian2DWith4TrisPerQuad(int nx, int ny, real_t sx=1.0, real_t sy=1.0)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*4 triangles.
void ReadInlineMesh(std::istream &input, bool generate_edges=false)
void SetPatchAttribute(int i, int attr)
Set the attribute of patch i, for a NURBS mesh.
void FinalizeTetMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a tetrahedral Mesh.
real_t GetLength(int i, int j) const
Return the length of the segment from node i to node j.
const FiniteElementSpace * GetNodalFESpace() const
void AddBdrQuadAsTriangles(const int *vi, int attr=1)
int AddPyramid(int v1, int v2, int v3, int v4, int v5, int attr=1)
Adds a pyramid to the mesh given by 5 vertices v1 through v5.
void Loader(std::istream &input, int generate_edges=0, std::string parse_tag="")
const Table & ElementToElementTable()
void ScaleElements(real_t sf)
void GenerateNCFaceInfo()
void ReadLineMesh(std::istream &input)
void ApplyLocalSlaveTransformation(FaceElementTransformations &FT, const FaceInfo &fi, bool is_ghost) const
real_t AggregateError(const Array< real_t > &elem_error, const int *fine, int nfine, int op)
Derefinement helper.
void CheckPartitioning(int *partitioning_)
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
Geometry::Type GetTypicalElementGeometry() const
If the local mesh is not empty, return GetElementGeometry(0); otherwise, return a typical Geometry pr...
void GetLocalPtToSegTransformation(IsoparametricTransformation &, int i) const
Used in GetFaceElementTransformations (...)
bool Nonconforming() const
Return a bool indicating whether this mesh is nonconforming.
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
virtual void MarkExternalBoundaries(Array< int > &bdr_marker, bool excl=true) const
Mark boundary attributes of external boundaries.
void PrintCharacteristics(Vector *Vh=NULL, Vector *Vk=NULL, std::ostream &os=mfem::out)
Compute and print mesh characteristics such as number of vertices, number of elements,...
static int ComposeQuadOrientations(int ori_a_b, int ori_b_c)
int AddTriangle(int v1, int v2, int v3, int attr=1)
Adds a triangle to the mesh given by 3 vertices v1 through v3.
static const int vtk_quadratic_wedge[18]
int EulerNumber2D() const
Equals 1 - num_holes.
AttributeSets bdr_attribute_sets
Named sets of boundary element attributes.
void AddBdrElements(Array< Element * > &bdr_elems, const Array< int > &be_to_face)
Add an array of boundary elements to the mesh, along with map from the elements to their faces.
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
void GetVertices(Vector &vert_coord) const
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
virtual int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type, does not count master nonconforming face...
void Make1D(int n, real_t sx=1.0)
void RefineNURBSFromFile(std::string ref_file)
const Element * GetElement(int i) const
Return pointer to the i'th element object.
int AddBdrPoint(int v, int attr=1)
void FinalizeWedgeMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a wedge Mesh.
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
void PrintTopo(std::ostream &os, const Array< int > &e_to_k, const int version, const std::string &comment="") const
Write the beginning of a NURBS mesh to os, specifying the NURBS patch topology. Optional file comment...
static Mesh MakeSimplicial(const Mesh &orig_mesh)
void SetPatchBdrAttribute(int i, int attr)
Set the attribute of patch boundary element i, for a NURBS mesh.
int GetNFaces() const
Return the number of faces in a 3D mesh.
int AddVertexAtMeanCenter(const int *vi, const int nverts, int dim=3)
static int GetTetOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
real_t GetGeckoElementOrdering(Array< int > &ordering, int iterations=4, int window=4, int period=2, int seed=0, bool verbose=false, real_t time_limit=0)
static int EncodeFaceInfo(int local_face_index, int orientation)
Given local_face_index and orientation, return the corresponding encoded "face info int".
bool FaceIsTrueInterior(int FaceNo) const
const CoarseFineTransformations & GetRefinementTransforms() const
void Make2D5QuadsFromQuad(int nx, int ny, real_t sx, real_t sy)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*5 quadrilaterals.
void GetLocalQuadToWdgTransformation(IsoparametricTransformation &loc, int i) const
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.
void SetAttribute(int i, int attr)
Set the attribute of element i.
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
static Mesh MakeCartesian2DWith5QuadsPerQuad(int nx, int ny, real_t sx=1.0, real_t sy=1.0)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*5 quadrilaterals.
const FaceGeometricFactors * GetFaceGeometricFactors(const IntegrationRule &ir, const int flags, FaceType type, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors for the faces corresponding to the given integration rule.
void PrintWithPartitioning(int *partitioning, std::ostream &os, int elem_attr=0) const
Prints the mesh with boundary elements given by the boundary of the subdomains, so that the boundary ...
void Clear()
Clear the contents of the Mesh.
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
void ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo, const std::string &xml_prefix="")
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
This function is not public anymore. Use GeneralRefinement instead.
int GetNE() const
Returns number of elements.
void Make3D(int nx, int ny, int nz, Element::Type type, real_t sx, real_t sy, real_t sz, bool sfc_ordering)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
virtual void Save(const std::string &fname, int precision=16) const
int AddTet(int v1, int v2, int v3, int v4, int attr=1)
Adds a tetrahedron to the mesh given by 4 vertices v1 through v4.
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const
int Dimension() const
Dimension of the reference space used within the elements.
ElementTransformation * GetTypicalElementTransformation()
If the local mesh is not empty return GetElementTransformation(0); otherwise, return the identity tra...
void RandomRefinement(real_t prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
void CheckDisplacements(const Vector &displacements, real_t &tmax)
friend class NURBSExtension
void AddTriangleFaceElement(int lf, int gf, int el, int v0, int v1, int v2)
void AddHexAs24TetsWithPoints(int *vi, std::map< std::array< int, 4 >, int > &hex_face_verts, int attr=1)
Adds 24 tetrahedrons to the mesh by splitting a hexahedron.
void GetNode(int i, real_t *coord) const
void ReorderElements(const Array< int > &ordering, bool reorder_vertices=true)
void GreenRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
void UpdateNodes()
Update the nodes of a curved mesh after the topological part of a Mesh::Operation,...
void PrintElementsWithPartitioning(int *partitioning, std::ostream &os, int interior_faces=0)
Mesh & operator=(Mesh &&mesh)
Move assignment operator.
static int InvertQuadOrientation(int ori)
Array< FaceGeometricFactors * > face_geom_factors
void GetLocalTriToPyrTransformation(IsoparametricTransformation &loc, int i) const
static Mesh MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
Create a refined (by any factor) version of orig_mesh.
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Element * ReadElementWithoutAttr(std::istream &input)
int AddBdrSegment(int v1, int v2, int attr=1)
bool DerefineByError(Array< real_t > &elem_error, real_t threshold, int nc_limit=0, int op=1)
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
int AddElement(Element *elem)
static int DecodeFaceInfoLocalIndex(int info)
Given a "face info int", return the local face index.
FaceInformation GetFaceInformation(int f) const
int GetNumFacesWithGhost() const
Return the number of faces (3D), edges (2D) or vertices (1D) including ghost faces.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
void RefineAtVertex(const Vertex &vert, real_t eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
static int InvertTriOrientation(int ori)
STable3D * GetElementToFaceTable(int ret_ftbl=0)
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
void Make3D24TetsFromHex(int nx, int ny, int nz, real_t sx, real_t sy, real_t sz)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz*24 tetrahedrons.
virtual bool NonconformingDerefinement(Array< real_t > &elem_error, real_t threshold, int nc_limit=0, int op=1)
NC version of GeneralDerefinement.
void AddVertexParents(int i, int p1, int p2)
Mark vertex i as nonconforming, with parent vertices p1 and p2.
MFEM_DEPRECATED void GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
Deprecated.
int GetPatchAttribute(int i) const
Return the attribute of patch i, for a NURBS mesh.
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
void Printer(std::ostream &os=mfem::out, std::string section_delimiter="", const std::string &comments="") const
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
ElementTransformation * GetBdrElementTransformation(int i)
Returns a pointer to the transformation defining the i-th boundary element.
IsoparametricTransformation Transformation
void Make2D4TrisFromQuad(int nx, int ny, real_t sx, real_t sy)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*4 triangles.
void GetLocalTriToWdgTransformation(IsoparametricTransformation &loc, int i) const
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
int SpaceDimension() const
Dimension of the physical space containing the mesh.
void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf, bool spacing=false)
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0, bool sfc_ordering=true)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
void GetCharacteristics(real_t &h_min, real_t &h_max, real_t &kappa_min, real_t &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
void ReadCubit(const std::string &filename, int &curved, int &read_gf)
Load a mesh from a Genesis file.
void SetNode(int i, const real_t *coord)
void GetNodes(Vector &node_coord) const
AttributeSets attribute_sets
Named sets of element attributes.
virtual void UniformRefinement3D()
Refine a mixed 3D mesh uniformly.
static int ComposeTriOrientations(int ori_a_b, int ori_b_c)
void PrintBdrVTU(std::string fname, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0)
void AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
int AddBdrTriangle(int v1, int v2, int v3, int attr=1)
void GetGeometricParametersFromJacobian(const DenseMatrix &J, real_t &volume, Vector &aspr, Vector &skew, Vector &ori) const
Computes geometric parameters associated with a Jacobian matrix in 2D/3D. These parameters are (1) Ar...
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
static int DecodeFaceInfoOrientation(int info)
Given a "face info int", return the face orientation.
void GetHilbertElementOrdering(Array< int > &ordering)
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
void AddQuadAs5QuadsWithPoints(int *vi, int attr=1)
Adds 5 quadrilaterals to the mesh by splitting a quadrilateral given by 4 vertices vi.
void ReadVTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo)
void PrintVTU(std::ostream &os, int ref=1, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0, bool bdr_elements=false)
virtual void UniformRefinement2D()
Refine a mixed 2D mesh uniformly.
static Mesh MakePeriodic(const Mesh &orig_mesh, const std::vector< int > &v2v)
Create a periodic mesh by identifying vertices of orig_mesh.
Element::Type GetFaceElementType(int Face) const
int CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i....
void AverageVertices(const int *indexes, int n, int result)
Averages the vertices with given indexes and saves the result in vertices[result].
real_t GetElementVolume(int i)
static Mesh LoadFromFile(const std::string &filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
static const int vtk_quadratic_hex[27]
void Swap(Mesh &other, bool non_geometry)
Array< Triple< int, int, int > > tmp_vertex_parents
virtual void GenerateBoundaryElements()
static void GetElementArrayEdgeTable(const Array< Element * > &elem_array, const DSTable &v_to_v, Table &el_to_edge)
std::vector< int > CreatePeriodicVertexMapping(const std::vector< Vector > &translations, real_t tol=1e-8) const
Creates a mapping v2v from the vertex indices of the mesh such that coincident vertices under the giv...
void Bisection(int i, const DSTable &, int *, int *, int *)
Bisect a triangle: element with index i is bisected.
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
void UniformRefinement2D_base(bool update_nodes=true)
IsoparametricTransformation BdrTransformation
void GetLocalSegToTriTransformation(IsoparametricTransformation &loc, int i) const
void DegreeElevate(int rel_degree, int degree=16)
int FindCoarseElement(int i)
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
void GetElementCenter(int i, Vector ¢er)
void AddHexAsPyramids(const int *vi, int attr=1)
Adds 6 pyramids to the mesh by splitting a hexahedron given by 8 vertices vi.
static void PrintElementsByGeometry(int dim, const Array< int > &num_elems_by_geom, std::ostream &os)
Auxiliary method used by PrintCharacteristics().
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Adds a hexahedron to the mesh given by 8 vertices v1 through v8.
static IntegrationPoint TransformBdrElementToFace(Geometry::Type geom, int o, const IntegrationPoint &ip)
For the vertex (1D), edge (2D), or face (3D) of a boundary element with the orientation o,...
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
int GetNBE() const
Returns number of boundary elements.
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
void AddQuadFaceElement(int lf, int gf, int el, int v0, int v1, int v2, int v3)
void ReadMFEMMesh(std::istream &input, int version, int &curved)
void KnotRemove(Array< Vector * > &kv)
For NURBS meshes, remove the knots in kv, for each direction.
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr) const
Builds the transformation defining the i-th edge element in EdTr. EdTr must be allocated in advance a...
void GetElementData(const Array< Element * > &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
void PrintSurfaces(const Table &Aface_face, std::ostream &os) const
Print set of disjoint surfaces:
bool IsSlaveFace(const FaceInfo &fi) const
void FreeElement(Element *E)
Array< Element * > boundary
virtual void MarkTetMeshForRefinement(const DSTable &v_to_v)
void GetPointMatrix(int i, DenseMatrix &pointmat) const
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
void GetLocalTriToTetTransformation(IsoparametricTransformation &loc, int i) const
virtual void PrintXG(std::ostream &os=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
NCMesh * ncmesh
Optional nonconforming mesh extension.
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
virtual void NURBSUniformRefinement(int rf=2, real_t tol=1.0e-12)
Refine NURBS mesh, with an optional refinement factor, generally anisotropic.
void DebugDump(std::ostream &os) const
Output an NCMesh-compatible debug dump.
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
bool RefineByError(const Array< real_t > &elem_error, real_t threshold, int nonconforming=-1, int nc_limit=0)
static Mesh MakeCartesian3DWith24TetsPerHex(int nx, int ny, int nz, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz*24 tetrahedrons.
STable3D * GetFacesTable()
Table * GetFaceEdgeTable() const
void EnsureNCMesh(bool simplices_nonconforming=false)
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
virtual MFEM_DEPRECATED void ReorientTetMesh()
void PrintVTK(std::ostream &os)
virtual void MarkNamedBoundaries(const std::string &set_name, Array< int > &bdr_marker) const
Mark boundary attributes in the named set.
void MoveNodes(const Vector &displacements)
Array< GeometricFactors * > geom_factors
Optional geometric factors.
void SetMeshGen()
Determine the mesh generator bitmask meshgen, see MeshGenerator().
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
void Make2D(int nx, int ny, Element::Type type, real_t sx, real_t sy, bool generate_edges, bool sfc_ordering)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
void AddHexAsTets(const int *vi, int attr=1)
Adds 6 tetrahedrons to the mesh by splitting a hexahedron given by 8 vertices vi.
FaceElementTransformations FaceElemTr
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
void GetLocalQuadToPyrTransformation(IsoparametricTransformation &loc, int i) const
void AddHexAsWedges(const int *vi, int attr=1)
Adds 2 wedges to the mesh by splitting a hexahedron given by 8 vertices vi.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Element * ReadElement(std::istream &input)
virtual bool HasBoundaryElements() const
Checks if the mesh has boundary elements.
void ScaleSubdomains(real_t sf)
void GetVertexToVertexTable(DSTable &) const
void GetLocalQuadToHexTransformation(IsoparametricTransformation &loc, int i) const
void Transform(void(*f)(const Vector &, Vector &))
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Geometry::Type GetTypicalFaceGeometry() const
If the local mesh is not empty, return GetFaceGeometry(0); otherwise return a typical face geometry p...
Geometry::Type GetElementBaseGeometry(int i) const
virtual void UnmarkInternalBoundaries(Array< int > &bdr_marker, bool excl=true) const
Unmark boundary attributes of internal boundaries.
@ LocalSlaveNonconforming
@ SharedSlaveNonconforming
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones.
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
void ChangeVertexDataOwnership(real_t *vertices, int len_vertices, bool zerocopy=false)
Set the internal Vertex array to point to the given vertices array without assuming ownership of the ...
Table * GetVertexToElementTable()
void InitRefinementTransforms()
Table * GetEdgeVertexTable() const
int * GeneratePartitioning(int nparts, int part_method=1)
Table * GetFaceToElementTable() const
void MarkTriMeshForRefinement()
void ReadNetgen2DMesh(std::istream &input, int &curved)
Array< Element * > elements
Array< int > attributes
A list of all unique element attributes used by the Mesh.
void AddPointFaceElement(int lf, int gf, int el)
Used in GenerateFaces()
void RemoveInternalBoundaries()
virtual void NonconformingRefinement(const Array< Refinement > &refinements, int nc_limit=0)
This function is not public anymore. Use GeneralRefinement instead.
void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
void MoveVertices(const Vector &displacements)
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
void DeleteGeometricFactors()
Destroy all GeometricFactors stored by the Mesh.
const Table & ElementToFaceTable() const
void AddQuadAs4TrisWithPoints(int *vi, int attr=1)
Adds 4 triangles to the mesh by splitting a quadrilateral given by 4 vertices vi.
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
void KnotInsert(Array< KnotVector * > &kv)
For NURBS meshes, insert the new knots in kv, for each direction.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
void OnMeshUpdated(Mesh *mesh)
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
const CoarseFineTransformations & GetRefinementTransforms() const
int Dimension() const
Return the dimension of the NCMesh.
BlockArray< Element > elements
static void GridSfcOrdering3D(int width, int height, int depth, Array< int > &coords)
Array< int > leaf_elements
finest elements, in Mesh ordering (+ ghosts)
virtual void LimitNCLevel(int max_nc_level)
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
virtual void Derefine(const Array< int > &derefs)
Array< real_t > coordinates
bool IsGhost(const Element &el) const
Return true if the Element el is a ghost element.
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
void Print(std::ostream &out, const std::string &comments="") const
int spaceDim
dimensions of the elements and the vertex coordinates
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
const Table & GetDerefinementTable()
void SetAttribute(int i, int attr)
Set the attribute of leaf element i, which is a Mesh element index.
int SpaceDimension() const
Return the space dimension of the NCMesh.
virtual void Refine(const Array< Refinement > &refinements)
static void GridSfcOrdering2D(int width, int height, Array< int > &coords)
NURBSExtension generally contains multiple NURBSPatch objects spanning an entire Mesh....
int GetNBE() const
Return the number of active boundary elements.
void GetCoarseningFactors(Array< int > &f) const
void SetPatchAttribute(int i, int attr)
Set the attribute for patch i, which is set to all elements in the patch.
const Array< int > & GetPatchElements(int patch)
Return the array of indices of all elements in patch patch.
void UniformRefinement(int rf=2)
Refine with optional refinement factor rf. Uniform means refinement is done everywhere by the same fa...
void Print(std::ostream &os, const std::string &comments="") const
Writes all patch data to the stream os.
void SetPatchBdrAttribute(int i, int attr)
Set the attribute for patch boundary element i to attr, which is set to all boundary elements in the ...
void SetCoordsFromPatches(Vector &Nodes)
Set FE coordinates in Nodes, using data from patches, and erase patches.
void Coarsen(int cf=2, real_t tol=1.0e-12)
int GetPatchAttribute(int i) const
Get the attribute for patch i, which is set to all elements in the patch.
void GetElementTopo(Array< Element * > &elements) const
Generate the active mesh elements and return them in elements.
const Array< int > & GetPatchBdrElements(int patch)
Return the array of indices of all boundary elements in patch patch.
int GetNKV() const
Return the number of KnotVectors.
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Get the local to global vertex index map lvert_vert.
bool HavePatches() const
Return true if at least 1 patch is defined, false otherwise.
void GetBdrElementTopo(Array< Element * > &boundary) const
Generate the active mesh boundary elements and return them in boundary.
void KnotRemove(Array< Vector * > &kv, real_t tol=1.0e-12)
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Get the local to global element index map lelem_elem.
void KnotInsert(Array< KnotVector * > &kv)
Insert knots from kv into all KnotVectors in all patches. The size of kv should be the same as knotVe...
int GetOrder() const
If all KnotVector orders are identical, return that number. Otherwise, return NURBSFECollection::Vari...
int GetNV() const
Return the local number of active vertices.
void ConvertToPatches(const Vector &Nodes)
Define patches in IKJ (B-net) format, using FE coordinates in Nodes.
int Dimension() const
Return the dimension of the reference space (not physical space).
void SetKnotsFromPatches()
Set KnotVectors from patches and construct mesh and space data.
int GetNE() const
Return the number of active elements.
int GetPatchBdrAttribute(int i) const
Get the attribute for boundary patch element i, which is set to all boundary elements in the patch.
void DegreeElevate(int rel_degree, int degree=16)
Call DegreeElevate for all KnotVectors of all patches. For each KnotVector, the new degree is max(old...
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Class for standard nodal finite elements.
Class used to extrude the nodes of a mesh.
void SetLayer(const int l)
NodeExtrudeCoefficient(const int dim, const int n_, const real_t s_)
void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip) override
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Class for parallel meshes.
Parallel version of NURBSExtension.
Data type Pyramid element.
Piecewise-(bi)quadratic continuous finite elements.
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
@ VALUES
Evaluate the values at quadrature points.
@ DERIVATIVES
Evaluate the derivatives at quadrature points.
@ DETERMINANTS
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
Data type quadrilateral element.
Symmetric 3D Table stored as an array of rows each of which has a stack of column,...
int Push(int r, int c, int f)
Check to see if this entry is in the table and add it to the table if it is not there....
int NumberOfElements()
Return the number of elements added to the table.
int Push4(int r, int c, int f, int t)
Check to see if this entry is in the table and add it to the table if it is not there....
Data type line segment element.
void GetVertices(Array< int > &v) const override
Get the indices defining the vertices.
void Start()
Start the stopwatch. The elapsed time is not cleared.
double UserTime()
Return the number of user seconds elapsed since the stopwatch was started.
void AddConnections(int r, const int *c, int nc)
void SetSize(int dim, int connections_per_row)
Set the size and the number of connections for the table.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
void AddConnection(int r, int c)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
int Size() const
Returns the number of TYPE I elements.
int Size_of_connections() const
void AddColumnsInRow(int r, int ncol)
void AddAColumnInRow(int r)
void SetDims(int rows, int nnz)
Data type tetrahedron element.
void Init(int ind1, int ind2, int ind3, int ind4, int attr=1, int ref_flag=0)
Initialize the vertex indices and the attribute of a Tetrahedron.
void PushTransform(int tr) override
Add 'tr' to the current chain of coarse-fine transformations.
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag) const
int GetRefinementFlag() const
void SetVertices(const Array< int > &v) override
Set the indices defining the vertices.
void ResetTransform(int tr) override
Set current coarse-fine transformation number.
unsigned GetTransform() const override
Return current coarse-fine transformation.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void GetMarkedFace(const int face, int *fv) const
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
void GetVertices(Array< int > &v) const override
Get the indices defining the vertices.
Data type triangle element.
void SetVertices(const Array< int > &v) override
Set the indices defining the vertices.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void PushTransform(int tr) override
Add 'tr' to the current chain of coarse-fine transformations.
void ResetTransform(int tr) override
Set current coarse-fine transformation number.
unsigned GetTransform() const override
Return current coarse-fine transformation.
void GetVertices(Array< int > &v) const override
Get the indices defining the vertices.
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
A general vector function coefficient.
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
real_t Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
void NewDataAndSize(real_t *d, int s)
Set the Vector data and size, deleting the old data, if owned.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
void cross3D(const Vector &vin, Vector &vout) const
real_t DistanceTo(const real_t *p) const
Compute the Euclidean distance to another vector.
void Print(const Mesh &mesh, const adios2stream::mode print_mode=mode::sync)
const std::string filename
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
int index(int i, int j, int nx, int ny)
void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void METIS_PartGraphVKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
Linear1DFiniteElement SegmentFE
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
PointFiniteElement PointFE
TriLinear3DFiniteElement HexahedronFE
void Swap(Array< T > &, Array< T > &)
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
GeometryRefiner GlobGeometryRefiner
int FindRoots(const Vector &z, Vector &x)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
real_t rand_real()
Generate a random real_t number in the interval [0,1) using rand().
void add(const Vector &v1, const Vector &v2, Vector &v)
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
void ShiftRight(int &a, int &b, int &c)
MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
Mesh * Extrude1D(Mesh *mesh, const int ny, const real_t sy, const bool closed)
Extrude a 1D mesh.
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
void WriteBase64WithSizeAndClear(std::ostream &os, std::vector< char > &buf, int compression_level)
Encode in base 64 (and potentially compress) the given data, write it to the output stream (with a he...
Mesh * Extrude2D(Mesh *mesh, const int nz, const real_t sz)
Extrude a 2D mesh.
VTKFormat
Data array format for VTK and VTU files.
@ ASCII
Data arrays will be written in ASCII format.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
void FindTMax(Vector &c, Vector &x, real_t &tmax, const real_t factor, const int Dim)
BiLinear2DFiniteElement QuadrilateralFE
void WriteBinaryOrASCII(std::ostream &os, std::vector< char > &buf, const T &val, const char *suffix, VTKFormat format)
Write either ASCII data to the stream or binary data to the buffer depending on the given format.
MFEM_EXPORT class LinearPyramidFiniteElement PyramidFE
const T & AsConst(const T &a)
Utility function similar to std::as_const in c++17.
double bisect(ElementTransformation &Tr, Coefficient *LvlSet)
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
void SortPairs(Pair< A, B > *pairs, int size)
Sort an array of Pairs with respect to the first element.
MemoryType
Memory types supported by MFEM.
@ HOST
Host memory; using new[] and delete[].
void CreateVTKElementConnectivity(Array< int > &con, Geometry::Type geom, int ref)
Create the VTK element connectivity array for a given element geometry and refinement level.
void FindPartitioningComponents(Table &elem_elem, const Array< int > &partitioning, Array< int > &component, Array< int > &num_comp)
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
@ NATIVE
Native ordering as defined by the FiniteElement.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
MFEM_EXPORT Linear2DFiniteElement TriangleFE
MPI_Comm GetGlobalMPI_Comm()
Get MFEM's "global" MPI communicator.
void XYZ_VectorFunction(const Vector &p, Vector &v)
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
real_t p(const Vector &x, real_t t)
Helper struct for defining a connectivity table, see Table::MakeFromList.
Defines the position of a fine element within a coarse element.
int parent
Coarse Element index in the coarse mesh.
static const int FaceVert[NumFaces][MaxFaceVert]
static const int Edges[NumEdges][2]
static const int FaceVert[NumFaces][MaxFaceVert]
static const int Edges[NumEdges][2]
static const int Edges[NumEdges][2]
static const int FaceVert[NumFaces][MaxFaceVert]
static const int Orient[NumOrient][NumVert]
static const int Orient[NumOrient][NumVert]
static const int Edges[NumEdges][2]
static const int Edges[NumEdges][2]
static const int FaceVert[NumFaces][MaxFaceVert]
static const int Orient[NumOrient][NumVert]
static const int Edges[NumEdges][2]
static const int Orient[NumOrient][NumVert]
EntityHelper(int dim_, const Array< int >(&entity_to_vertex_)[Geometry::NumGeom])
Entity FindEntity(int bytype_entity_id)
entity_to_vertex_type & entity_to_vertex
int geom_offsets[Geometry::NumGeom+1]
This structure stores the low level information necessary to interpret the configuration of elements ...
Lists all edges/faces in the nonconforming mesh.
Nonconforming edge/face within a bigger edge/face.
static const int HighOrderMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to arbitrary-order Lagrange VTK geometries.
static const int QuadraticMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to legacy quadratic VTK geometries/.
static const int * VertexPermutation[Geometry::NUM_GEOMETRIES]
Permutation from MFEM's vertex ordering to VTK's vertex ordering.
static const int Map[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to linear VTK geometries.
std::array< int, NCMesh::MaxFaceNodes > nodes