37#include <unordered_map>
38#include <unordered_set>
42#if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
47#if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
52 int*,
int*,
int*,
int*,
int*,
idxtype*);
54 int*,
int*,
int*,
int*,
int*,
idxtype*);
56 int*,
int*,
int*,
int*,
int*,
idxtype*);
121 return sqrt((d_hat * d_hat) / (dir * dir));
160 if (coord[d] < min(d)) { min(d) = coord[d]; }
161 if (coord[d] > max(d)) { max(d) = coord[d]; }
167 const bool use_boundary =
false;
176 for (
int i = 0; i < ne; i++)
193 for (
int j = 0; j < pointmat.
Width(); j++)
195 for (
int d = 0; d < pointmat.
Height(); d++)
197 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
198 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
220 h_max = kappa_max = -h_min;
221 if (
dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
229 if (Vh) { (*Vh)(i) = h; }
230 if (Vk) { (*Vk)(i) =
kappa; }
232 if (h < h_min) { h_min = h; }
233 if (h > h_max) { h_max = h; }
247 if (!num_elems_by_geom[g]) {
continue; }
248 if (!first) { os <<
" + "; }
256 real_t h_min, h_max, kappa_min, kappa_max;
258 os <<
"Mesh Characteristics:";
263 num_elems_by_geom = 0;
264 for (
int i = 0; i <
GetNE(); i++)
275 <<
"Number of vertices : " <<
GetNV() <<
'\n'
276 <<
"Number of elements : " <<
GetNE() <<
'\n'
277 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n';
282 <<
"Number of vertices : " <<
GetNV() <<
'\n'
283 <<
"Number of elements : " <<
GetNE() <<
'\n'
284 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
285 <<
"h_min : " << h_min <<
'\n'
286 <<
"h_max : " << h_max <<
'\n';
291 <<
"Number of vertices : " <<
GetNV() <<
'\n'
292 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
293 <<
"Number of elements : " <<
GetNE() <<
" -- ";
296 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
298 <<
"h_min : " << h_min <<
'\n'
299 <<
"h_max : " << h_max <<
'\n'
300 <<
"kappa_min : " << kappa_min <<
'\n'
301 <<
"kappa_max : " << kappa_max <<
'\n';
306 num_bdr_elems_by_geom = 0;
307 for (
int i = 0; i <
GetNBE(); i++)
312 num_faces_by_geom = 0;
319 <<
"Number of vertices : " <<
GetNV() <<
'\n'
320 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
321 <<
"Number of faces : " <<
GetNFaces() <<
" -- ";
324 <<
"Number of elements : " <<
GetNE() <<
" -- ";
327 <<
"Number of bdr elem : " <<
GetNBE() <<
" -- ";
331 <<
"h_min : " << h_min <<
'\n'
332 <<
"h_max : " << h_max <<
'\n'
333 <<
"kappa_min : " << kappa_min <<
'\n'
334 <<
"kappa_max : " << kappa_max <<
'\n';
336 os <<
'\n' << std::flush;
352 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
355 MFEM_ABORT(
"Unknown element type");
384 for (
int j = 0; j < n; j++)
386 pm(k,j) =
nodes(vdofs[n*k+j]);
435 int nv =
elements[i]->GetNVertices();
436 const int *v =
elements[i]->GetVertices();
441 for (
int j = 0; j < nv; j++)
443 pm(k, j) =
nodes(k*n+v[j]);
457 for (
int j = 0; j < n; j++)
459 pm(k,j) =
nodes(vdofs[n*k+j]);
493 for (
int j = 0; j < n; j++)
495 int idx = vdofs[n*k+j];
496 pm(k,j) =
nodes((idx<0)? -1-idx:idx);
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);
672 MFEM_ABORT(
"Not implemented.");
712 for (
int j = 0; j < 2; j++)
714 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
715 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
732 for (
int j = 0; j < 2; j++)
734 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
735 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
754 for (
int j = 0; j < 3; j++)
757 locpm(0, j) = vert.
x;
758 locpm(1, j) = vert.
y;
759 locpm(2, j) = vert.
z;
771 MFEM_VERIFY(i < 128,
"Local face index " << i/64
772 <<
" is not a triangular face of a wedge.");
780 for (
int j = 0; j < 3; j++)
783 locpm(0, j) = vert.
x;
784 locpm(1, j) = vert.
y;
785 locpm(2, j) = vert.
z;
796 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
797 <<
" is not a triangular face of a pyramid.");
805 for (
int j = 0; j < 3; j++)
808 locpm(0, j) = vert.
x;
809 locpm(1, j) = vert.
y;
810 locpm(2, j) = vert.
z;
827 for (
int j = 0; j < 4; j++)
830 locpm(0, j) = vert.
x;
831 locpm(1, j) = vert.
y;
832 locpm(2, j) = vert.
z;
844 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
845 <<
" is not a quadrilateral face of a wedge.");
851 for (
int j = 0; j < 4; j++)
854 locpm(0, j) = vert.
x;
855 locpm(1, j) = vert.
y;
856 locpm(2, j) = vert.
z;
867 MFEM_VERIFY(i < 64,
"Local face index " << i/64
868 <<
" is not a quadrilateral face of a pyramid.");
874 for (
int j = 0; j < 4; j++)
877 locpm(0, j) = vert.
x;
878 locpm(1, j) = vert.
y;
879 locpm(2, j) = vert.
z;
929 std::unordered_map<int, int> f_to_be;
930 for (
int i = 0; i <
GetNBE(); ++i)
940 for (
int f = 0;
f < nf; ++
f)
947 auto iter = f_to_be.find(
f);
948 if (iter != f_to_be.end())
950 const int be = iter->second;
972 for (
int i = 0; i <
GetNE(); ++i)
976 "Negative attribute on element " << i);
988 ifidcs.reserve(fidcs.Size());
1012const std::unordered_map<int, int> &
1075 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
1076 "face type " << face_type
1077 <<
" and element type " << elem_type <<
"\n");
1096 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
1097 "face type " << face_type
1098 <<
" and element type " << elem_type <<
"\n");
1130 FElTr.
Elem1 = &ElTr1;
1143 { MFEM_ABORT(
"NURBS mesh not supported!"); }
1146 FElTr.
Elem2 = &ElTr2;
1194 mfem::out <<
"\nInternal error: face id = " << FaceNo
1195 <<
", dist = " << dist <<
'\n';
1197 MFEM_ABORT(
"internal error");
1255 const FaceInfo &fi,
bool is_ghost)
const
1257#ifdef MFEM_THREAD_SAFE
1262 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1274 std::swap(composition(0,0), composition(0,1));
1275 std::swap(composition(1,0), composition(1,1));
1332 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1421 res.
Elem1No = element[0].index;
1422 res.Elem2No = element[1].index;
1423 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1424 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1425 res.NCFace = ncface;
1428 res.Elem1No = element[0].index;
1429 res.Elem2No = element[1].index;
1430 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1431 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1432 res.NCFace = ncface;
1435 res.Elem1No = element[0].index;
1436 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1439 res.Elem1No = element[0].index;
1440 res.Elem2No = -1 - element[1].index;
1441 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1442 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1445 res.Elem1No = element[0].index;
1446 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1449 res.Elem1No = element[0].index;
1450 res.Elem2No = -1 - element[1].index;
1451 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1452 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1457 res.Elem1No = element[0].index;
1458 res.Elem2No = -1 - element[1].index;
1459 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1460 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1468 os <<
"face topology=";
1478 os <<
"Non-conforming";
1485 os <<
"element[0].location=";
1499 os <<
"element[1].location=";
1513 os <<
"element[0].conformity=";
1530 os <<
"element[1].conformity=";
1547 os <<
"element[0].index=" << info.
element[0].
index <<
'\n'
1548 <<
"element[1].index=" << info.
element[1].
index <<
'\n'
1553 <<
"ncface=" << info.
ncface << std::endl;
1585 return faces[Face]->GetGeometryType();
1588 const int nc_face_id =
faces_info[Face].NCFace;
1590 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1651 "Could not determine a typical element Geometry!");
1660 face_marker.
SetSize(num_faces);
1662 for (
int f = 0;
f < num_faces;
f++)
1679 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1680 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1682 Array<bool> interior_bdr(max_bdr_attr); interior_bdr =
false;
1683 Array<bool> exterior_bdr(max_bdr_attr); exterior_bdr =
false;
1687 for (
int be = 0; be <
boundary.Size(); be++)
1689 const int bea =
boundary[be]->GetAttribute();
1691 if (bdr_marker[bea-1] != 0)
1697 interior_bdr[bea-1] =
true;
1701 exterior_bdr[bea-1] =
true;
1708 for (
int b = 0;
b < max_bdr_attr;
b++)
1710 if (bdr_marker[
b] != 0 && interior_bdr[
b])
1712 if (!excl || !exterior_bdr[
b])
1726 "Named set is not defined in this mesh!");
1728 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1732 for (
int b = 0;
b < max_bdr_attr;
b++)
1745 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1746 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1748 Array<bool> interior_bdr(max_bdr_attr); interior_bdr =
false;
1749 Array<bool> exterior_bdr(max_bdr_attr); exterior_bdr =
false;
1753 for (
int be = 0; be <
boundary.Size(); be++)
1755 const int bea =
boundary[be]->GetAttribute();
1761 interior_bdr[bea-1] =
true;
1765 exterior_bdr[bea-1] =
true;
1771 for (
int b = 0;
b < max_bdr_attr;
b++)
1773 if (bdr_marker[
b] == 0 && exterior_bdr[
b])
1775 if (!excl || !interior_bdr[
b])
1789 "Named set is not defined in this mesh!");
1790 MFEM_VERIFY(bdr_marker.
Size() >= max_bdr_attr,
1791 "bdr_marker must be at least bdr_attriburtes.Max() in length");
1795 for (
int b = 0;
b < max_bdr_attr;
b++)
1873 for (
int i = 0; i <
faces.Size(); i++)
1901#ifdef MFEM_USE_MEMALLOC
1939 if (bdr_face_attrs_changed)
1944 std::set<int> attribs;
1945 for (
int i = 0; i <
GetNBE(); i++)
1955 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1959 if (elem_attrs_changed)
1972 MFEM_WARNING(
"Non-positive attributes in the domain!");
1995static void CheckEnlarge(
Array<T> &array,
int size)
1997 if (size >= array.Size()) { array.SetSize(size + 1); }
2020 "invalid 'coords' size: " << coords.
Size());
2033 for (
int j = 0; j < 3; j++)
2035 vi[j] = (vp1[j] + vp2[j]) * 0.5;
2044 for (
int i = 0; i < nverts; i++)
2047 for (
int j = 0; j <
dim; j++)
2101 int vi[4] = {v1, v2, v3, v4};
2108#ifdef MFEM_USE_MEMALLOC
2148int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
2153 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
2166 static const int hex_to_tet[6][4] =
2168 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
2169 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
2173 for (
int i = 0; i < 6; i++)
2175 for (
int j = 0; j < 4; j++)
2177 ti[j] = vi[hex_to_tet[i][j]];
2185 static const int hex_to_wdg[2][6] =
2187 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
2191 for (
int i = 0; i < 2; i++)
2193 for (
int j = 0; j < 6; j++)
2195 ti[j] = vi[hex_to_wdg[i][j]];
2203 static const int hex_to_pyr[6][5] =
2205 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
2206 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
2210 for (
int i = 0; i < 6; i++)
2212 for (
int j = 0; j < 5; j++)
2214 ti[j] = vi[hex_to_pyr[i][j]];
2223 static const int quad_to_tri[4][2] =
2225 {0, 1}, {1, 2}, {2, 3}, {3, 0}
2231 ti[2] = elem_center_index;
2232 for (
int i = 0; i < num_faces; i++)
2234 for (
int j = 0; j < 2; j++)
2236 ti[j] = vi[quad_to_tri[i][j]];
2245 static const int quad_faces[4][2] =
2247 {0, 1}, {1, 2}, {2, 3}, {3, 0}
2251 for (
int i = 0; i < 4; i++)
2260 real_t r = 0.25, s = 0.25;
2261 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2262 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2267 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2268 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2273 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2274 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2279 vnew[0] = px(0)*(1-r)*(1-s) + px(1)*(r)*(1-s) + px(2)*r*s + px(3)*(1-r)*s;
2280 vnew[1] = py(0)*(1-r)*(1-s) + py(1)*(r)*(1-s) + py(2)*r*s + py(3)*(1-r)*s;
2284 static const int quad_faces_new[4][2] =
2286 { 1, 0}, { 2, 1}, { 3, 2}, { 0, 3}
2290 for (
int i = 0; i < num_faces; i++)
2292 for (
int j = 0; j < 2; j++)
2294 ti[j] = vi[quad_faces[i][j]];
2295 ti[j+2] = vnew_index[quad_faces_new[i][j]];
2303 std::map<std::array<int, 4>,
int> &hex_face_verts,
2309 return std::array<int, 4> {v[0], v[1], v[2], v[3]};
2313 static const int hex_to_tet[6][4] =
2315 { 0, 1, 2, 3 }, { 1, 2, 6, 5 }, { 5, 4, 7, 6},
2316 { 0, 1, 5, 4 }, { 2, 3, 7, 6 }, { 0,3, 7, 4}
2324 static const int tet_face[4][2] =
2326 {0, 1}, {1, 2}, {3, 2}, {3, 0}
2329 for (
int i = 0; i < num_faces; i++)
2331 for (
int j = 0; j < 4; j++)
2333 flist[j] = vi[hex_to_tet[i][j]];
2335 int face_center_index;
2337 auto t = get4arraysorted(flist);
2338 auto it = hex_face_verts.find(t);
2339 if (it == hex_face_verts.end())
2342 flist.
Size(), 3) - 1;
2343 hex_face_verts.insert({t, face_center_index});
2347 face_center_index = it->second;
2350 fti[2] = face_center_index;
2351 fti[3] = elem_center_index;
2352 for (
int j = 0; j < 4; j++)
2354 for (
int k = 0; k < 2; k++)
2356 fti[k] = flist[tet_face[j][k]];
2381 MFEM_ASSERT(bdr_elems.
Size() == new_be_to_face.
Size(),
"wrong size");
2382 for (
int i = 0; i < bdr_elems.
Size(); i++)
2433 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
2436 for (
int i = 0; i < 2; i++)
2438 for (
int j = 0; j < 3; j++)
2440 ti[j] = vi[quad_to_tri[i][j]];
2476 for (
int i = 0, j = 0; i <
faces_info.Size(); i++)
2492 "incorrect number of vertices: preallocated: " <<
vertices.Size()
2495 "incorrect number of elements: preallocated: " <<
elements.Size()
2498 "incorrect number of boundary elements: preallocated: "
2532 bool fix_orientation)
2535 if (fix_orientation)
2565 GeckoProgress(
real_t limit) : limit(limit) { sw.
Start(); }
2566 bool quit()
const override {
return limit > 0 && sw.
UserTime() > limit; }
2569class GeckoVerboseProgress :
public GeckoProgress
2575 GeckoVerboseProgress(
real_t limit) : GeckoProgress(limit) {}
2577 void beginorder(
const Graph* graph, Float cost)
const override
2578 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2579 void endorder(
const Graph* graph, Float cost)
const override
2580 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2582 void beginiter(
const Graph* graph,
2583 uint iter, uint maxiter, uint window)
const override
2585 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
2586 << window << std::flush;
2588 void enditer(
const Graph* graph, Float mincost, Float cost)
const override
2589 {
mfem::out <<
", cost = " << cost << endl; }
2594 int iterations,
int window,
2595 int period,
int seed,
bool verbose,
2601 GeckoProgress progress(time_limit);
2602 GeckoVerboseProgress vprogress(time_limit);
2605 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2613 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2615 const int *neighid = my_el_to_el.
GetRow(elemid);
2616 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2618 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2623 graph.
order(&functional, iterations, window, period, seed,
2624 verbose ? &vprogress : &progress);
2630 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2633 return graph.
cost();
2645 : coord(coord), dir(dir), points(points), mid(mid) {}
2647 bool operator()(
int i)
const
2649 return (points[3*i + coord] < mid) != dir;
2653static void HilbertSort2D(
int coord1,
2656 const Array<real_t> &points,
int *beg,
int *end,
2659 if (end - beg <= 1) {
return; }
2661 real_t xmid = (xmin + xmax)*0.5;
2662 real_t ymid = (ymin + ymax)*0.5;
2664 int coord2 = (coord1 + 1) % 2;
2667 int *p0 = beg, *p4 = end;
2668 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2669 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2670 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2674 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2675 ymin, xmin, ymid, xmid);
2677 if (p1 != p0 || p2 != p4)
2679 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2680 xmin, ymid, xmid, ymax);
2682 if (p2 != p0 || p3 != p4)
2684 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2685 xmid, ymid, xmax, ymax);
2689 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2690 ymid, xmax, ymin, xmid);
2694static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2695 const Array<real_t> &points,
int *beg,
int *end,
2699 if (end - beg <= 1) {
return; }
2701 real_t xmid = (xmin + xmax)*0.5;
2702 real_t ymid = (ymin + ymax)*0.5;
2703 real_t zmid = (zmin + zmax)*0.5;
2705 int coord2 = (coord1 + 1) % 3;
2706 int coord3 = (coord1 + 2) % 3;
2709 int *p0 = beg, *p8 = end;
2710 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2711 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2712 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2713 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2714 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2715 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2716 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2720 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2721 zmin, xmin, ymin, zmid, xmid, ymid);
2723 if (p1 != p0 || p2 != p8)
2725 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2726 ymin, zmid, xmin, ymid, zmax, xmid);
2728 if (p2 != p0 || p3 != p8)
2730 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2731 ymid, zmid, xmin, ymax, zmax, xmid);
2733 if (p3 != p0 || p4 != p8)
2735 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2736 xmin, ymax, zmid, xmid, ymid, zmin);
2738 if (p4 != p0 || p5 != p8)
2740 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2741 xmid, ymax, zmid, xmax, ymid, zmin);
2743 if (p5 != p0 || p6 != p8)
2745 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2746 ymax, zmid, xmax, ymid, zmax, xmid);
2748 if (p6 != p0 || p7 != p8)
2750 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2751 ymid, zmid, xmax, ymin, zmax, xmid);
2755 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2756 zmid, xmax, ymin, zmin, xmid, ymid);
2770 if (
spaceDim < 3) { points = 0.0; }
2773 for (
int i = 0; i <
GetNE(); i++)
2778 points[3*i + j] = center(j);
2785 indices.
Sort([&](
int a,
int b)
2786 {
return points[3*
a] < points[3*
b]; });
2791 HilbertSort2D(0,
false,
false,
2792 points, indices.
begin(), indices.
end(),
2793 min(0), min(1), max(0), max(1));
2798 HilbertSort3D(0,
false,
false,
false,
2799 points, indices.
begin(), indices.
end(),
2800 min(0), min(1), min(2), max(0), max(1), max(2));
2805 for (
int i = 0; i <
GetNE(); i++)
2807 ordering[indices[i]] = i;
2816 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2821 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
2825 MFEM_VERIFY(ordering.
Size() ==
GetNE(),
"invalid reordering array.")
2858 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2862 old_elem_node_vals[old_elid] =
new Vector(vals);
2868 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2870 int new_elid = ordering[old_elid];
2871 new_elements[new_elid] =
elements[old_elid];
2876 if (reorder_vertices)
2881 vertex_ordering = -1;
2883 int new_vertex_ind = 0;
2884 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2886 int *elem_vert =
elements[new_elid]->GetVertices();
2887 int nv =
elements[new_elid]->GetNVertices();
2888 for (
int vi = 0; vi < nv; ++vi)
2890 int old_vertex_ind = elem_vert[vi];
2891 if (vertex_ordering[old_vertex_ind] == -1)
2893 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2894 new_vertices[new_vertex_ind] =
vertices[old_vertex_ind];
2904 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2906 int *elem_vert =
elements[new_elid]->GetVertices();
2907 int nv =
elements[new_elid]->GetNVertices();
2908 for (
int vi = 0; vi < nv; ++vi)
2910 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2915 for (
int belid = 0; belid <
GetNBE(); ++belid)
2917 int *be_vert =
boundary[belid]->GetVertices();
2918 int nv =
boundary[belid]->GetNVertices();
2919 for (
int vi = 0; vi < nv; ++vi)
2921 be_vert[vi] = vertex_ordering[be_vert[vi]];
2950 nodes_fes->
Update(
false);
2953 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2955 int new_elid = ordering[old_elid];
2958 delete old_elem_node_vals[old_elid];
3007 length_idx[j].one =
GetLength(i, it.Column());
3008 length_idx[j].two = j;
3017 order[length_idx[i].two] = i;
3032 elements[i]->MarkEdge(v_to_v, order);
3039 boundary[i]->MarkEdge(v_to_v, order);
3046 if (*old_v_to_v && *old_elem_vert)
3053 if (*old_v_to_v == NULL)
3055 bool need_v_to_v =
false;
3063 if (dofs.
Size() > 0)
3075 if (*old_elem_vert == NULL)
3077 bool need_elem_vert =
false;
3079 for (
int i = 0; i <
GetNE(); i++)
3085 if (dofs.
Size() > 1)
3087 need_elem_vert =
true;
3093 *old_elem_vert =
new Table;
3095 for (
int i = 0; i <
GetNE(); i++)
3097 (*old_elem_vert)->AddColumnsInRow(i,
elements[i]->GetNVertices());
3099 (*old_elem_vert)->MakeJ();
3100 for (
int i = 0; i <
GetNE(); i++)
3105 (*old_elem_vert)->ShiftUpI();
3118 const int num_edge_dofs = old_dofs.
Size();
3130 if (num_edge_dofs > 0)
3139 const int old_i = (*old_v_to_v)(i, it.Column());
3140 const int new_i = it.Index();
3141 if (new_i == old_i) {
continue; }
3143 old_dofs.
SetSize(num_edge_dofs);
3144 new_dofs.
SetSize(num_edge_dofs);
3145 for (
int j = 0; j < num_edge_dofs; j++)
3147 old_dofs[j] = offset + old_i * num_edge_dofs + j;
3148 new_dofs[j] = offset + new_i * num_edge_dofs + j;
3152 for (
int j = 0; j < old_dofs.
Size(); j++)
3154 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3166 Table old_face_vertex;
3172 old_face_vertex.
MakeJ();
3175 faces[i]->GetNVertices());
3187 const int *old_v = old_face_vertex.
GetRow(i);
3189 switch (old_face_vertex.
RowSize(i))
3192 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
3196 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
3200 new_fdofs[new_i+1] = old_dofs.
Size();
3207 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
3210 switch (old_face_vertex.
RowSize(i))
3213 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
3214 new_v =
faces[new_i]->GetVertices();
3220 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
3221 new_v =
faces[new_i]->GetVertices();
3229 for (
int j = 0; j < old_dofs.
Size(); j++)
3232 const int old_j = dof_ord[j];
3233 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
3237 for (
int j = 0; j < old_dofs.
Size(); j++)
3239 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3261 for (
int i = 0; i <
GetNE(); i++)
3263 const int *old_v = old_elem_vert->
GetRow(i);
3264 const int *new_v =
elements[i]->GetVertices();
3271 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
3285 <<
" FE collection) are not supported yet!");
3289 MFEM_VERIFY(dof_ord != NULL,
3290 "FE collection '" << fec->
Name()
3291 <<
"' does not define reordering for "
3295 for (
int j = 0; j < new_dofs.
Size(); j++)
3298 const int old_j = dof_ord[j];
3299 new_dofs[old_j] = offset + j;
3301 offset += new_dofs.
Size();
3304 for (
int j = 0; j < old_dofs.
Size(); j++)
3306 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
3344 MFEM_ASSERT(
NURBSext,
"SetPatchAttribute is only for NURBS meshes");
3347 for (
auto e : elems)
3355 MFEM_ASSERT(
NURBSext,
"GetPatchAttribute is only for NURBS meshes");
3361 MFEM_ASSERT(
NURBSext,
"SetPatchBdrAttribute is only for NURBS meshes");
3365 for (
auto be : bdryelems)
3373 MFEM_ASSERT(
NURBSext,
"GetBdrPatchBdrAttribute is only for NURBS meshes");
3379 MFEM_VERIFY(
NURBSext,
"Must be a NURBS mesh");
3414 if (generate_edges == 1)
3432 bool fix_orientation)
3449 if (generate_edges == 1)
3514 bool generate_edges =
true;
3523 MFEM_VERIFY(
ncmesh == NULL,
"");
3558 if (
Dim > 1 && generate_edges)
3622 const bool check_orientation =
true;
3623 const bool curved = (
Nodes != NULL);
3624 const bool may_change_topology =
3626 ( check_orientation && fix_orientation &&
3630 Table *old_elem_vert = NULL;
3632 if (curved && may_change_topology)
3637 if (check_orientation)
3647 if (may_change_topology)
3652 delete old_elem_vert;
3673 for (
int i = 0; i < num_faces; i++)
3676 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
3677 " Interior face with incompatible orientations.");
3688 int NVert, NElem, NBdrElem;
3690 NVert = (nx+1) * (ny+1) * (nz+1);
3691 NElem = nx * ny * nz;
3692 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3701 NBdrElem += 2*nx*ny;
3706 NVert += nx * ny * nz;
3709 InitMesh(3, 3, NVert, NElem, NBdrElem);
3715 for (z = 0; z <= nz; z++)
3717 coord[2] = ((
real_t) z / nz) * sz;
3718 for (y = 0; y <= ny; y++)
3720 coord[1] = ((
real_t) y / ny) * sy;
3721 for (x = 0; x <= nx; x++)
3723 coord[0] = ((
real_t) x / nx) * sx;
3730 for (z = 0; z < nz; z++)
3732 coord[2] = (((
real_t) z + 0.5) / nz) * sz;
3733 for (y = 0; y < ny; y++)
3735 coord[1] = (((
real_t) y + 0.5) / ny) * sy;
3736 for (x = 0; x < nx; x++)
3738 coord[0] = (((
real_t) x + 0.5) / nx) * sx;
3745#define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3746#define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3753 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3755 for (
int k = 0; k < nx*ny*nz; k++)
3762 ind[0] = VTX(x , y , z );
3763 ind[1] = VTX(x+1, y , z );
3764 ind[2] = VTX(x+1, y+1, z );
3765 ind[3] = VTX(x , y+1, z );
3766 ind[4] = VTX(x , y , z+1);
3767 ind[5] = VTX(x+1, y , z+1);
3768 ind[6] = VTX(x+1, y+1, z+1);
3769 ind[7] = VTX(x , y+1, z+1);
3777 for (z = 0; z < nz; z++)
3779 for (y = 0; y < ny; y++)
3781 for (x = 0; x < nx; x++)
3784 ind[0] = VTX(x , y , z );
3785 ind[1] = VTX(x+1, y , z );
3786 ind[2] = VTX(x+1, y+1, z );
3787 ind[3] = VTX(x , y+1, z );
3788 ind[4] = VTX(x , y , z+1);
3789 ind[5] = VTX(x+1, y , z+1);
3790 ind[6] = VTX(x+1, y+1, z+1);
3791 ind[7] = VTX( x, y+1, z+1);
3803 ind[8] = VTXP(x, y, z);
3817 for (y = 0; y < ny; y++)
3819 for (x = 0; x < nx; x++)
3822 ind[0] = VTX(x , y , 0);
3823 ind[1] = VTX(x , y+1, 0);
3824 ind[2] = VTX(x+1, y+1, 0);
3825 ind[3] = VTX(x+1, y , 0);
3842 for (y = 0; y < ny; y++)
3844 for (x = 0; x < nx; x++)
3847 ind[0] = VTX(x , y , nz);
3848 ind[1] = VTX(x+1, y , nz);
3849 ind[2] = VTX(x+1, y+1, nz);
3850 ind[3] = VTX(x , y+1, nz);
3867 for (z = 0; z < nz; z++)
3869 for (y = 0; y < ny; y++)
3872 ind[0] = VTX(0 , y , z );
3873 ind[1] = VTX(0 , y , z+1);
3874 ind[2] = VTX(0 , y+1, z+1);
3875 ind[3] = VTX(0 , y+1, z );
3888 for (z = 0; z < nz; z++)
3890 for (y = 0; y < ny; y++)
3893 ind[0] = VTX(nx, y , z );
3894 ind[1] = VTX(nx, y+1, z );
3895 ind[2] = VTX(nx, y+1, z+1);
3896 ind[3] = VTX(nx, y , z+1);
3909 for (x = 0; x < nx; x++)
3911 for (z = 0; z < nz; z++)
3914 ind[0] = VTX(x , 0, z );
3915 ind[1] = VTX(x+1, 0, z );
3916 ind[2] = VTX(x+1, 0, z+1);
3917 ind[3] = VTX(x , 0, z+1);
3930 for (x = 0; x < nx; x++)
3932 for (z = 0; z < nz; z++)
3935 ind[0] = VTX(x , ny, z );
3936 ind[1] = VTX(x , ny, z+1);
3937 ind[2] = VTX(x+1, ny, z+1);
3938 ind[3] = VTX(x+1, ny, z );
3955 ofstream test_stream(
"debug.mesh");
3957 test_stream.close();
3985 for (
real_t j = 0; j < ny+1; j++)
3987 real_t cy = (j / ny) * sy;
3988 for (
real_t i = 0; i < nx+1; i++)
3990 real_t cx = (i / nx) * sx;
3997 for (
int y = 0; y < ny; y++)
3999 for (
int x = 0; x < nx; x++)
4001 ind[0] = x + y*(nx+1);
4002 ind[1] = x + 1 +y*(nx+1);
4003 ind[2] = x + 1 + (y+1)*(nx+1);
4004 ind[3] = x + (y+1)*(nx+1);
4010 for (
int i = 0; i < nx; i++)
4016 for (
int j = 0; j < ny; j++)
4059 for (
real_t j = 0; j < ny+1; j++)
4061 real_t cy = (j / ny) * sy;
4062 for (
real_t i = 0; i < nx+1; i++)
4064 real_t cx = (i / nx) * sx;
4071 for (
int y = 0; y < ny; y++)
4073 for (
int x = 0; x < nx; x++)
4075 ind[0] = x + y*(nx+1);
4076 ind[1] = x + 1 +y*(nx+1);
4077 ind[2] = x + 1 + (y+1)*(nx+1);
4078 ind[3] = x + (y+1)*(nx+1);
4084 for (
int i = 0; i < nx; i++)
4090 for (
int j = 0; j < ny; j++)
4116 const int NVert = (nx+1) * (ny+1) * (nz+1);
4117 const int NElem = nx * ny * nz * 24;
4118 const int NBdrElem = 2*(nx*ny+nx*nz+ny*nz)*4;
4120 InitMesh(3, 3, NVert, NElem, NBdrElem);
4125 for (
real_t z = 0; z <= nz; z++)
4127 coord[2] = ( z / nz) * sz;
4128 for (
real_t y = 0; y <= ny; y++)
4130 coord[1] = (y / ny) * sy;
4131 for (
real_t x = 0; x <= nx; x++)
4133 coord[0] = (x / nx) * sx;
4139 std::map<std::array<int, 4>,
int> hex_face_verts;
4140 auto VertexIndex = [nx, ny](
int xc,
int yc,
int zc)
4142 return xc + (yc + zc*(ny+1))*(nx+1);
4146 for (
int z = 0; z < nz; z++)
4148 for (
int y = 0; y < ny; y++)
4150 for (
int x = 0; x < nx; x++)
4153 ind[0] = VertexIndex(x , y , z );
4154 ind[1] = VertexIndex(x+1, y , z );
4155 ind[2] = VertexIndex(x+1, y+1, z );
4156 ind[3] = VertexIndex(x , y+1, z );
4157 ind[4] = VertexIndex(x , y , z+1);
4158 ind[5] = VertexIndex(x+1, y , z+1);
4159 ind[6] = VertexIndex(x+1, y+1, z+1);
4160 ind[7] = VertexIndex( x, y+1, z+1);
4168 hex_face_verts.clear();
4177 std::map<std::array<int, 3>,
int> tet_face_count;
4179 std::map<std::array<int, 3>,
int> face_count_map;
4184 return std::array<int, 3> {v[0], v[1], v[2]};
4193 for (
int j = 0; j < el_faces.
Size(); j++)
4196 auto t = get3array(vertidxs);
4197 auto it = tet_face_count.find(t);
4198 if (it == tet_face_count.end())
4200 tet_face_count.insert({t, 1});
4201 face_count_map.insert({t, el_faces[j]});
4210 for (
const auto &edge : tet_face_count)
4212 if (edge.second == 1)
4214 int facenum = (face_count_map.find(edge.first))->second;
4221 ofstream test_stream(
"debug.mesh");
4223 test_stream.close();
4232 bool generate_edges,
bool sfc_ordering)
4256 for (j = 0; j < ny+1; j++)
4258 cy = ((
real_t) j / ny) * sy;
4259 for (i = 0; i < nx+1; i++)
4261 cx = ((
real_t) i / nx) * sx;
4273 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
4275 for (k = 0; k < nx*ny; k++)
4279 ind[0] = i + j*(nx+1);
4280 ind[1] = i + 1 +j*(nx+1);
4281 ind[2] = i + 1 + (j+1)*(nx+1);
4282 ind[3] = i + (j+1)*(nx+1);
4289 for (j = 0; j < ny; j++)
4291 for (i = 0; i < nx; i++)
4293 ind[0] = i + j*(nx+1);
4294 ind[1] = i + 1 +j*(nx+1);
4295 ind[2] = i + 1 + (j+1)*(nx+1);
4296 ind[3] = i + (j+1)*(nx+1);
4305 for (i = 0; i < nx; i++)
4311 for (j = 0; j < ny; j++)
4333 for (j = 0; j < ny+1; j++)
4335 cy = ((
real_t) j / ny) * sy;
4336 for (i = 0; i < nx+1; i++)
4338 cx = ((
real_t) i / nx) * sx;
4347 for (j = 0; j < ny; j++)
4349 for (i = 0; i < nx; i++)
4351 ind[0] = i + j*(nx+1);
4352 ind[1] = i + 1 + (j+1)*(nx+1);
4353 ind[2] = i + (j+1)*(nx+1);
4356 ind[1] = i + 1 + j*(nx+1);
4357 ind[2] = i + 1 + (j+1)*(nx+1);
4365 for (i = 0; i < nx; i++)
4371 for (j = 0; j < ny; j++)
4381 MFEM_ABORT(
"Unsupported element type.");
4387 if (generate_edges == 1)
4425 for (j = 0; j < n+1; j++)
4431 for (j = 0; j < n; j++)
4458 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4510 for (
int i = 0; i <
faces.Size(); i++)
4552 if (
dynamic_cast<const ParMesh*
>(&mesh))
4564 if (mesh.
Nodes && copy_nodes)
4600 int refine,
bool fix_orientation)
4604 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
4605 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
4622 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
4632 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
4668 ref_factors = ref_factor;
4681Mesh::Mesh(
const std::string &filename,
int generate_edges,
int refine,
4682 bool fix_orientation)
4683 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4692 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
4696 Load(imesh, generate_edges, refine, fix_orientation);
4701 bool fix_orientation)
4702 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4705 Load(input, generate_edges, refine, fix_orientation);
4714 "Not enough vertices in external array : "
4715 "len_vertex_data = "<< len_vertex_data <<
", "
4720 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
4725 memcpy(vertex_data,
vertices.GetData(),
4734 int *element_attributes,
int num_elements,
4736 int *boundary_attributes,
int num_boundary_elements,
4738 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4740 if (space_dimension == -1)
4746 num_boundary_elements);
4749 int boundary_index_stride = num_boundary_elements > 0 ?
4753 vertices.MakeRef(
reinterpret_cast<Vertex*
>(vertices_), num_vertices);
4756 for (
int i = 0; i < num_elements; i++)
4759 elements[i]->SetVertices(element_indices + i * element_index_stride);
4760 elements[i]->SetAttribute(element_attributes[i]);
4764 for (
int i = 0; i < num_boundary_elements; i++)
4767 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
4768 boundary[i]->SetAttribute(boundary_attributes[i]);
4776 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4813 MFEM_ABORT(
"NURBS mesh has no patches.");
4827#ifdef MFEM_USE_MEMALLOC
4836 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
4849 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
4852 for (
int i = 0; i < nv; i++)
4865 for (
int j = 0; j < nv; j++)
4937 MFEM_ABORT(
"invalid element type: " << type);
4944 std::string parse_tag)
4946 int curved = 0, read_gf = 1;
4947 bool finalize_topo =
true;
4951 MFEM_ABORT(
"Input stream is not open");
4958 getline(input, mesh_type);
4962 int mfem_version = 0;
4963 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
4964 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
4965 else if (mesh_type ==
"MFEM mesh v1.3") { mfem_version = 13; }
4969 int mfem_nc_version = 0;
4970 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
4971 else if (mesh_type ==
"MFEM NC mesh v1.1") { mfem_nc_version = 11; }
4972 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
4980 if (mfem_version >= 12 && parse_tag.empty())
4982 parse_tag =
"mfem_mesh_end";
4986 else if (mfem_nc_version)
4988 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4995 MFEM_VERIFY(mfem_nc_version >= 10,
4996 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
4997 "used to load a parallel nonconforming mesh, sorry.");
5000 input, mfem_nc_version, curved, is_nc);
5005 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
5018 else if (mesh_type ==
"linemesh")
5022 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
5024 if (mesh_type ==
"curved_areamesh2")
5030 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
5034 else if (mesh_type ==
"TrueGrid")
5038 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
5040 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
5042 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
5043 "Unsupported VTK format");
5044 ReadVTKMesh(input, curved, read_gf, finalize_topo);
5046 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
5050 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
5054 else if (mesh_type ==
"MFEM NURBS NC-patch mesh v1.0")
5058 else if (mesh_type ==
"MFEM NURBS mesh v1.1")
5062 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
5067 else if (mesh_type ==
"$MeshFormat")
5072 ((mesh_type.size() > 2 &&
5073 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
5074 (mesh_type.size() > 3 &&
5075 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
5080#ifdef MFEM_USE_NETCDF
5083 MFEM_ABORT(
"NetCDF support requires configuration with"
5084 " MFEM_USE_NETCDF=YES");
5090 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
5091 " Use mfem::named_ifgzstream for input.");
5097 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
5123 bool generate_bdr =
false;
5128 if (curved && read_gf)
5142 if (mfem_version >= 12)
5148 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
5149 getline(input, line);
5155 if (line ==
"mfem_mesh_end") {
break; }
5157 while (line != parse_tag);
5159 else if (mfem_nc_version >= 10)
5164 MFEM_VERIFY(ident ==
"mfem_mesh_end",
5165 "invalid mesh: end of file tag not found");
5173 if (input.peek() ==
'p')
5176 MFEM_VERIFY(ident ==
"patch_cp",
"Invalid mesh format");
5185 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
5187 int i, j, ie, ib, iv, *v, nv;
5218 for (i = 0; i < num_pieces; i++)
5225 for (i = 0; i < num_pieces; i++)
5231 for (j = 0; j < m->
GetNE(); j++)
5236 for (j = 0; j < m->
GetNBE(); j++)
5241 for (
int k = 0; k < nv; k++)
5243 v[k] = lvert_vert[v[k]];
5248 for (j = 0; j < m->
GetNV(); j++)
5260 for (i = 0; i < num_pieces; i++)
5271 for (i = 0; i < num_pieces; i++)
5275 for (j = 0; j < m->
GetNE(); j++)
5280 for (
int k = 0; k < nv; k++)
5287 for (j = 0; j < m->
GetNBE(); j++)
5292 for (
int k = 0; k < nv; k++)
5299 for (j = 0; j < m->
GetNV(); j++)
5313 for (i = 0; i < num_pieces; i++)
5315 gf_array[i] = mesh_array[i]->
GetNodes();
5328 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
5331 ref_factors = ref_factor;
5342 int orig_ne = orig_mesh.
GetNE();
5343 MFEM_VERIFY(ref_factors.
Size() == orig_ne,
5344 "Number of refinement factors must equal number of elements")
5345 MFEM_VERIFY(orig_ne == 0 ||
5346 ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
5349 "Invalid refinement type. Must use closed basis type.");
5351 int min_ref = orig_ne > 0 ? ref_factors.
Min() : 1;
5352 int max_ref = orig_ne > 0 ? ref_factors.
Max() : 1;
5354 bool var_order = (min_ref != max_ref);
5367 for (
int i = 0; i < orig_ne; i++)
5383 for (
int el = 0; el < orig_ne; el++)
5395 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
5396 for (
int i = 0; i < phys_pts.
Width(); i++)
5405 for (
int k = 0; k < nvert; k++)
5408 v[k] = rdofs[c2h_map[cid]];
5421 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
5432 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
5438 for (
int k = 0; k < nvert; k++)
5441 v[k] = rdofs[c2h_map[cid]];
5463 for (
int iel = 0; iel < orig_ne; iel++)
5472 const int *node_map = NULL;
5475 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
5476 const int *vertex_map = vertex_fec.
GetDofMap(geom);
5477 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
5478 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
5481 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
5484 int iv = vertex_map[iv_lex];
5486 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
5488 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
5491 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
5502 using GeomRef = std::pair<Geometry::Type, int>;
5503 std::map<GeomRef, int> point_matrices_offsets;
5505 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5509 GeomRef id(geom, ref_factors[el_coarse]);
5510 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
5516 point_matrices_offsets[id] = n_point_matrices[geom];
5517 n_point_matrices[geom] += nref_el;
5524 int nmatrices = n_point_matrices[geom];
5531 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5534 int ref = ref_factors[el_coarse];
5535 int offset = point_matrices_offsets[GeomRef(geom, ref)];
5541 for (
int k = 0; k < nvert; k++)
5568 if (orig_mesh.
GetNodes() !=
nullptr)
5578 "Mesh::MakeSimplicial requires a properly oriented input mesh");
5580 "Mesh::MakeSimplicial does not support non-conforming meshes.")
5587 Mesh copy(orig_mesh);
5590 std::iota(parent_elements.
begin(), parent_elements.
end(), 0);
5591 return parent_elements;
5594 int nv = orig_mesh.
GetNV();
5595 int ne = orig_mesh.
GetNE();
5596 int nbe = orig_mesh.
GetNBE();
5609 int new_ne = 0, new_nbe = 0;
5610 for (
int i=0; i<ne; ++i)
5614 for (
int i=0; i<nbe; ++i)
5623 for (
int i=0; i<nv; ++i)
5633 if (vglobal ==
nullptr)
5636 std::iota(vglobal_id.
begin(), vglobal_id.
end(), 0);
5637 vglobal = vglobal_id.
GetData();
5641 constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
5642 constexpr int quad_ntris = 2;
5643 constexpr int prism_ntets = 3;
5647 static const int quad_trimap[2][nv_tri*quad_ntris] =
5659 static const int prism_rot[nv_prism*nv_prism] =
5668 static const int prism_f[nv_quad] = {1, 2, 5, 4};
5669 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
5683 static const int hex_rot[nv_hex*nv_hex] =
5685 0, 1, 2, 3, 4, 5, 6, 7,
5686 1, 0, 4, 5, 2, 3, 7, 6,
5687 2, 1, 5, 6, 3, 0, 4, 7,
5688 3, 0, 1, 2, 7, 4, 5, 6,
5689 4, 0, 3, 7, 5, 1, 2, 6,
5690 5, 1, 0, 4, 6, 2, 3, 7,
5691 6, 2, 1, 5, 7, 3, 0, 4,
5692 7, 3, 2, 6, 4, 0, 1, 5
5694 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
5695 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
5696 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
5697 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
5698 static const int hex_tetmap0[nv_tet*5] =
5705 static const int hex_tetmap1[nv_tet*6] =
5712 static const int hex_tetmap2[nv_tet*6] =
5719 static const int hex_tetmap3[nv_tet*6] =
5726 static const int *hex_tetmaps[4] =
5728 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
5731 auto find_min = [](
const int *
a,
int n) {
return std::min_element(
a,
a+n)-
a; };
5734 for (
int i=0; i<ne; ++i)
5736 const int *v = orig_mesh.
elements[i]->GetVertices();
5740 if (num_subdivisions[orig_geom] == 1)
5753 for (
int itri=0; itri<quad_ntris; ++itri)
5758 for (
int iv=0; iv<nv_tri; ++iv)
5760 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
5769 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
5772 int irot = find_min(vg, nv_prism);
5773 for (
int iv=0; iv<nv_prism; ++iv)
5775 int jv = prism_rot[iv + irot*nv_prism];
5780 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
5781 int j = find_min(q, nv_quad);
5782 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
5783 for (
int itet=0; itet<prism_ntets; ++itet)
5788 for (
int iv=0; iv<nv_tet; ++iv)
5790 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
5799 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
5803 int irot = find_min(vg, nv_hex);
5804 for (
int iv=0; iv<nv_hex; ++iv)
5806 int jv = hex_rot[iv + irot*nv_hex];
5816 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
5817 j = find_min(q, nv_quad);
5818 if (j == 0 || j == 2) { bitmask += 4; }
5820 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
5821 j = find_min(q, nv_quad);
5822 if (j == 1 || j == 3) { bitmask += 2; }
5824 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
5825 j = find_min(q, nv_quad);
5826 if (j == 0 || j == 2) { bitmask += 1; }
5829 int nrot = num_rot[bitmask];
5830 for (
int k=0; k<nrot; ++k)
5844 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
5845 int ntets = (ndiags == 0) ? 5 : 6;
5846 const int *tetmap = hex_tetmaps[ndiags];
5847 for (
int itet=0; itet<ntets; ++itet)
5852 for (
int iv=0; iv<nv_tet; ++iv)
5854 v2[iv] = vg[tetmap[itet + iv*ntets]];
5864 for (
int i=0; i<nbe; ++i)
5866 const int *v = orig_mesh.
boundary[i]->GetVertices();
5869 if (num_subdivisions[orig_geom] == 1)
5879 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
5881 int iv_min = find_min(vg, nv_quad);
5882 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
5883 for (
int itri=0; itri<quad_ntris; ++itri)
5888 for (
int iv=0; iv<nv_tri; ++iv)
5890 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
5897 MFEM_ABORT(
"Unreachable");
5908 return parent_elems;
5918 auto *orig_fespace = orig_mesh.
GetNodes()->FESpace();
5919 SetCurvature(orig_fespace->GetMaxElementOrder(), orig_fespace->IsDGSpace(),
5944 child_nodes_in_parent;
5945 for (
int i = 0; i < parent_elements.
Size(); i++)
5947 const int ip = parent_elements[i];
5949 orig_mesh.
GetNodes()->GetElementDofValues(ip, edofvals);
5971 for (
auto cv : child_vertices)
5972 for (
int ipv = 0; ipv < parent_vertices.
Size(); ipv++)
5973 if (cv == parent_vertices[ipv])
5984 child_nodes_in_parent.
SetSize(0);
5985 const auto *orig_FE = orig_mesh.
GetNodes()->FESpace()->GetFE(ip);
5986 for (
auto pn : node_map)
5988 child_nodes_in_parent.
Append(orig_FE->GetNodes()[pn]);
5991 shape.
SetSize(orig_FE->GetDof(),
5992 simplex_FE->GetDof());
5994 for (
int j = 0; j < simplex_FE->GetNodes().Size(); j++)
5996 const auto &simplex_node = simplex_FE->GetNodes()[j];
5999 simplex_node_in_orig.
Set3(
6000 child_nodes_in_parent[0].x +
6001 simplex_node.x * (child_nodes_in_parent[1].x - child_nodes_in_parent[0].x)
6002 + simplex_node.y * (child_nodes_in_parent[2].x - child_nodes_in_parent[0].x)
6003 + simplex_node.z * (child_nodes_in_parent[(
Dim > 2) ? 3 : 0].x -
6004 child_nodes_in_parent[0].x),
6005 child_nodes_in_parent[0].y +
6006 simplex_node.x * (child_nodes_in_parent[1].y - child_nodes_in_parent[0].y)
6007 + simplex_node.y * (child_nodes_in_parent[2].y - child_nodes_in_parent[0].y)
6008 + simplex_node.z * (child_nodes_in_parent[(
Dim > 2) ? 3 : 0].y -
6009 child_nodes_in_parent[0].y),
6010 child_nodes_in_parent[0].z +
6011 simplex_node.x * (child_nodes_in_parent[1].z - child_nodes_in_parent[0].z)
6012 + simplex_node.y * (child_nodes_in_parent[2].z - child_nodes_in_parent[0].z)
6013 + simplex_node.z * (child_nodes_in_parent[(
Dim > 2) ? 3 : 0].z -
6014 child_nodes_in_parent[0].z));
6016 orig_FE->CalcShape(simplex_node_in_orig, col);
6021 orig_mesh.
GetNodes()->GetElementDofValues(ip, edofvals);
6027 point_matrix.
SetSize(simplex_FE->GetDof(), sdim);
6028 MultAtB(shape, edofvals_mat, point_matrix);
6036 MFEM_ABORT(
"Internal Error!");
6044 Mesh periodic_mesh(orig_mesh,
true);
6050 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
6055 for (
int j = 0; j < nv; j++)
6061 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
6066 for (
int j = 0; j < nv; j++)
6073 return periodic_mesh;
6077 const std::vector<Vector> &translations,
real_t tol)
const
6081 Vector coord(sdim), at(sdim), dx(sdim);
6082 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
6083 xMax = xMin = xDiff = 0.0;
6086 unordered_set<int> bdr_v;
6087 for (
int be = 0; be <
GetNBE(); be++)
6092 for (
int i = 0; i < dofs.
Size(); i++)
6094 bdr_v.insert(dofs[i]);
6097 for (
int j = 0; j < sdim; j++)
6099 xMax[j] = max(xMax[j], coord[j]);
6100 xMin[j] = min(xMin[j], coord[j]);
6104 add(xMax, -1.0, xMin, xDiff);
6113 unordered_map<int, int> replica2primary;
6115 unordered_map<int, unordered_set<int>> primary2replicas;
6118 std::unique_ptr<KDTreeBase<int,real_t>> kdtree;
6119 if (sdim == 1) { kdtree.reset(
new KDTree1D); }
6120 else if (sdim == 2) { kdtree.reset(
new KDTree2D); }
6121 else if (sdim == 3) { kdtree.reset(
new KDTree3D); }
6122 else { MFEM_ABORT(
"Invalid space dimension."); }
6126 for (
const int v : bdr_v)
6128 primary2replicas[v];
6136 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
6138 if (r ==
p) {
return; }
6139 primary2replicas[
p].insert(r);
6140 replica2primary[r] =
p;
6141 for (
const int s : primary2replicas[r])
6143 primary2replicas[
p].insert(s);
6144 replica2primary[s] =
p;
6146 primary2replicas.erase(r);
6149 for (
unsigned int i = 0; i < translations.size(); i++)
6151 for (
int vi : bdr_v)
6154 add(coord, translations[i], at);
6156 const int vj = kdtree->FindClosestPoint(at.GetData());
6158 add(at, -1.0, coord, dx);
6160 if (dx.
Norml2() > dia*tol) {
continue; }
6165 const bool pi = primary2replicas.find(vi) != primary2replicas.end();
6166 const bool pj = primary2replicas.find(vj) != primary2replicas.end();
6172 make_replica(vj, vi);
6177 const int owner_of_vj = replica2primary[vj];
6179 make_replica(vi, owner_of_vj);
6185 const int owner_of_vi = replica2primary[vi];
6186 make_replica(vj, owner_of_vi);
6193 const int owner_of_vi = replica2primary[vi];
6194 const int owner_of_vj = replica2primary[vj];
6195 make_replica(owner_of_vj, owner_of_vi);
6200 std::vector<int> v2v(
GetNV());
6201 for (
size_t i = 0; i < v2v.size(); i++)
6203 v2v[i] =
static_cast<int>(i);
6205 for (
const auto &r2p : replica2primary)
6207 v2v[r2p.first] = r2p.second;
6214 MFEM_VERIFY(
NURBSext,
"Mesh::RefineNURBSFromFile: Not a NURBS mesh!");
6215 mfem::out<<
"Refining NURBS from refinement file: "<<ref_file<<endl;
6218 ifstream input(ref_file);
6225 mfem::out<<
"Knot vectors in ref_file: "<<nkv<<endl;
6227 MFEM_ABORT(
"Refine file does not have the correct number of knot vectors");
6232 for (
int kv = 0; kv < nkv; kv++)
6234 knotVec[kv] =
new Vector();
6235 knotVec[kv]->
Load(input);
6243 for (
int kv = 0; kv < nkv; kv++)
6253 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
6258 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
6275 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
6280 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
6297 mfem_error(
"Mesh::KnotRemove : Not a NURBS mesh!");
6302 mfem_error(
"Mesh::KnotRemove : KnotVector array size mismatch!");
6330 "Refinement factors must be defined for each dimension");
6336 const std::string &kvf)
6338 MFEM_VERIFY(
NURBSext,
"This type of refinement is only for NURBS meshes");
6347 cf1 = (cf1 &&
f == 1);
6357 MFEM_VERIFY(!usingKVF,
"This refinement type is not supported for this"
6358 " NURBS mesh type");
6366 for (
int i=0; i<cf.
Size(); ++i) { cf[i] *= rf[i]; }
6392 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
6416 for (
int i = 0; i <
elements.Size(); i++)
6426 for (
int i = 0; i <
boundary.Size(); i++)
6444 for (
int i = 0; i < vd; i++)
6511 input >> edge_to_ukv[j] >> v[0] >> v[1];
6514 edge_to_ukv[j] = -1 - edge_to_ukv[j];
6535 if (edge_to_ukv.
Size() == 0)
6547 const int NPKV = NP *
dim;
6548 constexpr int notset = -9999999;
6550 auto sign = [](
int i) {
return -1 - i; };
6551 auto unsign = [](
int i) {
return (i < 0) ? -1 - i : i; };
6553 auto edge_to_dim = [](
int i) {
return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };
6563 for (
int i = 0; i < NP; i++)
6567 edge_to_ukv[i] = (v[1] > v[0]) ? i : sign(i);
6578 edge_to_pkv = notset;
6584 for (
int i = 0; i < NPKV; i++)
6588 std::function<int(
int)> get_root;
6589 get_root = [&pkv_map, &get_root](
int i) ->
int
6591 return (pkv_map[i] == i) ? i : get_root(pkv_map[i]);
6593 auto unite = [&pkv_map, &get_root](
int i,
int j)
6595 const int ri = get_root(i);
6596 const int rj = get_root(j);
6597 if (ri == rj) {
return; }
6599 (ri < rj) ? pkv_map[rj] = ri : pkv_map[ri] = rj;
6603 for (
int p = 0;
p < NP;
p++)
6608 for (
int i = 0; i < edges.
Size(); i++)
6610 const int edge = edges[i];
6611 const int d = edge_to_dim(i);
6612 const int pkv =
p*
dim+d;
6615 if (edge_to_pkv[edge] != notset)
6617 const int pkv_other = unsign(edge_to_pkv[edge]);
6618 unite(pkv, pkv_other);
6624 edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : sign(pkv);
6632 for (
int i = 0; i < NPKV; i++)
6634 pkv_to_rpkv[i] = get_root(pkv_map[i]);
6635 ukv_to_rpkv[i] = pkv_to_rpkv[i];
6641 std::map<int, int> rpkv_to_ukv;
6642 for (
int i = 0; i < ukv_to_rpkv.
Size(); i++)
6644 rpkv_to_ukv[ukv_to_rpkv[i]] = i;
6651 const int pkv = unsign(edge_to_pkv[i]);
6652 const int rpkv = pkv_to_rpkv[pkv];
6653 const int ukv = rpkv_to_ukv[rpkv];
6654 edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? sign(ukv) : ukv;
6674 int inputNumOfEdges = -1;
6677 input >> inputNumOfEdges;
6679 MFEM_VERIFY(
NumOfEdges == inputNumOfEdges,
"");
6686 input >> ukv >> v[0] >> v[1];
6688 for (
int i=0; i<2; ++i)
6697 edge_to_ukv[j] = ukv;
6706 if (
p.Size() >= v.
Size())
6708 for (
int d = 0; d < v.
Size(); d++)
6716 for (d = 0; d <
p.Size(); d++)
6720 for ( ; d < v.
Size(); d++)
6733 nodes.ProjectCoefficient(xyz);
6763 const bool warn =
true;
6766 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
6770 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
6771 " H1-continuous mesh!\n "
6772 "If this is the desired behavior, you can silence"
6773 " this warning by converting\n "
6774 "the NURBS mesh to high-order mesh in advance by"
6775 " calling the method\n "
6776 "Mesh::SetCurvature().");
6807 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
6821 const int old_space_dim =
spaceDim;
6834 MFEM_ASSERT(
nodes != NULL,
"");
6838 nodes->GetNodalValues(vert_val, i+1);
6850 case 1:
return GetNV();
6886#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
6887static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
6892 int i, j, k, wo = 0, fo = 0;
6901 int *vi =
elements[i]->GetVertices();
6904 for (j = 0; j < 3; j++)
6908 for (j = 0; j < 2; j++)
6909 for (k = 0; k < 2; k++)
6911 J(j, k) = v[j+1][k] - v[0][k];
6932 MFEM_ABORT(
"Invalid 2D element type \""
6949 int *vi =
elements[i]->GetVertices();
6955 for (j = 0; j < 4; j++)
6959 for (j = 0; j < 3; j++)
6960 for (k = 0; k < 3; k++)
6962 J(j, k) = v[j+1][k] - v[0][k];
7022 MFEM_ABORT(
"Invalid 3D element type \""
7028#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
7031 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
7032 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
7036 MFEM_CONTRACT_VAR(fo);
7049 if (test[0] == base[0])
7050 if (test[1] == base[1])
7058 else if (test[0] == base[1])
7059 if (test[1] == base[0])
7068 if (test[1] == base[0])
7079 for (
int j = 0; j < 3; j++)
7080 if (test[aor[j]] != base[j])
7082 mfem::err <<
"Mesh::GetTriOrientation(...)" << endl;
7084 for (
int k = 0; k < 3; k++)
7089 for (
int k = 0; k < 3; k++)
7110 const int oo[6][6] =
7120 int ori_a_c = oo[ori_a_b][ori_b_c];
7126 const int inv_ori[6] = {0, 1, 4, 3, 2, 5};
7127 return inv_ori[ori];
7134 for (i = 0; i < 4; i++)
7135 if (test[i] == base[0])
7142 if (test[(i+1)%4] == base[1])
7151 for (
int j = 0; j < 4; j++)
7152 if (test[aor[j]] != base[j])
7154 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
7156 for (
int k = 0; k < 4; k++)
7161 for (
int k = 0; k < 4; k++)
7170 if (test[(i+1)%4] == base[1])
7187 const int oo[8][8] =
7189 {0, 1, 2, 3, 4, 5, 6, 7},
7190 {1, 0, 3, 2, 5, 4, 7, 6},
7191 {2, 7, 4, 1, 6, 3, 0, 5},
7192 {3, 6, 5, 0, 7, 2, 1, 4},
7193 {4, 5, 6, 7, 0, 1, 2, 3},
7194 {5, 4, 7, 6, 1, 0, 3, 2},
7195 {6, 3, 0, 5, 2, 7, 4, 1},
7196 {7, 2, 1, 4, 3, 6, 5, 0}
7199 int ori_a_c = oo[ori_a_b][ori_b_c];
7205 const int inv_ori[8] = {0, 1, 6, 3, 4, 5, 2, 7};
7206 return inv_ori[ori];
7217 if (test[0] == base[0])
7218 if (test[1] == base[1])
7219 if (test[2] == base[2])
7227 else if (test[2] == base[1])
7228 if (test[3] == base[2])
7237 if (test[1] == base[2])
7245 else if (test[1] == base[0])
7246 if (test[2] == base[1])
7247 if (test[0] == base[2])
7255 else if (test[3] == base[1])
7256 if (test[2] == base[2])
7265 if (test[3] == base[2])
7273 else if (test[2] == base[0])
7274 if (test[3] == base[1])
7275 if (test[0] == base[2])
7283 else if (test[0] == base[1])
7284 if (test[1] == base[2])
7293 if (test[3] == base[2])
7302 if (test[0] == base[1])
7303 if (test[2] == base[2])
7311 else if (test[1] == base[1])
7312 if (test[0] == base[2])
7321 if (test[1] == base[2])
7332 for (
int j = 0; j < 4; j++)
7333 if (test[aor[j]] != base[j])
7358 int *bv =
boundary[i]->GetVertices();
7378 if (
faces_info[fi].Elem2No >= 0) {
continue; }
7381 int *bv =
boundary[i]->GetVertices();
7383 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
7384 const int *fv =
faces[fi]->GetVertices();
7401 MFEM_ABORT(
"Invalid 2D boundary element type \""
7402 << bdr_type <<
"\"");
7407 if (orientation % 2 == 0) {
continue; }
7409 if (!fix_it) {
continue; }
7445 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
7463 MFEM_ASSERT(o >= 0 && o < 2,
"Invalid orientation for Geometry::SEGMENT!");
7475 MFEM_ASSERT(o >= 0 && o < 6,
"Invalid orientation for Geometry::TRIANGLE!");
7488 fip.
x = 1.0 - ip.
x - ip.
y;
7493 fip.
x = 1.0 - ip.
x - ip.
y;
7499 fip.
y = 1.0 - ip.
x - ip.
y;
7504 fip.
y = 1.0 - ip.
x - ip.
y;
7509 MFEM_ASSERT(o >= 0 && o < 8,
"Invalid orientation for Geometry::SQUARE!");
7553 MFEM_ABORT(
"Unsupported face geometry for TransformBdrElementToFace!");
7560 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
7571 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
7596 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
7597 "is not generated.");
7600 const int *v =
elements[i]->GetVertices();
7601 const int ne =
elements[i]->GetNEdges();
7603 for (
int j = 0; j < ne; j++)
7605 const int *e =
elements[i]->GetEdgeVertices(j);
7606 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7617 const int *v =
boundary[i]->GetVertices();
7618 cor[0] = (v[0] < v[1]) ? (1) : (-1);
7631 const int *v =
boundary[i]->GetVertices();
7632 const int ne =
boundary[i]->GetNEdges();
7634 for (
int j = 0; j < ne; j++)
7636 const int *e =
boundary[i]->GetEdgeVertices(j);
7637 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7649 const int *v =
faces[i]->GetVertices();
7650 o[0] = (v[0] < v[1]) ? (1) : (-1);
7662 const int *v =
faces[i]->GetVertices();
7663 const int ne =
faces[i]->GetNEdges();
7665 for (
int j = 0; j < ne; j++)
7667 const int *e =
faces[i]->GetEdgeVertices(j);
7668 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7696 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
7743 const int nv =
elements[i]->GetNVertices();
7744 const int *v =
elements[i]->GetVertices();
7745 for (
int j = 0; j < nv; j++)
7755 const int nv =
elements[i]->GetNVertices();
7756 const int *v =
elements[i]->GetVertices();
7757 for (
int j = 0; j < nv; j++)
7776 const int nv =
boundary[i]->GetNVertices();
7777 const int *v =
boundary[i]->GetVertices();
7778 for (
int j = 0; j < nv; j++)
7784 vert_bdr_elem->
MakeJ();
7788 const int nv =
boundary[i]->GetNVertices();
7789 const int *v =
boundary[i]->GetVertices();
7790 for (
int j = 0; j < nv; j++)
7798 return vert_bdr_elem;
7837 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
7841 int n = el_faces.
Size();
7844 for (
int j = 0; j < n; j++)
7848 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
7852 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
7853 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
7870 for (
auto f : elem_faces)
7891 const int *bv =
boundary[i]->GetVertices();
7901 default: MFEM_ABORT(
"invalid geometry");
7910 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
7913 const int *bv =
boundary[bdr_el]->GetVertices();
7921 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
7928 int bdr_el,
int &el,
int &info)
const
7933 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
7936 const int *bv =
boundary[bdr_el]->GetVertices();
7944 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
7983 for (j = 0; j < nv; j++)
7985 pointmat(k, j) =
vertices[v[j]](k);
8000 for (j = 0; j < nv; j++)
8002 pointmat(k, j) =
vertices[v[j]](k);
8014 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
8017 return sqrt(length);
8025 for (
int i = 0; i < elem_array.
Size(); i++)
8030 for (
int i = 0; i < elem_array.
Size(); i++)
8032 const int *v = elem_array[i]->GetVertices();
8033 const int ne = elem_array[i]->GetNEdges();
8034 for (
int j = 0; j < ne; j++)
8036 const int *e = elem_array[i]->GetEdgeVertices(j);
8050 v_to_v.
Push(v[0], v[1]);
8057 const int *v =
elements[i]->GetVertices();
8058 const int ne =
elements[i]->GetNEdges();
8059 for (
int j = 0; j < ne; j++)
8061 const int *e =
elements[i]->GetEdgeVertices(j);
8062 v_to_v.
Push(v[e[0]], v[e[1]]);
8070 int i, NumberOfEdges;
8086 const int *v =
boundary[i]->GetVertices();
8100 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
8104 return NumberOfEdges;
8163 if (
faces[gf] == NULL)
8195 if (
faces[gf] == NULL)
8205 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
8206 "Interior edge found between 2D elements "
8208 <<
" and " << el <<
".");
8209 int *v =
faces[gf]->GetVertices();
8211 if (v[1] == v0 && v[0] == v1)
8215 else if (v[0] == v0 && v[1] == v1)
8225 MFEM_ABORT(
"internal error");
8231 int v0,
int v1,
int v2)
8233 if (
faces[gf] == NULL)
8243 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
8244 "Interior triangular face found connecting elements "
8246 <<
" and " << el <<
".");
8247 int orientation, vv[3] = { v0, v1, v2 };
8254 faces_info[gf].Elem2Inf = 64 * lf + orientation;
8259 int v0,
int v1,
int v2,
int v3)
8271 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
8272 "Interior quadrilateral face found connecting elements "
8274 <<
" and " << el <<
".");
8275 int vv[4] = { v0, v1, v2, v3 };
8300 faces.SetSize(nfaces);
8302 for (
int i = 0; i < nfaces; ++i)
8321 const int ne =
elements[i]->GetNEdges();
8322 for (
int j = 0; j < ne; j++)
8324 const int *e =
elements[i]->GetEdgeVertices(j);
8335 for (
int j = 0; j < 4; j++)
8339 v[fv[0]], v[fv[1]], v[fv[2]]);
8345 for (
int j = 0; j < 2; j++)
8349 v[fv[0]], v[fv[1]], v[fv[2]]);
8351 for (
int j = 2; j < 5; j++)
8355 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8361 for (
int j = 0; j < 1; j++)
8365 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8367 for (
int j = 1; j < 5; j++)
8371 v[fv[0]], v[fv[1]], v[fv[2]]);
8377 for (
int j = 0; j < 6; j++)
8381 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8386 MFEM_ABORT(
"Unexpected type of Element.");
8394 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
8405 nc_faces_info.Reserve(list.masters.Size() + list.slaves.Size());
8412 if (master.index >= nfaces) {
continue; }
8418 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
8419 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
8425 if (slave.index < 0 ||
8426 slave.index >= nfaces ||
8427 slave.master >= nfaces)
8446 list.point_matrices[slave.geom][slave.matrix]));
8455 const int *v =
elements[i]->GetVertices();
8460 for (
int j = 0; j < 4; j++)
8463 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8469 for (
int j = 0; j < 1; j++)
8472 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8474 for (
int j = 1; j < 5; j++)
8477 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8483 for (
int j = 0; j < 2; j++)
8486 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
8488 for (
int j = 2; j < 5; j++)
8491 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8499 for (
int j = 0; j < 6; j++)
8502 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
8531 for (
int j = 0; j < 4; j++)
8535 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
8541 for (
int j = 0; j < 2; j++)
8545 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
8547 for (
int j = 2; j < 5; j++)
8551 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
8557 for (
int j = 0; j < 1; j++)
8561 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
8563 for (
int j = 1; j < 5; j++)
8567 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
8575 for (
int j = 0; j < 6; j++)
8579 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
8584 MFEM_ABORT(
"Unexpected type of Element.");
8598 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
8603 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
8607 MFEM_ABORT(
"Unexpected type of boundary Element.");
8621void Rotate3(
int &
a,
int &
b,
int &c)
8653 Table *old_elem_vert = NULL;
8664 int *v =
elements[i]->GetVertices();
8666 Rotate3(v[0], v[1], v[2]);
8669 Rotate3(v[1], v[2], v[3]);
8682 int *v =
boundary[i]->GetVertices();
8684 Rotate3(v[0], v[1], v[2]);
8700 delete old_elem_vert;
8716 if (
p[i] < pmin[i]) { pmin[i] =
p[i]; }
8717 if (
p[i] > pmax[i]) { pmax[i] =
p[i]; }
8731 for (
int i =
spaceDim-1; i >= 0; i--)
8733 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
8734 if (idx < 0) { idx = 0; }
8735 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
8736 part = part * nxyz[i] + idx;
8738 partitioning[el] = part;
8741 return partitioning;
8751#ifdef MFEM_USE_METIS
8753 int print_messages = 1;
8756 int init_flag, fin_flag;
8757 MPI_Initialized(&init_flag);
8758 MPI_Finalized(&fin_flag);
8759 if (init_flag && !fin_flag)
8763 if (rank != 0) { print_messages = 0; }
8767 int i, *partitioning;
8777 partitioning[i] = 0;
8784 partitioning[i] = i;
8790#ifndef MFEM_USE_METIS_5
8802 bool freedata =
false;
8804 idx_t *mpartitioning;
8807 if (
sizeof(
idx_t) ==
sizeof(int))
8811 mpartitioning = (
idx_t*) partitioning;
8820 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
8821 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
8822 mpartitioning =
new idx_t[n];
8825#ifndef MFEM_USE_METIS_5
8828 METIS_SetDefaultOptions(options);
8829 options[METIS_OPTION_CONTIG] = 1;
8837 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
8842 if (part_method >= 0 && part_method <= 2)
8844 for (i = 0; i < n; i++)
8850 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
8856 if (part_method == 0 || part_method == 3)
8858#ifndef MFEM_USE_METIS_5
8887 " error in METIS_PartGraphRecursive!");
8894 if (part_method == 1 || part_method == 4)
8896#ifndef MFEM_USE_METIS_5
8925 " error in METIS_PartGraphKway!");
8932 if (part_method == 2 || part_method == 5)
8934#ifndef MFEM_USE_METIS_5
8947 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
8964 " error in METIS_PartGraphKway!");
8972 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
8976 nparts = (int) mparts;
8977 if (mpartitioning != (
idx_t*)partitioning)
8981 partitioning[k] = mpartitioning[k];
8988 delete[] mpartitioning;
9004 auto count_partition_elements = [&]()
9006 for (i = 0; i < nparts; i++)
9014 psize[partitioning[i]].one++;
9018 for (i = 0; i < nparts; i++)
9020 if (psize[i].one == 0) { empty_parts++; }
9024 count_partition_elements();
9032 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
9033 << empty_parts <<
" empty parts!"
9034 <<
" Applying a simple fix ..." << endl;
9039 for (i = nparts-1; i > nparts-1-empty_parts; i--)
9046 for (i = nparts-1; i > nparts-1-empty_parts; i--)
9048 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
9054 partitioning[j] = psize[nparts-1-i].two;
9061 count_partition_elements();
9065 return partitioning;
9069 mfem_error(
"Mesh::GeneratePartitioning(...): "
9070 "MFEM was compiled without Metis.");
9084 int num_elem, *i_elem_elem, *j_elem_elem;
9086 num_elem = elem_elem.
Size();
9087 i_elem_elem = elem_elem.
GetI();
9088 j_elem_elem = elem_elem.
GetJ();
9093 int stack_p, stack_top_p, elem;
9097 for (i = 0; i < num_elem; i++)
9099 if (partitioning[i] > num_part)
9101 num_part = partitioning[i];
9108 for (i = 0; i < num_part; i++)
9115 for (elem = 0; elem < num_elem; elem++)
9117 if (component[elem] >= 0)
9122 component[elem] = num_comp[partitioning[elem]]++;
9124 elem_stack[stack_top_p++] = elem;
9126 for ( ; stack_p < stack_top_p; stack_p++)
9128 i = elem_stack[stack_p];
9129 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
9132 if (partitioning[k] == partitioning[i])
9134 if (component[k] < 0)
9136 component[k] = component[i];
9137 elem_stack[stack_top_p++] = k;
9139 else if (component[k] != component[i])
9151 int i, n_empty, n_mcomp;
9159 n_empty = n_mcomp = 0;
9160 for (i = 0; i < num_comp.
Size(); i++)
9161 if (num_comp[i] == 0)
9165 else if (num_comp[i] > 1)
9172 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
9173 <<
"The following subdomains are empty :\n";
9174 for (i = 0; i < num_comp.
Size(); i++)
9175 if (num_comp[i] == 0)
9183 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
9184 <<
"The following subdomains are NOT connected :\n";
9185 for (i = 0; i < num_comp.
Size(); i++)
9186 if (num_comp[i] > 1)
9192 if (n_empty == 0 && n_mcomp == 0)
9193 mfem::out <<
"Mesh::CheckPartitioning(...) : "
9194 "All subdomains are connected." << endl;
9218 c(0) =
a[0]*
a[3]-
a[1]*
a[2];
9219 c(1) =
a[0]*
b[3]-
a[1]*
b[2]+
b[0]*
a[3]-
b[1]*
a[2];
9220 c(2) =
b[0]*
b[3]-
b[1]*
b[2];
9241 c(0) = (
a[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
9242 a[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
9243 a[2] * (
a[3] *
a[7] -
a[4] *
a[6]));
9245 c(1) = (
b[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
9246 b[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
9247 b[2] * (
a[3] *
a[7] -
a[4] *
a[6]) +
9249 a[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
9250 a[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
9251 a[2] * (
b[3] *
a[7] -
b[4] *
a[6]) +
9253 a[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
9254 a[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
9255 a[2] * (
a[3] *
b[7] -
a[4] *
b[6]));
9257 c(2) = (
a[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
9258 a[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
9259 a[2] * (
b[3] *
b[7] -
b[4] *
b[6]) +
9261 b[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
9262 b[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
9263 b[2] * (
a[3] *
b[7] -
a[4] *
b[6]) +
9265 b[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
9266 b[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
9267 b[2] * (
b[3] *
a[7] -
b[4] *
a[6]));
9269 c(3) = (
b[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
9270 b[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
9271 b[2] * (
b[3] *
b[7] -
b[4] *
b[6]));
9317 real_t a = z(2),
b = z(1), c = z(0);
9325 x(0) = x(1) = -0.5 *
b /
a;
9330 x(0) = -(x(1) = fabs(0.5 * sqrt(D) /
a));
9338 t = -0.5 * (
b + sqrt(D));
9342 t = -0.5 * (
b - sqrt(D));
9356 real_t a = z(2)/z(3),
b = z(1)/z(3), c = z(0)/z(3);
9360 real_t R = (2 *
a *
a *
a - 9 *
a *
b + 27 * c) / 54;
9368 x(0) = x(1) = x(2) = -
a / 3;
9376 x(0) = -2 * sqrtQ -
a / 3;
9377 x(1) = x(2) = sqrtQ -
a / 3;
9381 x(0) = x(1) = - sqrtQ -
a / 3;
9382 x(2) = 2 * sqrtQ -
a / 3;
9389 real_t theta = acos(R / sqrt(Q3));
9392 x0 = A * cos(theta / 3) -
a / 3;
9393 x1 = A * cos((theta + 2.0 * M_PI) / 3) -
a / 3;
9394 x2 = A * cos((theta - 2.0 * M_PI) / 3) -
a / 3;
9419 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
9423 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
9425 x(0) = A + Q / A -
a / 3;
9434 const real_t factor,
const int Dim)
9437 c(0) = c0 * (1.0 - pow(factor, -Dim));
9439 for (
int j = 0; j < nr; j++)
9451 c(0) = c0 * (1.0 - pow(factor, Dim));
9453 for (
int j = 0; j < nr; j++)
9472 const real_t factor = 2.0;
9487 for (
int k = 0; k < nv; k++)
9490 V(j, k) = displacements(v[k]+j*nvs);
9520 for (
int j = 0; j < nv; j++)
9546 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
9549 vertices[i](j) += displacements(j*nv+i);
9557 for (
int i = 0; i < nv; i++)
9560 vert_coord(j*nv+i) =
vertices[i](j);
9566 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
9569 vertices[i](j) = vert_coord(j*nv+i);
9599 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
9616 (*Nodes) += displacements;
9631 node_coord = (*Nodes);
9643 (*Nodes) = node_coord;
9706 for (j = 1; j < n; j++)
9743 int quad_counter = 0;
9765 const int attr =
elements[i]->GetAttribute();
9766 int *v =
elements[i]->GetVertices();
9772 for (
int ei = 0; ei < 3; ei++)
9774 for (
int k = 0; k < 2; k++)
9782 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
9784 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
9786 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
9788 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
9792 const int qe = quad_counter;
9796 for (
int ei = 0; ei < 4; ei++)
9798 for (
int k = 0; k < 2; k++)
9806 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
9808 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
9810 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
9812 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
9816 MFEM_ABORT(
"unknown element type: " << el_type);
9826 const int attr =
boundary[i]->GetAttribute();
9827 int *v =
boundary[i]->GetVertices();
9836 static const real_t A = 0.0, B = 0.5, C = 1.0;
9837 static real_t tri_children[2*3*4] =
9844 static real_t quad_children[2*4*4] =
9858 for (
int i = 0; i <
elements.Size(); i++)
9879 if (!
Nodes || update_nodes)
9909 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
9912 int NumOfQuadFaces = 0;
9918 for (
int i = 0; i <
faces.Size(); i++)
9922 f2qf[i] = NumOfQuadFaces;
9933 int hex_counter = 0;
9936 for (
int i = 0; i <
elements.Size(); i++)
9945 int pyr_counter = 0;
9948 for (
int i = 0; i <
elements.Size(); i++)
9966 DSTable *v_to_v_ptr = v_to_v_p;
9982 std::sort(row_start, J_v2v.
end());
9985 for (
int i = 0; i < J_v2v.
Size(); i++)
9987 e2v[J_v2v[i].two] = i;
10000 it.SetIndex(e2v[it.Index()]);
10010 const int oelem = oface + NumOfQuadFaces;
10023 const int attr =
elements[i]->GetAttribute();
10024 int *v =
elements[i]->GetVertices();
10031 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
10039 for (
int ei = 0; ei < 6; ei++)
10041 for (
int k = 0; k < 2; k++)
10051 const int rt_algo = 1;
10062 real_t len_sqr, min_len;
10064 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
10065 sqr(J(1,0)-J(1,1)-J(1,2)) +
10066 sqr(J(2,0)-J(2,1)-J(2,2));
10069 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
10070 sqr(J(1,1)-J(1,0)-J(1,2)) +
10071 sqr(J(2,1)-J(2,0)-J(2,2));
10072 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
10074 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
10075 sqr(J(1,2)-J(1,0)-J(1,1)) +
10076 sqr(J(2,2)-J(2,0)-J(2,1));
10077 if (len_sqr < min_len) { rt = 2; }
10082 real_t Em_data[18], Js_data[9], Jp_data[9];
10084 DenseMatrix Js(Js_data, 3, 3), Jp(Jp_data, 3, 3);
10087 for (
int s = 0; s < 3; s++)
10089 for (
int t = 0; t < 3; t++)
10091 Em(t,s) = 0.5*J(t,s);
10094 for (
int t = 0; t < 3; t++)
10096 Em(t,3) = 0.5*(J(t,0)+J(t,1));
10097 Em(t,4) = 0.5*(J(t,0)+J(t,2));
10098 Em(t,5) = 0.5*(J(t,1)+J(t,2));
10102 for (
int t = 0; t < 3; t++)
10104 Js(t,0) = Em(t,5)-Em(t,0);
10105 Js(t,1) = Em(t,1)-Em(t,0);
10106 Js(t,2) = Em(t,2)-Em(t,0);
10110 for (
int t = 0; t < 3; t++)
10112 Js(t,0) = Em(t,5)-Em(t,0);
10113 Js(t,1) = Em(t,2)-Em(t,0);
10114 Js(t,2) = Em(t,4)-Em(t,0);
10118 kappa_min = std::max(ar1, ar2);
10122 for (
int t = 0; t < 3; t++)
10124 Js(t,0) = Em(t,0)-Em(t,1);
10125 Js(t,1) = Em(t,4)-Em(t,1);
10126 Js(t,2) = Em(t,2)-Em(t,1);
10130 for (
int t = 0; t < 3; t++)
10132 Js(t,0) = Em(t,2)-Em(t,1);
10133 Js(t,1) = Em(t,4)-Em(t,1);
10134 Js(t,2) = Em(t,5)-Em(t,1);
10138 kappa = std::max(ar1, ar2);
10139 if (
kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
10142 for (
int t = 0; t < 3; t++)
10144 Js(t,0) = Em(t,0)-Em(t,2);
10145 Js(t,1) = Em(t,1)-Em(t,2);
10146 Js(t,2) = Em(t,3)-Em(t,2);
10150 for (
int t = 0; t < 3; t++)
10152 Js(t,0) = Em(t,1)-Em(t,2);
10153 Js(t,1) = Em(t,5)-Em(t,2);
10154 Js(t,2) = Em(t,3)-Em(t,2);
10158 kappa = std::max(ar1, ar2);
10159 if (
kappa < kappa_min) { rt = 2; }
10162 static const int mv_all[3][4][4] =
10164 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
10165 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
10166 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
10168 const int (&mv)[4][4] = mv_all[rt];
10170#ifndef MFEM_USE_MEMALLOC
10171 new_elements[j+0] =
10172 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
10173 new_elements[j+1] =
10174 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
10175 new_elements[j+2] =
10176 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
10177 new_elements[j+3] =
10178 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
10180 for (
int k = 0; k < 4; k++)
10182 new_elements[j+4+k] =
10183 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
10184 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
10188 new_elements[j+0] = tet =
TetMemory.Alloc();
10189 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
10191 new_elements[j+1] = tet =
TetMemory.Alloc();
10192 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
10194 new_elements[j+2] = tet =
TetMemory.Alloc();
10195 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
10197 new_elements[j+3] = tet =
TetMemory.Alloc();
10198 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
10200 for (
int k = 0; k < 4; k++)
10202 new_elements[j+4+k] = tet =
TetMemory.Alloc();
10203 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
10204 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
10207 for (
int k = 0; k < 4; k++)
10212 for (
int k = 0; k < 4; k++)
10226 for (
int fi = 2; fi < 5; fi++)
10228 for (
int k = 0; k < 4; k++)
10235 for (
int ei = 0; ei < 9; ei++)
10237 for (
int k = 0; k < 2; k++)
10244 const int qf2 = f2qf[
f[2]];
10245 const int qf3 = f2qf[
f[3]];
10246 const int qf4 = f2qf[
f[4]];
10248 new_elements[j++] =
10249 new Wedge(v[0], oedge+e[0], oedge+e[2],
10250 oedge+e[6], oface+qf2, oface+qf4, attr);
10252 new_elements[j++] =
10253 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
10254 oface+qf3, oface+qf4, oface+qf2, attr);
10256 new_elements[j++] =
10257 new Wedge(oedge+e[0], v[1], oedge+e[1],
10258 oface+qf2, oedge+e[7], oface+qf3, attr);
10260 new_elements[j++] =
10261 new Wedge(oedge+e[2], oedge+e[1], v[2],
10262 oface+qf4, oface+qf3, oedge+e[8], attr);
10264 new_elements[j++] =
10265 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
10266 v[3], oedge+e[3], oedge+e[5], attr);
10268 new_elements[j++] =
10269 new Wedge(oface+qf3, oface+qf4, oface+qf2,
10270 oedge+e[4], oedge+e[5], oedge+e[3], attr);
10272 new_elements[j++] =
10273 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
10274 oedge+e[3], v[4], oedge+e[4], attr);
10276 new_elements[j++] =
10277 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
10278 oedge+e[5], oedge+e[4], v[5], attr);
10287 for (
int fi = 0; fi < 1; fi++)
10289 for (
int k = 0; k < 4; k++)
10296 for (
int ei = 0; ei < 8; ei++)
10298 for (
int k = 0; k < 2; k++)
10305 const int qf0 = f2qf[
f[0]];
10307 new_elements[j++] =
10308 new Pyramid(v[0], oedge+e[0], oface+qf0,
10309 oedge+e[3], oedge+e[4], attr);
10311 new_elements[j++] =
10312 new Pyramid(oedge+e[0], v[1], oedge+e[1],
10313 oface+qf0, oedge+e[5], attr);
10315 new_elements[j++] =
10316 new Pyramid(oface+qf0, oedge+e[1], v[2],
10317 oedge+e[2], oedge+e[6], attr);
10319 new_elements[j++] =
10320 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
10321 v[3], oedge+e[7], attr);
10323 new_elements[j++] =
10324 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
10325 oedge+e[7], v[4], attr);
10327 new_elements[j++] =
10328 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
10329 oedge+e[4], oface+qf0, attr);
10331#ifndef MFEM_USE_MEMALLOC
10332 new_elements[j++] =
10333 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
10336 new_elements[j++] =
10337 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
10340 new_elements[j++] =
10341 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
10344 new_elements[j++] =
10345 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
10349 new_elements[j++] = tet =
TetMemory.Alloc();
10350 tet->
Init(oedge+e[0], oedge+e[4], oedge+e[5],
10353 new_elements[j++] = tet =
TetMemory.Alloc();
10354 tet->
Init(oedge+e[1], oedge+e[5], oedge+e[6],
10357 new_elements[j++] = tet =
TetMemory.Alloc();
10358 tet->
Init(oedge+e[2], oedge+e[6], oedge+e[7],
10361 new_elements[j++] = tet =
TetMemory.Alloc();
10362 tet->
Init(oedge+e[3], oedge+e[7], oedge+e[4],
10375 const int he = hex_counter;
10380 if (f2qf.
Size() == 0)
10386 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[
f[k]]; }
10392 for (
int fi = 0; fi < 6; fi++)
10394 for (
int k = 0; k < 4; k++)
10401 for (
int ei = 0; ei < 12; ei++)
10403 for (
int k = 0; k < 2; k++)
10410 new_elements[j++] =
10411 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
10412 oedge+e[3], oedge+e[8], oface+qf[1],
10413 oelem+he, oface+qf[4], attr);
10414 new_elements[j++] =
10415 new Hexahedron(oedge+e[0], v[1], oedge+e[1],
10416 oface+qf[0], oface+qf[1], oedge+e[9],
10417 oface+qf[2], oelem+he, attr);
10418 new_elements[j++] =
10419 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
10420 oedge+e[2], oelem+he, oface+qf[2],
10421 oedge+e[10], oface+qf[3], attr);
10422 new_elements[j++] =
10423 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
10424 v[3], oface+qf[4], oelem+he,
10425 oface+qf[3], oedge+e[11], attr);
10426 new_elements[j++] =
10427 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
10428 oface+qf[4], v[4], oedge+e[4],
10429 oface+qf[5], oedge+e[7], attr);
10430 new_elements[j++] =
10431 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
10432 oelem+he, oedge+e[4], v[5],
10433 oedge+e[5], oface+qf[5], attr);
10434 new_elements[j++] =
10435 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
10436 oface+qf[3], oface+qf[5], oedge+e[5],
10437 v[6], oedge+e[6], attr);
10438 new_elements[j++] =
10439 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
10440 oedge+e[11], oedge+e[7], oface+qf[5],
10441 oedge+e[6], v[7], attr);
10446 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
10458 const int attr =
boundary[i]->GetAttribute();
10459 int *v =
boundary[i]->GetVertices();
10466 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
10472 new_boundary[j++] =
10473 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
10474 new_boundary[j++] =
10475 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
10476 new_boundary[j++] =
10477 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
10478 new_boundary[j++] =
10479 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
10486 new_boundary[j++] =
10487 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
10488 new_boundary[j++] =
10489 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
10490 new_boundary[j++] =
10491 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
10492 new_boundary[j++] =
10493 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
10497 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
10503 static const real_t A = 0.0, B = 0.5, C = 1.0, D = -1.0;
10504 static real_t tet_children[3*4*16] =
10506 A,A,A, B,A,A, A,B,A, A,A,B,
10507 B,A,A, C,A,A, B,B,A, B,A,B,
10508 A,B,A, B,B,A, A,C,A, A,B,B,
10509 A,A,B, B,A,B, A,B,B, A,A,C,
10514 B,A,A, A,B,B, A,B,A, A,A,B,
10515 B,A,A, A,B,B, A,A,B, B,A,B,
10516 B,A,A, A,B,B, B,A,B, B,B,A,
10517 B,A,A, A,B,B, B,B,A, A,B,A,
10519 A,B,A, B,A,A, B,A,B, A,A,B,
10520 A,B,A, A,A,B, B,A,B, A,B,B,
10521 A,B,A, A,B,B, B,A,B, B,B,A,
10522 A,B,A, B,B,A, B,A,B, B,A,A,
10524 A,A,B, B,A,A, A,B,A, B,B,A,
10525 A,A,B, A,B,A, A,B,B, B,B,A,
10526 A,A,B, A,B,B, B,A,B, B,B,A,
10527 A,A,B, B,A,B, B,A,A, B,B,A
10529 static real_t pyr_children[3*5*10] =
10531 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
10532 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
10533 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
10534 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
10535 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
10536 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
10537 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
10538 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
10539 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
10540 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
10542 static real_t pri_children[3*6*8] =
10544 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
10545 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
10546 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
10547 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
10548 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
10549 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
10550 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
10551 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
10553 static real_t hex_children[3*8*8] =
10555 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
10556 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
10557 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
10558 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
10559 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
10560 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
10561 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
10562 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
10574 for (
int i = 0; i <
elements.Size(); i++)
10605 int i, j, ind, nedges;
10612 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
10626 for (j = 0; j < marked_el.
Size(); j++)
10631 int new_v = cnv + j, new_e = cne + j;
10640 static real_t seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
10642 UseExternalData(seg_children, 1, 2, 3);
10655 int *edge1 =
new int[nedges];
10656 int *edge2 =
new int[nedges];
10657 int *middle =
new int[nedges];
10659 for (i = 0; i < nedges; i++)
10661 edge1[i] = edge2[i] = middle[i] = -1;
10667 for (j = 1; j < v.
Size(); j++)
10669 ind = v_to_v(v[j-1], v[j]);
10670 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
10672 ind = v_to_v(v[0], v[v.
Size()-1]);
10673 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
10677 for (i = 0; i < marked_el.
Size(); i++)
10683 int need_refinement;
10686 need_refinement = 0;
10687 for (i = 0; i < nedges; i++)
10689 if (middle[i] != -1 && edge1[i] != -1)
10691 need_refinement = 1;
10696 while (need_refinement == 1);
10699 int v1[2], v2[2],
bisect, temp;
10701 for (i = 0; i < temp; i++)
10704 bisect = v_to_v(v[0], v[1]);
10705 if (middle[
bisect] != -1)
10709 v1[0] = v[0]; v1[1] = middle[
bisect];
10710 v2[0] = middle[
bisect]; v2[1] = v[1];
10716 mfem_error(
"Only bisection of segment is implemented"
10740 MFEM_VERIFY(
GetNE() == 0 ||
10742 "tetrahedral mesh is not marked for refinement:"
10743 " call Finalize(true)");
10750 for (i = 0; i < marked_el.
Size(); i++)
10756 for (i = 0; i < marked_el.
Size(); i++)
10765 for (i = 0; i < marked_el.
Size(); i++)
10782 int need_refinement;
10787 need_refinement = 0;
10795 if (
elements[i]->NeedRefinement(v_to_v))
10797 need_refinement = 1;
10802 while (need_refinement == 1);
10809 need_refinement = 0;
10811 if (
boundary[i]->NeedRefinement(v_to_v))
10813 need_refinement = 1;
10817 while (need_refinement == 1);
10848 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
10849 "not supported. Project the NURBS to Nodes first.");
10859 if (!refinements.
Size())
10880 Swap(*mesh2,
false);
10892 const int *fine,
int nfine,
int op)
10894 real_t error = (op == 3) ? std::pow(elem_error[fine[0]],
10895 2.0) : elem_error[fine[0]];
10897 for (
int i = 1; i < nfine; i++)
10899 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
10901 real_t err_fine = elem_error[fine[i]];
10904 case 0: error = std::min(error, err_fine);
break;
10905 case 1: error += err_fine;
break;
10906 case 2: error = std::max(error, err_fine);
break;
10907 case 3: error += std::pow(err_fine, 2.0);
break;
10908 default: MFEM_ABORT(
"Invalid operation.");
10911 return (op == 3) ? std::sqrt(error) : error;
10915 real_t threshold,
int nc_limit,
int op)
10917 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
10918 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
10919 "Project the NURBS to Nodes first.");
10932 for (
int i = 0; i < dt.
Size(); i++)
10934 if (nc_limit > 0 && !level_ok[i]) {
continue; }
10939 if (error < threshold) { derefs.
Append(i); }
10942 if (!derefs.
Size()) {
return false; }
10949 Swap(*mesh2,
false);
10963 int nc_limit,
int op)
10973 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
10980 int nc_limit,
int op)
10983 for (
int i = 0; i < tmp.
Size(); i++)
10985 tmp[i] = elem_error(i);
11028 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
11074#ifdef MFEM_USE_MEMALLOC
11111 for (
int i = 0; i < elem_array.
Size(); i++)
11113 if (elem_array[i]->GetGeometryType() == geom)
11118 elem_vtx.
SetSize(nv*num_elems);
11122 for (
int i = 0; i < elem_array.
Size(); i++)
11128 elem_vtx.
Append(loc_vtx);
11136 for (
int i = 0; i < nelem; i++) { list[i] = i; }
11152 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
11164 default: MFEM_ABORT(
"internal error");
11178 bool noInitialCoarsening =
true;
11179 for (
auto f : initialCoarsening)
11181 noInitialCoarsening = (noInitialCoarsening &&
f == 1);
11184 if (noInitialCoarsening)
11204 bool divisible =
true;
11205 for (
int i=0; i<rf.
Size(); ++i)
11208 divisible = divisible && cf * rf[i] == initialCoarsening[i];
11211 MFEM_VERIFY(divisible,
"Invalid coarsening");
11226 int nonconforming,
int nc_limit)
11236 else if (nonconforming < 0)
11257 for (
int i = 0; i < refinements.
Size(); i++)
11259 el_to_refine[i] = refinements[i].index;
11263 int type, rt = (refinements.
Size() ? refinements[0].GetType() : 7);
11264 if (rt == 1 || rt == 2 || rt == 4)
11268 else if (rt == 3 || rt == 5 || rt == 6)
11286 for (
int i = 0; i < el_to_refine.
Size(); i++)
11288 refinements[i] =
Refinement(el_to_refine[i]);
11295 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
11296 "Please project the NURBS to Nodes first, with SetCurvature().");
11299 MFEM_VERIFY(
ncmesh != NULL ||
dynamic_cast<const ParMesh*
>(
this) == NULL,
11300 "Sorry, converting a conforming ParMesh to an NC mesh is "
11308 (simplices_nonconforming && (
meshgen & 0x1)) )
11321 for (
int i = 0; i <
GetNE(); i++)
11328 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
11340 for (
int i = 0; i <
GetNE(); i++)
11343 bool refine =
false;
11344 for (
int j = 0; j < v.
Size(); j++)
11347 for (
int l = 0; l <
spaceDim; l++)
11352 if (dist <= eps*eps) { refine =
true;
break; }
11363 int nonconforming,
int nc_limit)
11365 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
11367 for (
int i = 0; i <
GetNE(); i++)
11369 if (elem_error[i] > threshold)
11383 int nonconforming,
int nc_limit)
11386 elem_error.
Size());
11387 return RefineByError(tmp, threshold, nonconforming, nc_limit);
11392 int *edge1,
int *edge2,
int *middle)
11395 int v[2][4], v_new,
bisect, t;
11407 bisect = v_to_v(vert[0], vert[1]);
11408 MFEM_ASSERT(
bisect >= 0,
"");
11410 if (middle[
bisect] == -1)
11413 for (
int d = 0; d <
spaceDim; d++)
11438 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
11439 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
11460 bisect = v_to_v(v[1][0], v[1][1]);
11461 MFEM_ASSERT(
bisect >= 0,
"");
11467 else if (edge2[
bisect] == i)
11476 MFEM_ABORT(
"Bisection for now works only for triangles.");
11483 int v[2][4], v_new,
bisect, t;
11493 "TETRAHEDRON element is not marked for refinement.");
11502 for (
int j = 0; j < 3; j++)
11515 int type, old_redges[2], flag;
11518 int new_type, new_redges[2][2];
11521 new_redges[0][0] = 2;
11522 new_redges[0][1] = 1;
11523 new_redges[1][0] = 2;
11524 new_redges[1][1] = 1;
11525 int tr1 = -1, tr2 = -1;
11526 switch (old_redges[0])
11529 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
11534 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
11538 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
11541 switch (old_redges[1])
11544 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
11549 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
11553 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
11560#ifdef MFEM_USE_MEMALLOC
11596 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
11603 int v[2][3], v_new,
bisect, t;
11615 MFEM_ASSERT(
bisect >= 0,
"");
11617 MFEM_ASSERT(v_new != -1,
"");
11621 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
11622 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
11632 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
11638 int *edge1,
int *edge2,
int *middle)
11641 int j, v1[3], v2[3], v3[3], v4[3], v_new[3],
bisect[3];
11650 bisect[0] = v_to_v(v[0],v[1]);
11651 bisect[1] = v_to_v(v[1],v[2]);
11652 bisect[2] = v_to_v(v[0],v[2]);
11655 for (j = 0; j < 3; j++)
11657 if (middle[
bisect[j]] == -1)
11660 for (
int d = 0; d <
spaceDim; d++)
11668 if (edge1[
bisect[j]] == i)
11673 middle[
bisect[j]] = v_new[j];
11677 v_new[j] = middle[
bisect[j]];
11686 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
11687 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
11688 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
11689 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
11723 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
11759 for (
int i = 0; i < elem_geoms.
Size(); i++)
11767 std::map<unsigned, int> mat_no;
11771 for (
int j = 0; j <
elements.Size(); j++)
11774 unsigned code =
elements[j]->GetTransform();
11777 int &matrix = mat_no[code];
11778 if (!matrix) { matrix =
static_cast<int>(mat_no.size()); }
11785 pmats.
SetSize(
Dim,
Dim+1,
static_cast<int>((mat_no.size())));
11788 std::map<unsigned, int>::iterator it;
11789 for (it = mat_no.begin(); it != mat_no.end(); ++it)
11803 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
11804 " geom = " << geom);
11814 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
11823 os <<
"areamesh2\n\n";
11827 os <<
"curved_areamesh2\n\n";
11836 os <<
boundary[i]->GetAttribute();
11837 for (j = 0; j < v.
Size(); j++)
11839 os <<
' ' << v[j] + 1;
11851 for (j = 0; j < v.
Size(); j++)
11853 os <<
' ' << v[j] + 1;
11865 for (j = 1; j <
Dim; j++)
11882 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
11890 os <<
"NETGEN_Neutral_Format\n";
11895 for (j = 0; j <
Dim; j++)
11908 os <<
elements[i]->GetAttribute();
11909 for (j = 0; j < nv; j++)
11911 os <<
' ' << ind[j]+1;
11922 os <<
boundary[i]->GetAttribute();
11923 for (j = 0; j < nv; j++)
11925 os <<
' ' << ind[j]+1;
11937 <<
" 0 0 0 0 0 0 0\n"
11938 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
11940 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
11941 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11945 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
11951 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
11952 for (j = 0; j < nv; j++)
11954 os <<
' ' << ind[j]+1;
11963 os <<
boundary[i]->GetAttribute();
11964 for (j = 0; j < nv; j++)
11966 os <<
' ' << ind[j]+1;
11968 os <<
" 1.0 1.0 1.0 1.0\n";
11977 const std::string &comments)
const
12012 os <<
"\n# mesh curvature GridFunction";
12017 os <<
"\nmfem_mesh_end" << endl;
12024 os << (!set_names && section_delimiter.empty()
12025 ?
"MFEM mesh v1.0\n" :
12026 (!set_names ?
"MFEM mesh v1.2\n" :
"MFEM mesh v1.3\n"));
12028 if (set_names && section_delimiter.empty())
12030 section_delimiter =
"mfem_mesh_end";
12034 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
12037 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12042 "# TETRAHEDRON = 4\n"
12048 os <<
"\ndimension\n" <<
Dim;
12058 os <<
"\nattribute_sets\n";
12070 os <<
"\nbdr_attribute_sets\n";
12095 if (!section_delimiter.empty())
12098 << section_delimiter << endl;
12103 const int version,
const std::string &comments)
const
12105 MFEM_VERIFY(version == 10 || version == 11,
"Invalid NURBS mesh version");
12110 os <<
"MFEM NURBS mesh v" << int(version / 10) <<
"." << version % 10 <<
"\n";
12113 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
12116 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12122 os <<
"\ndimension\n" <<
Dim
12147 int ki = e_to_k[i];
12155 for (
int j=0; j<2; ++j)
12163 const int s = vert[0];
12169 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
12180 ofstream ofs(fname);
12181 ofs.precision(precision);
12185#ifdef MFEM_USE_ADIOS2
12195 "# vtk DataFile Version 3.0\n"
12196 "Generated by MFEM\n"
12198 "DATASET UNSTRUCTURED_GRID\n";
12211 for ( ; j < 3; j++)
12222 for (
int i = 0; i <
Nodes->
FESpace()->GetNDofs(); i++)
12227 os << (*Nodes)(vdofs[0]);
12231 os <<
' ' << (*Nodes)(vdofs[j]);
12233 for ( ; j < 3; j++)
12247 size +=
elements[i]->GetNVertices() + 1;
12252 const int *v =
elements[i]->GetVertices();
12253 const int nv =
elements[i]->GetNVertices();
12257 for (
int j = 0; j < nv; j++)
12259 os <<
' ' << v[perm ? perm[j] : j];
12272 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
12273 "Point meshes should have a single dof per element");
12274 size += dofs.
Size() + 1;
12279 if (!strcmp(fec_name,
"Linear") ||
12280 !strcmp(fec_name,
"H1_0D_P1") ||
12281 !strcmp(fec_name,
"H1_1D_P1") ||
12282 !strcmp(fec_name,
"H1_2D_P1") ||
12283 !strcmp(fec_name,
"H1_3D_P1"))
12287 else if (!strcmp(fec_name,
"Quadratic") ||
12288 !strcmp(fec_name,
"H1_1D_P2") ||
12289 !strcmp(fec_name,
"H1_2D_P2") ||
12290 !strcmp(fec_name,
"H1_3D_P2"))
12296 mfem::err <<
"Mesh::PrintVTK : can not save '"
12297 << fec_name <<
"' elements!" << endl;
12306 for (
int j = 0; j < dofs.
Size(); j++)
12308 os <<
' ' << dofs[j];
12311 else if (order == 2)
12313 const int *vtk_mfem;
12314 switch (
elements[i]->GetGeometryType())
12328 for (
int j = 0; j < dofs.
Size(); j++)
12330 os <<
' ' << dofs[vtk_mfem[j]];
12340 int vtk_cell_type = 5;
12344 os << vtk_cell_type <<
'\n';
12349 <<
"SCALARS material int\n"
12350 <<
"LOOKUP_TABLE default\n";
12353 os <<
elements[i]->GetAttribute() <<
'\n';
12360 bool high_order_output,
12361 int compression_level,
12364 int ref = (high_order_output &&
Nodes)
12367 fname = fname +
".vtu";
12368 std::fstream os(fname.c_str(),std::ios::out);
12369 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"2.2\"";
12370 if (compression_level != 0)
12372 os <<
" compressor=\"vtkZLibDataCompressor\"";
12375 os <<
"<UnstructuredGrid>\n";
12376 PrintVTU(os, ref, format, high_order_output, compression_level, bdr_elements);
12377 os <<
"</Piece>\n";
12378 os <<
"</UnstructuredGrid>\n";
12379 os <<
"</VTKFile>" << std::endl;
12386 bool high_order_output,
12387 int compression_level)
12389 PrintVTU(fname, format, high_order_output, compression_level,
true);
12393 bool high_order_output,
int compression_level,
12401 std::vector<char> buf;
12403 auto get_geom = [&](
int i)
12411 int np = 0, nc_ref = 0;
12412 for (
int i = 0; i < ne; i++)
12421 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
12422 << (high_order_output ? ne : nc_ref) <<
"\">\n";
12425 os <<
"<Points>\n";
12426 os <<
"<DataArray type=\"" << type_str
12427 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
12428 for (
int i = 0; i < ne; i++)
12441 for (
int j = 0; j < pmat.
Width(); j++)
12467 os <<
"</DataArray>" << std::endl;
12468 os <<
"</Points>" << std::endl;
12470 os <<
"<Cells>" << std::endl;
12471 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
12472 << fmt_str <<
"\">" << std::endl;
12474 std::vector<int> offset;
12477 if (high_order_output)
12480 for (
int iel = 0; iel < ne; iel++)
12484 int nnodes = local_connectivity.
Size();
12485 for (
int i=0; i<nnodes; ++i)
12492 offset.push_back(np);
12498 for (
int i = 0; i < ne; i++)
12504 for (
int j = 0; j < RG.
Size(); )
12507 offset.push_back(coff);
12509 for (
int k = 0; k < nv; k++, j++)
12523 os <<
"</DataArray>" << std::endl;
12525 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
12526 << fmt_str <<
"\">" << std::endl;
12528 for (
size_t ii=0; ii<offset.size(); ii++)
12536 os <<
"</DataArray>" << std::endl;
12537 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
12538 << fmt_str <<
"\">" << std::endl;
12540 const int *vtk_geom_map =
12542 for (
int i = 0; i < ne; i++)
12545 uint8_t vtk_cell_type = 5;
12547 vtk_cell_type = vtk_geom_map[geom];
12549 if (high_order_output)
12558 for (
int j = 0; j < RG.
Size(); j += nv)
12568 os <<
"</DataArray>" << std::endl;
12569 os <<
"</Cells>" << std::endl;
12571 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
12572 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
12573 << fmt_str <<
"\">" << std::endl;
12574 for (
int i = 0; i < ne; i++)
12577 if (high_order_output)
12596 os <<
"</DataArray>" << std::endl;
12597 os <<
"</CellData>" << std::endl;
12608 "# vtk DataFile Version 3.0\n"
12609 "Generated by MFEM\n"
12611 "DATASET UNSTRUCTURED_GRID\n";
12616 os <<
"FIELD FieldData 1\n"
12626 np = nc = size = 0;
12627 for (
int i = 0; i <
GetNE(); i++)
12636 os <<
"POINTS " << np <<
" double\n";
12638 for (
int i = 0; i <
GetNE(); i++)
12645 for (
int j = 0; j < pmat.
Width(); j++)
12647 os << pmat(0, j) <<
' ';
12650 os << pmat(1, j) <<
' ';
12662 os << 0.0 <<
' ' << 0.0;
12669 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
12671 for (
int i = 0; i <
GetNE(); i++)
12678 for (
int j = 0; j < RG.
Size(); )
12681 for (
int k = 0; k < nv; k++, j++)
12683 os <<
' ' << np + RG[j];
12689 os <<
"CELL_TYPES " << nc <<
'\n';
12690 for (
int i = 0; i <
GetNE(); i++)
12698 for (
int j = 0; j < RG.
Size(); j += nv)
12700 os << vtk_cell_type <<
'\n';
12704 os <<
"CELL_DATA " << nc <<
'\n'
12705 <<
"SCALARS material int\n"
12706 <<
"LOOKUP_TABLE default\n";
12707 for (
int i = 0; i <
GetNE(); i++)
12715 os << attr <<
'\n';
12722 srand((
unsigned)time(0));
12724 int el0 = (int)floor(
a *
GetNE());
12726 os <<
"SCALARS element_coloring int\n"
12727 <<
"LOOKUP_TABLE default\n";
12728 for (
int i = 0; i <
GetNE(); i++)
12735 os << coloring[i] + 1 <<
'\n';
12741 os <<
"POINT_DATA " << np <<
'\n' << flush;
12744#ifdef MFEM_USE_HDF5
12751#ifdef MFEM_PARALLEL_HDF5
12752 VTKHDF vtkhdf(fname, pmesh->GetComm());
12753 vtkhdf.
SaveMesh(*
this, high_order);
12756 MFEM_ABORT(
"Requires HDF5 library with parallel support enabled");
12761 vtkhdf.
SaveMesh(*
this, high_order);
12768 int delete_el_to_el = (
el_to_el) ? (0) : (1);
12770 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
12773 const int *i_el_el = el_el.
GetI();
12774 const int *j_el_el = el_el.
GetJ();
12779 stack_p = stack_top_p = 0;
12780 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
12782 if (colors[el] != -2)
12788 el_stack[stack_top_p++] = el;
12790 for ( ; stack_p < stack_top_p; stack_p++)
12792 int i = el_stack[stack_p];
12793 int num_nb = i_el_el[i+1] - i_el_el[i];
12794 if (max_num_col < num_nb + 1)
12796 max_num_col = num_nb + 1;
12798 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
12800 int k = j_el_el[j];
12801 if (colors[k] == -2)
12804 el_stack[stack_top_p++] = k;
12812 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
12814 int i = el_stack[stack_p], col;
12816 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
12818 col = colors[j_el_el[j]];
12821 col_marker[col] = 1;
12825 for (col = 0; col < max_num_col; col++)
12826 if (col_marker[col] == 0)
12834 if (delete_el_to_el)
12842 int elem_attr)
const
12844 if (
Dim != 3 &&
Dim != 2) {
return; }
12846 int i, j, k, l, nv, nbe, *v;
12848 os <<
"MFEM mesh v1.0\n";
12852 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12857 "# TETRAHEDRON = 4\n"
12862 os <<
"\ndimension\n" <<
Dim
12866 os << int((elem_attr) ? partitioning[i]+1 :
elements[i]->GetAttribute())
12867 <<
' ' <<
elements[i]->GetGeometryType();
12870 for (j = 0; j < nv; j++)
12882 l = partitioning[l];
12897 os <<
"\nboundary\n" << nbe <<
'\n';
12903 l = partitioning[l];
12906 nv =
faces[i]->GetNVertices();
12907 v =
faces[i]->GetVertices();
12908 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
12909 for (j = 0; j < nv; j++)
12916 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
12917 for (j = nv-1; j >= 0; j--)
12928 nv =
faces[i]->GetNVertices();
12929 v =
faces[i]->GetVertices();
12930 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
12931 for (j = 0; j < nv; j++)
12962 int interior_faces)
12964 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
12965 if (
Dim != 3 &&
Dim != 2) {
return; }
12974 int nv =
elements[i]->GetNVertices();
12975 const int *ind =
elements[i]->GetVertices();
12976 for (
int j = 0; j < nv; j++)
12986 voff[i] = vcount[i-1] + voff[i-1];
12992 vown[i] =
new int[vcount[i]];
13004 int nv =
elements[i]->GetNVertices();
13005 const int *ind =
elements[i]->GetVertices();
13006 for (
int j = 0; j < nv; j++)
13009 vown[ind[j]][vcount[ind[j]]] = i;
13015 vcount[i] = voff[i+1] - voff[i];
13019 for (
int i = 0; i < edge_el.
Size(); i++)
13021 const int *el = edge_el.
GetRow(i);
13024 int k = partitioning[el[0]];
13025 int l = partitioning[el[1]];
13026 if (interior_faces || k != l)
13038 os <<
"areamesh2\n\n" << nbe <<
'\n';
13040 for (
int i = 0; i < edge_el.
Size(); i++)
13042 const int *el = edge_el.
GetRow(i);
13045 int k = partitioning[el[0]];
13046 int l = partitioning[el[1]];
13047 if (interior_faces || k != l)
13052 for (
int j = 0; j < 2; j++)
13053 for (
int s = 0; s < vcount[ev[j]]; s++)
13054 if (vown[ev[j]][s] == el[0])
13056 os <<
' ' << voff[ev[j]]+s+1;
13060 for (
int j = 1; j >= 0; j--)
13061 for (
int s = 0; s < vcount[ev[j]]; s++)
13062 if (vown[ev[j]][s] == el[1])
13064 os <<
' ' << voff[ev[j]]+s+1;
13071 int k = partitioning[el[0]];
13075 for (
int j = 0; j < 2; j++)
13076 for (
int s = 0; s < vcount[ev[j]]; s++)
13077 if (vown[ev[j]][s] == el[0])
13079 os <<
' ' << voff[ev[j]]+s+1;
13089 int nv =
elements[i]->GetNVertices();
13090 const int *ind =
elements[i]->GetVertices();
13091 os << partitioning[i]+1 <<
' ';
13093 for (
int j = 0; j < nv; j++)
13095 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
13096 vown[ind[j]][vcount[ind[j]]] = i;
13103 vcount[i] = voff[i+1] - voff[i];
13109 for (
int k = 0; k < vcount[i]; k++)
13111 for (
int j = 0; j <
Dim; j++)
13121 os <<
"NETGEN_Neutral_Format\n";
13125 for (
int k = 0; k < vcount[i]; k++)
13127 for (
int j = 0; j <
Dim; j++)
13138 int nv =
elements[i]->GetNVertices();
13139 const int *ind =
elements[i]->GetVertices();
13140 os << partitioning[i]+1;
13141 for (
int j = 0; j < nv; j++)
13143 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
13144 vown[ind[j]][vcount[ind[j]]] = i;
13151 vcount[i] = voff[i+1] - voff[i];
13161 int k = partitioning[
faces_info[i].Elem1No];
13162 l = partitioning[l];
13163 if (interior_faces || k != l)
13180 int k = partitioning[
faces_info[i].Elem1No];
13181 l = partitioning[l];
13182 if (interior_faces || k != l)
13184 int nv =
faces[i]->GetNVertices();
13185 const int *ind =
faces[i]->GetVertices();
13187 for (
int j = 0; j < nv; j++)
13188 for (
int s = 0; s < vcount[ind[j]]; s++)
13189 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13191 os <<
' ' << voff[ind[j]]+s+1;
13195 for (
int j = nv-1; j >= 0; j--)
13196 for (
int s = 0; s < vcount[ind[j]]; s++)
13197 if (vown[ind[j]][s] ==
faces_info[i].Elem2No)
13199 os <<
' ' << voff[ind[j]]+s+1;
13206 int k = partitioning[
faces_info[i].Elem1No];
13207 int nv =
faces[i]->GetNVertices();
13208 const int *ind =
faces[i]->GetVertices();
13210 for (
int j = 0; j < nv; j++)
13211 for (
int s = 0; s < vcount[ind[j]]; s++)
13212 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13214 os <<
' ' << voff[ind[j]]+s+1;
13230 int k = partitioning[
faces_info[i].Elem1No];
13231 l = partitioning[l];
13232 if (interior_faces || k != l)
13245 <<
" 0 0 0 0 0 0 0\n"
13246 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
13247 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
13248 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
13249 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
13252 for (
int k = 0; k < vcount[i]; k++)
13253 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
13258 int nv =
elements[i]->GetNVertices();
13259 const int *ind =
elements[i]->GetVertices();
13260 os << i+1 <<
' ' << partitioning[i]+1;
13261 for (
int j = 0; j < nv; j++)
13263 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
13264 vown[ind[j]][vcount[ind[j]]] = i;
13271 vcount[i] = voff[i+1] - voff[i];
13280 int k = partitioning[
faces_info[i].Elem1No];
13281 l = partitioning[l];
13282 if (interior_faces || k != l)
13284 int nv =
faces[i]->GetNVertices();
13285 const int *ind =
faces[i]->GetVertices();
13287 for (
int j = 0; j < nv; j++)
13288 for (
int s = 0; s < vcount[ind[j]]; s++)
13289 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13291 os <<
' ' << voff[ind[j]]+s+1;
13293 os <<
" 1.0 1.0 1.0 1.0\n";
13295 for (
int j = nv-1; j >= 0; j--)
13296 for (
int s = 0; s < vcount[ind[j]]; s++)
13297 if (vown[ind[j]][s] ==
faces_info[i].Elem2No)
13299 os <<
' ' << voff[ind[j]]+s+1;
13301 os <<
" 1.0 1.0 1.0 1.0\n";
13306 int k = partitioning[
faces_info[i].Elem1No];
13307 int nv =
faces[i]->GetNVertices();
13308 const int *ind =
faces[i]->GetVertices();
13310 for (
int j = 0; j < nv; j++)
13311 for (
int s = 0; s < vcount[ind[j]]; s++)
13312 if (vown[ind[j]][s] ==
faces_info[i].Elem1No)
13314 os <<
' ' << voff[ind[j]]+s+1;
13316 os <<
" 1.0 1.0 1.0 1.0\n";
13340 " NURBS mesh is not supported!");
13344 os <<
"MFEM mesh v1.0\n";
13348 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
13353 "# TETRAHEDRON = 4\n"
13358 os <<
"\ndimension\n" <<
Dim
13366 const int *
const i_AF_f = Aface_face.
GetI();
13367 const int *
const j_AF_f = Aface_face.
GetJ();
13369 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
13370 for (
const int * iface = j_AF_f + i_AF_f[iAF];
13371 iface < j_AF_f + i_AF_f[iAF+1];
13374 os << iAF+1 <<
' ';
13407 int *nbea =
new int[na];
13414 for (i = 0; i < na; i++)
13426 for (k = 0; k < vert.
Size(); k++)
13438 for (k = 0; k < vert.
Size(); k++)
13439 if (vn[vert[k]] == 1)
13444 cg[bea*
spaceDim+j] += pointmat(j,k);
13455 for (k = 0; k < vert.
Size(); k++)
13460 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
13477 int *nbea =
new int[na];
13484 for (i = 0; i < na; i++)
13496 for (k = 0; k < vert.
Size(); k++)
13508 for (k = 0; k < vert.
Size(); k++)
13509 if (vn[vert[k]] == 1)
13514 cg[bea*
spaceDim+j] += pointmat(j,k);
13525 for (k = 0; k < vert.
Size(); k++)
13530 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
13546 for (
int i = 0; i <
vertices.Size(); i++)
13548 for (
int j = 0; j <
spaceDim; j++)
13569 "incompatible vector dimensions");
13577 for (
int d = 0; d <
spaceDim; d++)
13597 for (
int i = 0; i <
GetNE(); i++)
13602 for (
int j = 0; j < nv; j++)
13607 for (
int i = 0; i <
GetNBE(); i++)
13612 for (
int j = 0; j < nv; j++)
13618 for (
int i = 0; i < v2v.
Size(); i++)
13623 v2v[i] = num_vert++;
13627 if (num_vert == v2v.
Size()) {
return; }
13629 Vector nodes_by_element;
13634 for (
int i = 0; i <
GetNE(); i++)
13641 for (
int i = 0; i <
GetNE(); i++)
13650 for (
int i = 0; i <
GetNE(); i++)
13655 for (
int j = 0; j < nv; j++)
13660 for (
int i = 0; i <
GetNBE(); i++)
13665 for (
int j = 0; j < nv; j++)
13689 for (
int i = 0; i <
GetNE(); i++)
13702 int num_bdr_elem = 0;
13703 int new_bel_to_edge_nnz = 0;
13704 for (
int i = 0; i <
GetNBE(); i++)
13720 if (num_bdr_elem ==
GetNBE()) {
return; }
13724 Table *new_bel_to_edge = NULL;
13726 new_be_to_face.
Reserve(num_bdr_elem);
13729 new_bel_to_edge =
new Table;
13730 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
13732 for (
int i = 0; i <
GetNBE(); i++)
13737 int row = new_be_to_face.
Size();
13743 int *new_e = new_bel_to_edge->
GetRow(row);
13744 for (
int j = 0; j < ne; j++)
13748 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
13765 for (
int i = 0; i < attribs.
Size(); i++)
13777#ifdef MFEM_USE_MEMALLOC
13804 const int npts = point_mat.
Width();
13805 if (!npts) {
return 0; }
13806 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
13810 if (!
GetNE()) {
return 0; }
13819 min_dist = std::numeric_limits<real_t>::max();
13823 for (
int i = 0; i <
GetNE(); i++)
13827 for (
int k = 0; k < npts; k++)
13830 if (dist < min_dist(k))
13832 min_dist(k) = dist;
13841 for (
int k = 0; k < npts; k++)
13845 int res = inv_tr->
Transform(pt, ips[k]);
13848 elem_ids[k] = e_idx[k];
13852 if (pts_found != npts)
13856 for (
int k = 0; k < npts; k++)
13858 if (elem_ids[k] != -1) {
continue; }
13862 for (
int v = 0; v < elvertices.
Size(); v++)
13864 int vv = elvertices[v];
13866 const int* els = vtoel->
GetRow(vv);
13867 for (
int e = 0; e < ne; e++)
13869 if (els[e] == e_idx[k]) {
continue; }
13871 int res = inv_tr->
Transform(pt, ips[k]);
13874 elem_ids[k] = els[e];
13886 for (
int e = 0; e < neigh.
Size(); e++)
13892 int res = inv_tr->
Transform(pt, ips[k]);
13905 if (inv_trans == NULL) {
delete inv_tr; }
13907 if (warn && pts_found != npts)
13909 MFEM_WARNING((npts-pts_found) <<
" points were not found");
13924 MFEM_VERIFY(
Dim == 2 ||
Dim == 3,
"Only 2D/3D meshes supported right now.");
13925 MFEM_VERIFY(
Dim ==
spaceDim,
"Surface meshes not currently supported.");
13942 skew(0) = std::atan2(J.
Det(), col1 * col2);
13945 ori(0) = std::atan2(J(1,0), J(0,0));
13952 Vector col1, col2, col3;
13963 col1unit *= 1.0/len1;
13964 col2unit *= 1.0/len2;
13965 col3unit *= 1.0/len3;
13971 aspr(0) = len1/std::sqrt(len2*len3),
13972 aspr(1) = len2/std::sqrt(len1*len3);
13975 aspr(2) = std::sqrt(len1/(len2*len3)),
13976 aspr(3) = std::sqrt(len2/(len1*len3));
13979 Vector crosscol12, crosscol13;
13980 col1.
cross3D(col2, crosscol12);
13981 col1.
cross3D(col3, crosscol13);
13982 skew(0) = std::acos(col1unit*col2unit);
13983 skew(1) = std::acos(col1unit*col3unit);
13984 skew(2) = std::atan(len1*volume/(crosscol12*crosscol13));
13990 for (
int d=0; d<
Dim; d++) { rot(d, 0) = col1unit(d); }
13994 rot1 *= col1unit*col2unit;
13996 col1unit.
cross3D(col2unit, rot1);
13998 for (
int d=0; d <
Dim; d++) { rot(d, 1) = rot2(d); }
14001 for (
int d=0; d <
Dim; d++) { rot(d, 2) = rot1(d); }
14002 real_t delta = sqrt(pow(rot(2,1)-rot(1,2), 2.0) +
14003 pow(rot(0,2)-rot(2,0), 2.0) +
14004 pow(rot(1,0)-rot(0,1), 2.0));
14009 for (
int d = 0; d <
Dim; d++) { Iden(d, d) = 1.0; };
14015 MFEM_ABORT(
"Invalid rotation matrix. Contact TMOP Developers.");
14020 ori(0) = (1./
delta)*(rot(2,1)-rot(1,2));
14021 ori(1) = (1./
delta)*(rot(0,2)-rot(2,0));
14022 ori(2) = (1./
delta)*(rot(1,0)-rot(0,1));
14023 ori(3) = std::acos(0.5*(rot.
Trace()-1.0));
14032 entity_to_vertex(entity_to_vertex_)
14034 int geom_offset = 0;
14048 while (geom_offsets[geom+1] <= bytype_entity_id) { geom++; }
14052 const int geom_elem_id = bytype_entity_id - geom_offsets[geom];
14054 return { geom, nv, v };
14059 os <<
"MFEM mesh v1.2\n";
14063 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
14068 "# TETRAHEDRON = 4\n"
14075 os <<
"\ndimension\n" <<
dim;
14081 "invalid MeshPart state");
14084 "invalid MeshPart state");
14085 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
14087 const int bytype_elem_id = have_element_map ?
14092 for (
int i = 0; i < ent.
num_verts; i++)
14094 os <<
' ' << ent.
verts[i];
14104 "invalid MeshPart state");
14107 "invalid MeshPart state");
14110 const int bytype_bdr_id = have_boundary_map ?
14115 for (
int i = 0; i < ent.
num_verts; i++)
14117 os <<
' ' << ent.
verts[i];
14127 os << sdim <<
'\n';
14131 for (
int d = 1; d < sdim; d++)
14144 os <<
"\nmfem_serial_mesh_end\n";
14148 os <<
"\ncommunication_groups\n";
14149 os <<
"number_of_groups " << num_groups <<
"\n\n";
14151 os <<
"# number of entities in each group, followed by ranks in group\n";
14152 for (
int group_id = 0; group_id < num_groups; ++group_id)
14157 for (
int group_member_index = 0; group_member_index < group_size;
14158 ++group_member_index)
14160 os <<
' ' << group_ptr[group_member_index];
14171 MFEM_VERIFY(g2v.
RowSize(0) == 0,
"internal erroor");
14175 MFEM_VERIFY(g2ev.
RowSize(0) == 0,
"internal erroor");
14180 MFEM_VERIFY(g2tv.
RowSize(0) == 0,
"internal erroor");
14181 MFEM_VERIFY(g2qv.
RowSize(0) == 0,
"internal erroor");
14182 const int total_shared_faces =
14184 os <<
"total_shared_faces " << total_shared_faces <<
'\n';
14186 os <<
"\n# group 0 has no shared entities\n";
14187 for (
int gr = 1; gr < num_groups; gr++)
14190 const int nv = g2v.
RowSize(gr);
14191 const int *sv = g2v.
GetRow(gr);
14192 os <<
"\n# group " << gr <<
"\nshared_vertices " << nv <<
'\n';
14193 for (
int i = 0; i < nv; i++)
14195 os << sv[i] <<
'\n';
14200 const int ne = g2ev.
RowSize(gr)/2;
14201 const int *se = g2ev.
GetRow(gr);
14202 os <<
"\nshared_edges " << ne <<
'\n';
14203 for (
int i = 0; i < ne; i++)
14205 const int *v = se + 2*i;
14206 os << v[0] <<
' ' << v[1] <<
'\n';
14211 const int nt = g2tv.
RowSize(gr)/3;
14212 const int *st = g2tv.
GetRow(gr);
14213 const int nq = g2qv.
RowSize(gr)/4;
14215 os <<
"\nshared_faces " << nt+nq <<
'\n';
14216 for (
int i = 0; i < nt; i++)
14219 const int *v = st + 3*i;
14220 for (
int j = 0; j < 3; j++) { os <<
' ' << v[j]; }
14223 for (
int i = 0; i < nq; i++)
14226 const int *v =
sq + 4*i;
14227 for (
int j = 0; j < 4; j++) { os <<
' ' << v[j]; }
14234 os <<
"\nmfem_mesh_end" << endl;
14251 "invalid MeshPart state");
14254 "invalid MeshPart state");
14256 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
14258 const int bytype_elem_id = have_element_map ?
14269 static_cast<Tetrahedron*
>(el)->SetRefinementFlag(ref_flag);
14271 mesh->AddElement(el);
14279 "invalid MeshPart state");
14282 "invalid MeshPart state");
14285 const int bytype_bdr_id = have_boundary_map ?
14291 mesh->AddBdrElement(bdr);
14298 MFEM_ASSERT(!
nodes,
"invalid MeshPart state");
14299 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
14307 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
14309 mesh->AddVertex(0., 0., 0.);
14314 mesh->FinalizeTopology(
false);
14322 const int *partitioning_,
14354 for (
int i = 0; i < boundary_to_part.
Size(); i++)
14356 int face, o, el1, el2;
14359 boundary_to_part[i] =
14365 for (
int i = 0; i < boundary_to_part.
Size(); i++)
14374 for (
int i = 0; i < boundary_to_part.
Size(); i++)
14389 delete vert_element;
14396 MFEM_VERIFY(0 <= part_id && part_id < num_parts,
14397 "invalid part_id = " << part_id
14398 <<
", num_parts = " << num_parts);
14431 mesh_part.
nodes.reset(
nullptr);
14433 mesh_part.
mesh.reset(
nullptr);
14441 int geom_marker = 0, num_geom = 0;
14442 for (
int i = 0; i < num_elems; i++)
14448 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
14450 "overflow in 'entity_to_vertex[geom]', geom: "
14475 if ((geom_marker & (1 << geom)) == 0)
14477 geom_marker |= (1 << geom);
14492 offsets[g] = offset;
14496 for (
int i = 0; i < num_elems; i++)
14507 geom_marker = 0; num_geom = 0;
14508 for (
int i = 0; i < num_bdr_elems; i++)
14514 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
14516 "overflow in 'entity_to_vertex[geom]', geom: "
14520 if ((geom_marker & (1 << geom)) == 0)
14522 geom_marker |= (1 << geom);
14533 offsets[g] = offset;
14537 for (
int i = 0; i < num_bdr_elems; i++)
14548 std::unordered_set<int> vertex_set;
14549 for (
int i = 0; i < num_elems; i++)
14555 vertex_set.insert(v, v + nv);
14557 vertex_loc_to_glob.
SetSize(
static_cast<int>(vertex_set.size()));
14558 std::copy(vertex_set.begin(), vertex_set.end(),
14559 vertex_loc_to_glob.
begin());
14561 vertex_loc_to_glob.
Sort();
14571 for (
int i = 0; i < vert_array.
Size(); i++)
14573 const int glob_id = vert_array[i];
14574 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
14575 MFEM_ASSERT(loc_id >= 0,
"internal error: global vertex id not found");
14576 vert_array[i] = loc_id;
14583 MFEM_VERIFY(numeric_limits<int>::max()/sdim >= vertex_loc_to_glob.
Size(),
14584 "overflow in 'vertex_coordinates', num_vertices = "
14585 << vertex_loc_to_glob.
Size() <<
", sdim = " << sdim);
14587 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
14590 for (
int d = 0; d < sdim; d++)
14607 mesh_part.
mesh->NewNodes(*mesh_part.
nodes,
false);
14629 std::unordered_set<int> face_set;
14632 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14634 const int glob_elem_id = elem_list[loc_elem_id];
14635 const int nfaces = elem_to_face.
RowSize(glob_elem_id);
14636 const int *faces = elem_to_face.
GetRow(glob_elem_id);
14637 face_set.insert(faces, faces + nfaces);
14641 for (
int glob_face_id : face_set)
14645 if (el[1] < 0) {
continue; }
14648 MFEM_ASSERT(el[0] == part_id || el[1] == part_id,
"internal error");
14649 if (el[0] != part_id || el[1] != part_id)
14652 const int group_id = groups.
Insert(group);
14656 shared_faces.
Sort();
14667 std::unordered_set<int> edge_set;
14670 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14672 const int glob_elem_id = elem_list[loc_elem_id];
14673 const int nedges = elem_to_edge.
RowSize(glob_elem_id);
14674 const int *edges = elem_to_edge.
GetRow(glob_elem_id);
14675 edge_set.insert(edges, edges + nedges);
14679 for (
int glob_edge_id : edge_set)
14685 for (
int j = 0; j < nelem; j++)
14691 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
14692 if (group.
Size() > 1)
14694 const int group_id = groups.
Insert(group);
14698 shared_edges.
Sort();
14709 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
14712 const int glob_vertex_id = vertex_loc_to_glob[i];
14717 for (
int j = 0; j < nelem; j++)
14723 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
14724 if (group.
Size() > 1)
14726 const int group_id = groups.
Insert(group);
14733 const int num_groups = groups.
Size();
14739 Table &group__shared_vertex_to_vertex =
14741 group__shared_vertex_to_vertex.
MakeI(num_groups);
14742 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
14744 const int group_id = shared_verts[sv].two;
14747 group__shared_vertex_to_vertex.
MakeJ();
14748 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
14750 const int glob_vertex_id = shared_verts[sv].one;
14751 const int group_id = shared_verts[sv].two;
14752 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(glob_vertex_id);
14753 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
14754 group__shared_vertex_to_vertex.
AddConnection(group_id, loc_vertex_id);
14756 group__shared_vertex_to_vertex.
ShiftUpI();
14761 Table &group__shared_edge_to_vertex =
14763 group__shared_edge_to_vertex.
MakeI(num_groups);
14764 for (
int se = 0; se < shared_edges.
Size(); se++)
14766 const int group_id = shared_edges[se].two;
14769 group__shared_edge_to_vertex.
MakeJ();
14771 for (
int se = 0; se < shared_edges.
Size(); se++)
14773 const int glob_edge_id = shared_edges[se].one;
14774 const int group_id = shared_edges[se].two;
14775 const int *v = edge_to_vertex.
GetRow(glob_edge_id);
14776 for (
int i = 0; i < 2; i++)
14778 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(v[i]);
14779 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
14780 group__shared_edge_to_vertex.
AddConnection(group_id, loc_vertex_id);
14783 group__shared_edge_to_vertex.
ShiftUpI();
14790 Table &group__shared_tria_to_vertex =
14792 Table &group__shared_quad_to_vertex =
14795 group__shared_tria_to_vertex.
MakeI(num_groups);
14796 group__shared_quad_to_vertex.
MakeI(num_groups);
14797 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
14799 const int glob_face_id = shared_faces[sf].one;
14800 const int group_id = shared_faces[sf].two;
14805 group__shared_tria_to_vertex.
MakeJ();
14806 group__shared_quad_to_vertex.
MakeJ();
14807 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
14809 const int glob_face_id = shared_faces[sf].one;
14810 const int group_id = shared_faces[sf].two;
14846 for (
int i = 0; i < vertex_ids.
Size(); i++)
14848 const int glob_id = vertex_ids[i];
14849 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
14850 MFEM_ASSERT(loc_id >= 0,
"internal error");
14851 vertex_ids[i] = loc_id;
14854 AddConnections(group_id, vertex_ids, vertex_ids.
Size());
14856 group__shared_tria_to_vertex.
ShiftUpI();
14857 group__shared_quad_to_vertex.
ShiftUpI();
14861std::unique_ptr<FiniteElementSpace>
14869 return std::unique_ptr<FiniteElementSpace>(
14871 global_fespace.
FEColl(),
14876std::unique_ptr<GridFunction>
14881 std::unique_ptr<GridFunction> local_gf(
new GridFunction(&local_fespace));
14889 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14891 const int glob_elem_id = elem_list[loc_elem_id];
14898 local_gf->SetSubVector(lvdofs, loc_vals);
14912 "mixed meshes are not supported!");
14913 MFEM_ASSERT(
mesh->
GetNodes(),
"meshes without nodes are not supported!");
14926 Compute(
nodes, d_mt);
14936 const int vdim = fespace->
GetVDim();
14937 const int NE = fespace->
GetNE();
14938 const int ND = fe->
GetDof();
14941 unsigned eval_flags = 0;
14943 Device::GetDeviceMemoryType();
14966 qi->DisableTensorProducts(!use_tensor_products);
14974 Vector Enodes(vdim*ND*NE, my_d_mt);
14975 elem_restr->Mult(
nodes, Enodes);
14976 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
14996 const int vdim = fespace->
GetVDim();
15012 unsigned eval_flags = 0;
15059 V(1) = s * ((ip.
y + layer) / n);
15064 V(2) = s * ((ip.
z + layer) / n);
15073 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
15077 int nvy = (closed) ? (ny) : (ny + 1);
15078 int nvt = mesh->
GetNV() * nvy;
15087 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
15092 for (
int i = 0; i < mesh->
GetNV(); i++)
15095 for (
int j = 0; j < nvy; j++)
15097 vc[1] = sy * (
real_t(j) / ny);
15103 for (
int i = 0; i < mesh->
GetNE(); i++)
15108 for (
int j = 0; j < ny; j++)
15111 qv[0] = vert[0] * nvy + j;
15112 qv[1] = vert[1] * nvy + j;
15113 qv[2] = vert[1] * nvy + (j + 1) % nvy;
15114 qv[3] = vert[0] * nvy + (j + 1) % nvy;
15120 for (
int i = 0; i < mesh->
GetNBE(); i++)
15125 for (
int j = 0; j < ny; j++)
15128 sv[0] = vert[0] * nvy + j;
15129 sv[1] = vert[0] * nvy + (j + 1) % nvy;
15145 for (
int i = 0; i < mesh->
GetNE(); i++)
15151 sv[0] = vert[0] * nvy;
15152 sv[1] = vert[1] * nvy;
15156 sv[0] = vert[1] * nvy + ny;
15157 sv[1] = vert[0] * nvy + ny;
15173 string cname = name;
15174 if (cname ==
"Linear")
15178 else if (cname ==
"Quadratic")
15182 else if (cname ==
"Cubic")
15186 else if (!strncmp(name,
"H1_", 3))
15190 else if (!strncmp(name,
"L2_T", 4))
15194 else if (!strncmp(name,
"L2_", 3))
15201 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
15213 for (
int i = 0; i < mesh->
GetNE(); i++)
15216 for (
int j = ny-1; j >= 0; j--)
15233 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
15238 int nvt = mesh->
GetNV() * nvz;
15243 bool wdgMesh =
false;
15244 bool hexMesh =
false;
15248 for (
int i = 0; i < mesh->
GetNV(); i++)
15252 for (
int j = 0; j < nvz; j++)
15254 vc[2] = sz * (
real_t(j) / nz);
15260 for (
int i = 0; i < mesh->
GetNE(); i++)
15270 for (
int j = 0; j < nz; j++)
15273 pv[0] = vert[0] * nvz + j;
15274 pv[1] = vert[1] * nvz + j;
15275 pv[2] = vert[2] * nvz + j;
15276 pv[3] = vert[0] * nvz + (j + 1) % nvz;
15277 pv[4] = vert[1] * nvz + (j + 1) % nvz;
15278 pv[5] = vert[2] * nvz + (j + 1) % nvz;
15285 for (
int j = 0; j < nz; j++)
15288 hv[0] = vert[0] * nvz + j;
15289 hv[1] = vert[1] * nvz + j;
15290 hv[2] = vert[2] * nvz + j;
15291 hv[3] = vert[3] * nvz + j;
15292 hv[4] = vert[0] * nvz + (j + 1) % nvz;
15293 hv[5] = vert[1] * nvz + (j + 1) % nvz;
15294 hv[6] = vert[2] * nvz + (j + 1) % nvz;
15295 hv[7] = vert[3] * nvz + (j + 1) % nvz;
15297 mesh3d->
AddHex(hv, attr);
15301 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
15302 << geom <<
"\'" << endl;
15308 for (
int i = 0; i < mesh->
GetNBE(); i++)
15313 for (
int j = 0; j < nz; j++)
15316 qv[0] = vert[0] * nvz + j;
15317 qv[1] = vert[1] * nvz + j;
15318 qv[2] = vert[1] * nvz + (j + 1) % nvz;
15319 qv[3] = vert[0] * nvz + (j + 1) % nvz;
15328 for (
int i = 0; i < mesh->
GetNE(); i++)
15339 tv[0] = vert[0] * nvz;
15340 tv[1] = vert[2] * nvz;
15341 tv[2] = vert[1] * nvz;
15345 tv[0] = vert[0] * nvz + nz;
15346 tv[1] = vert[1] * nvz + nz;
15347 tv[2] = vert[2] * nvz + nz;
15355 qv[0] = vert[0] * nvz;
15356 qv[1] = vert[3] * nvz;
15357 qv[2] = vert[2] * nvz;
15358 qv[3] = vert[1] * nvz;
15362 qv[0] = vert[0] * nvz + nz;
15363 qv[1] = vert[1] * nvz + nz;
15364 qv[2] = vert[2] * nvz + nz;
15365 qv[3] = vert[3] * nvz + nz;
15371 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
15372 << geom <<
"\'" << endl;
15378 if ( hexMesh && wdgMesh )
15382 else if ( hexMesh )
15386 else if ( wdgMesh )
15399 string cname = name;
15400 if (cname ==
"Linear")
15404 else if (cname ==
"Quadratic")
15408 else if (cname ==
"Cubic")
15412 else if (!strncmp(name,
"H1_", 3))
15416 else if (!strncmp(name,
"L2_T", 4))
15420 else if (!strncmp(name,
"L2_", 3))
15427 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
15439 for (
int i = 0; i < mesh->
GetNE(); i++)
15442 for (
int j = nz-1; j >= 0; j--)
15477 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
15478 <<
" 0 0 " << i <<
" -1 0\n";
15485 real_t mid[3] = {0, 0, 0};
15486 for (
int j = 0; j < 2; j++)
15488 for (
int k = 0; k <
spaceDim; k++)
15494 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
15495 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
Node::Index insert_node(Float length=1)
Arc::Index insert_arc(Node::Index i, Node::Index j, Float w=1, Float b=1)
uint rank(Node::Index i) const
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int FindSorted(const T &el) const
Do bisection search for 'el' in a sorted array; return -1 if not found.
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
void DeleteAll()
Delete the whole array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
T * GetData()
Returns the data.
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
T * 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 is the named attribute set is present.
Array< int > GetAttributeSetMarker(const std::string &set_name) const
Return a marker array corresponding to a named attribute set.
void Print(std::ostream &out=mfem::out, int width=-1) const
Print the contents of the container to an output stream.
void Copy(AttributeSets ©) const
Create a copy of the internal data to the provided copy.
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Piecewise-(bi)cubic continuous finite elements.
int NumberOfEntries() const
Data type dense matrix using column-major storage.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
const real_t * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
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.
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
const FiniteElementCollection * FEColl() const
Mesh * GetMesh() const
Returns the mesh.
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
int GetVDim() const
Returns the vector dimension of the finite element space.
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified face.
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
virtual int GetMaxElementOrder() const
Return the maximum polynomial order over all elements.
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified element.
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering f_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Abstract class for all finite elements.
int GetDim() const
Returns the reference space dimension for the finite element.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Vector X
Mapped (physical) coordinates of all quadrature points.
const IntegrationRule * IntRule
Vector detJ
Determinants of the Jacobians at all quadrature points.
Vector J
Jacobians of the element transformations at all quadrature points.
GeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, MemoryType d_mt=MemoryType::DEFAULT)
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
static const int Dimension[NumGeom]
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
static const char * Name[NumGeom]
static const int NumVerts[NumGeom]
const IntegrationRule * GetVertices(int GeomType) const
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType.
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
static int GetInverseOrientation(Type geom_type, int orientation)
Return the inverse of the given orientation for the specified geometry type.
static const int DimStart[MaxDim+2]
Class for grid function - Vector with associated FE space.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec_owned and fes.
FiniteElementSpace * FESpace()
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int VectorDim() const
Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying fes.
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values 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.
void ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
int GetPatchBdrAttribute(int i) const
Return the attribute of patch boundary element i, for a NURBS mesh.
int GetElementToEdgeTable(Table &)
void 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
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 (...)
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".
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)
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
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)
For NURBS meshes, remove the knots in kv, for each direction.
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr) const
Builds the transformation defining the i-th edge element in EdTr. EdTr must be allocated in advance a...
void GetElementData(const Array< Element * > &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
void PrintSurfaces(const Table &Aface_face, std::ostream &os) const
Print set of disjoint surfaces:
bool IsSlaveFace(const FaceInfo &fi) const
void FreeElement(Element *E)
Array< Element * > boundary
virtual void MarkTetMeshForRefinement(const DSTable &v_to_v)
void GetPointMatrix(int i, DenseMatrix &pointmat) const
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
void GetLocalTriToTetTransformation(IsoparametricTransformation &loc, int i) const
virtual void PrintXG(std::ostream &os=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
NCMesh * ncmesh
Optional nonconforming mesh extension.
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
virtual void NURBSUniformRefinement(int rf=2, real_t tol=1.0e-12)
Refine NURBS mesh, with an optional refinement factor, generally anisotropic.
void DebugDump(std::ostream &os) const
Output an NCMesh-compatible debug dump.
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
bool RefineByError(const Array< real_t > &elem_error, real_t threshold, int nonconforming=-1, int nc_limit=0)
static Mesh MakeCartesian3DWith24TetsPerHex(int nx, int ny, int nz, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz*24 tetrahedrons.
STable3D * GetFacesTable()
Table * GetFaceEdgeTable() const
void EnsureNCMesh(bool simplices_nonconforming=false)
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
virtual MFEM_DEPRECATED void ReorientTetMesh()
void PrintVTK(std::ostream &os)
virtual void MarkNamedBoundaries(const std::string &set_name, Array< int > &bdr_marker) const
Mark boundary attributes in the named set.
void 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.
void GetLocalQuadToPyrTransformation(IsoparametricTransformation &loc, int i) const
void AddHexAsWedges(const int *vi, int attr=1)
Adds 2 wedges to the mesh by splitting a hexahedron given by 8 vertices vi.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Element * ReadElement(std::istream &input)
virtual bool HasBoundaryElements() const
Checks if the mesh has boundary elements.
void ScaleSubdomains(real_t sf)
void GetVertexToVertexTable(DSTable &) const
void GetLocalQuadToHexTransformation(IsoparametricTransformation &loc, int i) const
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Geometry::Type GetTypicalFaceGeometry() const
If the local mesh is not empty, return GetFaceGeometry(0); otherwise return a typical face geometry p...
Geometry::Type GetElementBaseGeometry(int i) const
virtual void UnmarkInternalBoundaries(Array< int > &bdr_marker, bool excl=true) const
Unmark boundary attributes of internal boundaries.
@ LocalSlaveNonconforming
@ SharedSlaveNonconforming
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones.
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()
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)
For NURBS meshes, insert the new knots in kv, for each direction.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
void OnMeshUpdated(Mesh *mesh)
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
const CoarseFineTransformations & GetRefinementTransforms() const
int Dimension() const
Return the dimension of the NCMesh.
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 ...
void SetCoordsFromPatches(Vector &Nodes)
Set FE coordinates in Nodes, using data from patches, and erase patches.
void Coarsen(int cf=2, real_t tol=1.0e-12)
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 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
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
GeometryRefiner GlobGeometryRefiner
int FindRoots(const Vector &z, Vector &x)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
real_t rand_real()
Generate a random real_t number in the interval [0,1) using rand().
void add(const Vector &v1, const Vector &v2, Vector &v)
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
void ShiftRight(int &a, int &b, int &c)
MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
Mesh * Extrude1D(Mesh *mesh, const int ny, const real_t sy, const bool closed)
Extrude a 1D mesh.
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
void 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.
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