37#include <unordered_map>
38#include <unordered_set>
43#if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
48#if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
53 int*,
int*,
int*,
int*,
int*,
idxtype*);
55 int*,
int*,
int*,
int*,
int*,
idxtype*);
57 int*,
int*,
int*,
int*,
int*,
idxtype*);
122 return sqrt((d_hat * d_hat) / (dir * dir));
161 if (coord[d] < min(d)) { min(d) = coord[d]; }
162 if (coord[d] > max(d)) { max(d) = coord[d]; }
168 const bool use_boundary =
false;
177 for (
int i = 0; i < ne; i++)
194 for (
int j = 0; j < pointmat.
Width(); j++)
196 for (
int d = 0; d < pointmat.
Height(); d++)
198 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
199 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
221 h_max = kappa_max = -h_min;
222 if (
dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
230 if (Vh) { (*Vh)(i) = h; }
231 if (Vk) { (*Vk)(i) =
kappa; }
233 if (h < h_min) { h_min = h; }
234 if (h > h_max) { h_max = h; }
248 if (!num_elems_by_geom[g]) {
continue; }
249 if (!first) { os <<
" + "; }
257 real_t h_min, h_max, kappa_min, kappa_max;
259 os <<
"Mesh Characteristics:";
264 num_elems_by_geom = 0;
265 for (
int i = 0; i <
GetNE(); i++)
276 <<
"Number of vertices : " <<
GetNV() <<
'\n'
277 <<
"Number of elements : " <<
GetNE() <<
'\n'
278 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n';
283 <<
"Number of vertices : " <<
GetNV() <<
'\n'
284 <<
"Number of elements : " <<
GetNE() <<
'\n'
285 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
286 <<
"h_min : " << h_min <<
'\n'
287 <<
"h_max : " << h_max <<
'\n';
292 <<
"Number of vertices : " <<
GetNV() <<
'\n'
293 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
294 <<
"Number of elements : " <<
GetNE() <<
" -- ";
297 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
299 <<
"h_min : " << h_min <<
'\n'
300 <<
"h_max : " << h_max <<
'\n'
301 <<
"kappa_min : " << kappa_min <<
'\n'
302 <<
"kappa_max : " << kappa_max <<
'\n';
307 num_bdr_elems_by_geom = 0;
308 for (
int i = 0; i <
GetNBE(); i++)
313 num_faces_by_geom = 0;
320 <<
"Number of vertices : " <<
GetNV() <<
'\n'
321 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
322 <<
"Number of faces : " <<
GetNFaces() <<
" -- ";
325 <<
"Number of elements : " <<
GetNE() <<
" -- ";
328 <<
"Number of bdr elem : " <<
GetNBE() <<
" -- ";
332 <<
"h_min : " << h_min <<
'\n'
333 <<
"h_max : " << h_max <<
'\n'
334 <<
"kappa_min : " << kappa_min <<
'\n'
335 <<
"kappa_max : " << kappa_max <<
'\n';
337 os <<
'\n' << std::flush;
353 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
356 MFEM_ABORT(
"Unknown element type");
385 for (
int j = 0; j < n; j++)
387 pm(k,j) =
nodes(vdofs[n*k+j]);
436 int nv =
elements[i]->GetNVertices();
437 const int *v =
elements[i]->GetVertices();
442 for (
int j = 0; j < nv; j++)
444 pm(k, j) =
nodes(k*n+v[j]);
458 for (
int j = 0; j < n; j++)
460 pm(k,j) =
nodes(vdofs[n*k+j]);
494 for (
int j = 0; j < n; j++)
503 int elem_id, face_info;
519 "Mesh requires nodal Finite Element.");
528 ElTr->
SetFE(face_el);
550 const int *v = (
Dim == 1) ? &FaceNo :
faces[FaceNo]->GetVertices();
551 const int nv = (
Dim == 1) ? 1 :
faces[FaceNo]->GetNVertices();
555 for (
int j = 0; j < nv; j++)
575 for (
int j = 0; j < n; j++)
577 pm(i, j) =
nodes(vdofs[n*i+j]);
596 "Mesh requires nodal Finite Element.");
626 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
643 for (
int j = 0; j < nv; j++)
663 for (
int j = 0; j < n; j++)
665 pm(i, j) =
nodes(vdofs[n*i+j]);
668 EdTr->
SetFE(edge_el);
678 MFEM_VERIFY(faces_e.
Size() > 0,
"Edge not found in any face!");
679 const int face_no = faces_e[0];
684 const int local_idx = edges_f.
Find(EdgeNo);
685 MFEM_ASSERT(local_idx >= 0,
"Edge not found on the face!");
686 const int edge_ori = oris_f[local_idx] > 0 ? 0 : 1;
705 MFEM_ABORT(
"Unsupported face type for edge transformation!");
725 MFEM_ABORT(
"Unsupported finite element collection.");
747 EdTr->
SetFE(edge_el);
787 for (
int j = 0; j < 2; j++)
789 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
790 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
807 for (
int j = 0; j < 2; j++)
809 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
810 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
829 for (
int j = 0; j < 3; j++)
832 locpm(0, j) = vert.
x;
833 locpm(1, j) = vert.
y;
834 locpm(2, j) = vert.
z;
846 MFEM_VERIFY(i < 128,
"Local face index " << i/64
847 <<
" is not a triangular face of a wedge.");
855 for (
int j = 0; j < 3; j++)
858 locpm(0, j) = vert.
x;
859 locpm(1, j) = vert.
y;
860 locpm(2, j) = vert.
z;
871 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
872 <<
" is not a triangular face of a pyramid.");
880 for (
int j = 0; j < 3; j++)
883 locpm(0, j) = vert.
x;
884 locpm(1, j) = vert.
y;
885 locpm(2, j) = vert.
z;
902 for (
int j = 0; j < 4; j++)
905 locpm(0, j) = vert.
x;
906 locpm(1, j) = vert.
y;
907 locpm(2, j) = vert.
z;
919 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
920 <<
" is not a quadrilateral face of a wedge.");
926 for (
int j = 0; j < 4; j++)
929 locpm(0, j) = vert.
x;
930 locpm(1, j) = vert.
y;
931 locpm(2, j) = vert.
z;
942 MFEM_VERIFY(i < 64,
"Local face index " << i/64
943 <<
" is not a quadrilateral face of a pyramid.");
949 for (
int j = 0; j < 4; j++)
952 locpm(0, j) = vert.
x;
953 locpm(1, j) = vert.
y;
954 locpm(2, j) = vert.
z;
1004 std::unordered_map<int, int> f_to_be;
1005 for (
int i = 0; i <
GetNBE(); ++i)
1015 for (
int f = 0;
f < nf; ++
f)
1022 auto iter = f_to_be.find(
f);
1023 if (iter != f_to_be.end())
1025 const int be = iter->second;
1047 for (
int i = 0; i <
GetNE(); ++i)
1051 "Negative attribute on element " << i);
1063 ifidcs.reserve(fidcs.Size());
1087const std::unordered_map<int, int> &
1150 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
1151 "face type " << face_type
1152 <<
" and element type " << elem_type <<
"\n");
1171 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
1172 "face type " << face_type
1173 <<
" and element type " << elem_type <<
"\n");
1205 FElTr.
Elem1 = &ElTr1;
1218 { MFEM_ABORT(
"NURBS mesh not supported!"); }
1221 FElTr.
Elem2 = &ElTr2;
1269 mfem::out <<
"\nInternal error: face id = " << FaceNo
1270 <<
", dist = " << dist <<
'\n';
1272 MFEM_ABORT(
"internal error");
1330 const FaceInfo &fi,
bool is_ghost)
const
1332#ifdef MFEM_THREAD_SAFE
1337 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1349 std::swap(composition(0,0), composition(0,1));
1350 std::swap(composition(1,0), composition(1,1));
1407 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1496 res.
Elem1No = element[0].index;
1497 res.Elem2No = element[1].index;
1498 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1499 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1500 res.NCFace = ncface;
1503 res.Elem1No = element[0].index;
1504 res.Elem2No = element[1].index;
1505 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1506 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1507 res.NCFace = ncface;
1510 res.Elem1No = element[0].index;
1511 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1514 res.Elem1No = element[0].index;
1516 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1517 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1520 res.Elem1No = element[0].index;
1521 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1524 res.Elem1No = element[0].index;
1526 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1527 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1532 res.Elem1No = element[0].index;
1534 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1535 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1543 os <<
"face topology=";
1553 os <<
"Non-conforming";
1560 os <<
"element[0].location=";
1574 os <<
"element[1].location=";
1588 os <<
"element[0].conformity=";
1605 os <<
"element[1].conformity=";
1622 os <<
"element[0].index=" << info.
element[0].
index <<
'\n'
1623 <<
"element[1].index=" << info.
element[1].
index <<
'\n'
1628 <<
"ncface=" << info.
ncface << std::endl;
1660 return faces[Face]->GetGeometryType();
1663 const int nc_face_id =
faces_info[Face].NCFace;
1665 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1728 "Could not determine a typical element Geometry!");
1737 face_marker.
SetSize(num_faces);
1739 for (
int f = 0;
f < num_faces;
f++)
1756 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1757 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1759 Array<bool> interior_bdr(max_bdr_attr); interior_bdr =
false;
1760 Array<bool> exterior_bdr(max_bdr_attr); exterior_bdr =
false;
1764 for (
int be = 0; be <
boundary.Size(); be++)
1766 const int bea =
boundary[be]->GetAttribute();
1768 if (bdr_marker[bea-1] != 0)
1774 interior_bdr[bea-1] =
true;
1778 exterior_bdr[bea-1] =
true;
1785 for (
int b = 0;
b < max_bdr_attr;
b++)
1787 if (bdr_marker[
b] != 0 && interior_bdr[
b])
1789 if (!excl || !exterior_bdr[
b])
1803 "Named set is not defined in this mesh!");
1805 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1809 for (
int b = 0;
b < max_bdr_attr;
b++)
1822 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1823 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1825 Array<bool> interior_bdr(max_bdr_attr); interior_bdr =
false;
1826 Array<bool> exterior_bdr(max_bdr_attr); exterior_bdr =
false;
1830 for (
int be = 0; be <
boundary.Size(); be++)
1832 const int bea =
boundary[be]->GetAttribute();
1838 interior_bdr[bea-1] =
true;
1842 exterior_bdr[bea-1] =
true;
1848 for (
int b = 0;
b < max_bdr_attr;
b++)
1850 if (bdr_marker[
b] == 0 && exterior_bdr[
b])
1852 if (!excl || !interior_bdr[
b])
1866 "Named set is not defined in this mesh!");
1867 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1868 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1872 for (
int b = 0;
b < max_bdr_attr;
b++)
1951 for (
int i = 0; i <
faces.Size(); i++)
1979#ifdef MFEM_USE_MEMALLOC
2018 if (bdr_face_attrs_changed)
2023 std::set<int> attribs;
2024 for (
int i = 0; i <
GetNBE(); i++)
2034 MFEM_WARNING(
"Non-positive attributes on the boundary!");
2038 if (elem_attrs_changed)
2051 MFEM_WARNING(
"Non-positive attributes in the domain!");
2074static void CheckEnlarge(
Array<T> &array,
int size)
2076 if (size >= array.Size()) { array.SetSize(size + 1); }
2099 "invalid 'coords' size: " << coords.
Size());
2112 for (
int j = 0; j < 3; j++)
2114 vi[j] = (vp1[j] + vp2[j]) * 0.5;
2123 for (
int i = 0; i < nverts; i++)
2126 for (
int j = 0; j <
dim; j++)
2180 int vi[4] = {v1, v2, v3, v4};
2187#ifdef MFEM_USE_MEMALLOC
2227int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
2232 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
2245 static const int hex_to_tet[6][4] =
2247 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
2248 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
2252 for (
int i = 0; i < 6; i++)
2254 for (
int j = 0; j < 4; j++)
2256 ti[j] = vi[hex_to_tet[i][j]];
2264 static const int hex_to_wdg[2][6] =
2266 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
2270 for (
int i = 0; i < 2; i++)
2272 for (
int j = 0; j < 6; j++)
2274 ti[j] = vi[hex_to_wdg[i][j]];
2282 static const int hex_to_pyr[6][5] =
2284 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
2285 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
2289 for (
int i = 0; i < 6; i++)
2291 for (
int j = 0; j < 5; j++)
2293 ti[j] = vi[hex_to_pyr[i][j]];
2302 static const int quad_to_tri[4][2] =
2304 {0, 1}, {1, 2}, {2, 3}, {3, 0}
2310 ti[2] = elem_center_index;
2311 for (
int i = 0; i < num_faces; i++)
2313 for (
int j = 0; j < 2; j++)
2315 ti[j] = vi[quad_to_tri[i][j]];
2324 static const int quad_faces[4][2] =
2326 {0, 1}, {1, 2}, {2, 3}, {3, 0}
2330 for (
int i = 0; i < 4; i++)
2339 real_t r = 0.25, s = 0.25;
2340 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2341 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2346 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2347 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2352 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2353 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2358 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2359 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2363 static const int quad_faces_new[4][2] =
2365 { 1, 0}, { 2, 1}, { 3, 2}, { 0, 3}
2369 for (
int i = 0; i < num_faces; i++)
2371 for (
int j = 0; j < 2; j++)
2373 ti[j] = vi[quad_faces[i][j]];
2374 ti[j+2] = vnew_index[quad_faces_new[i][j]];
2382 std::map<std::array<int, 4>,
int> &hex_face_verts,
2388 return std::array<int, 4> {v[0], v[1], v[2], v[3]};
2392 static const int hex_to_tet[6][4] =
2394 { 0, 1, 2, 3 }, { 1, 2, 6, 5 }, { 5, 4, 7, 6},
2395 { 0, 1, 5, 4 }, { 2, 3, 7, 6 }, { 0,3, 7, 4}
2403 static const int tet_face[4][2] =
2405 {0, 1}, {1, 2}, {3, 2}, {3, 0}
2408 for (
int i = 0; i < num_faces; i++)
2410 for (
int j = 0; j < 4; j++)
2412 flist[j] = vi[hex_to_tet[i][j]];
2414 int face_center_index;
2416 auto t = get4arraysorted(flist);
2417 auto it = hex_face_verts.find(t);
2418 if (it == hex_face_verts.end())
2421 flist.
Size(), 3) - 1;
2422 hex_face_verts.insert({t, face_center_index});
2426 face_center_index = it->second;
2429 fti[2] = face_center_index;
2430 fti[3] = elem_center_index;
2431 for (
int j = 0; j < 4; j++)
2433 for (
int k = 0; k < 2; k++)
2435 fti[k] = flist[tet_face[j][k]];
2460 MFEM_ASSERT(bdr_elems.
Size() == new_be_to_face.
Size(),
"wrong size");
2461 for (
int i = 0; i < bdr_elems.
Size(); i++)
2512 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
2515 for (
int i = 0; i < 2; i++)
2517 for (
int j = 0; j < 3; j++)
2519 ti[j] = vi[quad_to_tri[i][j]];
2555 for (
int i = 0, j = 0; i <
faces_info.Size(); i++)
2571 "incorrect number of vertices: preallocated: " <<
vertices.Size()
2574 "incorrect number of elements: preallocated: " <<
elements.Size()
2577 "incorrect number of boundary elements: preallocated: "
2611 bool fix_orientation)
2614 if (fix_orientation)
2644 GeckoProgress(
real_t limit) : limit(limit) { sw.
Start(); }
2645 bool quit()
const override {
return limit > 0 && sw.
UserTime() > limit; }
2648class GeckoVerboseProgress :
public GeckoProgress
2654 GeckoVerboseProgress(
real_t limit) : GeckoProgress(limit) {}
2656 void beginorder(
const Graph* graph, Float cost)
const override
2657 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2658 void endorder(
const Graph* graph, Float cost)
const override
2659 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2661 void beginiter(
const Graph* graph,
2662 uint iter, uint maxiter, uint window)
const override
2664 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
2665 << window << std::flush;
2667 void enditer(
const Graph* graph, Float mincost, Float cost)
const override
2668 {
mfem::out <<
", cost = " << cost << endl; }
2673 int iterations,
int window,
2674 int period,
int seed,
bool verbose,
2680 GeckoProgress progress(time_limit);
2681 GeckoVerboseProgress vprogress(time_limit);
2684 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2692 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2694 const int *neighid = my_el_to_el.
GetRow(elemid);
2695 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2697 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2702 graph.
order(&functional, iterations, window, period, seed,
2703 verbose ? &vprogress : &progress);
2709 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2712 return graph.
cost();
2724 : coord(coord), dir(dir), points(points), mid(mid) {}
2726 bool operator()(
int i)
const
2728 return (points[3*i + coord] < mid) != dir;
2732static void HilbertSort2D(
int coord1,
2735 const Array<real_t> &points,
int *beg,
int *end,
2738 if (end - beg <= 1) {
return; }
2740 real_t xmid = (xmin + xmax)*0.5;
2741 real_t ymid = (ymin + ymax)*0.5;
2743 int coord2 = (coord1 + 1) % 2;
2746 int *p0 = beg, *p4 = end;
2747 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2748 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2749 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2753 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2754 ymin, xmin, ymid, xmid);
2756 if (p1 != p0 || p2 != p4)
2758 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2759 xmin, ymid, xmid, ymax);
2761 if (p2 != p0 || p3 != p4)
2763 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2764 xmid, ymid, xmax, ymax);
2768 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2769 ymid, xmax, ymin, xmid);
2773static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2774 const Array<real_t> &points,
int *beg,
int *end,
2778 if (end - beg <= 1) {
return; }
2780 real_t xmid = (xmin + xmax)*0.5;
2781 real_t ymid = (ymin + ymax)*0.5;
2782 real_t zmid = (zmin + zmax)*0.5;
2784 int coord2 = (coord1 + 1) % 3;
2785 int coord3 = (coord1 + 2) % 3;
2788 int *p0 = beg, *p8 = end;
2789 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2790 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2791 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2792 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2793 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2794 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2795 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2799 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2800 zmin, xmin, ymin, zmid, xmid, ymid);
2802 if (p1 != p0 || p2 != p8)
2804 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2805 ymin, zmid, xmin, ymid, zmax, xmid);
2807 if (p2 != p0 || p3 != p8)
2809 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2810 ymid, zmid, xmin, ymax, zmax, xmid);
2812 if (p3 != p0 || p4 != p8)
2814 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2815 xmin, ymax, zmid, xmid, ymid, zmin);
2817 if (p4 != p0 || p5 != p8)
2819 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2820 xmid, ymax, zmid, xmax, ymid, zmin);
2822 if (p5 != p0 || p6 != p8)
2824 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2825 ymax, zmid, xmax, ymid, zmax, xmid);
2827 if (p6 != p0 || p7 != p8)
2829 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2830 ymid, zmid, xmax, ymin, zmax, xmid);
2834 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2835 zmid, xmax, ymin, zmin, xmid, ymid);
2849 if (
spaceDim < 3) { points = 0.0; }
2852 for (
int i = 0; i <
GetNE(); i++)
2857 points[3*i + j] = center(j);
2864 indices.
Sort([&](
int a,
int b)
2865 {
return points[3*
a] < points[3*
b]; });
2870 HilbertSort2D(0,
false,
false,
2871 points, indices.
begin(), indices.
end(),
2872 min(0), min(1), max(0), max(1));
2877 HilbertSort3D(0,
false,
false,
false,
2878 points, indices.
begin(), indices.
end(),
2879 min(0), min(1), min(2), max(0), max(1), max(2));
2884 for (
int i = 0; i <
GetNE(); i++)
2886 ordering[indices[i]] = i;
2895 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2900 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
2904 MFEM_VERIFY(ordering.
Size() ==
GetNE(),
"invalid reordering array.")
2943 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2947 old_elem_node_vals[old_elid] =
new Vector(vals);
2953 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2955 int new_elid = ordering[old_elid];
2956 new_elements[new_elid] =
elements[old_elid];
2961 if (reorder_vertices)
2966 vertex_ordering = -1;
2968 int new_vertex_ind = 0;
2969 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2971 int *elem_vert =
elements[new_elid]->GetVertices();
2972 int nv =
elements[new_elid]->GetNVertices();
2973 for (
int vi = 0; vi < nv; ++vi)
2975 int old_vertex_ind = elem_vert[vi];
2976 if (vertex_ordering[old_vertex_ind] == -1)
2978 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2979 new_vertices[new_vertex_ind] =
vertices[old_vertex_ind];
2989 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2991 int *elem_vert =
elements[new_elid]->GetVertices();
2992 int nv =
elements[new_elid]->GetNVertices();
2993 for (
int vi = 0; vi < nv; ++vi)
2995 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
3000 for (
int belid = 0; belid <
GetNBE(); ++belid)
3002 int *be_vert =
boundary[belid]->GetVertices();
3003 int nv =
boundary[belid]->GetNVertices();
3004 for (
int vi = 0; vi < nv; ++vi)
3006 be_vert[vi] = vertex_ordering[be_vert[vi]];
3040 bdr_perm.
Sort([
this](
int a,
int b)
3050 new_boundary[new_i] =
boundary[bdr_perm[new_i]];
3051 new_be_to_face[new_i] =
be_to_face[bdr_perm[new_i]];
3069 int *new_I = new_bel_to_edge->
GetI();
3070 int *new_J = new_bel_to_edge->
GetJ();
3074 const int old_i = bdr_perm[new_i];
3077 for (
int k = 0; k < nrow; ++k)
3079 new_J[new_I[new_i] + k] = old_J[k];
3081 new_I[new_i + 1] = new_I[new_i] + nrow;
3095 nodes_fes->
Update(
false);
3098 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
3100 int new_elid = ordering[old_elid];
3103 delete old_elem_node_vals[old_elid];
3152 length_idx[j].one =
GetLength(i, it.Column());
3153 length_idx[j].two = j;
3162 order[length_idx[i].two] = i;
3177 elements[i]->MarkEdge(v_to_v, order);
3184 boundary[i]->MarkEdge(v_to_v, order);
3191 if (*old_v_to_v && *old_elem_vert)
3198 if (*old_v_to_v == NULL)
3200 bool need_v_to_v =
false;
3208 if (dofs.
Size() > 0)
3220 if (*old_elem_vert == NULL)
3222 bool need_elem_vert =
false;
3224 for (
int i = 0; i <
GetNE(); i++)
3230 if (dofs.
Size() > 1)
3232 need_elem_vert =
true;
3238 *old_elem_vert =
new Table;
3240 for (
int i = 0; i <
GetNE(); i++)
3242 (*old_elem_vert)->AddColumnsInRow(i,
elements[i]->GetNVertices());
3244 (*old_elem_vert)->MakeJ();
3245 for (
int i = 0; i <
GetNE(); i++)
3250 (*old_elem_vert)->ShiftUpI();
3263 const int num_edge_dofs = old_dofs.
Size();
3275 if (num_edge_dofs > 0)
3284 const int old_i = (*old_v_to_v)(i, it.Column());
3285 const int new_i = it.Index();
3286 if (new_i == old_i) {
continue; }
3288 old_dofs.
SetSize(num_edge_dofs);
3289 new_dofs.
SetSize(num_edge_dofs);
3290 for (
int j = 0; j < num_edge_dofs; j++)
3292 old_dofs[j] = offset + old_i * num_edge_dofs + j;
3293 new_dofs[j] = offset + new_i * num_edge_dofs + j;
3297 for (
int j = 0; j < old_dofs.
Size(); j++)
3299 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3311 Table old_face_vertex;
3317 old_face_vertex.
MakeJ();
3320 faces[i]->GetNVertices());
3332 const int *old_v = old_face_vertex.
GetRow(i);
3334 switch (old_face_vertex.
RowSize(i))
3337 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
3341 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
3345 new_fdofs[new_i+1] = old_dofs.
Size();
3352 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
3355 switch (old_face_vertex.
RowSize(i))
3358 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
3359 new_v =
faces[new_i]->GetVertices();
3365 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
3366 new_v =
faces[new_i]->GetVertices();
3374 for (
int j = 0; j < old_dofs.
Size(); j++)
3377 const int old_j = dof_ord[j];
3378 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
3382 for (
int j = 0; j < old_dofs.
Size(); j++)
3384 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3406 for (
int i = 0; i <
GetNE(); i++)
3410 if (old_dofs.
Size() < 2)
3412 offset += old_dofs.
Size();
3416 const int *old_v = old_elem_vert->
GetRow(i);
3417 const int *new_v =
elements[i]->GetVertices();
3424 offset += old_dofs.
Size();
3430 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
3444 <<
" FE collection) are not supported yet!");
3448 MFEM_VERIFY(dof_ord != NULL,
3449 "FE collection '" << fec->
Name()
3450 <<
"' does not define reordering (" << new_or <<
") for "
3453 for (
int j = 0; j < new_dofs.
Size(); j++)
3456 const int old_j = dof_ord[j];
3457 new_dofs[old_j] = offset + j;
3459 offset += new_dofs.
Size();
3462 for (
int j = 0; j < old_dofs.
Size(); j++)
3464 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3502 MFEM_ASSERT(
NURBSext,
"SetPatchAttribute is only for NURBS meshes");
3505 for (
auto e : elems)
3513 MFEM_ASSERT(
NURBSext,
"GetPatchAttribute is only for NURBS meshes");
3519 MFEM_ASSERT(
NURBSext,
"SetPatchBdrAttribute is only for NURBS meshes");
3523 for (
auto be : bdryelems)
3531 MFEM_ASSERT(
NURBSext,
"GetBdrPatchBdrAttribute is only for NURBS meshes");
3537 MFEM_VERIFY(
NURBSext,
"Must be a NURBS mesh");
3572 if (generate_edges == 1)
3590 bool fix_orientation)
3607 if (generate_edges == 1)
3672 bool generate_edges =
true;
3681 MFEM_VERIFY(
ncmesh == NULL,
"");
3716 if (
Dim > 1 && generate_edges)
3780 const bool check_orientation =
true;
3781 const bool curved = (
Nodes != NULL);
3782 const bool may_change_topology =
3784 ( check_orientation && fix_orientation &&
3788 Table *old_elem_vert = NULL;
3790 if (curved && may_change_topology)
3795 if (check_orientation)
3805 if (may_change_topology)
3810 delete old_elem_vert;
3831 for (
int i = 0; i < num_faces; i++)
3834 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
3835 " Interior face with incompatible orientations.");
3846 int NVert, NElem, NBdrElem;
3848 NVert = (nx+1) * (ny+1) * (nz+1);
3849 NElem = nx * ny * nz;
3850 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3859 NBdrElem += 2*nx*ny;
3864 NVert += nx * ny * nz;
3867 InitMesh(3, 3, NVert, NElem, NBdrElem);
3873 for (z = 0; z <= nz; z++)
3875 coord[2] = ((
real_t) z / nz) * sz;
3876 for (y = 0; y <= ny; y++)
3878 coord[1] = ((
real_t) y / ny) * sy;
3879 for (x = 0; x <= nx; x++)
3881 coord[0] = ((
real_t) x / nx) * sx;
3888 for (z = 0; z < nz; z++)
3890 coord[2] = (((
real_t) z + 0.5) / nz) * sz;
3891 for (y = 0; y < ny; y++)
3893 coord[1] = (((
real_t) y + 0.5) / ny) * sy;
3894 for (x = 0; x < nx; x++)
3896 coord[0] = (((
real_t) x + 0.5) / nx) * sx;
3903#define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3904#define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3911 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3913 for (
int k = 0; k < nx*ny*nz; k++)
3920 ind[0] = VTX(x , y , z );
3921 ind[1] = VTX(x+1, y , z );
3922 ind[2] = VTX(x+1, y+1, z );
3923 ind[3] = VTX(x , y+1, z );
3924 ind[4] = VTX(x , y , z+1);
3925 ind[5] = VTX(x+1, y , z+1);
3926 ind[6] = VTX(x+1, y+1, z+1);
3927 ind[7] = VTX(x , y+1, z+1);
3935 for (z = 0; z < nz; z++)
3937 for (y = 0; y < ny; y++)
3939 for (x = 0; x < nx; x++)
3942 ind[0] = VTX(x , y , z );
3943 ind[1] = VTX(x+1, y , z );
3944 ind[2] = VTX(x+1, y+1, z );
3945 ind[3] = VTX(x , y+1, z );
3946 ind[4] = VTX(x , y , z+1);
3947 ind[5] = VTX(x+1, y , z+1);
3948 ind[6] = VTX(x+1, y+1, z+1);
3949 ind[7] = VTX( x, y+1, z+1);
3961 ind[8] = VTXP(x, y, z);
3975 for (y = 0; y < ny; y++)
3977 for (x = 0; x < nx; x++)
3980 ind[0] = VTX(x , y , 0);
3981 ind[1] = VTX(x , y+1, 0);
3982 ind[2] = VTX(x+1, y+1, 0);
3983 ind[3] = VTX(x+1, y , 0);
4000 for (y = 0; y < ny; y++)
4002 for (x = 0; x < nx; x++)
4005 ind[0] = VTX(x , y , nz);
4006 ind[1] = VTX(x+1, y , nz);
4007 ind[2] = VTX(x+1, y+1, nz);
4008 ind[3] = VTX(x , y+1, nz);
4025 for (z = 0; z < nz; z++)
4027 for (y = 0; y < ny; y++)
4030 ind[0] = VTX(0 , y , z );
4031 ind[1] = VTX(0 , y , z+1);
4032 ind[2] = VTX(0 , y+1, z+1);
4033 ind[3] = VTX(0 , y+1, z );
4046 for (z = 0; z < nz; z++)
4048 for (y = 0; y < ny; y++)
4051 ind[0] = VTX(nx, y , z );
4052 ind[1] = VTX(nx, y+1, z );
4053 ind[2] = VTX(nx, y+1, z+1);
4054 ind[3] = VTX(nx, y , z+1);
4067 for (x = 0; x < nx; x++)
4069 for (z = 0; z < nz; z++)
4072 ind[0] = VTX(x , 0, z );
4073 ind[1] = VTX(x+1, 0, z );
4074 ind[2] = VTX(x+1, 0, z+1);
4075 ind[3] = VTX(x , 0, z+1);
4088 for (x = 0; x < nx; x++)
4090 for (z = 0; z < nz; z++)
4093 ind[0] = VTX(x , ny, z );
4094 ind[1] = VTX(x , ny, z+1);
4095 ind[2] = VTX(x+1, ny, z+1);
4096 ind[3] = VTX(x+1, ny, z );
4113 ofstream test_stream(
"debug.mesh");
4115 test_stream.close();
4143 for (
real_t j = 0; j < ny+1; j++)
4145 real_t cy = (j / ny) * sy;
4146 for (
real_t i = 0; i < nx+1; i++)
4148 real_t cx = (i / nx) * sx;
4155 for (
int y = 0; y < ny; y++)
4157 for (
int x = 0; x < nx; x++)
4159 ind[0] = x + y*(nx+1);
4160 ind[1] = x + 1 +y*(nx+1);
4161 ind[2] = x + 1 + (y+1)*(nx+1);
4162 ind[3] = x + (y+1)*(nx+1);
4168 for (
int i = 0; i < nx; i++)
4174 for (
int j = 0; j < ny; j++)
4217 for (
real_t j = 0; j < ny+1; j++)
4219 real_t cy = (j / ny) * sy;
4220 for (
real_t i = 0; i < nx+1; i++)
4222 real_t cx = (i / nx) * sx;
4229 for (
int y = 0; y < ny; y++)
4231 for (
int x = 0; x < nx; x++)
4233 ind[0] = x + y*(nx+1);
4234 ind[1] = x + 1 +y*(nx+1);
4235 ind[2] = x + 1 + (y+1)*(nx+1);
4236 ind[3] = x + (y+1)*(nx+1);
4242 for (
int i = 0; i < nx; i++)
4248 for (
int j = 0; j < ny; j++)
4274 const int NVert = (nx+1) * (ny+1) * (nz+1);
4275 const int NElem = nx * ny * nz * 24;
4276 const int NBdrElem = 2*(nx*ny+nx*nz+ny*nz)*4;
4278 InitMesh(3, 3, NVert, NElem, NBdrElem);
4283 for (
real_t z = 0; z <= nz; z++)
4285 coord[2] = ( z / nz) * sz;
4286 for (
real_t y = 0; y <= ny; y++)
4288 coord[1] = (y / ny) * sy;
4289 for (
real_t x = 0; x <= nx; x++)
4291 coord[0] = (x / nx) * sx;
4297 std::map<std::array<int, 4>,
int> hex_face_verts;
4298 auto VertexIndex = [nx, ny](
int xc,
int yc,
int zc)
4300 return xc + (yc + zc*(ny+1))*(nx+1);
4304 for (
int z = 0; z < nz; z++)
4306 for (
int y = 0; y < ny; y++)
4308 for (
int x = 0; x < nx; x++)
4311 ind[0] = VertexIndex(x , y , z );
4312 ind[1] = VertexIndex(x+1, y , z );
4313 ind[2] = VertexIndex(x+1, y+1, z );
4314 ind[3] = VertexIndex(x , y+1, z );
4315 ind[4] = VertexIndex(x , y , z+1);
4316 ind[5] = VertexIndex(x+1, y , z+1);
4317 ind[6] = VertexIndex(x+1, y+1, z+1);
4318 ind[7] = VertexIndex( x, y+1, z+1);
4326 hex_face_verts.clear();
4335 std::map<std::array<int, 3>,
int> tet_face_count;
4337 std::map<std::array<int, 3>,
int> face_count_map;
4342 return std::array<int, 3> {v[0], v[1], v[2]};
4351 for (
int j = 0; j < el_faces.
Size(); j++)
4354 auto t = get3array(vertidxs);
4355 auto it = tet_face_count.find(t);
4356 if (it == tet_face_count.end())
4358 tet_face_count.insert({t, 1});
4359 face_count_map.insert({t, el_faces[j]});
4368 for (
const auto &edge : tet_face_count)
4370 if (edge.second == 1)
4372 int facenum = (face_count_map.find(edge.first))->second;
4379 ofstream test_stream(
"debug.mesh");
4381 test_stream.close();
4390 bool generate_edges,
bool sfc_ordering)
4414 for (j = 0; j < ny+1; j++)
4416 cy = ((
real_t) j / ny) * sy;
4417 for (i = 0; i < nx+1; i++)
4419 cx = ((
real_t) i / nx) * sx;
4431 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
4433 for (k = 0; k < nx*ny; k++)
4437 ind[0] = i + j*(nx+1);
4438 ind[1] = i + 1 +j*(nx+1);
4439 ind[2] = i + 1 + (j+1)*(nx+1);
4440 ind[3] = i + (j+1)*(nx+1);
4447 for (j = 0; j < ny; j++)
4449 for (i = 0; i < nx; i++)
4451 ind[0] = i + j*(nx+1);
4452 ind[1] = i + 1 +j*(nx+1);
4453 ind[2] = i + 1 + (j+1)*(nx+1);
4454 ind[3] = i + (j+1)*(nx+1);
4463 for (i = 0; i < nx; i++)
4469 for (j = 0; j < ny; j++)
4491 for (j = 0; j < ny+1; j++)
4493 cy = ((
real_t) j / ny) * sy;
4494 for (i = 0; i < nx+1; i++)
4496 cx = ((
real_t) i / nx) * sx;
4505 for (j = 0; j < ny; j++)
4507 for (i = 0; i < nx; i++)
4509 ind[0] = i + j*(nx+1);
4510 ind[1] = i + 1 + (j+1)*(nx+1);
4511 ind[2] = i + (j+1)*(nx+1);
4514 ind[1] = i + 1 + j*(nx+1);
4515 ind[2] = i + 1 + (j+1)*(nx+1);
4523 for (i = 0; i < nx; i++)
4529 for (j = 0; j < ny; j++)
4539 MFEM_ABORT(
"Unsupported element type.");
4545 if (generate_edges == 1)
4583 for (j = 0; j < n+1; j++)
4589 for (j = 0; j < n; j++)
4616 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4668 for (
int i = 0; i <
faces.Size(); i++)
4711 if (
dynamic_cast<const ParMesh*
>(&mesh))
4723 if (mesh.
Nodes && copy_nodes)
4759 int refine,
bool fix_orientation)
4763 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
4764 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
4781 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
4791 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
4827 ref_factors = ref_factor;
4840Mesh::Mesh(
const std::string &filename,
int generate_edges,
int refine,
4841 bool fix_orientation)
4842 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4851 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
4855 Load(imesh, generate_edges, refine, fix_orientation);
4860 bool fix_orientation)
4861 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4864 Load(input, generate_edges, refine, fix_orientation);
4873 "Not enough vertices in external array : "
4874 "len_vertex_data = "<< len_vertex_data <<
", "
4879 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
4884 memcpy(vertex_data,
vertices.GetData(),
4893 int *element_attributes,
int num_elements,
4895 int *boundary_attributes,
int num_boundary_elements,
4897 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4899 if (space_dimension == -1)
4905 num_boundary_elements);
4908 int boundary_index_stride = num_boundary_elements > 0 ?
4912 vertices.MakeRef(
reinterpret_cast<Vertex*
>(vertices_), num_vertices);
4915 for (
int i = 0; i < num_elements; i++)
4918 elements[i]->SetVertices(element_indices + i * element_index_stride);
4919 elements[i]->SetAttribute(element_attributes[i]);
4923 for (
int i = 0; i < num_boundary_elements; i++)
4926 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
4927 boundary[i]->SetAttribute(boundary_attributes[i]);
4935 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4973 MFEM_ABORT(
"NURBS mesh has no patches.");
4987#ifdef MFEM_USE_MEMALLOC
4996 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
5009 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
5012 for (
int i = 0; i < nv; i++)
5025 for (
int j = 0; j < nv; j++)
5097 MFEM_ABORT(
"invalid element type: " << type);
5104 std::string parse_tag)
5106 int curved = 0, read_gf = 1;
5107 bool finalize_topo =
true;
5111 MFEM_ABORT(
"Input stream is not open");
5118 getline(input, mesh_type);
5122 int mfem_version = 0;
5123 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
5124 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
5125 else if (mesh_type ==
"MFEM mesh v1.3") { mfem_version = 13; }
5129 int mfem_nc_version = 0;
5130 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
5131 else if (mesh_type ==
"MFEM NC mesh v1.1") { mfem_nc_version = 11; }
5132 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
5140 if (mfem_version >= 12 && parse_tag.empty())
5142 parse_tag =
"mfem_mesh_end";
5146 else if (mfem_nc_version)
5148 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
5155 MFEM_VERIFY(mfem_nc_version >= 10,
5156 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
5157 "used to load a parallel nonconforming mesh, sorry.");
5160 input, mfem_nc_version, curved, is_nc);
5165 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
5178 else if (mesh_type ==
"linemesh")
5182 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
5184 if (mesh_type ==
"curved_areamesh2")
5190 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
5194 else if (mesh_type ==
"TrueGrid")
5198 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
5200 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
5202 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
5203 "Unsupported VTK format");
5204 ReadVTKMesh(input, curved, read_gf, finalize_topo);
5206 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
5210 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
5214 else if (mesh_type ==
"MFEM NURBS NC-patch mesh v1.0")
5218 else if (mesh_type ==
"MFEM NURBS mesh v1.1")
5222 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
5227 else if (mesh_type ==
"$MeshFormat")
5230 finalize_topo =
false;
5231 curved =
Nodes !=
nullptr;
5235 ((mesh_type.size() > 2 &&
5236 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
5237 (mesh_type.size() > 3 &&
5238 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
5243#ifdef MFEM_USE_NETCDF
5246 MFEM_ABORT(
"NetCDF support requires configuration with"
5247 " MFEM_USE_NETCDF=YES");
5253 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
5254 " Use mfem::named_ifgzstream for input.");
5260 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
5286 bool generate_bdr =
false;
5291 if (curved && read_gf)
5305 if (mfem_version >= 12)
5311 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
5312 getline(input, line);
5318 if (line ==
"mfem_mesh_end") {
break; }
5320 while (line != parse_tag);
5322 else if (mfem_nc_version >= 10)
5327 MFEM_VERIFY(ident ==
"mfem_mesh_end",
5328 "invalid mesh: end of file tag not found");
5336 if (input.peek() ==
'p')
5339 MFEM_VERIFY(ident ==
"patch_cp",
"Invalid mesh format");
5348 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
5350 int i, j, ie, ib, iv, *v, nv;
5381 for (i = 0; i < num_pieces; i++)
5388 for (i = 0; i < num_pieces; i++)
5394 for (j = 0; j < m->
GetNE(); j++)
5399 for (j = 0; j < m->
GetNBE(); j++)
5404 for (
int k = 0; k < nv; k++)
5406 v[k] = lvert_vert[v[k]];
5411 for (j = 0; j < m->
GetNV(); j++)
5423 for (i = 0; i < num_pieces; i++)
5434 for (i = 0; i < num_pieces; i++)
5438 for (j = 0; j < m->
GetNE(); j++)
5443 for (
int k = 0; k < nv; k++)
5450 for (j = 0; j < m->
GetNBE(); j++)
5455 for (
int k = 0; k < nv; k++)
5462 for (j = 0; j < m->
GetNV(); j++)
5476 for (i = 0; i < num_pieces; i++)
5478 gf_array[i] = mesh_array[i]->
GetNodes();
5491 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
5494 ref_factors = ref_factor;
5505 int orig_ne = orig_mesh.
GetNE();
5506 MFEM_VERIFY(ref_factors.
Size() == orig_ne,
5507 "Number of refinement factors must equal number of elements")
5508 MFEM_VERIFY(orig_ne == 0 ||
5509 ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
5512 "Invalid refinement type. Must use closed basis type.");
5514 int min_ref = orig_ne > 0 ? ref_factors.
Min() : 1;
5515 int max_ref = orig_ne > 0 ? ref_factors.
Max() : 1;
5517 bool var_order = (min_ref != max_ref);
5530 for (
int i = 0; i < orig_ne; i++)
5546 for (
int el = 0; el < orig_ne; el++)
5558 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
5559 for (
int i = 0; i < phys_pts.
Width(); i++)
5568 for (
int k = 0; k < nvert; k++)
5571 v[k] = rdofs[c2h_map[cid]];
5584 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
5595 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
5601 for (
int k = 0; k < nvert; k++)
5604 v[k] = rdofs[c2h_map[cid]];
5626 for (
int iel = 0; iel < orig_ne; iel++)
5635 const int *node_map = NULL;
5638 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
5639 const int *vertex_map = vertex_fec.
GetDofMap(geom);
5640 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
5641 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
5644 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
5647 int iv = vertex_map[iv_lex];
5649 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
5651 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
5654 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
5665 using GeomRef = std::pair<Geometry::Type, int>;
5666 std::map<GeomRef, int> point_matrices_offsets;
5668 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5672 GeomRef id(geom, ref_factors[el_coarse]);
5673 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
5679 point_matrices_offsets[id] = n_point_matrices[geom];
5680 n_point_matrices[geom] += nref_el;
5687 int nmatrices = n_point_matrices[geom];
5694 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5697 int ref = ref_factors[el_coarse];
5698 int offset = point_matrices_offsets[GeomRef(geom, ref)];
5704 for (
int k = 0; k < nvert; k++)
5731 if (orig_mesh.
GetNodes() !=
nullptr)
5741 "Mesh::MakeSimplicial requires a properly oriented input mesh");
5743 "Mesh::MakeSimplicial does not support non-conforming meshes.")
5750 Mesh copy(orig_mesh);
5753 std::iota(parent_elements.
begin(), parent_elements.
end(), 0);
5754 return parent_elements;
5757 int nv = orig_mesh.
GetNV();
5758 int ne = orig_mesh.
GetNE();
5759 int nbe = orig_mesh.
GetNBE();
5772 int new_ne = 0, new_nbe = 0;
5773 for (
int i=0; i<ne; ++i)
5777 for (
int i=0; i<nbe; ++i)
5786 for (
int i=0; i<nv; ++i)
5796 if (vglobal ==
nullptr)
5799 std::iota(vglobal_id.
begin(), vglobal_id.
end(), 0);
5800 vglobal = vglobal_id.
GetData();
5804 constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
5805 constexpr int quad_ntris = 2;
5806 constexpr int prism_ntets = 3;
5810 static const int quad_trimap[2][nv_tri*quad_ntris] =
5822 static const int prism_rot[nv_prism*nv_prism] =
5831 static const int prism_f[nv_quad] = {1, 2, 5, 4};
5832 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
5846 static const int hex_rot[nv_hex*nv_hex] =
5848 0, 1, 2, 3, 4, 5, 6, 7,
5849 1, 0, 4, 5, 2, 3, 7, 6,
5850 2, 1, 5, 6, 3, 0, 4, 7,
5851 3, 0, 1, 2, 7, 4, 5, 6,
5852 4, 0, 3, 7, 5, 1, 2, 6,
5853 5, 1, 0, 4, 6, 2, 3, 7,
5854 6, 2, 1, 5, 7, 3, 0, 4,
5855 7, 3, 2, 6, 4, 0, 1, 5
5857 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
5858 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
5859 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
5860 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
5861 static const int hex_tetmap0[nv_tet*5] =
5868 static const int hex_tetmap1[nv_tet*6] =
5875 static const int hex_tetmap2[nv_tet*6] =
5882 static const int hex_tetmap3[nv_tet*6] =
5889 static const int *hex_tetmaps[4] =
5891 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
5894 auto find_min = [](
const int *
a,
int n) {
return std::min_element(
a,
a+n)-
a; };
5897 for (
int i=0; i<ne; ++i)
5899 const int *v = orig_mesh.
elements[i]->GetVertices();
5903 if (num_subdivisions[orig_geom] == 1)
5916 for (
int itri=0; itri<quad_ntris; ++itri)
5921 for (
int iv=0; iv<nv_tri; ++iv)
5923 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
5932 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
5935 int irot = find_min(vg, nv_prism);
5936 for (
int iv=0; iv<nv_prism; ++iv)
5938 int jv = prism_rot[iv + irot*nv_prism];
5943 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
5944 int j = find_min(q, nv_quad);
5945 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
5946 for (
int itet=0; itet<prism_ntets; ++itet)
5951 for (
int iv=0; iv<nv_tet; ++iv)
5953 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
5962 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
5966 int irot = find_min(vg, nv_hex);
5967 for (
int iv=0; iv<nv_hex; ++iv)
5969 int jv = hex_rot[iv + irot*nv_hex];
5979 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
5980 j = find_min(q, nv_quad);
5981 if (j == 0 || j == 2) { bitmask += 4; }
5983 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
5984 j = find_min(q, nv_quad);
5985 if (j == 1 || j == 3) { bitmask += 2; }
5987 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
5988 j = find_min(q, nv_quad);
5989 if (j == 0 || j == 2) { bitmask += 1; }
5992 int nrot = num_rot[bitmask];
5993 for (
int k=0; k<nrot; ++k)
6007 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
6008 int ntets = (ndiags == 0) ? 5 : 6;
6009 const int *tetmap = hex_tetmaps[ndiags];
6010 for (
int itet=0; itet<ntets; ++itet)
6015 for (
int iv=0; iv<nv_tet; ++iv)
6017 v2[iv] = vg[tetmap[itet + iv*ntets]];
6027 for (
int i=0; i<nbe; ++i)
6029 const int *v = orig_mesh.
boundary[i]->GetVertices();
6032 if (num_subdivisions[orig_geom] == 1)
6042 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
6044 int iv_min = find_min(vg, nv_quad);
6045 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
6046 for (
int itri=0; itri<quad_ntris; ++itri)
6051 for (
int iv=0; iv<nv_tri; ++iv)
6053 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
6060 MFEM_ABORT(
"Unreachable");
6071 return parent_elems;
6081 auto *orig_fespace = orig_mesh.
GetNodes()->FESpace();
6082 SetCurvature(orig_fespace->GetMaxElementOrder(), orig_fespace->IsDGSpace(),
6107 child_nodes_in_parent;
6108 for (
int i = 0; i < parent_elements.
Size(); i++)
6110 const int ip = parent_elements[i];
6112 orig_mesh.
GetNodes()->GetElementDofValues(ip, edofvals);
6134 for (
auto cv : child_vertices)
6135 for (
int ipv = 0; ipv < parent_vertices.
Size(); ipv++)
6136 if (cv == parent_vertices[ipv])
6147 child_nodes_in_parent.
SetSize(0);
6148 const auto *orig_FE = orig_mesh.
GetNodes()->FESpace()->GetFE(ip);
6149 for (
auto pn : node_map)
6151 child_nodes_in_parent.
Append(orig_FE->GetNodes()[pn]);
6154 shape.
SetSize(orig_FE->GetDof(),
6155 simplex_FE->GetDof());
6157 for (
int j = 0; j < simplex_FE->GetNodes().Size(); j++)
6159 const auto &simplex_node = simplex_FE->GetNodes()[j];
6162 simplex_node_in_orig.
Set3(
6163 child_nodes_in_parent[0].x +
6164 simplex_node.x * (child_nodes_in_parent[1].x - child_nodes_in_parent[0].x)
6165 + simplex_node.y * (child_nodes_in_parent[2].x - child_nodes_in_parent[0].x)
6166 + simplex_node.z * (child_nodes_in_parent[(
Dim > 2) ? 3 : 0].x -
6167 child_nodes_in_parent[0].x),
6168 child_nodes_in_parent[0].y +
6169 simplex_node.x * (child_nodes_in_parent[1].y - child_nodes_in_parent[0].y)
6170 + simplex_node.y * (child_nodes_in_parent[2].y - child_nodes_in_parent[0].y)
6171 + simplex_node.z * (child_nodes_in_parent[(
Dim > 2) ? 3 : 0].y -
6172 child_nodes_in_parent[0].y),
6173 child_nodes_in_parent[0].z +
6174 simplex_node.x * (child_nodes_in_parent[1].z - child_nodes_in_parent[0].z)
6175 + simplex_node.y * (child_nodes_in_parent[2].z - child_nodes_in_parent[0].z)
6176 + simplex_node.z * (child_nodes_in_parent[(
Dim > 2) ? 3 : 0].z -
6177 child_nodes_in_parent[0].z));
6179 orig_FE->CalcShape(simplex_node_in_orig, col);
6184 orig_mesh.
GetNodes()->GetElementDofValues(ip, edofvals);
6190 point_matrix.
SetSize(simplex_FE->GetDof(), sdim);
6191 MultAtB(shape, edofvals_mat, point_matrix);
6199 MFEM_ABORT(
"Internal Error!");
6207 Mesh periodic_mesh(orig_mesh,
true);
6213 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
6218 for (
int j = 0; j < nv; j++)
6224 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
6229 for (
int j = 0; j < nv; j++)
6236 return periodic_mesh;
6240 const std::vector<Vector> &translations,
real_t tol)
const
6244 Vector coord(sdim), at(sdim), dx(sdim);
6245 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
6246 xMax = xMin = xDiff = 0.0;
6249 unordered_set<int> bdr_v;
6250 for (
int be = 0; be <
GetNBE(); be++)
6255 for (
int i = 0; i < dofs.
Size(); i++)
6257 bdr_v.insert(dofs[i]);
6260 for (
int j = 0; j < sdim; j++)
6262 xMax[j] = max(xMax[j], coord[j]);
6263 xMin[j] = min(xMin[j], coord[j]);
6267 add(xMax, -1.0, xMin, xDiff);
6276 unordered_map<int, int> replica2primary;
6278 unordered_map<int, unordered_set<int>> primary2replicas;
6281 std::unique_ptr<KDTreeBase<int,real_t>> kdtree;
6282 if (sdim == 1) { kdtree.reset(
new KDTree1D); }
6283 else if (sdim == 2) { kdtree.reset(
new KDTree2D); }
6284 else if (sdim == 3) { kdtree.reset(
new KDTree3D); }
6285 else { MFEM_ABORT(
"Invalid space dimension."); }
6289 for (
const int v : bdr_v)
6291 primary2replicas[v];
6299 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
6301 if (r ==
p) {
return; }
6302 primary2replicas[
p].insert(r);
6303 replica2primary[r] =
p;
6304 for (
const int s : primary2replicas[r])
6306 primary2replicas[
p].insert(s);
6307 replica2primary[s] =
p;
6309 primary2replicas.erase(r);
6312 for (
unsigned int i = 0; i < translations.size(); i++)
6314 for (
int vi : bdr_v)
6317 add(coord, translations[i], at);
6319 const int vj = kdtree->FindClosestPoint(at.GetData());
6321 add(at, -1.0, coord, dx);
6323 if (dx.
Norml2() > dia*tol) {
continue; }
6328 const bool pi = primary2replicas.find(vi) != primary2replicas.end();
6329 const bool pj = primary2replicas.find(vj) != primary2replicas.end();
6335 make_replica(vj, vi);
6340 const int owner_of_vj = replica2primary[vj];
6342 make_replica(vi, owner_of_vj);
6348 const int owner_of_vi = replica2primary[vi];
6349 make_replica(vj, owner_of_vi);
6356 const int owner_of_vi = replica2primary[vi];
6357 const int owner_of_vj = replica2primary[vj];
6358 make_replica(owner_of_vj, owner_of_vi);
6363 std::vector<int> v2v(
GetNV());
6364 for (
size_t i = 0; i < v2v.size(); i++)
6366 v2v[i] =
static_cast<int>(i);
6368 for (
const auto &r2p : replica2primary)
6370 v2v[r2p.first] = r2p.second;
6377 MFEM_VERIFY(
NURBSext,
"Mesh::RefineNURBSFromFile: Not a NURBS mesh!");
6378 mfem::out<<
"Refining NURBS from refinement file: "<<ref_file<<endl;
6381 ifstream input(ref_file);
6388 mfem::out<<
"Knot vectors in ref_file: "<<nkv<<endl;
6390 MFEM_ABORT(
"Refine file does not have the correct number of knot vectors");
6395 for (
int kv = 0; kv < nkv; kv++)
6397 knotVec[kv] =
new Vector();
6398 knotVec[kv]->
Load(input);
6406 for (
int kv = 0; kv < nkv; kv++)
6416 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
6421 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
6438 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
6443 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
6460 mfem_error(
"Mesh::KnotRemove : Not a NURBS mesh!");
6465 mfem_error(
"Mesh::KnotRemove : KnotVector array size mismatch!");
6493 "Refinement factors must be defined for each dimension");
6499 const std::string &kvf)
6501 MFEM_VERIFY(
NURBSext,
"This type of refinement is only for NURBS meshes");
6510 cf1 = (cf1 &&
f == 1);
6520 MFEM_VERIFY(!usingKVF,
"This refinement type is not supported for this"
6521 " NURBS mesh type");
6529 for (
int i=0; i<cf.
Size(); ++i) { cf[i] *= rf[i]; }
6555 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
6579 for (
int i = 0; i <
elements.Size(); i++)
6589 for (
int i = 0; i <
boundary.Size(); i++)
6608 for (
int i = 0; i < vd; i++)
6675 input >> edge_to_ukv[j] >> v[0] >> v[1];
6699 if (edge_to_ukv.
Size() == 0)
6713 const int NPKV = NP *
dim;
6714 constexpr int notset = -9999999;
6716 auto edge_to_dim = [](
int i) {
return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };
6726 for (
int i = 0; i < NP; i++)
6741 edge_to_pkv = notset;
6747 for (
int i = 0; i < NPKV; i++)
6751 std::function<int(
int)> get_root;
6752 get_root = [&pkv_map, &get_root](
int i) ->
int
6754 return (pkv_map[i] == i) ? i : get_root(pkv_map[i]);
6756 auto unite = [&pkv_map, &get_root](
int i,
int j)
6758 const int ri = get_root(i);
6759 const int rj = get_root(j);
6760 if (ri == rj) {
return; }
6762 (ri < rj) ? pkv_map[rj] = ri : pkv_map[ri] = rj;
6766 for (
int p = 0;
p < NP;
p++)
6771 for (
int i = 0; i < edges.
Size(); i++)
6773 const int edge = edges[i];
6774 const int d = edge_to_dim(i);
6775 const int pkv =
p*
dim+d;
6778 if (edge_to_pkv[edge] != notset)
6780 const int pkv_other =
UnsignIndex(edge_to_pkv[edge]);
6781 unite(pkv, pkv_other);
6787 edge_to_pkv[edge] = (v[1] > v[0]) ? pkv :
FlipIndexSign(pkv);
6795 for (
int i = 0; i < NPKV; i++)
6797 pkv_to_rpkv[i] = get_root(pkv_map[i]);
6798 ukv_to_rpkv[i] = pkv_to_rpkv[i];
6804 std::map<int, int> rpkv_to_ukv;
6805 for (
int i = 0; i < ukv_to_rpkv.
Size(); i++)
6807 rpkv_to_ukv[ukv_to_rpkv[i]] = i;
6815 const int rpkv = pkv_to_rpkv[pkv];
6816 const int ukv = rpkv_to_ukv[rpkv];
6817 edge_to_ukv[i] = (edge_to_pkv[i] < 0) ?
FlipIndexSign(ukv) : ukv;
6826 if (
dim == 1) {
return; }
6833 auto faceNeighbors = [&](
int p,
int kv, std::unordered_set<int> &nghb)
6838 for (
auto face : pfaces)
6845 const int skv = edge_to_ukv[e];
6846 if (skv == kv ||
FlipIndexSign(skv) == kv) { hasKV =
true; }
6851 face2elem->
GetRow(face, row);
6852 for (
auto elem : row) { nghb.insert(elem); }
6857 std::vector<std::vector<int>> dir_edges;
6878 bool initKV =
false;
6880 auto setPatchDirections = [&](
int p,
int kv,
Array<bool> &edgeSet,
6881 std::unordered_set<int> &visited)
6887 for (
int i = 0; i < pe.
Size(); i++)
6889 ukvs[i] = edge_to_ukv[pe[i]];
6895 for (
int d=0; d<
dim; ++d)
6897 const int skv = edge_to_ukv[pe[dir_edges[d][0]]];
6900 for (
auto e : dir_edges[d])
6901 if (!edgeSet[pe[e]])
6914 int ref_edge0 = dir_edges[thisDir][0];
6915 for (
auto ref_edge : dir_edges[thisDir])
6917 const int edge = pe[ref_edge];
6920 ref_edge0 = ref_edge;
6924 if (initKV && !edgeSet[pe[ref_edge0]])
6933 edgeSet[pe[ref_edge0]] =
true;
6934 for (
auto i : dir_edges[thisDir])
6941 const int edge = pe[i];
6946 MFEM_ASSERT(!edgeSet[edge],
"");
6950 edgeSet[edge] =
true;
6959 std::unordered_set<int> unset;
6963 for (
int iter=0; iter<max_iter; ++iter)
6974 std::list<int> nextPatches;
6975 std::unordered_set<int> nextSet;
6976 std::unordered_set<int> visited;
6978 if (unset.size() == 0)
6983 const int p0 = *unset.begin();
6984 nextPatches.push_back(p0);
6990 for (
int d=0; d<
dim; ++d)
6992 for (
auto e : dir_edges[d])
6993 if (!edgeSet[pe[e]])
7005 const int kv =
UnsignIndex(edge_to_ukv[pe[dir_edges[unsetDim][0]]]);
7009 while (nextPatches.size() > 0)
7011 const int p = nextPatches.front();
7012 nextPatches.pop_front();
7016 const bool somethingSet = setPatchDirections(
p, kv, edgeSet, visited);
7023 std::unordered_set<int> neighbors;
7024 faceNeighbors(
p, kv, neighbors);
7028 for (
auto edge : pe)
7041 for (
auto n : neighbors)
7043 if (n !=
p && visited.count(n) == 0 && unset.count(n) > 0)
7045 if (nextSet.count(n) == 0)
7047 nextPatches.push_back(n);
7057 for (
auto eset : edgeSet)
7064 MFEM_ASSERT(allSet && unset.size() == 0,
"Some edge is not set");
7086 int inputNumOfEdges = -1;
7089 input >> inputNumOfEdges;
7091 MFEM_VERIFY(
NumOfEdges == inputNumOfEdges,
"");
7098 input >> ukv >> v[0] >> v[1];
7100 for (
int i=0; i<2; ++i)
7109 edge_to_ukv[j] = ukv;
7118 if (
p.Size() >= v.
Size())
7120 for (
int d = 0; d < v.
Size(); d++)
7128 for (d = 0; d <
p.Size(); d++)
7132 for ( ; d < v.
Size(); d++)
7145 nodes.ProjectCoefficient(xyz);
7175 const bool warn =
true;
7178 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
7182 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
7183 " H1-continuous mesh!\n "
7184 "If this is the desired behavior, you can silence"
7185 " this warning by converting\n "
7186 "the NURBS mesh to high-order mesh in advance by"
7187 " calling the method\n "
7188 "Mesh::SetCurvature().");
7220 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
7235 const int old_space_dim =
spaceDim;
7248 MFEM_ASSERT(
nodes != NULL,
"");
7252 nodes->GetNodalValues(vert_val, i+1);
7265 for (
int e = 0; e <
GetNE(); e++)
7278 detvals(q) = Jac.
Weight();
7291 int det_order =
Dim*mesh_poly_deg-1;
7296 auto detgf = std::make_unique<GridFunction>(fespace_det);
7297 detgf->MakeOwner(fec_det);
7306 case 1:
return GetNV();
7342#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
7343static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
7348 int i, j, k, wo = 0, fo = 0;
7357 int *vi =
elements[i]->GetVertices();
7360 for (j = 0; j < 3; j++)
7364 for (j = 0; j < 2; j++)
7365 for (k = 0; k < 2; k++)
7367 J(j, k) = v[j+1][k] - v[0][k];
7388 MFEM_ABORT(
"Invalid 2D element type \""
7405 int *vi =
elements[i]->GetVertices();
7411 for (j = 0; j < 4; j++)
7415 for (j = 0; j < 3; j++)
7416 for (k = 0; k < 3; k++)
7418 J(j, k) = v[j+1][k] - v[0][k];
7478 MFEM_ABORT(
"Invalid 3D element type \""
7484#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
7487 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
7488 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
7492 MFEM_CONTRACT_VAR(fo);
7505 if (test[0] == base[0])
7506 if (test[1] == base[1])
7514 else if (test[0] == base[1])
7515 if (test[1] == base[0])
7524 if (test[1] == base[0])
7535 for (
int j = 0; j < 3; j++)
7536 if (test[aor[j]] != base[j])
7538 mfem::err <<
"Mesh::GetTriOrientation(...)" << endl;
7540 for (
int k = 0; k < 3; k++)
7545 for (
int k = 0; k < 3; k++)
7566 const int oo[6][6] =
7576 int ori_a_c = oo[ori_a_b][ori_b_c];
7582 const int inv_ori[6] = {0, 1, 4, 3, 2, 5};
7583 return inv_ori[ori];
7590 for (i = 0; i < 4; i++)
7591 if (test[i] == base[0])
7598 if (test[(i+1)%4] == base[1])
7607 for (
int j = 0; j < 4; j++)
7608 if (test[aor[j]] != base[j])
7610 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
7612 for (
int k = 0; k < 4; k++)
7617 for (
int k = 0; k < 4; k++)
7626 if (test[(i+1)%4] == base[1])
7643 const int oo[8][8] =
7645 {0, 1, 2, 3, 4, 5, 6, 7},
7646 {1, 0, 3, 2, 5, 4, 7, 6},
7647 {2, 7, 4, 1, 6, 3, 0, 5},
7648 {3, 6, 5, 0, 7, 2, 1, 4},
7649 {4, 5, 6, 7, 0, 1, 2, 3},
7650 {5, 4, 7, 6, 1, 0, 3, 2},
7651 {6, 3, 0, 5, 2, 7, 4, 1},
7652 {7, 2, 1, 4, 3, 6, 5, 0}
7655 int ori_a_c = oo[ori_a_b][ori_b_c];
7661 const int inv_ori[8] = {0, 1, 6, 3, 4, 5, 2, 7};
7662 return inv_ori[ori];
7673 if (test[0] == base[0])
7674 if (test[1] == base[1])
7675 if (test[2] == base[2])
7683 else if (test[2] == base[1])
7684 if (test[3] == base[2])
7693 if (test[1] == base[2])
7701 else if (test[1] == base[0])
7702 if (test[2] == base[1])
7703 if (test[0] == base[2])
7711 else if (test[3] == base[1])
7712 if (test[2] == base[2])
7721 if (test[3] == base[2])
7729 else if (test[2] == base[0])
7730 if (test[3] == base[1])
7731 if (test[0] == base[2])
7739 else if (test[0] == base[1])
7740 if (test[1] == base[2])
7749 if (test[3] == base[2])
7758 if (test[0] == base[1])
7759 if (test[2] == base[2])
7767 else if (test[1] == base[1])
7768 if (test[0] == base[2])
7777 if (test[1] == base[2])
7788 for (
int j = 0; j < 4; j++)
7789 if (test[aor[j]] != base[j])
7814 int *bv =
boundary[i]->GetVertices();
7834 if (
faces_info[fi].Elem2No >= 0) {
continue; }
7837 int *bv =
boundary[i]->GetVertices();
7839 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
7840 const int *fv =
faces[fi]->GetVertices();
7857 MFEM_ABORT(
"Invalid 2D boundary element type \""
7858 << bdr_type <<
"\"");
7863 if (orientation % 2 == 0) {
continue; }
7865 if (!fix_it) {
continue; }
7901 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
7919 MFEM_ASSERT(o >= 0 && o < 2,
"Invalid orientation for Geometry::SEGMENT!");
7931 MFEM_ASSERT(o >= 0 && o < 6,
"Invalid orientation for Geometry::TRIANGLE!");
7944 fip.
x = 1.0 - ip.
x - ip.
y;
7949 fip.
x = 1.0 - ip.
x - ip.
y;
7955 fip.
y = 1.0 - ip.
x - ip.
y;
7960 fip.
y = 1.0 - ip.
x - ip.
y;
7965 MFEM_ASSERT(o >= 0 && o < 8,
"Invalid orientation for Geometry::SQUARE!");
8009 MFEM_ABORT(
"Unsupported face geometry for TransformBdrElementToFace!");
8016 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
8027 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
8052 const int *v =
elements[i]->GetVertices();
8053 cor[0] = (v[0] < v[1]) ? (1) : (-1);
8063 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
8064 "is not generated.");
8067 const int *v =
elements[i]->GetVertices();
8068 const int ne =
elements[i]->GetNEdges();
8070 for (
int j = 0; j < ne; j++)
8072 const int *e =
elements[i]->GetEdgeVertices(j);
8073 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
8084 const int *v =
boundary[i]->GetVertices();
8085 cor[0] = (v[0] < v[1]) ? (1) : (-1);
8098 const int *v =
boundary[i]->GetVertices();
8099 const int ne =
boundary[i]->GetNEdges();
8101 for (
int j = 0; j < ne; j++)
8103 const int *e =
boundary[i]->GetEdgeVertices(j);
8104 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
8116 const int *v =
faces[i]->GetVertices();
8117 o[0] = (v[0] < v[1]) ? (1) : (-1);
8129 const int *v =
faces[i]->GetVertices();
8130 const int ne =
faces[i]->GetNEdges();
8132 for (
int j = 0; j < ne; j++)
8134 const int *e =
faces[i]->GetEdgeVertices(j);
8135 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
8163 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
8226 const int nv =
elements[i]->GetNVertices();
8227 const int *v =
elements[i]->GetVertices();
8228 for (
int j = 0; j < nv; j++)
8238 const int nv =
elements[i]->GetNVertices();
8239 const int *v =
elements[i]->GetVertices();
8240 for (
int j = 0; j < nv; j++)
8259 const int nv =
boundary[i]->GetNVertices();
8260 const int *v =
boundary[i]->GetVertices();
8261 for (
int j = 0; j < nv; j++)
8267 vert_bdr_elem->
MakeJ();
8271 const int nv =
boundary[i]->GetNVertices();
8272 const int *v =
boundary[i]->GetVertices();
8273 for (
int j = 0; j < nv; j++)
8281 return vert_bdr_elem;
8320 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
8324 int n = el_faces.
Size();
8327 for (
int j = 0; j < n; j++)
8331 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
8335 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
8336 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
8353 for (
auto f : elem_faces)
8374 const int *bv =
boundary[i]->GetVertices();
8384 default: MFEM_ABORT(
"invalid geometry");
8393 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
8396 const int *bv =
boundary[bdr_el]->GetVertices();
8404 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
8411 int bdr_el,
int &el,
int &info)
const
8416 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
8419 const int *bv =
boundary[bdr_el]->GetVertices();
8427 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
8466 for (j = 0; j < nv; j++)
8468 pointmat(k, j) =
vertices[v[j]](k);
8483 for (j = 0; j < nv; j++)
8485 pointmat(k, j) =
vertices[v[j]](k);
8497 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
8500 return sqrt(length);
8508 for (
int i = 0; i < elem_array.
Size(); i++)
8513 for (
int i = 0; i < elem_array.
Size(); i++)
8515 const int *v = elem_array[i]->GetVertices();
8516 const int ne = elem_array[i]->GetNEdges();
8517 for (
int j = 0; j < ne; j++)
8519 const int *e = elem_array[i]->GetEdgeVertices(j);
8533 v_to_v.
Push(v[0], v[1]);
8540 const int *v =
elements[i]->GetVertices();
8541 const int ne =
elements[i]->GetNEdges();
8542 for (
int j = 0; j < ne; j++)
8544 const int *e =
elements[i]->GetEdgeVertices(j);
8545 v_to_v.
Push(v[e[0]], v[e[1]]);
8553 int i, NumberOfEdges;
8569 const int *v =
boundary[i]->GetVertices();
8583 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
8587 return NumberOfEdges;
8646 if (
faces[gf] == NULL)
8678 if (
faces[gf] == NULL)
8688 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
8689 "Interior edge found between 2D elements "
8691 <<
" and " << el <<
".");
8692 int *v =
faces[gf]->GetVertices();
8694 if (v[1] == v0 && v[0] == v1)
8698 else if (v[0] == v0 && v[1] == v1)
8708 MFEM_ABORT(
"internal error");
8714 int v0,
int v1,
int v2)
8716 if (
faces[gf] == NULL)
8726 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
8727 "Interior triangular face found connecting elements "
8729 <<
" and " << el <<
".");
8730 int orientation, vv[3] = { v0, v1, v2 };
8737 faces_info[gf].Elem2Inf = 64 * lf + orientation;
8742 int v0,
int v1,
int v2,
int v3)
8754 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
8755 "Interior quadrilateral face found connecting elements "
8757 <<
" and " << el <<
".");
8758 int vv[4] = { v0, v1, v2, v3 };
8783 faces.SetSize(nfaces);
8785 for (
int i = 0; i < nfaces; ++i)
8804 const int ne =
elements[i]->GetNEdges();
8805 for (
int j = 0; j < ne; j++)
8807 const int *e =
elements[i]->GetEdgeVertices(j);
8818 for (
int j = 0; j < 4; j++)
8822 v[fv[0]], v[fv[1]], v[fv[2]]);
8828 for (
int j = 0; j < 2; j++)
8832 v[fv[0]], v[fv[1]], v[fv[2]]);
8834 for (
int j = 2; j < 5; j++)
8838 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8844 for (
int j = 0; j < 1; j++)
8848 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8850 for (
int j = 1; j < 5; j++)
8854 v[fv[0]], v[fv[1]], v[fv[2]]);
8860 for (
int j = 0; j < 6; j++)
8864 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8869 MFEM_ABORT(
"Unexpected type of Element.");
8877 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
8888 nc_faces_info.Reserve(list.masters.Size() + list.slaves.Size());
8895 if (master.index >= nfaces) {
continue; }
8901 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
8902 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
8908 if (slave.index < 0 ||
8909 slave.index >= nfaces ||
8910 slave.master >= nfaces)
8929 list.point_matrices[slave.geom][slave.matrix]));
8938 const int *v =
elements[i]->GetVertices();
8943 for (
int j = 0; j < 4; j++)
8946 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8952 for (
int j = 0; j < 1; j++)
8955 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8957 for (
int j = 1; j < 5; j++)
8960 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8966 for (
int j = 0; j < 2; j++)
8969 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8971 for (
int j = 2; j < 5; j++)
8974 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8982 for (
int j = 0; j < 6; j++)
8985 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
9014 for (
int j = 0; j < 4; j++)
9018 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
9024 for (
int j = 0; j < 2; j++)
9028 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
9030 for (
int j = 2; j < 5; j++)
9034 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
9040 for (
int j = 0; j < 1; j++)
9044 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
9046 for (
int j = 1; j < 5; j++)
9050 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
9058 for (
int j = 0; j < 6; j++)
9062 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
9067 MFEM_ABORT(
"Unexpected type of Element.");
9081 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
9086 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
9090 MFEM_ABORT(
"Unexpected type of boundary Element.");
9104void Rotate3(
int &
a,
int &
b,
int &c)
9136 Table *old_elem_vert = NULL;
9147 int *v =
elements[i]->GetVertices();
9149 Rotate3(v[0], v[1], v[2]);
9152 Rotate3(v[1], v[2], v[3]);
9165 int *v =
boundary[i]->GetVertices();
9167 Rotate3(v[0], v[1], v[2]);
9183 delete old_elem_vert;
9199 if (
p[i] < pmin[i]) { pmin[i] =
p[i]; }
9200 if (
p[i] > pmax[i]) { pmax[i] =
p[i]; }
9214 for (
int i =
spaceDim-1; i >= 0; i--)
9216 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
9217 if (idx < 0) { idx = 0; }
9218 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
9219 part = part * nxyz[i] + idx;
9221 partitioning[el] = part;
9224 return partitioning;
9234#ifdef MFEM_USE_METIS
9236 int print_messages = 1;
9239 int init_flag, fin_flag;
9240 MPI_Initialized(&init_flag);
9241 MPI_Finalized(&fin_flag);
9242 if (init_flag && !fin_flag)
9246 if (rank != 0) { print_messages = 0; }
9250 int i, *partitioning;
9260 partitioning[i] = 0;
9267 partitioning[i] = i;
9273#ifndef MFEM_USE_METIS_5
9285 bool freedata =
false;
9287 idx_t *mpartitioning;
9290 if (
sizeof(
idx_t) ==
sizeof(int))
9294 mpartitioning = (
idx_t*) partitioning;
9303 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
9304 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
9305 mpartitioning =
new idx_t[n];
9308#ifndef MFEM_USE_METIS_5
9311 METIS_SetDefaultOptions(options);
9312 options[METIS_OPTION_CONTIG] = 1;
9320 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
9325 if (part_method >= 0 && part_method <= 2)
9327 for (i = 0; i < n; i++)
9333 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
9339 if (part_method == 0 || part_method == 3)
9341#ifndef MFEM_USE_METIS_5
9370 " error in METIS_PartGraphRecursive!");
9377 if (part_method == 1 || part_method == 4)
9379#ifndef MFEM_USE_METIS_5
9408 " error in METIS_PartGraphKway!");
9415 if (part_method == 2 || part_method == 5)
9417#ifndef MFEM_USE_METIS_5
9430 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
9447 " error in METIS_PartGraphKway!");
9455 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
9459 nparts = (int) mparts;
9460 if (mpartitioning != (
idx_t*)partitioning)
9464 partitioning[k] = mpartitioning[k];
9471 delete[] mpartitioning;
9487 auto count_partition_elements = [&]()
9489 for (i = 0; i < nparts; i++)
9497 psize[partitioning[i]].one++;
9501 for (i = 0; i < nparts; i++)
9503 if (psize[i].one == 0) { empty_parts++; }
9507 count_partition_elements();
9515 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
9516 << empty_parts <<
" empty parts!"
9517 <<
" Applying a simple fix ..." << endl;
9522 for (i = nparts-1; i > nparts-1-empty_parts; i--)
9529 for (i = nparts-1; i > nparts-1-empty_parts; i--)
9531 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
9537 partitioning[j] = psize[nparts-1-i].two;
9544 count_partition_elements();
9548 return partitioning;
9552 mfem_error(
"Mesh::GeneratePartitioning(...): "
9553 "MFEM was compiled without Metis.");
9567 int num_elem, *i_elem_elem, *j_elem_elem;
9569 num_elem = elem_elem.
Size();
9570 i_elem_elem = elem_elem.
GetI();
9571 j_elem_elem = elem_elem.
GetJ();
9576 int stack_p, stack_top_p, elem;
9580 for (i = 0; i < num_elem; i++)
9582 if (partitioning[i] > num_part)
9584 num_part = partitioning[i];
9591 for (i = 0; i < num_part; i++)
9598 for (elem = 0; elem < num_elem; elem++)
9600 if (component[elem] >= 0)
9605 component[elem] = num_comp[partitioning[elem]]++;
9607 elem_stack[stack_top_p++] = elem;
9609 for ( ; stack_p < stack_top_p; stack_p++)
9611 i = elem_stack[stack_p];
9612 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
9615 if (partitioning[k] == partitioning[i])
9617 if (component[k] < 0)
9619 component[k] = component[i];
9620 elem_stack[stack_top_p++] = k;
9622 else if (component[k] != component[i])
9634 int i, n_empty, n_mcomp;
9642 n_empty = n_mcomp = 0;
9643 for (i = 0; i < num_comp.
Size(); i++)
9644 if (num_comp[i] == 0)
9648 else if (num_comp[i] > 1)
9655 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
9656 <<
"The following subdomains are empty :\n";
9657 for (i = 0; i < num_comp.
Size(); i++)
9658 if (num_comp[i] == 0)
9666 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
9667 <<
"The following subdomains are NOT connected :\n";
9668 for (i = 0; i < num_comp.
Size(); i++)
9669 if (num_comp[i] > 1)
9675 if (n_empty == 0 && n_mcomp == 0)
9676 mfem::out <<
"Mesh::CheckPartitioning(...) : "
9677 "All subdomains are connected." << endl;
9701 c(0) =
a[0]*
a[3]-
a[1]*
a[2];
9702 c(1) =
a[0]*
b[3]-
a[1]*
b[2]+
b[0]*
a[3]-
b[1]*
a[2];
9703 c(2) =
b[0]*
b[3]-
b[1]*
b[2];
9724 c(0) = (
a[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
9725 a[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
9726 a[2] * (
a[3] *
a[7] -
a[4] *
a[6]));
9728 c(1) = (
b[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
9729 b[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
9730 b[2] * (
a[3] *
a[7] -
a[4] *
a[6]) +
9732 a[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
9733 a[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
9734 a[2] * (
b[3] *
a[7] -
b[4] *
a[6]) +
9736 a[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
9737 a[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
9738 a[2] * (
a[3] *
b[7] -
a[4] *
b[6]));
9740 c(2) = (
a[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
9741 a[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
9742 a[2] * (
b[3] *
b[7] -
b[4] *
b[6]) +
9744 b[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
9745 b[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
9746 b[2] * (
a[3] *
b[7] -
a[4] *
b[6]) +
9748 b[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
9749 b[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
9750 b[2] * (
b[3] *
a[7] -
b[4] *
a[6]));
9752 c(3) = (
b[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
9753 b[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
9754 b[2] * (
b[3] *
b[7] -
b[4] *
b[6]));
9800 real_t a = z(2),
b = z(1), c = z(0);
9808 x(0) = x(1) = -0.5 *
b /
a;
9813 x(0) = -(x(1) = fabs(0.5 * sqrt(D) /
a));
9821 t = -0.5 * (
b + sqrt(D));
9825 t = -0.5 * (
b - sqrt(D));
9839 real_t a = z(2)/z(3),
b = z(1)/z(3), c = z(0)/z(3);
9843 real_t R = (2 *
a *
a *
a - 9 *
a *
b + 27 * c) / 54;
9851 x(0) = x(1) = x(2) = -
a / 3;
9859 x(0) = -2 * sqrtQ -
a / 3;
9860 x(1) = x(2) = sqrtQ -
a / 3;
9864 x(0) = x(1) = - sqrtQ -
a / 3;
9865 x(2) = 2 * sqrtQ -
a / 3;
9872 real_t theta = acos(R / sqrt(Q3));
9875 x0 = A * cos(theta / 3) -
a / 3;
9876 x1 = A * cos((theta + 2.0 * M_PI) / 3) -
a / 3;
9877 x2 = A * cos((theta - 2.0 * M_PI) / 3) -
a / 3;
9902 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
9906 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
9908 x(0) = A + Q / A -
a / 3;
9917 const real_t factor,
const int Dim)
9920 c(0) = c0 * (1.0 - pow(factor, -Dim));
9922 for (
int j = 0; j < nr; j++)
9934 c(0) = c0 * (1.0 - pow(factor, Dim));
9936 for (
int j = 0; j < nr; j++)
9955 const real_t factor = 2.0;
9970 for (
int k = 0; k < nv; k++)
9973 V(j, k) = displacements(v[k]+j*nvs);
10003 for (
int j = 0; j < nv; j++)
10022 mfem_error(
"Mesh::CheckDisplacements(...)");
10029 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
10030 for (
int j = 0; j <
spaceDim; j++)
10032 vertices[i](j) += displacements(j*nv+i);
10040 for (
int i = 0; i < nv; i++)
10041 for (
int j = 0; j <
spaceDim; j++)
10043 vert_coord(j*nv+i) =
vertices[i](j);
10051 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
10052 for (
int j = 0; j <
spaceDim; j++)
10054 vertices[i](j) = vert_coord(j*nv+i);
10063 for (
int j = 0; j <
spaceDim; j++)
10070 for (
int j = 0; j <
spaceDim; j++)
10082 for (
int j = 0; j <
spaceDim; j++)
10084 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
10089 for (
int j = 0; j <
spaceDim; j++)
10101 (*Nodes) += displacements;
10116 node_coord = (*Nodes);
10128 (*Nodes) = node_coord;
10191 for (j = 1; j < n; j++)
10228 int quad_counter = 0;
10250 const int attr =
elements[i]->GetAttribute();
10251 int *v =
elements[i]->GetVertices();
10257 for (
int ei = 0; ei < 3; ei++)
10259 for (
int k = 0; k < 2; k++)
10266 new_elements[j++] =
10267 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
10268 new_elements[j++] =
10269 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
10270 new_elements[j++] =
10271 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
10272 new_elements[j++] =
10273 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
10277 const int qe = quad_counter;
10281 for (
int ei = 0; ei < 4; ei++)
10283 for (
int k = 0; k < 2; k++)
10290 new_elements[j++] =
10291 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
10292 new_elements[j++] =
10293 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
10294 new_elements[j++] =
10295 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
10296 new_elements[j++] =
10297 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
10301 MFEM_ABORT(
"unknown element type: " << el_type);
10311 const int attr =
boundary[i]->GetAttribute();
10312 int *v =
boundary[i]->GetVertices();
10321 static const real_t A = 0.0, B = 0.5, C = 1.0;
10322 static real_t tri_children[2*3*4] =
10329 static real_t quad_children[2*4*4] =
10331 A,A, B,A, B,B, A,B,
10332 B,A, C,A, C,B, B,B,
10333 B,B, C,B, C,C, B,C,
10343 for (
int i = 0; i <
elements.Size(); i++)
10364 if (!
Nodes || update_nodes)
10394 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
10397 int NumOfQuadFaces = 0;
10403 for (
int i = 0; i <
faces.Size(); i++)
10407 f2qf[i] = NumOfQuadFaces;
10418 int hex_counter = 0;
10421 for (
int i = 0; i <
elements.Size(); i++)
10430 int pyr_counter = 0;
10433 for (
int i = 0; i <
elements.Size(); i++)
10451 DSTable *v_to_v_ptr = v_to_v_p;
10467 std::sort(row_start, J_v2v.
end());
10470 for (
int i = 0; i < J_v2v.
Size(); i++)
10472 e2v[J_v2v[i].two] = i;
10485 it.SetIndex(e2v[it.Index()]);
10495 const int oelem = oface + NumOfQuadFaces;
10508 const int attr =
elements[i]->GetAttribute();
10509 int *v =
elements[i]->GetVertices();
10516 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
10524 for (
int ei = 0; ei < 6; ei++)
10526 for (
int k = 0; k < 2; k++)
10536 const int rt_algo = 1;
10547 real_t len_sqr, min_len;
10549 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
10550 sqr(J(1,0)-J(1,1)-J(1,2)) +
10551 sqr(J(2,0)-J(2,1)-J(2,2));
10554 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
10555 sqr(J(1,1)-J(1,0)-J(1,2)) +
10556 sqr(J(2,1)-J(2,0)-J(2,2));
10557 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
10559 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
10560 sqr(J(1,2)-J(1,0)-J(1,1)) +
10561 sqr(J(2,2)-J(2,0)-J(2,1));
10562 if (len_sqr < min_len) { rt = 2; }
10567 real_t Em_data[18], Js_data[9], Jp_data[9];
10569 DenseMatrix Js(Js_data, 3, 3), Jp(Jp_data, 3, 3);
10572 for (
int s = 0; s < 3; s++)
10574 for (
int t = 0; t < 3; t++)
10576 Em(t,s) = 0.5*J(t,s);
10579 for (
int t = 0; t < 3; t++)
10581 Em(t,3) = 0.5*(J(t,0)+J(t,1));
10582 Em(t,4) = 0.5*(J(t,0)+J(t,2));
10583 Em(t,5) = 0.5*(J(t,1)+J(t,2));
10587 for (
int t = 0; t < 3; t++)
10589 Js(t,0) = Em(t,5)-Em(t,0);
10590 Js(t,1) = Em(t,1)-Em(t,0);
10591 Js(t,2) = Em(t,2)-Em(t,0);
10595 for (
int t = 0; t < 3; t++)
10597 Js(t,0) = Em(t,5)-Em(t,0);
10598 Js(t,1) = Em(t,2)-Em(t,0);
10599 Js(t,2) = Em(t,4)-Em(t,0);
10603 kappa_min = std::max(ar1, ar2);
10607 for (
int t = 0; t < 3; t++)
10609 Js(t,0) = Em(t,0)-Em(t,1);
10610 Js(t,1) = Em(t,4)-Em(t,1);
10611 Js(t,2) = Em(t,2)-Em(t,1);
10615 for (
int t = 0; t < 3; t++)
10617 Js(t,0) = Em(t,2)-Em(t,1);
10618 Js(t,1) = Em(t,4)-Em(t,1);
10619 Js(t,2) = Em(t,5)-Em(t,1);
10623 kappa = std::max(ar1, ar2);
10624 if (
kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
10627 for (
int t = 0; t < 3; t++)
10629 Js(t,0) = Em(t,0)-Em(t,2);
10630 Js(t,1) = Em(t,1)-Em(t,2);
10631 Js(t,2) = Em(t,3)-Em(t,2);
10635 for (
int t = 0; t < 3; t++)
10637 Js(t,0) = Em(t,1)-Em(t,2);
10638 Js(t,1) = Em(t,5)-Em(t,2);
10639 Js(t,2) = Em(t,3)-Em(t,2);
10643 kappa = std::max(ar1, ar2);
10644 if (
kappa < kappa_min) { rt = 2; }
10647 static const int mv_all[3][4][4] =
10649 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
10650 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
10651 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
10653 const int (&mv)[4][4] = mv_all[rt];
10655#ifndef MFEM_USE_MEMALLOC
10656 new_elements[j+0] =
10657 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
10658 new_elements[j+1] =
10659 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
10660 new_elements[j+2] =
10661 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
10662 new_elements[j+3] =
10663 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
10665 for (
int k = 0; k < 4; k++)
10667 new_elements[j+4+k] =
10668 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
10669 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
10673 new_elements[j+0] = tet =
TetMemory.Alloc();
10674 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
10676 new_elements[j+1] = tet =
TetMemory.Alloc();
10677 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
10679 new_elements[j+2] = tet =
TetMemory.Alloc();
10680 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
10682 new_elements[j+3] = tet =
TetMemory.Alloc();
10683 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
10685 for (
int k = 0; k < 4; k++)
10687 new_elements[j+4+k] = tet =
TetMemory.Alloc();
10688 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
10689 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
10692 for (
int k = 0; k < 4; k++)
10697 for (
int k = 0; k < 4; k++)
10711 for (
int fi = 2; fi < 5; fi++)
10713 for (
int k = 0; k < 4; k++)
10720 for (
int ei = 0; ei < 9; ei++)
10722 for (
int k = 0; k < 2; k++)
10729 const int qf2 = f2qf[
f[2]];
10730 const int qf3 = f2qf[
f[3]];
10731 const int qf4 = f2qf[
f[4]];
10733 new_elements[j++] =
10734 new Wedge(v[0], oedge+e[0], oedge+e[2],
10735 oedge+e[6], oface+qf2, oface+qf4, attr);
10737 new_elements[j++] =
10738 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
10739 oface+qf3, oface+qf4, oface+qf2, attr);
10741 new_elements[j++] =
10742 new Wedge(oedge+e[0], v[1], oedge+e[1],
10743 oface+qf2, oedge+e[7], oface+qf3, attr);
10745 new_elements[j++] =
10746 new Wedge(oedge+e[2], oedge+e[1], v[2],
10747 oface+qf4, oface+qf3, oedge+e[8], attr);
10749 new_elements[j++] =
10750 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
10751 v[3], oedge+e[3], oedge+e[5], attr);
10753 new_elements[j++] =
10754 new Wedge(oface+qf3, oface+qf4, oface+qf2,
10755 oedge+e[4], oedge+e[5], oedge+e[3], attr);
10757 new_elements[j++] =
10758 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
10759 oedge+e[3], v[4], oedge+e[4], attr);
10761 new_elements[j++] =
10762 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
10763 oedge+e[5], oedge+e[4], v[5], attr);
10772 for (
int fi = 0; fi < 1; fi++)
10774 for (
int k = 0; k < 4; k++)
10781 for (
int ei = 0; ei < 8; ei++)
10783 for (
int k = 0; k < 2; k++)
10790 const int qf0 = f2qf[
f[0]];
10792 new_elements[j++] =
10793 new Pyramid(v[0], oedge+e[0], oface+qf0,
10794 oedge+e[3], oedge+e[4], attr);
10796 new_elements[j++] =
10797 new Pyramid(oedge+e[0], v[1], oedge+e[1],
10798 oface+qf0, oedge+e[5], attr);
10800 new_elements[j++] =
10801 new Pyramid(oface+qf0, oedge+e[1], v[2],
10802 oedge+e[2], oedge+e[6], attr);
10804 new_elements[j++] =
10805 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
10806 v[3], oedge+e[7], attr);
10808 new_elements[j++] =
10809 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
10810 oedge+e[7], v[4], attr);
10812 new_elements[j++] =
10813 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
10814 oedge+e[4], oface+qf0, attr);
10816#ifndef MFEM_USE_MEMALLOC
10817 new_elements[j++] =
10818 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
10821 new_elements[j++] =
10822 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
10825 new_elements[j++] =
10826 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
10829 new_elements[j++] =
10830 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
10834 new_elements[j++] = tet =
TetMemory.Alloc();
10835 tet->
Init(oedge+e[0], oedge+e[4], oedge+e[5],
10838 new_elements[j++] = tet =
TetMemory.Alloc();
10839 tet->
Init(oedge+e[1], oedge+e[5], oedge+e[6],
10842 new_elements[j++] = tet =
TetMemory.Alloc();
10843 tet->
Init(oedge+e[2], oedge+e[6], oedge+e[7],
10846 new_elements[j++] = tet =
TetMemory.Alloc();
10847 tet->
Init(oedge+e[3], oedge+e[7], oedge+e[4],
10860 const int he = hex_counter;
10865 if (f2qf.
Size() == 0)
10871 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[
f[k]]; }
10877 for (
int fi = 0; fi < 6; fi++)
10879 for (
int k = 0; k < 4; k++)
10886 for (
int ei = 0; ei < 12; ei++)
10888 for (
int k = 0; k < 2; k++)
10895 new_elements[j++] =
10896 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
10897 oedge+e[3], oedge+e[8], oface+qf[1],
10898 oelem+he, oface+qf[4], attr);
10899 new_elements[j++] =
10900 new Hexahedron(oedge+e[0], v[1], oedge+e[1],
10901 oface+qf[0], oface+qf[1], oedge+e[9],
10902 oface+qf[2], oelem+he, attr);
10903 new_elements[j++] =
10904 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
10905 oedge+e[2], oelem+he, oface+qf[2],
10906 oedge+e[10], oface+qf[3], attr);
10907 new_elements[j++] =
10908 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
10909 v[3], oface+qf[4], oelem+he,
10910 oface+qf[3], oedge+e[11], attr);
10911 new_elements[j++] =
10912 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
10913 oface+qf[4], v[4], oedge+e[4],
10914 oface+qf[5], oedge+e[7], attr);
10915 new_elements[j++] =
10916 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
10917 oelem+he, oedge+e[4], v[5],
10918 oedge+e[5], oface+qf[5], attr);
10919 new_elements[j++] =
10920 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
10921 oface+qf[3], oface+qf[5], oedge+e[5],
10922 v[6], oedge+e[6], attr);
10923 new_elements[j++] =
10924 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
10925 oedge+e[11], oedge+e[7], oface+qf[5],
10926 oedge+e[6], v[7], attr);
10931 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
10943 const int attr =
boundary[i]->GetAttribute();
10944 int *v =
boundary[i]->GetVertices();
10951 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
10957 new_boundary[j++] =
10958 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
10959 new_boundary[j++] =
10960 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
10961 new_boundary[j++] =
10962 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
10963 new_boundary[j++] =
10964 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
10971 new_boundary[j++] =
10972 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
10973 new_boundary[j++] =
10974 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
10975 new_boundary[j++] =
10976 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
10977 new_boundary[j++] =
10978 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
10982 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
10988 static const real_t A = 0.0, B = 0.5, C = 1.0, D = -1.0;
10989 static real_t tet_children[3*4*16] =
10991 A,A,A, B,A,A, A,B,A, A,A,B,
10992 B,A,A, C,A,A, B,B,A, B,A,B,
10993 A,B,A, B,B,A, A,C,A, A,B,B,
10994 A,A,B, B,A,B, A,B,B, A,A,C,
10999 B,A,A, A,B,B, A,B,A, A,A,B,
11000 B,A,A, A,B,B, A,A,B, B,A,B,
11001 B,A,A, A,B,B, B,A,B, B,B,A,
11002 B,A,A, A,B,B, B,B,A, A,B,A,
11004 A,B,A, B,A,A, B,A,B, A,A,B,
11005 A,B,A, A,A,B, B,A,B, A,B,B,
11006 A,B,A, A,B,B, B,A,B, B,B,A,
11007 A,B,A, B,B,A, B,A,B, B,A,A,
11009 A,A,B, B,A,A, A,B,A, B,B,A,
11010 A,A,B, A,B,A, A,B,B, B,B,A,
11011 A,A,B, A,B,B, B,A,B, B,B,A,
11012 A,A,B, B,A,B, B,A,A, B,B,A
11014 static real_t pyr_children[3*5*10] =
11016 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
11017 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
11018 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
11019 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
11020 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
11021 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
11022 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
11023 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
11024 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
11025 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
11027 static real_t pri_children[3*6*8] =
11029 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
11030 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
11031 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
11032 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
11033 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
11034 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
11035 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
11036 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
11038 static real_t hex_children[3*8*8] =
11040 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
11041 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
11042 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
11043 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
11044 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
11045 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
11046 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
11047 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
11059 for (
int i = 0; i <
elements.Size(); i++)
11090 int i, j, ind, nedges;
11097 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
11111 for (j = 0; j < marked_el.
Size(); j++)
11116 int new_v = cnv + j, new_e = cne + j;
11125 static real_t seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
11127 UseExternalData(seg_children, 1, 2, 3);
11140 int *edge1 =
new int[nedges];
11141 int *edge2 =
new int[nedges];
11142 int *middle =
new int[nedges];
11144 for (i = 0; i < nedges; i++)
11146 edge1[i] = edge2[i] = middle[i] = -1;
11152 for (j = 1; j < v.
Size(); j++)
11154 ind = v_to_v(v[j-1], v[j]);
11155 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
11157 ind = v_to_v(v[0], v[v.
Size()-1]);
11158 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
11162 for (i = 0; i < marked_el.
Size(); i++)
11168 int need_refinement;
11171 need_refinement = 0;
11172 for (i = 0; i < nedges; i++)
11174 if (middle[i] != -1 && edge1[i] != -1)
11176 need_refinement = 1;
11181 while (need_refinement == 1);
11184 int v1[2], v2[2],
bisect, temp;
11186 for (i = 0; i < temp; i++)
11189 bisect = v_to_v(v[0], v[1]);
11190 if (middle[
bisect] != -1)
11194 v1[0] = v[0]; v1[1] = middle[
bisect];
11195 v2[0] = middle[
bisect]; v2[1] = v[1];
11201 mfem_error(
"Only bisection of segment is implemented"
11225 MFEM_VERIFY(
GetNE() == 0 ||
11227 "tetrahedral mesh is not marked for refinement:"
11228 " call Finalize(true)");
11235 for (i = 0; i < marked_el.
Size(); i++)
11241 for (i = 0; i < marked_el.
Size(); i++)
11250 for (i = 0; i < marked_el.
Size(); i++)
11267 int need_refinement;
11272 need_refinement = 0;
11280 if (
elements[i]->NeedRefinement(v_to_v))
11282 need_refinement = 1;
11287 while (need_refinement == 1);
11294 need_refinement = 0;
11296 if (
boundary[i]->NeedRefinement(v_to_v))
11298 need_refinement = 1;
11302 while (need_refinement == 1);
11333 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
11334 "not supported. Project the NURBS to Nodes first.");
11344 if (!refinements.
Size())
11365 Swap(*mesh2,
false);
11377 const int *fine,
int nfine,
int op)
11379 real_t error = (op == 3) ? std::pow(elem_error[fine[0]],
11380 2.0) : elem_error[fine[0]];
11382 for (
int i = 1; i < nfine; i++)
11384 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
11386 real_t err_fine = elem_error[fine[i]];
11389 case 0: error = std::min(error, err_fine);
break;
11390 case 1: error += err_fine;
break;
11391 case 2: error = std::max(error, err_fine);
break;
11392 case 3: error += std::pow(err_fine, 2.0);
break;
11393 default: MFEM_ABORT(
"Invalid operation.");
11396 return (op == 3) ? std::sqrt(error) : error;
11400 real_t threshold,
int nc_limit,
int op)
11402 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
11403 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
11404 "Project the NURBS to Nodes first.");
11417 for (
int i = 0; i < dt.
Size(); i++)
11419 if (nc_limit > 0 && !level_ok[i]) {
continue; }
11424 if (error < threshold) { derefs.
Append(i); }
11427 if (!derefs.
Size()) {
return false; }
11434 Swap(*mesh2,
false);
11448 int nc_limit,
int op)
11458 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
11465 int nc_limit,
int op)
11468 for (
int i = 0; i < tmp.
Size(); i++)
11470 tmp[i] = elem_error(i);
11513 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
11562#ifdef MFEM_USE_MEMALLOC
11599 for (
int i = 0; i < elem_array.
Size(); i++)
11601 if (elem_array[i]->GetGeometryType() == geom)
11606 elem_vtx.
SetSize(nv*num_elems);
11610 for (
int i = 0; i < elem_array.
Size(); i++)
11616 elem_vtx.
Append(loc_vtx);
11624 for (
int i = 0; i < nelem; i++) { list[i] = i; }
11640 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
11652 default: MFEM_ABORT(
"internal error");
11666 bool noInitialCoarsening =
true;
11667 for (
auto f : initialCoarsening)
11669 noInitialCoarsening = (noInitialCoarsening &&
f == 1);
11672 if (noInitialCoarsening)
11692 bool divisible =
true;
11693 for (
int i=0; i<rf.
Size(); ++i)
11696 divisible = divisible && cf * rf[i] == initialCoarsening[i];
11699 MFEM_VERIFY(divisible,
"Invalid coarsening");
11714 int nonconforming,
int nc_limit)
11724 else if (nonconforming < 0)
11745 for (
int i = 0; i < refinements.
Size(); i++)
11747 el_to_refine[i] = refinements[i].index;
11751 int type, rt = (refinements.
Size() ? refinements[0].GetType() : 7);
11752 if (rt == 1 || rt == 2 || rt == 4)
11756 else if (rt == 3 || rt == 5 || rt == 6)
11774 for (
int i = 0; i < el_to_refine.
Size(); i++)
11776 refinements[i] =
Refinement(el_to_refine[i]);
11783 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
11784 "Please project the NURBS to Nodes first, with SetCurvature().");
11787 MFEM_VERIFY(
ncmesh != NULL ||
dynamic_cast<const ParMesh*
>(
this) == NULL,
11788 "Sorry, converting a conforming ParMesh to an NC mesh is "
11796 (simplices_nonconforming && (
meshgen & 0x1)) )
11809 for (
int i = 0; i <
GetNE(); i++)
11816 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
11828 for (
int i = 0; i <
GetNE(); i++)
11831 bool refine =
false;
11832 for (
int j = 0; j < v.
Size(); j++)
11835 for (
int l = 0; l <
spaceDim; l++)
11840 if (dist <= eps*eps) { refine =
true;
break; }
11851 int nonconforming,
int nc_limit)
11853 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
11855 for (
int i = 0; i <
GetNE(); i++)
11857 if (elem_error[i] > threshold)
11871 int nonconforming,
int nc_limit)
11874 elem_error.
Size());
11875 return RefineByError(tmp, threshold, nonconforming, nc_limit);
11880 int *edge1,
int *edge2,
int *middle)
11883 int v[2][4], v_new,
bisect, t;
11895 bisect = v_to_v(vert[0], vert[1]);
11896 MFEM_ASSERT(
bisect >= 0,
"");
11898 if (middle[
bisect] == -1)
11901 for (
int d = 0; d <
spaceDim; d++)
11926 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
11927 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
11948 bisect = v_to_v(v[1][0], v[1][1]);
11949 MFEM_ASSERT(
bisect >= 0,
"");
11955 else if (edge2[
bisect] == i)
11964 MFEM_ABORT(
"Bisection for now works only for triangles.");
11971 int v[2][4], v_new,
bisect, t;
11981 "TETRAHEDRON element is not marked for refinement.");
11990 for (
int j = 0; j < 3; j++)
12003 int type, old_redges[2], flag;
12006 int new_type, new_redges[2][2];
12009 new_redges[0][0] = 2;
12010 new_redges[0][1] = 1;
12011 new_redges[1][0] = 2;
12012 new_redges[1][1] = 1;
12013 int tr1 = -1, tr2 = -1;
12014 switch (old_redges[0])
12017 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
12022 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
12026 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
12029 switch (old_redges[1])
12032 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
12037 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
12041 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
12048#ifdef MFEM_USE_MEMALLOC
12084 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
12091 int v[2][3], v_new,
bisect, t;
12103 MFEM_ASSERT(
bisect >= 0,
"");
12105 MFEM_ASSERT(v_new != -1,
"");
12109 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
12110 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
12120 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
12126 int *edge1,
int *edge2,
int *middle)
12129 int j, v1[3], v2[3], v3[3], v4[3], v_new[3],
bisect[3];
12138 bisect[0] = v_to_v(v[0],v[1]);
12139 bisect[1] = v_to_v(v[1],v[2]);
12140 bisect[2] = v_to_v(v[0],v[2]);
12143 for (j = 0; j < 3; j++)
12145 if (middle[
bisect[j]] == -1)
12148 for (
int d = 0; d <
spaceDim; d++)
12156 if (edge1[
bisect[j]] == i)
12161 middle[
bisect[j]] = v_new[j];
12165 v_new[j] = middle[
bisect[j]];
12174 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
12175 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
12176 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
12177 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
12211 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
12247 for (
int i = 0; i < elem_geoms.
Size(); i++)
12255 std::map<unsigned, int> mat_no;
12259 for (
int j = 0; j <
elements.Size(); j++)
12262 unsigned code =
elements[j]->GetTransform();
12265 int &matrix = mat_no[code];
12266 if (!matrix) { matrix =
static_cast<int>(mat_no.size()); }
12273 pmats.
SetSize(
Dim,
Dim+1,
static_cast<int>((mat_no.size())));
12276 std::map<unsigned, int>::iterator it;
12277 for (it = mat_no.begin(); it != mat_no.end(); ++it)
12291 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
12292 " geom = " << geom);
12302 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
12311 os <<
"areamesh2\n\n";
12315 os <<
"curved_areamesh2\n\n";
12324 os <<
boundary[i]->GetAttribute();
12325 for (j = 0; j < v.
Size(); j++)
12327 os <<
' ' << v[j] + 1;
12339 for (j = 0; j < v.
Size(); j++)
12341 os <<
' ' << v[j] + 1;
12353 for (j = 1; j <
Dim; j++)
12370 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
12378 os <<
"NETGEN_Neutral_Format\n";
12383 for (j = 0; j <
Dim; j++)
12396 os <<
elements[i]->GetAttribute();
12397 for (j = 0; j < nv; j++)
12399 os <<
' ' << ind[j]+1;
12410 os <<
boundary[i]->GetAttribute();
12411 for (j = 0; j < nv; j++)
12413 os <<
' ' << ind[j]+1;
12425 <<
" 0 0 0 0 0 0 0\n"
12426 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
12428 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
12429 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
12433 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
12439 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
12440 for (j = 0; j < nv; j++)
12442 os <<
' ' << ind[j]+1;
12451 os <<
boundary[i]->GetAttribute();
12452 for (j = 0; j < nv; j++)
12454 os <<
' ' << ind[j]+1;
12456 os <<
" 1.0 1.0 1.0 1.0\n";
12465 const std::string &comments)
const
12500 os <<
"\n# mesh curvature GridFunction";
12505 os <<
"\nmfem_mesh_end" << endl;
12512 os << (!set_names && section_delimiter.empty()
12513 ?
"MFEM mesh v1.0\n" :
12514 (!set_names ?
"MFEM mesh v1.2\n" :
"MFEM mesh v1.3\n"));
12516 if (set_names && section_delimiter.empty())
12518 section_delimiter =
"mfem_mesh_end";
12522 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
12525 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12530 "# TETRAHEDRON = 4\n"
12536 os <<
"\ndimension\n" <<
Dim;
12546 os <<
"\nattribute_sets\n";
12558 os <<
"\nbdr_attribute_sets\n";
12583 if (!section_delimiter.empty())
12586 << section_delimiter << endl;
12591 const int version,
const std::string &comments)
const
12593 MFEM_VERIFY(version == 10 || version == 11,
"Invalid NURBS mesh version");
12598 os <<
"MFEM NURBS mesh v" << int(version / 10) <<
"." << version % 10 <<
"\n";
12601 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
12604 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12610 os <<
"\ndimension\n" <<
Dim
12639 os <<
"\nedges\n" << ne <<
'\n';
12640 for (
int i = 0; i < ne; i++)
12642 const int *v =
elements[i]->GetVertices();
12643 int v0 = v[0], v1 = v[1];
12645 int ki = e_to_k[i];
12646 const bool flip = (ki < 0);
12647 if (flip) { ki = -1 - ki; }
12651 if ((v0 > v1) != flip) { std::swap(v0, v1); }
12653 os << ki <<
' ' << v0 <<
' ' << v1 <<
'\n';
12671 for (
int j=0; j<2; ++j)
12679 const int s = vert[0];
12685 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
12696 ofstream ofs(fname);
12697 ofs.precision(precision);
12701#ifdef MFEM_USE_ADIOS2
12711 "# vtk DataFile Version 3.0\n"
12712 "Generated by MFEM\n"
12714 "DATASET UNSTRUCTURED_GRID\n";
12727 for ( ; j < 3; j++)
12738 for (
int i = 0; i <
Nodes->
FESpace()->GetNDofs(); i++)
12743 os << (*Nodes)(vdofs[0]);
12747 os <<
' ' << (*Nodes)(vdofs[j]);
12749 for ( ; j < 3; j++)
12763 size +=
elements[i]->GetNVertices() + 1;
12768 const int *v =
elements[i]->GetVertices();
12769 const int nv =
elements[i]->GetNVertices();
12773 for (
int j = 0; j < nv; j++)
12775 os <<
' ' << v[perm ? perm[j] : j];
12788 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
12789 "Point meshes should have a single dof per element");
12790 size += dofs.
Size() + 1;
12795 if (!strcmp(fec_name,
"Linear") ||
12796 !strcmp(fec_name,
"H1_0D_P1") ||
12797 !strcmp(fec_name,
"H1_1D_P1") ||
12798 !strcmp(fec_name,
"H1_2D_P1") ||
12799 !strcmp(fec_name,
"H1_3D_P1"))
12803 else if (!strcmp(fec_name,
"Quadratic") ||
12804 !strcmp(fec_name,
"H1_1D_P2") ||
12805 !strcmp(fec_name,
"H1_2D_P2") ||
12806 !strcmp(fec_name,
"H1_3D_P2"))
12812 mfem::err <<
"Mesh::PrintVTK : can not save '"
12813 << fec_name <<
"' elements!" << endl;
12822 for (
int j = 0; j < dofs.
Size(); j++)
12824 os <<
' ' << dofs[j];
12827 else if (order == 2)
12829 const int *vtk_mfem;
12830 switch (
elements[i]->GetGeometryType())
12844 for (
int j = 0; j < dofs.
Size(); j++)
12846 os <<
' ' << dofs[vtk_mfem[j]];
12856 int vtk_cell_type = 5;
12860 os << vtk_cell_type <<
'\n';
12865 <<
"SCALARS material int\n"
12866 <<
"LOOKUP_TABLE default\n";
12869 os <<
elements[i]->GetAttribute() <<
'\n';
12876 bool high_order_output,
12877 int compression_level,
12880 int ref = (high_order_output &&
Nodes)
12883 fname = fname +
".vtu";
12884 std::fstream os(fname.c_str(),std::ios::out);
12885 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"2.2\"";
12886 if (compression_level != 0)
12888 os <<
" compressor=\"vtkZLibDataCompressor\"";
12891 os <<
"<UnstructuredGrid>\n";
12892 PrintVTU(os, ref, format, high_order_output, compression_level, bdr_elements);
12893 os <<
"</Piece>\n";
12894 os <<
"</UnstructuredGrid>\n";
12895 os <<
"</VTKFile>" << std::endl;
12902 bool high_order_output,
12903 int compression_level)
12905 PrintVTU(fname, format, high_order_output, compression_level,
true);
12909 bool high_order_output,
int compression_level,
12917 std::vector<char> buf;
12919 auto get_geom = [&](
int i)
12927 int np = 0, nc_ref = 0;
12928 for (
int i = 0; i < ne; i++)
12937 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
12938 << (high_order_output ? ne : nc_ref) <<
"\">\n";
12941 os <<
"<Points>\n";
12942 os <<
"<DataArray type=\"" << type_str
12943 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
12944 for (
int i = 0; i < ne; i++)
12957 for (
int j = 0; j < pmat.
Width(); j++)
12983 os <<
"</DataArray>" << std::endl;
12984 os <<
"</Points>" << std::endl;
12986 os <<
"<Cells>" << std::endl;
12987 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
12988 << fmt_str <<
"\">" << std::endl;
12990 std::vector<int> offset;
12993 if (high_order_output)
12996 for (
int iel = 0; iel < ne; iel++)
13000 int nnodes = local_connectivity.
Size();
13001 for (
int i=0; i<nnodes; ++i)
13008 offset.push_back(np);
13014 for (
int i = 0; i < ne; i++)
13020 for (
int j = 0; j < RG.
Size(); )
13023 offset.push_back(coff);
13025 for (
int k = 0; k < nv; k++, j++)
13039 os <<
"</DataArray>" << std::endl;
13041 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
13042 << fmt_str <<
"\">" << std::endl;
13044 for (
size_t ii=0; ii<offset.size(); ii++)
13052 os <<
"</DataArray>" << std::endl;
13053 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
13054 << fmt_str <<
"\">" << std::endl;
13056 const int *vtk_geom_map =
13058 for (
int i = 0; i < ne; i++)
13061 uint8_t vtk_cell_type = 5;
13063 vtk_cell_type = vtk_geom_map[geom];
13065 if (high_order_output)
13074 for (
int j = 0; j < RG.
Size(); j += nv)
13084 os <<
"</DataArray>" << std::endl;
13085 os <<
"</Cells>" << std::endl;
13087 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
13088 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
13089 << fmt_str <<
"\">" << std::endl;
13090 for (
int i = 0; i < ne; i++)
13093 if (high_order_output)
13112 os <<
"</DataArray>" << std::endl;
13113 os <<
"</CellData>" << std::endl;
13124 "# vtk DataFile Version 3.0\n"
13125 "Generated by MFEM\n"
13127 "DATASET UNSTRUCTURED_GRID\n";
13132 os <<
"FIELD FieldData 1\n"
13142 np = nc = size = 0;
13143 for (
int i = 0; i <
GetNE(); i++)
13152 os <<
"POINTS " << np <<
" double\n";
13154 for (
int i = 0; i <
GetNE(); i++)
13161 for (
int j = 0; j < pmat.
Width(); j++)
13163 os << pmat(0, j) <<
' ';
13166 os << pmat(1, j) <<
' ';
13178 os << 0.0 <<
' ' << 0.0;
13185 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
13187 for (
int i = 0; i <
GetNE(); i++)
13194 for (
int j = 0; j < RG.
Size(); )
13197 for (
int k = 0; k < nv; k++, j++)
13199 os <<
' ' << np + RG[j];
13205 os <<
"CELL_TYPES " << nc <<
'\n';
13206 for (
int i = 0; i <
GetNE(); i++)
13214 for (
int j = 0; j < RG.
Size(); j += nv)
13216 os << vtk_cell_type <<
'\n';
13220 os <<
"CELL_DATA " << nc <<
'\n'
13221 <<
"SCALARS material int\n"
13222 <<
"LOOKUP_TABLE default\n";
13223 for (
int i = 0; i <
GetNE(); i++)
13231 os << attr <<
'\n';
13238 srand((
unsigned)time(0));
13240 int el0 = (int)floor(
a *
GetNE());
13242 os <<
"SCALARS element_coloring int\n"
13243 <<
"LOOKUP_TABLE default\n";
13244 for (
int i = 0; i <
GetNE(); i++)
13251 os << coloring[i] + 1 <<
'\n';
13257 os <<
"POINT_DATA " << np <<
'\n' << flush;
13260#ifdef MFEM_USE_HDF5
13267#ifdef MFEM_PARALLEL_HDF5
13268 VTKHDF vtkhdf(fname, pmesh->GetComm());
13269 vtkhdf.
SaveMesh(*
this, high_order);
13272 MFEM_ABORT(
"Requires HDF5 library with parallel support enabled");
13277 vtkhdf.
SaveMesh(*
this, high_order);
13284 int delete_el_to_el = (
el_to_el) ? (0) : (1);
13286 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
13289 const int *i_el_el = el_el.
GetI();
13290 const int *j_el_el = el_el.
GetJ();
13295 stack_p = stack_top_p = 0;
13296 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
13298 if (colors[el] != -2)
13304 el_stack[stack_top_p++] = el;
13306 for ( ; stack_p < stack_top_p; stack_p++)
13308 int i = el_stack[stack_p];
13309 int num_nb = i_el_el[i+1] - i_el_el[i];
13310 if (max_num_col < num_nb + 1)
13312 max_num_col = num_nb + 1;
13314 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
13316 int k = j_el_el[j];
13317 if (colors[k] == -2)
13320 el_stack[stack_top_p++] = k;
13328 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
13330 int i = el_stack[stack_p], col;
13332 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
13334 col = colors[j_el_el[j]];
13337 col_marker[col] = 1;
13341 for (col = 0; col < max_num_col; col++)
13342 if (col_marker[col] == 0)
13350 if (delete_el_to_el)
13358 int elem_attr)
const
13360 if (
Dim != 3 &&
Dim != 2) {
return; }
13362 int i, j, k, l, nv, nbe, *v;
13364 os <<
"MFEM mesh v1.0\n";
13368 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
13373 "# TETRAHEDRON = 4\n"
13378 os <<
"\ndimension\n" <<
Dim
13382 os << int((elem_attr) ? partitioning[i]+1 :
elements[i]->GetAttribute())
13383 <<
' ' <<
elements[i]->GetGeometryType();
13386 for (j = 0; j < nv; j++)
13398 l = partitioning[l];
13413 os <<
"\nboundary\n" << nbe <<
'\n';
13419 l = partitioning[l];
13422 nv =
faces[i]->GetNVertices();
13423 v =
faces[i]->GetVertices();
13424 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
13425 for (j = 0; j < nv; j++)
13432 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
13433 for (j = nv-1; j >= 0; j--)
13444 nv =
faces[i]->GetNVertices();
13445 v =
faces[i]->GetVertices();
13446 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
13447 for (j = 0; j < nv; j++)
13478 int interior_faces)
13480 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
13481 if (
Dim != 3 &&
Dim != 2) {
return; }
13490 int nv =
elements[i]->GetNVertices();
13491 const int *ind =
elements[i]->GetVertices();
13492 for (
int j = 0; j < nv; j++)
13502 voff[i] = vcount[i-1] + voff[i-1];
13508 vown[i] =
new int[vcount[i]];
13520 int nv =
elements[i]->GetNVertices();
13521 const int *ind =
elements[i]->GetVertices();
13522 for (
int j = 0; j < nv; j++)
13525 vown[ind[j]][vcount[ind[j]]] = i;
13531 vcount[i] = voff[i+1] - voff[i];
13535 for (
int i = 0; i < edge_el.
Size(); i++)
13537 const int *el = edge_el.
GetRow(i);
13540 int k = partitioning[el[0]];
13541 int l = partitioning[el[1]];
13542 if (interior_faces || k != l)
13554 os <<
"areamesh2\n\n" << nbe <<
'\n';
13556 for (
int i = 0; i < edge_el.
Size(); i++)
13558 const int *el = edge_el.
GetRow(i);
13561 int k = partitioning[el[0]];
13562 int l = partitioning[el[1]];
13563 if (interior_faces || k != l)
13568 for (
int j = 0; j < 2; j++)
13569 for (
int s = 0; s < vcount[ev[j]]; s++)
13570 if (vown[ev[j]][s] == el[0])
13572 os <<
' ' << voff[ev[j]]+s+1;
13576 for (
int j = 1; j >= 0; j--)
13577 for (
int s = 0; s < vcount[ev[j]]; s++)
13578 if (vown[ev[j]][s] == el[1])
13580 os <<
' ' << voff[ev[j]]+s+1;
13587 int k = partitioning[el[0]];
13591 for (
int j = 0; j < 2; j++)
13592 for (
int s = 0; s < vcount[ev[j]]; s++)
13593 if (vown[ev[j]][s] == el[0])
13595 os <<
' ' << voff[ev[j]]+s+1;
13605 int nv =
elements[i]->GetNVertices();
13606 const int *ind =
elements[i]->GetVertices();
13607 os << partitioning[i]+1 <<
' ';
13609 for (
int j = 0; j < nv; j++)
13611 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
13612 vown[ind[j]][vcount[ind[j]]] = i;
13619 vcount[i] = voff[i+1] - voff[i];
13625 for (
int k = 0; k < vcount[i]; k++)
13627 for (
int j = 0; j <
Dim; j++)
13637 os <<
"NETGEN_Neutral_Format\n";
13641 for (
int k = 0; k < vcount[i]; k++)
13643 for (
int j = 0; j <
Dim; j++)
13654 int nv =
elements[i]->GetNVertices();
13655 const int *ind =
elements[i]->GetVertices();
13656 os << partitioning[i]+1;
13657 for (
int j = 0; j < nv; j++)
13659 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
13660 vown[ind[j]][vcount[ind[j]]] = i;
13667 vcount[i] = voff[i+1] - voff[i];
13677 int k = partitioning[
faces_info[i].Elem1No];
13678 l = partitioning[l];
13679 if (interior_faces || k != l)
13696 int k = partitioning[
faces_info[i].Elem1No];
13697 l = partitioning[l];
13698 if (interior_faces || k != l)
13700 int nv =
faces[i]->GetNVertices();
13701 const int *ind =
faces[i]->GetVertices();
13703 for (
int j = 0; j < nv; j++)
13704 for (
int s = 0; s < vcount[ind[j]]; s++)
13705 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13707 os <<
' ' << voff[ind[j]]+s+1;
13711 for (
int j = nv-1; j >= 0; j--)
13712 for (
int s = 0; s < vcount[ind[j]]; s++)
13713 if (vown[ind[j]][s] ==
faces_info[i].Elem2No)
13715 os <<
' ' << voff[ind[j]]+s+1;
13722 int k = partitioning[
faces_info[i].Elem1No];
13723 int nv =
faces[i]->GetNVertices();
13724 const int *ind =
faces[i]->GetVertices();
13726 for (
int j = 0; j < nv; j++)
13727 for (
int s = 0; s < vcount[ind[j]]; s++)
13728 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13730 os <<
' ' << voff[ind[j]]+s+1;
13746 int k = partitioning[
faces_info[i].Elem1No];
13747 l = partitioning[l];
13748 if (interior_faces || k != l)
13761 <<
" 0 0 0 0 0 0 0\n"
13762 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
13763 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
13764 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
13765 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
13768 for (
int k = 0; k < vcount[i]; k++)
13769 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
13774 int nv =
elements[i]->GetNVertices();
13775 const int *ind =
elements[i]->GetVertices();
13776 os << i+1 <<
' ' << partitioning[i]+1;
13777 for (
int j = 0; j < nv; j++)
13779 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
13780 vown[ind[j]][vcount[ind[j]]] = i;
13787 vcount[i] = voff[i+1] - voff[i];
13796 int k = partitioning[
faces_info[i].Elem1No];
13797 l = partitioning[l];
13798 if (interior_faces || k != l)
13800 int nv =
faces[i]->GetNVertices();
13801 const int *ind =
faces[i]->GetVertices();
13803 for (
int j = 0; j < nv; j++)
13804 for (
int s = 0; s < vcount[ind[j]]; s++)
13805 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13807 os <<
' ' << voff[ind[j]]+s+1;
13809 os <<
" 1.0 1.0 1.0 1.0\n";
13811 for (
int j = nv-1; j >= 0; j--)
13812 for (
int s = 0; s < vcount[ind[j]]; s++)
13813 if (vown[ind[j]][s] ==
faces_info[i].Elem2No)
13815 os <<
' ' << voff[ind[j]]+s+1;
13817 os <<
" 1.0 1.0 1.0 1.0\n";
13822 int k = partitioning[
faces_info[i].Elem1No];
13823 int nv =
faces[i]->GetNVertices();
13824 const int *ind =
faces[i]->GetVertices();
13826 for (
int j = 0; j < nv; j++)
13827 for (
int s = 0; s < vcount[ind[j]]; s++)
13828 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13830 os <<
' ' << voff[ind[j]]+s+1;
13832 os <<
" 1.0 1.0 1.0 1.0\n";
13856 " NURBS mesh is not supported!");
13860 os <<
"MFEM mesh v1.0\n";
13864 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
13869 "# TETRAHEDRON = 4\n"
13874 os <<
"\ndimension\n" <<
Dim
13882 const int *
const i_AF_f = Aface_face.
GetI();
13883 const int *
const j_AF_f = Aface_face.
GetJ();
13885 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
13886 for (
const int * iface = j_AF_f + i_AF_f[iAF];
13887 iface < j_AF_f + i_AF_f[iAF+1];
13890 os << iAF+1 <<
' ';
13923 int *nbea =
new int[na];
13930 for (i = 0; i < na; i++)
13942 for (k = 0; k < vert.
Size(); k++)
13954 for (k = 0; k < vert.
Size(); k++)
13955 if (vn[vert[k]] == 1)
13960 cg[bea*
spaceDim+j] += pointmat(j,k);
13971 for (k = 0; k < vert.
Size(); k++)
13976 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
13993 int *nbea =
new int[na];
14000 for (i = 0; i < na; i++)
14012 for (k = 0; k < vert.
Size(); k++)
14024 for (k = 0; k < vert.
Size(); k++)
14025 if (vn[vert[k]] == 1)
14030 cg[bea*
spaceDim+j] += pointmat(j,k);
14041 for (k = 0; k < vert.
Size(); k++)
14046 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
14062 for (
int i = 0; i <
vertices.Size(); i++)
14064 for (
int j = 0; j <
spaceDim; j++)
14085 "incompatible vector dimensions");
14093 for (
int d = 0; d <
spaceDim; d++)
14113 for (
int i = 0; i <
GetNE(); i++)
14118 for (
int j = 0; j < nv; j++)
14123 for (
int i = 0; i <
GetNBE(); i++)
14128 for (
int j = 0; j < nv; j++)
14134 for (
int i = 0; i < v2v.
Size(); i++)
14139 v2v[i] = num_vert++;
14143 if (num_vert == v2v.
Size()) {
return; }
14145 Vector nodes_by_element;
14150 for (
int i = 0; i <
GetNE(); i++)
14157 for (
int i = 0; i <
GetNE(); i++)
14166 for (
int i = 0; i <
GetNE(); i++)
14171 for (
int j = 0; j < nv; j++)
14176 for (
int i = 0; i <
GetNBE(); i++)
14181 for (
int j = 0; j < nv; j++)
14205 for (
int i = 0; i <
GetNE(); i++)
14218 int num_bdr_elem = 0;
14219 int new_bel_to_edge_nnz = 0;
14220 for (
int i = 0; i <
GetNBE(); i++)
14236 if (num_bdr_elem ==
GetNBE()) {
return; }
14240 Table *new_bel_to_edge = NULL;
14242 new_be_to_face.
Reserve(num_bdr_elem);
14245 new_bel_to_edge =
new Table;
14246 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
14248 for (
int i = 0; i <
GetNBE(); i++)
14253 int row = new_be_to_face.
Size();
14259 int *new_e = new_bel_to_edge->
GetRow(row);
14260 for (
int j = 0; j < ne; j++)
14264 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
14281 for (
int i = 0; i < attribs.
Size(); i++)
14293#ifdef MFEM_USE_MEMALLOC
14320 const int npts = point_mat.
Width();
14321 if (!npts) {
return 0; }
14322 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
14326 if (!
GetNE()) {
return 0; }
14335 min_dist = std::numeric_limits<real_t>::max();
14339 for (
int i = 0; i <
GetNE(); i++)
14343 for (
int k = 0; k < npts; k++)
14346 if (dist < min_dist(k))
14348 min_dist(k) = dist;
14357 for (
int k = 0; k < npts; k++)
14361 int res = inv_tr->
Transform(pt, ips[k]);
14364 elem_ids[k] = e_idx[k];
14368 if (pts_found != npts)
14372 for (
int k = 0; k < npts; k++)
14374 if (elem_ids[k] != -1) {
continue; }
14378 for (
int v = 0; v < elvertices.
Size(); v++)
14380 int vv = elvertices[v];
14382 const int* els = vtoel->
GetRow(vv);
14383 for (
int e = 0; e < ne; e++)
14385 if (els[e] == e_idx[k]) {
continue; }
14387 int res = inv_tr->
Transform(pt, ips[k]);
14390 elem_ids[k] = els[e];
14402 for (
int e = 0; e < neigh.
Size(); e++)
14408 int res = inv_tr->
Transform(pt, ips[k]);
14421 if (inv_trans == NULL) {
delete inv_tr; }
14423 if (warn && pts_found != npts)
14425 MFEM_WARNING((npts-pts_found) <<
" points were not found");
14440 MFEM_VERIFY(
Dim == 2 ||
Dim == 3,
"Only 2D/3D meshes supported right now.");
14441 MFEM_VERIFY(
Dim ==
spaceDim,
"Surface meshes not currently supported.");
14458 skew(0) = std::atan2(J.
Det(), col1 * col2);
14461 ori(0) = std::atan2(J(1,0), J(0,0));
14468 Vector col1, col2, col3;
14479 col1unit *= 1.0/len1;
14480 col2unit *= 1.0/len2;
14481 col3unit *= 1.0/len3;
14487 aspr(0) = len1/std::sqrt(len2*len3),
14488 aspr(1) = len2/std::sqrt(len1*len3);
14491 aspr(2) = std::sqrt(len1/(len2*len3)),
14492 aspr(3) = std::sqrt(len2/(len1*len3));
14495 Vector crosscol12, crosscol13;
14496 col1.
cross3D(col2, crosscol12);
14497 col1.
cross3D(col3, crosscol13);
14498 skew(0) = std::acos(col1unit*col2unit);
14499 skew(1) = std::acos(col1unit*col3unit);
14500 skew(2) = std::atan(len1*volume/(crosscol12*crosscol13));
14506 for (
int d=0; d<
Dim; d++) { rot(d, 0) = col1unit(d); }
14510 rot1 *= col1unit*col2unit;
14512 col1unit.
cross3D(col2unit, rot1);
14514 for (
int d=0; d <
Dim; d++) { rot(d, 1) = rot2(d); }
14517 for (
int d=0; d <
Dim; d++) { rot(d, 2) = rot1(d); }
14518 real_t delta = sqrt(pow(rot(2,1)-rot(1,2), 2.0) +
14519 pow(rot(0,2)-rot(2,0), 2.0) +
14520 pow(rot(1,0)-rot(0,1), 2.0));
14525 for (
int d = 0; d <
Dim; d++) { Iden(d, d) = 1.0; };
14531 MFEM_ABORT(
"Invalid rotation matrix. Contact TMOP Developers.");
14536 ori(0) = (1./
delta)*(rot(2,1)-rot(1,2));
14537 ori(1) = (1./
delta)*(rot(0,2)-rot(2,0));
14538 ori(2) = (1./
delta)*(rot(1,0)-rot(0,1));
14539 ori(3) = std::acos(0.5*(rot.
Trace()-1.0));
14548 entity_to_vertex(entity_to_vertex_)
14550 int geom_offset = 0;
14564 while (geom_offsets[geom+1] <= bytype_entity_id) { geom++; }
14568 const int geom_elem_id = bytype_entity_id - geom_offsets[geom];
14570 return { geom, nv, v };
14575 os <<
"MFEM mesh v1.2\n";
14579 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
14584 "# TETRAHEDRON = 4\n"
14591 os <<
"\ndimension\n" <<
dim;
14597 "invalid MeshPart state");
14600 "invalid MeshPart state");
14601 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
14603 const int bytype_elem_id = have_element_map ?
14608 for (
int i = 0; i < ent.
num_verts; i++)
14610 os <<
' ' << ent.
verts[i];
14620 "invalid MeshPart state");
14623 "invalid MeshPart state");
14626 const int bytype_bdr_id = have_boundary_map ?
14631 for (
int i = 0; i < ent.
num_verts; i++)
14633 os <<
' ' << ent.
verts[i];
14643 os << sdim <<
'\n';
14647 for (
int d = 1; d < sdim; d++)
14660 os <<
"\nmfem_serial_mesh_end\n";
14664 os <<
"\ncommunication_groups\n";
14665 os <<
"number_of_groups " << num_groups <<
"\n\n";
14667 os <<
"# number of entities in each group, followed by ranks in group\n";
14668 for (
int group_id = 0; group_id < num_groups; ++group_id)
14673 for (
int group_member_index = 0; group_member_index < group_size;
14674 ++group_member_index)
14676 os <<
' ' << group_ptr[group_member_index];
14687 MFEM_VERIFY(g2v.
RowSize(0) == 0,
"internal erroor");
14691 MFEM_VERIFY(g2ev.
RowSize(0) == 0,
"internal erroor");
14696 MFEM_VERIFY(g2tv.
RowSize(0) == 0,
"internal erroor");
14697 MFEM_VERIFY(g2qv.
RowSize(0) == 0,
"internal erroor");
14698 const int total_shared_faces =
14700 os <<
"total_shared_faces " << total_shared_faces <<
'\n';
14702 os <<
"\n# group 0 has no shared entities\n";
14703 for (
int gr = 1; gr < num_groups; gr++)
14706 const int nv = g2v.
RowSize(gr);
14707 const int *sv = g2v.
GetRow(gr);
14708 os <<
"\n# group " << gr <<
"\nshared_vertices " << nv <<
'\n';
14709 for (
int i = 0; i < nv; i++)
14711 os << sv[i] <<
'\n';
14716 const int ne = g2ev.
RowSize(gr)/2;
14717 const int *se = g2ev.
GetRow(gr);
14718 os <<
"\nshared_edges " << ne <<
'\n';
14719 for (
int i = 0; i < ne; i++)
14721 const int *v = se + 2*i;
14722 os << v[0] <<
' ' << v[1] <<
'\n';
14727 const int nt = g2tv.
RowSize(gr)/3;
14728 const int *st = g2tv.
GetRow(gr);
14729 const int nq = g2qv.
RowSize(gr)/4;
14731 os <<
"\nshared_faces " << nt+nq <<
'\n';
14732 for (
int i = 0; i < nt; i++)
14735 const int *v = st + 3*i;
14736 for (
int j = 0; j < 3; j++) { os <<
' ' << v[j]; }
14739 for (
int i = 0; i < nq; i++)
14742 const int *v =
sq + 4*i;
14743 for (
int j = 0; j < 4; j++) { os <<
' ' << v[j]; }
14750 os <<
"\nmfem_mesh_end" << endl;
14767 "invalid MeshPart state");
14770 "invalid MeshPart state");
14772 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
14774 const int bytype_elem_id = have_element_map ?
14785 static_cast<Tetrahedron*
>(el)->SetRefinementFlag(ref_flag);
14787 mesh->AddElement(el);
14795 "invalid MeshPart state");
14798 "invalid MeshPart state");
14801 const int bytype_bdr_id = have_boundary_map ?
14807 mesh->AddBdrElement(bdr);
14814 MFEM_ASSERT(!
nodes,
"invalid MeshPart state");
14815 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
14823 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
14825 mesh->AddVertex(0., 0., 0.);
14830 mesh->FinalizeTopology(
false);
14838 const int *partitioning_,
14870 for (
int i = 0; i < boundary_to_part.
Size(); i++)
14872 int face, o, el1, el2;
14875 boundary_to_part[i] =
14881 for (
int i = 0; i < boundary_to_part.
Size(); i++)
14890 for (
int i = 0; i < boundary_to_part.
Size(); i++)
14905 delete vert_element;
14912 MFEM_VERIFY(0 <= part_id && part_id < num_parts,
14913 "invalid part_id = " << part_id
14914 <<
", num_parts = " << num_parts);
14947 mesh_part.
nodes.reset(
nullptr);
14949 mesh_part.
mesh.reset(
nullptr);
14957 int geom_marker = 0, num_geom = 0;
14958 for (
int i = 0; i < num_elems; i++)
14964 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
14966 "overflow in 'entity_to_vertex[geom]', geom: "
14991 if ((geom_marker & (1 << geom)) == 0)
14993 geom_marker |= (1 << geom);
15008 offsets[g] = offset;
15012 for (
int i = 0; i < num_elems; i++)
15023 geom_marker = 0; num_geom = 0;
15024 for (
int i = 0; i < num_bdr_elems; i++)
15030 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
15032 "overflow in 'entity_to_vertex[geom]', geom: "
15036 if ((geom_marker & (1 << geom)) == 0)
15038 geom_marker |= (1 << geom);
15049 offsets[g] = offset;
15053 for (
int i = 0; i < num_bdr_elems; i++)
15064 std::unordered_set<int> vertex_set;
15065 for (
int i = 0; i < num_elems; i++)
15071 vertex_set.insert(v, v + nv);
15073 vertex_loc_to_glob.
SetSize(
static_cast<int>(vertex_set.size()));
15074 std::copy(vertex_set.begin(), vertex_set.end(),
15075 vertex_loc_to_glob.
begin());
15077 vertex_loc_to_glob.
Sort();
15087 for (
int i = 0; i < vert_array.
Size(); i++)
15089 const int glob_id = vert_array[i];
15090 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
15091 MFEM_ASSERT(loc_id >= 0,
"internal error: global vertex id not found");
15092 vert_array[i] = loc_id;
15099 MFEM_VERIFY(numeric_limits<int>::max()/sdim >= vertex_loc_to_glob.
Size(),
15100 "overflow in 'vertex_coordinates', num_vertices = "
15101 << vertex_loc_to_glob.
Size() <<
", sdim = " << sdim);
15103 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
15106 for (
int d = 0; d < sdim; d++)
15123 mesh_part.
mesh->NewNodes(*mesh_part.
nodes,
false);
15145 std::unordered_set<int> face_set;
15148 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
15150 const int glob_elem_id = elem_list[loc_elem_id];
15151 const int nfaces = elem_to_face.
RowSize(glob_elem_id);
15152 const int *faces = elem_to_face.
GetRow(glob_elem_id);
15153 face_set.insert(faces, faces + nfaces);
15157 for (
int glob_face_id : face_set)
15161 if (el[1] < 0) {
continue; }
15164 MFEM_ASSERT(el[0] == part_id || el[1] == part_id,
"internal error");
15165 if (el[0] != part_id || el[1] != part_id)
15168 const int group_id = groups.
Insert(group);
15172 shared_faces.
Sort();
15183 std::unordered_set<int> edge_set;
15186 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
15188 const int glob_elem_id = elem_list[loc_elem_id];
15189 const int nedges = elem_to_edge.
RowSize(glob_elem_id);
15190 const int *edges = elem_to_edge.
GetRow(glob_elem_id);
15191 edge_set.insert(edges, edges + nedges);
15195 for (
int glob_edge_id : edge_set)
15201 for (
int j = 0; j < nelem; j++)
15207 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
15208 if (group.
Size() > 1)
15210 const int group_id = groups.
Insert(group);
15214 shared_edges.
Sort();
15225 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
15228 const int glob_vertex_id = vertex_loc_to_glob[i];
15233 for (
int j = 0; j < nelem; j++)
15239 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
15240 if (group.
Size() > 1)
15242 const int group_id = groups.
Insert(group);
15249 const int num_groups = groups.
Size();
15255 Table &group__shared_vertex_to_vertex =
15257 group__shared_vertex_to_vertex.
MakeI(num_groups);
15258 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
15260 const int group_id = shared_verts[sv].two;
15263 group__shared_vertex_to_vertex.
MakeJ();
15264 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
15266 const int glob_vertex_id = shared_verts[sv].one;
15267 const int group_id = shared_verts[sv].two;
15268 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(glob_vertex_id);
15269 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
15270 group__shared_vertex_to_vertex.
AddConnection(group_id, loc_vertex_id);
15272 group__shared_vertex_to_vertex.
ShiftUpI();
15277 Table &group__shared_edge_to_vertex =
15279 group__shared_edge_to_vertex.
MakeI(num_groups);
15280 for (
int se = 0; se < shared_edges.
Size(); se++)
15282 const int group_id = shared_edges[se].two;
15285 group__shared_edge_to_vertex.
MakeJ();
15287 for (
int se = 0; se < shared_edges.
Size(); se++)
15289 const int glob_edge_id = shared_edges[se].one;
15290 const int group_id = shared_edges[se].two;
15291 const int *v = edge_to_vertex.
GetRow(glob_edge_id);
15292 for (
int i = 0; i < 2; i++)
15294 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(v[i]);
15295 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
15296 group__shared_edge_to_vertex.
AddConnection(group_id, loc_vertex_id);
15299 group__shared_edge_to_vertex.
ShiftUpI();
15306 Table &group__shared_tria_to_vertex =
15308 Table &group__shared_quad_to_vertex =
15311 group__shared_tria_to_vertex.
MakeI(num_groups);
15312 group__shared_quad_to_vertex.
MakeI(num_groups);
15313 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
15315 const int glob_face_id = shared_faces[sf].one;
15316 const int group_id = shared_faces[sf].two;
15321 group__shared_tria_to_vertex.
MakeJ();
15322 group__shared_quad_to_vertex.
MakeJ();
15323 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
15325 const int glob_face_id = shared_faces[sf].one;
15326 const int group_id = shared_faces[sf].two;
15362 for (
int i = 0; i < vertex_ids.
Size(); i++)
15364 const int glob_id = vertex_ids[i];
15365 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
15366 MFEM_ASSERT(loc_id >= 0,
"internal error");
15367 vertex_ids[i] = loc_id;
15370 AddConnections(group_id, vertex_ids, vertex_ids.
Size());
15372 group__shared_tria_to_vertex.
ShiftUpI();
15373 group__shared_quad_to_vertex.
ShiftUpI();
15377std::unique_ptr<FiniteElementSpace>
15385 return std::unique_ptr<FiniteElementSpace>(
15387 global_fespace.
FEColl(),
15392std::unique_ptr<GridFunction>
15397 std::unique_ptr<GridFunction> local_gf(
new GridFunction(&local_fespace));
15405 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
15407 const int glob_elem_id = elem_list[loc_elem_id];
15414 local_gf->SetSubVector(lvdofs, loc_vals);
15428 "mixed meshes are not supported!");
15429 MFEM_ASSERT(
mesh->
GetNodes(),
"meshes without nodes are not supported!");
15442 Compute(
nodes, d_mt);
15452 const int vdim = fespace->
GetVDim();
15453 const int NE = fespace->
GetNE();
15454 const int ND = fe->
GetDof();
15457 unsigned eval_flags = 0;
15459 Device::GetDeviceMemoryType();
15482 qi->DisableTensorProducts(!use_tensor_products);
15490 Vector Enodes(vdim*ND*NE, my_d_mt);
15491 elem_restr->Mult(
nodes, Enodes);
15492 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
15512 const int vdim = fespace->
GetVDim();
15528 unsigned eval_flags = 0;
15575 V(1) = s * ((ip.
y + layer) / n);
15580 V(2) = s * ((ip.
z + layer) / n);
15589 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
15593 int nvy = (closed) ? (ny) : (ny + 1);
15594 int nvt = mesh->
GetNV() * nvy;
15603 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
15608 for (
int i = 0; i < mesh->
GetNV(); i++)
15611 for (
int j = 0; j < nvy; j++)
15613 vc[1] = sy * (
real_t(j) / ny);
15619 for (
int i = 0; i < mesh->
GetNE(); i++)
15624 for (
int j = 0; j < ny; j++)
15627 qv[0] = vert[0] * nvy + j;
15628 qv[1] = vert[1] * nvy + j;
15629 qv[2] = vert[1] * nvy + (j + 1) % nvy;
15630 qv[3] = vert[0] * nvy + (j + 1) % nvy;
15636 for (
int i = 0; i < mesh->
GetNBE(); i++)
15641 for (
int j = 0; j < ny; j++)
15644 sv[0] = vert[0] * nvy + j;
15645 sv[1] = vert[0] * nvy + (j + 1) % nvy;
15661 for (
int i = 0; i < mesh->
GetNE(); i++)
15667 sv[0] = vert[0] * nvy;
15668 sv[1] = vert[1] * nvy;
15672 sv[0] = vert[1] * nvy + ny;
15673 sv[1] = vert[0] * nvy + ny;
15689 string cname = name;
15690 if (cname ==
"Linear")
15694 else if (cname ==
"Quadratic")
15698 else if (cname ==
"Cubic")
15702 else if (!strncmp(name,
"H1_", 3))
15706 else if (!strncmp(name,
"L2_T", 4))
15710 else if (!strncmp(name,
"L2_", 3))
15717 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
15729 for (
int i = 0; i < mesh->
GetNE(); i++)
15732 for (
int j = ny-1; j >= 0; j--)
15749 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
15754 int nvt = mesh->
GetNV() * nvz;
15759 bool wdgMesh =
false;
15760 bool hexMesh =
false;
15764 for (
int i = 0; i < mesh->
GetNV(); i++)
15768 for (
int j = 0; j < nvz; j++)
15770 vc[2] = sz * (
real_t(j) / nz);
15776 for (
int i = 0; i < mesh->
GetNE(); i++)
15786 for (
int j = 0; j < nz; j++)
15789 pv[0] = vert[0] * nvz + j;
15790 pv[1] = vert[1] * nvz + j;
15791 pv[2] = vert[2] * nvz + j;
15792 pv[3] = vert[0] * nvz + (j + 1) % nvz;
15793 pv[4] = vert[1] * nvz + (j + 1) % nvz;
15794 pv[5] = vert[2] * nvz + (j + 1) % nvz;
15801 for (
int j = 0; j < nz; j++)
15804 hv[0] = vert[0] * nvz + j;
15805 hv[1] = vert[1] * nvz + j;
15806 hv[2] = vert[2] * nvz + j;
15807 hv[3] = vert[3] * nvz + j;
15808 hv[4] = vert[0] * nvz + (j + 1) % nvz;
15809 hv[5] = vert[1] * nvz + (j + 1) % nvz;
15810 hv[6] = vert[2] * nvz + (j + 1) % nvz;
15811 hv[7] = vert[3] * nvz + (j + 1) % nvz;
15813 mesh3d->
AddHex(hv, attr);
15817 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
15818 << geom <<
"\'" << endl;
15824 for (
int i = 0; i < mesh->
GetNBE(); i++)
15829 for (
int j = 0; j < nz; j++)
15832 qv[0] = vert[0] * nvz + j;
15833 qv[1] = vert[1] * nvz + j;
15834 qv[2] = vert[1] * nvz + (j + 1) % nvz;
15835 qv[3] = vert[0] * nvz + (j + 1) % nvz;
15844 for (
int i = 0; i < mesh->
GetNE(); i++)
15855 tv[0] = vert[0] * nvz;
15856 tv[1] = vert[2] * nvz;
15857 tv[2] = vert[1] * nvz;
15861 tv[0] = vert[0] * nvz + nz;
15862 tv[1] = vert[1] * nvz + nz;
15863 tv[2] = vert[2] * nvz + nz;
15871 qv[0] = vert[0] * nvz;
15872 qv[1] = vert[3] * nvz;
15873 qv[2] = vert[2] * nvz;
15874 qv[3] = vert[1] * nvz;
15878 qv[0] = vert[0] * nvz + nz;
15879 qv[1] = vert[1] * nvz + nz;
15880 qv[2] = vert[2] * nvz + nz;
15881 qv[3] = vert[3] * nvz + nz;
15887 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
15888 << geom <<
"\'" << endl;
15894 if ( hexMesh && wdgMesh )
15898 else if ( hexMesh )
15902 else if ( wdgMesh )
15915 string cname = name;
15916 if (cname ==
"Linear")
15920 else if (cname ==
"Quadratic")
15924 else if (cname ==
"Cubic")
15928 else if (!strncmp(name,
"H1_", 3))
15932 else if (!strncmp(name,
"L2_T", 4))
15936 else if (!strncmp(name,
"L2_", 3))
15943 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
15955 for (
int i = 0; i < mesh->
GetNE(); i++)
15958 for (
int j = nz-1; j >= 0; j--)
15974 MFEM_VERIFY(
dim > 1,
"Not implemented for 1D meshes.");
15977 auto factor3 = [](
int N)
15979 for (
int i =
static_cast<int>(round(cbrt(N))); i > 0; i--)
15980 {
if (N % i == 0) {
return i; } }
15985 auto factor2 = [](
int N)
15987 for (
int i =
static_cast<int>(round(sqrt(N))); i > 0; i--)
15988 {
if (N % i == 0) {
return i; } }
15993 const int ref_factor = (
dim == 2) ? 4 : 8;
15997 int el0 = elem_per_mpi;
15998 while (el0 % ref_factor == 0)
16007 int mpi_x, mpi_y, mpi_z;
16008 int el0_x, el0_y, el0_z;
16011 mpi_x = factor2(mpi_cnt);
16012 mpi_y = mpi_cnt / mpi_x;
16015 el0_y = factor2(el0);
16016 el0_x = el0 / el0_y;
16020 mpi_x = factor3(mpi_cnt);
16021 mpi_y = factor2(mpi_cnt / mpi_x);
16022 mpi_z = mpi_cnt / mpi_x / mpi_y;
16025 el0_z = factor3(el0);
16026 el0_y = factor2(el0 / el0_z);
16027 el0_x = el0 / el0_y / el0_z;
16030 if (print &&
dim == 2)
16032 int elem_par_x = mpi_x * el0_x * pow(2, par_ref),
16033 elem_par_y = mpi_y * el0_y * pow(2, par_ref);
16035 mfem::out <<
"--- Mesh generation: \n";
16036 mfem::out <<
"Par mesh: " << elem_par_x <<
" x " << elem_par_y
16037 <<
" (" << elem_par_x * elem_par_y <<
" elements)\n"
16039 << el0_x * pow(2, par_ref) <<
" x "
16040 << el0_y * pow(2, par_ref)
16041 <<
" (" << el0_x * pow(2, 2*par_ref) * el0_y <<
" elements)\n"
16042 <<
"MPI blocks: " << mpi_x <<
" x " << mpi_y
16043 <<
" (" << mpi_x * mpi_y <<
" mpi tasks)\n" <<
"-\n"
16045 << mpi_x * el0_x <<
" x " << mpi_y * el0_y
16046 <<
" (" << mpi_x * el0_x * mpi_y * el0_y <<
" elements)\n"
16047 <<
"Elem / task: " << el0_x <<
" x " << el0_y << std::endl
16048 <<
"Par refine: " << par_ref << std::endl;
16052 if (print &&
dim == 3)
16054 int elem_par_x = mpi_x * el0_x * pow(2, par_ref),
16055 elem_par_y = mpi_y * el0_y * pow(2, par_ref),
16056 elem_par_z = mpi_z * el0_z * pow(2, par_ref);
16058 mfem::out <<
"--- Mesh generation: \n";
16060 << elem_par_x <<
" x " << elem_par_y <<
" x " << elem_par_z
16061 <<
" (" << elem_par_x*elem_par_y*elem_par_z <<
" elements)\n"
16063 << el0_x * pow(2, par_ref) <<
" x "
16064 << el0_y * pow(2, par_ref) <<
" x "
16065 << el0_z * pow(2, par_ref)
16066 <<
" (" << el0_x*pow(2, 3*par_ref)*el0_y*el0_z <<
" elements)\n"
16067 <<
"MPI blocks: " << mpi_x <<
" x " << mpi_y <<
" x " << mpi_z
16068 <<
" (" << mpi_x * mpi_y * mpi_z <<
" mpi tasks)\n" <<
"-\n"
16070 << mpi_x*el0_x <<
" x " << mpi_y*el0_y <<
" x " << mpi_z*el0_z
16071 <<
" (" << mpi_x*el0_x*mpi_y*el0_y*mpi_z*el0_z <<
" elements)\n"
16073 << el0_x <<
" x " << el0_y <<
" x " << el0_z << std::endl
16074 <<
"Par refine: " << par_ref << std::endl;
16084 nxyz[0] = mpi_x; nxyz[1] = mpi_y;
16091 nxyz[0] = mpi_x; nxyz[1] = mpi_y; nxyz[2] = mpi_z;
16094 const int NE = mesh.
GetNE();
16097 std::copy(p_raw.get(), p_raw.get() + NE, partitioning.
GetData());
16124 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
16125 <<
" 0 0 " << i <<
" -1 0\n";
16132 real_t mid[3] = {0, 0, 0};
16133 for (
int j = 0; j < 2; j++)
16135 for (
int k = 0; k <
spaceDim; k++)
16141 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
16142 << 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 Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
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 * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
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 * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
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 if the named attribute set is present.
ArraysByName< int > attr_sets
Named sets of attributes.
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...
@ GaussLobatto
Closed type.
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).
void GetColumnReference(int c, Vector &col)
real_t * Data() const
Returns the matrix data array. Warning: this method casts away constness.
real_t * GetData() const
Returns the matrix data array. Warning: this method casts away constness.
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.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
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()
int VectorDim() const
Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying fes.
virtual void ProjectCoefficient(Coefficient &coeff, ProjectType type=ProjectType::DEFAULT)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values at the vertices of element i for the 1-based 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 the "id" of the item whose parents are p1, p2, this "id" corresponding to the index of the item i...
int FindId(int p1, int p2) const
Find the "id" of an item whose parents are p1, p2. Return -1 if it does not 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
void Set3(const real_t x1, const real_t x2, const real_t x3)
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.
int GetPatchBdrAttribute(int i) const
Return the attribute of patch boundary element i, for a NURBS mesh.
int GetElementToEdgeTable(Table &)
void LoadNonconformingPatchTopo(std::istream &input, Array< int > &edge_to_ukv)
Read NURBS patch/macro-element mesh (MFEM NURBS NC-patch mesh format)
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
Table * GetEdgeFaceTable() const
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).
const Array< int > & GetFaceIndices(FaceType ftype) const
Map from boundary or interior face indices to mesh face indices.
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.
std::unordered_map< int, int > inv_face_indices[2]
cache for FaceIndices(ftype)
const Table & ElementToEdgeTable() const
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.
Array< int > face_indices[2]
cache for FaceIndices(ftype)
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)
Array< int > bdr_face_attrs_cache
internal cache for boundary element attributes
Geometry::Type GetElementGeometry(int i) const
Geometry::Type GetBdrElementGeometry(int i) const
void MakeHigherOrderSimplicial_(const Mesh &orig_mesh, const Array< int > &parent_elements)
Helper function for constructing higher order nodes from a mesh transformed into simplices....
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)
void ComputeFaceInfo(FaceType ftype) const
compute face_indices[ftype] and inv_face_indices[type]
IsoparametricTransformation FaceTransformation
Array< NCFaceInfo > nc_faces_info
Array< int > MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
Internal helper user in MakeSimplicial (and ParMesh::MakeSimplicial). Optional return is used in asse...
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 (...)
void CorrectPatchTopoOrientations(Array< int > &edge_to_ukv) const
Set signs to ensure knotvectors are pointed in the same direction.
bool Nonconforming() const
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".
std::unique_ptr< GridFunction > GetJacobianDeterminantGF() const
Create a GridFunction representing the Jacobian determinant.
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.
bool IsMixedMesh() const
Returns true if the mesh is a mixed mesh, false otherwise.
const Array< int > & GetElementAttributes() const
Returns the attributes for all elements in this mesh. The i'th entry of the array is the attribute of...
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
Print the mesh to the given stream using the default MFEM mesh format.
void RefineNURBS(bool usingKVF, real_t tol, const Array< int > &rf, const std::string &kvf)
Refine the NURBS mesh with default refinement factors in rf for each dimension.
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.
void Transform(std::function< void(const Vector &, Vector &)> f)
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)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1, int pyr_type=1)
Set the curvature of the mesh nodes using the given polynomial degree.
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 SaveVTKHDF(const std::string &fname, bool high_order=true)
Save the Mesh in VTKHDF format.
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
Return the indices of the elements sharing face Face.
void Printer(std::ostream &os=mfem::out, std::string section_delimiter="", const std::string &comments="") const
If NURBS mesh, write NURBS format. If NCMesh, write mfem v1.1 format. If section_delimiter is empty,...
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 GetNURBSPatches(Array< NURBSPatch * > &patches)
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 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.
virtual void SetAttributes(bool elem_attrs_changed=true, bool bdr_face_attrs_changed=true)
Determine the sets of unique attribute values in domain if elem_attrs_changed and boundary elements i...
void SetNode(int i, const real_t *coord)
void GetNodes(Vector &node_coord) const
const Array< int > & GetBdrFaceAttributes() const
Returns the attributes for all boundary elements in this mesh.
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)
const std::unordered_map< int, int > & GetInvFaceIndices(FaceType ftype) const
Inverse of the map FaceIndices(ftype)
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 GetEdgeToUniqueKnotvector(Array< int > &edge_to_ukv, Array< int > &ukv_to_rpkv) const
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....
Array< int > elem_attrs_cache
internal cache for element attributes
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 PrintTopoEdges(std::ostream &out, const Array< int > &e_to_k, bool vmap=false) const
Write the patch topology edges of a NURBS mesh (see PrintTopo()).
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)
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 ReadNURBSMesh(std::istream &input, int &curved, int &read_gf, bool spacing=false, bool nc=false)
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.
std::unique_ptr< L2_SegmentElement > EdgeTransfElement
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.
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 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 ReadGmshMesh(std::istream &input)
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones.
virtual void RefineNURBSWithKVFactors(int rf, const std::string &kvf)
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_ukv)
Read NURBS patch/macro-element mesh.
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()
void UpdateJacobianDeterminantGF(GridFunction &detgf) const
Update Jacobian determinant values in a given gridfunction.
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 MoveVertices(const Vector &displacements)
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)
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.
void Print(std::ostream &out, const std::string &comments="", bool nurbs=false) const
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)
Array< int > vertex_nodeId
vertex-index to node-id map, see UpdateVertices
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.
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 GetPatches(Array< NURBSPatch * > &patches)
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.
virtual void ReadCoarsePatchCP(std::istream &input)
Read the control points for coarse patches.
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 ...
int GetPatchSpaceDimension() const
Return the physical dimension of the NURBS geometry.
void Coarsen(int cf=2, real_t tol=1.0e-12)
Coarsen with optional coarsening factor cf.
int GetPatchAttribute(int i) const
Get the attribute for patch i, which is set to all elements in the patch.
void SetCoordsFromPatches(Vector &Nodes, int vdim)
Set FE coordinates in Nodes, using data from patches, with physical vector dimension vdim,...
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.
bool NonconformingPatches() const
Return true if the patch topology mesh is nonconforming.
virtual void RefineWithKVFactors(int rf, const std::string &kvf_filename, bool coarsened)
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.
virtual void PrintCoarsePatches(std::ostream &os)
Print control points for coarse patches.
void FullyCoarsen()
Fully coarsen all structured patches, for non-nested refinement of a mesh with a nonconforming patch ...
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.
Table stores the connectivity of elements of TYPE I to elements of TYPE II. For example,...
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)
int Push(int i, int j)
Establish connection between element i and element j in the table.
void AddConnection(int r, int c)
void Finalize()
Finalize the table initialization.
int Size() const
Returns the number of TYPE I elements.
int Size_of_connections() const
Returns the number of connections in the table.
void AddColumnsInRow(int r, int ncol)
void AddAColumnInRow(int r)
void SetDims(int rows, int nnz)
Set the rows and the number of all connections for the table.
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.
Low-level class for writing VTKHDF data (for use in ParaView).
void SaveMesh(const Mesh &mesh, bool high_order=true, int ref=-1)
Save the mesh, appending as a new time step.
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.
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
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
Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print, int &par_ref, Array< int > &partitioning)
Constructs the smallest possible [0,1]^dim serial mesh that can be used later to obtain a ParMesh wit...
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_HOST_DEVICE int FlipIndexSign(int i)
Signed indices i -> -1 - i are used as a convention to encode orientation.
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_HOST_DEVICE int UnsignIndex(int i)
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
void Swap(T &a, T &b)
Swap objects of type T. The operation is performed using the most specialized swap function from the ...
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
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
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.
constexpr real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
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