34#include <unordered_map>
35#include <unordered_set>
39#if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
44#if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
49 int*,
int*,
int*,
int*,
int*,
idxtype*);
51 int*,
int*,
int*,
int*,
int*,
idxtype*);
53 int*,
int*,
int*,
int*,
int*,
idxtype*);
118 return sqrt((d_hat * d_hat) / (dir * dir));
157 if (coord[d] < min(d)) { min(d) = coord[d]; }
158 if (coord[d] > max(d)) { max(d) = coord[d]; }
164 const bool use_boundary =
false;
173 for (
int i = 0; i < ne; i++)
190 for (
int j = 0; j < pointmat.
Width(); j++)
192 for (
int d = 0; d < pointmat.
Height(); d++)
194 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
195 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
217 h_max = kappa_max = -h_min;
218 if (
dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
226 if (Vh) { (*Vh)(i) = h; }
227 if (Vk) { (*Vk)(i) =
kappa; }
229 if (h < h_min) { h_min = h; }
230 if (h > h_max) { h_max = h; }
244 if (!num_elems_by_geom[g]) {
continue; }
245 if (!first) { os <<
" + "; }
253 real_t h_min, h_max, kappa_min, kappa_max;
255 os <<
"Mesh Characteristics:";
260 num_elems_by_geom = 0;
261 for (
int i = 0; i <
GetNE(); i++)
272 <<
"Number of vertices : " <<
GetNV() <<
'\n'
273 <<
"Number of elements : " <<
GetNE() <<
'\n'
274 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n';
279 <<
"Number of vertices : " <<
GetNV() <<
'\n'
280 <<
"Number of elements : " <<
GetNE() <<
'\n'
281 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
282 <<
"h_min : " << h_min <<
'\n'
283 <<
"h_max : " << h_max <<
'\n';
288 <<
"Number of vertices : " <<
GetNV() <<
'\n'
289 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
290 <<
"Number of elements : " <<
GetNE() <<
" -- ";
293 <<
"Number of bdr elem : " <<
GetNBE() <<
'\n'
295 <<
"h_min : " << h_min <<
'\n'
296 <<
"h_max : " << h_max <<
'\n'
297 <<
"kappa_min : " << kappa_min <<
'\n'
298 <<
"kappa_max : " << kappa_max <<
'\n';
303 num_bdr_elems_by_geom = 0;
304 for (
int i = 0; i <
GetNBE(); i++)
309 num_faces_by_geom = 0;
316 <<
"Number of vertices : " <<
GetNV() <<
'\n'
317 <<
"Number of edges : " <<
GetNEdges() <<
'\n'
318 <<
"Number of faces : " <<
GetNFaces() <<
" -- ";
321 <<
"Number of elements : " <<
GetNE() <<
" -- ";
324 <<
"Number of bdr elem : " <<
GetNBE() <<
" -- ";
328 <<
"h_min : " << h_min <<
'\n'
329 <<
"h_max : " << h_max <<
'\n'
330 <<
"kappa_min : " << kappa_min <<
'\n'
331 <<
"kappa_max : " << kappa_max <<
'\n';
333 os <<
'\n' << std::flush;
349 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
352 MFEM_ABORT(
"Unknown element type");
381 for (
int j = 0; j < n; j++)
383 pm(k,j) = nodes(vdofs[n*k+j]);
409 int nv =
elements[i]->GetNVertices();
410 const int *v =
elements[i]->GetVertices();
415 for (
int j = 0; j < nv; j++)
417 pm(k, j) = nodes(k*n+v[j]);
431 for (
int j = 0; j < n; j++)
433 pm(k,j) = nodes(vdofs[n*k+j]);
467 for (
int j = 0; j < n; j++)
469 pm(k,j) = nodes(vdofs[n*k+j]);
476 int elem_id, face_info;
492 "Mesh requires nodal Finite Element.");
501 ElTr->
SetFE(face_el);
523 const int *v = (
Dim == 1) ? &FaceNo :
faces[FaceNo]->GetVertices();
524 const int nv = (
Dim == 1) ? 1 :
faces[FaceNo]->GetNVertices();
528 for (
int j = 0; j < nv; j++)
548 for (
int j = 0; j < n; j++)
550 pm(i, j) = nodes(vdofs[n*i+j]);
569 "Mesh requires nodal Finite Element.");
599 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
616 for (
int j = 0; j < nv; j++)
636 for (
int j = 0; j < n; j++)
638 pm(i, j) = nodes(vdofs[n*i+j]);
641 EdTr->
SetFE(edge_el);
645 MFEM_ABORT(
"Not implemented.");
685 for (
int j = 0; j < 2; j++)
687 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
688 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
705 for (
int j = 0; j < 2; j++)
707 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
708 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
727 for (
int j = 0; j < 3; j++)
730 locpm(0, j) = vert.
x;
731 locpm(1, j) = vert.
y;
732 locpm(2, j) = vert.
z;
744 MFEM_VERIFY(i < 128,
"Local face index " << i/64
745 <<
" is not a triangular face of a wedge.");
753 for (
int j = 0; j < 3; j++)
756 locpm(0, j) = vert.
x;
757 locpm(1, j) = vert.
y;
758 locpm(2, j) = vert.
z;
769 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
770 <<
" is not a triangular face of a pyramid.");
778 for (
int j = 0; j < 3; j++)
781 locpm(0, j) = vert.
x;
782 locpm(1, j) = vert.
y;
783 locpm(2, j) = vert.
z;
800 for (
int j = 0; j < 4; j++)
803 locpm(0, j) = vert.
x;
804 locpm(1, j) = vert.
y;
805 locpm(2, j) = vert.
z;
817 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
818 <<
" is not a quadrilateral face of a wedge.");
824 for (
int j = 0; j < 4; j++)
827 locpm(0, j) = vert.
x;
828 locpm(1, j) = vert.
y;
829 locpm(2, j) = vert.
z;
840 MFEM_VERIFY(i < 64,
"Local face index " << i/64
841 <<
" is not a quadrilateral face of a pyramid.");
847 for (
int j = 0; j < 4; j++)
850 locpm(0, j) = vert.
x;
851 locpm(1, j) = vert.
y;
852 locpm(2, j) = vert.
z;
951 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
952 "face type " << face_type
953 <<
" and element type " << elem_type <<
"\n");
972 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
973 "face type " << face_type
974 <<
" and element type " << elem_type <<
"\n");
1006 FElTr.
Elem1 = &ElTr1;
1019 { MFEM_ABORT(
"NURBS mesh not supported!"); }
1022 FElTr.
Elem2 = &ElTr2;
1070 mfem::out <<
"\nInternal error: face id = " << FaceNo
1071 <<
", dist = " << dist <<
'\n';
1073 MFEM_ABORT(
"internal error");
1131 const FaceInfo &fi,
bool is_ghost)
const
1133#ifdef MFEM_THREAD_SAFE
1138 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1150 std::swap(composition(0,0), composition(0,1));
1151 std::swap(composition(1,0), composition(1,1));
1208 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1297 res.
Elem1No = element[0].index;
1298 res.Elem2No = element[1].index;
1299 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1300 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1301 res.NCFace = ncface;
1304 res.Elem1No = element[0].index;
1305 res.Elem2No = element[1].index;
1306 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1307 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1308 res.NCFace = ncface;
1311 res.Elem1No = element[0].index;
1312 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1315 res.Elem1No = element[0].index;
1316 res.Elem2No = -1 - element[1].index;
1317 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1318 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1321 res.Elem1No = element[0].index;
1322 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1325 res.Elem1No = element[0].index;
1326 res.Elem2No = -1 - element[1].index;
1327 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1328 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1333 res.Elem1No = element[0].index;
1334 res.Elem2No = -1 - element[1].index;
1335 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1336 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1344 os <<
"face topology=";
1354 os <<
"Non-conforming";
1361 os <<
"element[0].location=";
1375 os <<
"element[1].location=";
1389 os <<
"element[0].conformity=";
1406 os <<
"element[1].conformity=";
1423 os <<
"element[0].index=" << info.
element[0].
index <<
'\n'
1424 <<
"element[1].index=" << info.
element[1].
index <<
'\n'
1429 <<
"ncface=" << info.
ncface << std::endl;
1461 return faces[Face]->GetGeometryType();
1464 const int nc_face_id =
faces_info[Face].NCFace;
1466 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1557 for (
int i = 0; i <
faces.Size(); i++)
1585#ifdef MFEM_USE_MEMALLOC
1609 for (
int i = 0; i < attribs.
Size(); i++)
1618 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1622 for (
int i = 0; i < attribs.
Size(); i++)
1631 MFEM_WARNING(
"Non-positive attributes in the domain!");
1653static void CheckEnlarge(
Array<T> &array,
int size)
1655 if (size >= array.Size()) { array.SetSize(size + 1); }
1678 "invalid 'coords' size: " << coords.
Size());
1691 for (
int j = 0; j < 3; j++)
1693 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1702 for (
int i = 0; i < nverts; i++)
1705 for (
int j = 0; j <
dim; j++)
1759 int vi[4] = {v1, v2, v3, v4};
1766#ifdef MFEM_USE_MEMALLOC
1806int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1811 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1824 static const int hex_to_tet[6][4] =
1826 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1827 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1831 for (
int i = 0; i < 6; i++)
1833 for (
int j = 0; j < 4; j++)
1835 ti[j] = vi[hex_to_tet[i][j]];
1843 static const int hex_to_wdg[2][6] =
1845 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1849 for (
int i = 0; i < 2; i++)
1851 for (
int j = 0; j < 6; j++)
1853 ti[j] = vi[hex_to_wdg[i][j]];
1861 static const int hex_to_pyr[6][5] =
1863 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
1864 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
1868 for (
int i = 0; i < 6; i++)
1870 for (
int j = 0; j < 5; j++)
1872 ti[j] = vi[hex_to_pyr[i][j]];
1881 static const int quad_to_tri[4][2] =
1883 {0, 1}, {1, 2}, {2, 3}, {3, 0}
1889 ti[2] = elem_center_index;
1890 for (
int i = 0; i < num_faces; i++)
1892 for (
int j = 0; j < 2; j++)
1894 ti[j] = vi[quad_to_tri[i][j]];
1903 static const int quad_faces[4][2] =
1905 {0, 1}, {1, 2}, {2, 3}, {3, 0}
1909 for (
int i = 0; i < 4; i++)
1919 vnew[0] = px(0)*(1-r)*(1-
s) + px(1)*(r)*(1-
s) + px(2)*r*
s + px(3)*(1-r)*
s;
1920 vnew[1] = py(0)*(1-r)*(1-
s) + py(1)*(r)*(1-
s) + py(2)*r*
s + py(3)*(1-r)*
s;
1925 vnew[0] = px(0)*(1-r)*(1-
s) + px(1)*(r)*(1-
s) + px(2)*r*
s + px(3)*(1-r)*
s;
1926 vnew[1] = py(0)*(1-r)*(1-
s) + py(1)*(r)*(1-
s) + py(2)*r*
s + py(3)*(1-r)*
s;
1931 vnew[0] = px(0)*(1-r)*(1-
s) + px(1)*(r)*(1-
s) + px(2)*r*
s + px(3)*(1-r)*
s;
1932 vnew[1] = py(0)*(1-r)*(1-
s) + py(1)*(r)*(1-
s) + py(2)*r*
s + py(3)*(1-r)*
s;
1937 vnew[0] = px(0)*(1-r)*(1-
s) + px(1)*(r)*(1-
s) + px(2)*r*
s + px(3)*(1-r)*
s;
1938 vnew[1] = py(0)*(1-r)*(1-
s) + py(1)*(r)*(1-
s) + py(2)*r*
s + py(3)*(1-r)*
s;
1942 static const int quad_faces_new[4][2] =
1944 { 1, 0}, { 2, 1}, { 3, 2}, { 0, 3}
1948 for (
int i = 0; i < num_faces; i++)
1950 for (
int j = 0; j < 2; j++)
1952 ti[j] = vi[quad_faces[i][j]];
1953 ti[j+2] = vnew_index[quad_faces_new[i][j]];
1961 std::map<std::array<int, 4>,
int> &hex_face_verts,
1967 return std::array<int, 4> {v[0], v[1], v[2], v[3]};
1971 static const int hex_to_tet[6][4] =
1973 { 0, 1, 2, 3 }, { 1, 2, 6, 5 }, { 5, 4, 7, 6},
1974 { 0, 1, 5, 4 }, { 2, 3, 7, 6 }, { 0,3, 7, 4}
1982 static const int tet_face[4][2] =
1984 {0, 1}, {1, 2}, {3, 2}, {3, 0}
1987 for (
int i = 0; i < num_faces; i++)
1989 for (
int j = 0; j < 4; j++)
1991 flist[j] = vi[hex_to_tet[i][j]];
1993 int face_center_index;
1995 auto t = get4arraysorted(flist);
1996 auto it = hex_face_verts.find(
t);
1997 if (it == hex_face_verts.end())
2000 flist.
Size(), 3) - 1;
2001 hex_face_verts.insert({
t, face_center_index});
2005 face_center_index = it->second;
2008 fti[2] = face_center_index;
2009 fti[3] = elem_center_index;
2010 for (
int j = 0; j < 4; j++)
2012 for (
int k = 0; k < 2; k++)
2014 fti[k] = flist[tet_face[j][k]];
2079 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
2082 for (
int i = 0; i < 2; i++)
2084 for (
int j = 0; j < 3; j++)
2086 ti[j] = vi[quad_to_tri[i][j]];
2122 for (
int i = 0, j = 0; i <
faces_info.Size(); i++)
2138 "incorrect number of vertices: preallocated: " <<
vertices.Size()
2141 "incorrect number of elements: preallocated: " <<
elements.Size()
2144 "incorrect number of boundary elements: preallocated: "
2178 bool fix_orientation)
2181 if (fix_orientation)
2211 GeckoProgress(
real_t limit) : limit(limit) { sw.
Start(); }
2212 virtual bool quit()
const {
return limit > 0 && sw.
UserTime() > limit; }
2215class GeckoVerboseProgress :
public GeckoProgress
2221 GeckoVerboseProgress(
real_t limit) : GeckoProgress(limit) {}
2223 virtual void beginorder(
const Graph* graph, Float cost)
const
2224 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2225 virtual void endorder(
const Graph* graph, Float cost)
const
2226 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2228 virtual void beginiter(
const Graph* graph,
2229 uint iter, uint maxiter, uint window)
const
2231 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
2232 << window << std::flush;
2234 virtual void enditer(
const Graph* graph, Float mincost, Float cost)
const
2235 {
mfem::out <<
", cost = " << cost << endl; }
2240 int iterations,
int window,
2241 int period,
int seed,
bool verbose,
2247 GeckoProgress progress(time_limit);
2248 GeckoVerboseProgress vprogress(time_limit);
2251 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2259 for (
int elemid = 0; elemid <
GetNE(); ++elemid)
2261 const int *neighid = my_el_to_el.
GetRow(elemid);
2262 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2264 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2269 graph.
order(&functional, iterations, window, period, seed,
2270 verbose ? &vprogress : &progress);
2276 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2279 return graph.
cost();
2291 : coord(coord), dir(dir), points(points), mid(mid) {}
2293 bool operator()(
int i)
const
2295 return (points[3*i + coord] < mid) != dir;
2299static void HilbertSort2D(
int coord1,
2302 const Array<real_t> &points,
int *beg,
int *end,
2305 if (end - beg <= 1) {
return; }
2307 real_t xmid = (xmin + xmax)*0.5;
2308 real_t ymid = (ymin + ymax)*0.5;
2310 int coord2 = (coord1 + 1) % 2;
2313 int *p0 = beg, *p4 = end;
2314 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2315 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2316 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2320 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2321 ymin, xmin, ymid, xmid);
2323 if (p1 != p0 || p2 != p4)
2325 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2326 xmin, ymid, xmid, ymax);
2328 if (p2 != p0 || p3 != p4)
2330 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2331 xmid, ymid, xmax, ymax);
2335 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2336 ymid, xmax, ymin, xmid);
2340static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2341 const Array<real_t> &points,
int *beg,
int *end,
2345 if (end - beg <= 1) {
return; }
2347 real_t xmid = (xmin + xmax)*0.5;
2348 real_t ymid = (ymin + ymax)*0.5;
2349 real_t zmid = (zmin + zmax)*0.5;
2351 int coord2 = (coord1 + 1) % 3;
2352 int coord3 = (coord1 + 2) % 3;
2355 int *p0 = beg, *p8 = end;
2356 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2357 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2358 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2359 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2360 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2361 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2362 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2366 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2367 zmin, xmin, ymin, zmid, xmid, ymid);
2369 if (p1 != p0 || p2 != p8)
2371 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2372 ymin, zmid, xmin, ymid, zmax, xmid);
2374 if (p2 != p0 || p3 != p8)
2376 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2377 ymid, zmid, xmin, ymax, zmax, xmid);
2379 if (p3 != p0 || p4 != p8)
2381 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2382 xmin, ymax, zmid, xmid, ymid, zmin);
2384 if (p4 != p0 || p5 != p8)
2386 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2387 xmid, ymax, zmid, xmax, ymid, zmin);
2389 if (p5 != p0 || p6 != p8)
2391 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2392 ymax, zmid, xmax, ymid, zmax, xmid);
2394 if (p6 != p0 || p7 != p8)
2396 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2397 ymid, zmid, xmax, ymin, zmax, xmid);
2401 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2402 zmid, xmax, ymin, zmin, xmid, ymid);
2416 if (
spaceDim < 3) { points = 0.0; }
2419 for (
int i = 0; i <
GetNE(); i++)
2424 points[3*i + j] = center(j);
2431 indices.
Sort([&](
int a,
int b)
2432 {
return points[3*
a] < points[3*
b]; });
2437 HilbertSort2D(0,
false,
false,
2438 points, indices.
begin(), indices.
end(),
2439 min(0), min(1), max(0), max(1));
2444 HilbertSort3D(0,
false,
false,
false,
2445 points, indices.
begin(), indices.
end(),
2446 min(0), min(1), min(2), max(0), max(1), max(2));
2451 for (
int i = 0; i <
GetNE(); i++)
2453 ordering[indices[i]] = i;
2462 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2467 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
2471 MFEM_VERIFY(ordering.
Size() ==
GetNE(),
"invalid reordering array.")
2504 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2508 old_elem_node_vals[old_elid] =
new Vector(vals);
2514 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2516 int new_elid = ordering[old_elid];
2517 new_elements[new_elid] =
elements[old_elid];
2522 if (reorder_vertices)
2527 vertex_ordering = -1;
2529 int new_vertex_ind = 0;
2530 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2532 int *elem_vert =
elements[new_elid]->GetVertices();
2533 int nv =
elements[new_elid]->GetNVertices();
2534 for (
int vi = 0; vi < nv; ++vi)
2536 int old_vertex_ind = elem_vert[vi];
2537 if (vertex_ordering[old_vertex_ind] == -1)
2539 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2540 new_vertices[new_vertex_ind] =
vertices[old_vertex_ind];
2550 for (
int new_elid = 0; new_elid <
GetNE(); ++new_elid)
2552 int *elem_vert =
elements[new_elid]->GetVertices();
2553 int nv =
elements[new_elid]->GetNVertices();
2554 for (
int vi = 0; vi < nv; ++vi)
2556 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2561 for (
int belid = 0; belid <
GetNBE(); ++belid)
2563 int *be_vert =
boundary[belid]->GetVertices();
2564 int nv =
boundary[belid]->GetNVertices();
2565 for (
int vi = 0; vi < nv; ++vi)
2567 be_vert[vi] = vertex_ordering[be_vert[vi]];
2596 nodes_fes->
Update(
false);
2599 for (
int old_elid = 0; old_elid <
GetNE(); ++old_elid)
2601 int new_elid = ordering[old_elid];
2604 delete old_elem_node_vals[old_elid];
2653 length_idx[j].one =
GetLength(i, it.Column());
2654 length_idx[j].two = j;
2663 order[length_idx[i].two] = i;
2678 elements[i]->MarkEdge(v_to_v, order);
2685 boundary[i]->MarkEdge(v_to_v, order);
2692 if (*old_v_to_v && *old_elem_vert)
2699 if (*old_v_to_v == NULL)
2701 bool need_v_to_v =
false;
2709 if (dofs.
Size() > 0)
2721 if (*old_elem_vert == NULL)
2723 bool need_elem_vert =
false;
2725 for (
int i = 0; i <
GetNE(); i++)
2731 if (dofs.
Size() > 1)
2733 need_elem_vert =
true;
2739 *old_elem_vert =
new Table;
2741 for (
int i = 0; i <
GetNE(); i++)
2743 (*old_elem_vert)->AddColumnsInRow(i,
elements[i]->GetNVertices());
2745 (*old_elem_vert)->MakeJ();
2746 for (
int i = 0; i <
GetNE(); i++)
2751 (*old_elem_vert)->ShiftUpI();
2764 const int num_edge_dofs = old_dofs.
Size();
2775 if (num_edge_dofs > 0)
2784 const int old_i = (*old_v_to_v)(i, it.Column());
2785 const int new_i = it.Index();
2786 if (new_i == old_i) {
continue; }
2788 old_dofs.
SetSize(num_edge_dofs);
2789 new_dofs.
SetSize(num_edge_dofs);
2790 for (
int j = 0; j < num_edge_dofs; j++)
2792 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2793 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2797 for (
int j = 0; j < old_dofs.
Size(); j++)
2799 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2811 Table old_face_vertex;
2817 old_face_vertex.
MakeJ();
2820 faces[i]->GetNVertices());
2832 const int *old_v = old_face_vertex.
GetRow(i);
2834 switch (old_face_vertex.
RowSize(i))
2837 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2841 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2845 new_fdofs[new_i+1] = old_dofs.
Size();
2852 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2855 switch (old_face_vertex.
RowSize(i))
2858 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2859 new_v =
faces[new_i]->GetVertices();
2865 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2866 new_v =
faces[new_i]->GetVertices();
2874 for (
int j = 0; j < old_dofs.
Size(); j++)
2877 const int old_j = dof_ord[j];
2878 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2882 for (
int j = 0; j < old_dofs.
Size(); j++)
2884 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2906 for (
int i = 0; i <
GetNE(); i++)
2908 const int *old_v = old_elem_vert->
GetRow(i);
2909 const int *new_v =
elements[i]->GetVertices();
2916 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2930 <<
" FE collection) are not supported yet!");
2934 MFEM_VERIFY(dof_ord != NULL,
2935 "FE collection '" << fec->
Name()
2936 <<
"' does not define reordering for "
2940 for (
int j = 0; j < new_dofs.
Size(); j++)
2943 const int old_j = dof_ord[j];
2944 new_dofs[old_j] = offset + j;
2946 offset += new_dofs.
Size();
2949 for (
int j = 0; j < old_dofs.
Size(); j++)
2951 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2989 MFEM_ASSERT(
NURBSext,
"SetPatchAttribute is only for NURBS meshes");
2992 for (
auto e : elems)
3000 MFEM_ASSERT(
NURBSext,
"GetPatchAttribute is only for NURBS meshes");
3006 MFEM_ASSERT(
NURBSext,
"SetPatchBdrAttribute is only for NURBS meshes");
3010 for (
auto be : bdryelems)
3018 MFEM_ASSERT(
NURBSext,
"GetBdrPatchBdrAttribute is only for NURBS meshes");
3046 if (generate_edges == 1)
3064 bool fix_orientation)
3081 if (generate_edges == 1)
3147 bool generate_edges =
true;
3156 MFEM_VERIFY(
ncmesh == NULL,
"");
3191 if (
Dim > 1 && generate_edges)
3255 const bool check_orientation =
true;
3256 const bool curved = (
Nodes != NULL);
3257 const bool may_change_topology =
3259 ( check_orientation && fix_orientation &&
3263 Table *old_elem_vert = NULL;
3265 if (curved && may_change_topology)
3270 if (check_orientation)
3280 if (may_change_topology)
3285 delete old_elem_vert;
3306 for (
int i = 0; i < num_faces; i++)
3309 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
3310 " Interior face with incompatible orientations.");
3321 int NVert, NElem, NBdrElem;
3323 NVert = (nx+1) * (ny+1) * (nz+1);
3324 NElem = nx * ny * nz;
3325 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3334 NBdrElem += 2*nx*ny;
3339 NVert += nx * ny * nz;
3342 InitMesh(3, 3, NVert, NElem, NBdrElem);
3348 for (z = 0; z <= nz; z++)
3350 coord[2] = ((
real_t) z / nz) * sz;
3351 for (y = 0; y <= ny; y++)
3353 coord[1] = ((
real_t) y / ny) * sy;
3354 for (x = 0; x <= nx; x++)
3356 coord[0] = ((
real_t) x / nx) * sx;
3363 for (z = 0; z < nz; z++)
3365 coord[2] = (((
real_t) z + 0.5) / nz) * sz;
3366 for (y = 0; y < ny; y++)
3368 coord[1] = (((
real_t) y + 0.5) / ny) * sy;
3369 for (x = 0; x < nx; x++)
3371 coord[0] = (((
real_t) x + 0.5) / nx) * sx;
3378#define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3379#define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3386 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3388 for (
int k = 0; k < nx*ny*nz; k++)
3395 ind[0] = VTX(x , y , z );
3396 ind[1] = VTX(x+1, y , z );
3397 ind[2] = VTX(x+1, y+1, z );
3398 ind[3] = VTX(x , y+1, z );
3399 ind[4] = VTX(x , y , z+1);
3400 ind[5] = VTX(x+1, y , z+1);
3401 ind[6] = VTX(x+1, y+1, z+1);
3402 ind[7] = VTX(x , y+1, z+1);
3410 for (z = 0; z < nz; z++)
3412 for (y = 0; y < ny; y++)
3414 for (x = 0; x < nx; x++)
3417 ind[0] = VTX(x , y , z );
3418 ind[1] = VTX(x+1, y , z );
3419 ind[2] = VTX(x+1, y+1, z );
3420 ind[3] = VTX(x , y+1, z );
3421 ind[4] = VTX(x , y , z+1);
3422 ind[5] = VTX(x+1, y , z+1);
3423 ind[6] = VTX(x+1, y+1, z+1);
3424 ind[7] = VTX( x, y+1, z+1);
3436 ind[8] = VTXP(x, y, z);
3450 for (y = 0; y < ny; y++)
3452 for (x = 0; x < nx; x++)
3455 ind[0] = VTX(x , y , 0);
3456 ind[1] = VTX(x , y+1, 0);
3457 ind[2] = VTX(x+1, y+1, 0);
3458 ind[3] = VTX(x+1, y , 0);
3475 for (y = 0; y < ny; y++)
3477 for (x = 0; x < nx; x++)
3480 ind[0] = VTX(x , y , nz);
3481 ind[1] = VTX(x+1, y , nz);
3482 ind[2] = VTX(x+1, y+1, nz);
3483 ind[3] = VTX(x , y+1, nz);
3500 for (z = 0; z < nz; z++)
3502 for (y = 0; y < ny; y++)
3505 ind[0] = VTX(0 , y , z );
3506 ind[1] = VTX(0 , y , z+1);
3507 ind[2] = VTX(0 , y+1, z+1);
3508 ind[3] = VTX(0 , y+1, z );
3521 for (z = 0; z < nz; z++)
3523 for (y = 0; y < ny; y++)
3526 ind[0] = VTX(nx, y , z );
3527 ind[1] = VTX(nx, y+1, z );
3528 ind[2] = VTX(nx, y+1, z+1);
3529 ind[3] = VTX(nx, y , z+1);
3542 for (x = 0; x < nx; x++)
3544 for (z = 0; z < nz; z++)
3547 ind[0] = VTX(x , 0, z );
3548 ind[1] = VTX(x+1, 0, z );
3549 ind[2] = VTX(x+1, 0, z+1);
3550 ind[3] = VTX(x , 0, z+1);
3563 for (x = 0; x < nx; x++)
3565 for (z = 0; z < nz; z++)
3568 ind[0] = VTX(x , ny, z );
3569 ind[1] = VTX(x , ny, z+1);
3570 ind[2] = VTX(x+1, ny, z+1);
3571 ind[3] = VTX(x+1, ny, z );
3587 ofstream test_stream(
"debug.mesh");
3589 test_stream.close();
3617 for (
real_t j = 0; j < ny+1; j++)
3619 real_t cy = (j / ny) * sy;
3620 for (
real_t i = 0; i < nx+1; i++)
3622 real_t cx = (i / nx) * sx;
3629 for (
int y = 0; y < ny; y++)
3631 for (
int x = 0; x < nx; x++)
3633 ind[0] = x + y*(nx+1);
3634 ind[1] = x + 1 +y*(nx+1);
3635 ind[2] = x + 1 + (y+1)*(nx+1);
3636 ind[3] = x + (y+1)*(nx+1);
3642 for (
int i = 0; i < nx; i++)
3648 for (
int j = 0; j < ny; j++)
3691 for (
real_t j = 0; j < ny+1; j++)
3693 real_t cy = (j / ny) * sy;
3694 for (
real_t i = 0; i < nx+1; i++)
3696 real_t cx = (i / nx) * sx;
3703 for (
int y = 0; y < ny; y++)
3705 for (
int x = 0; x < nx; x++)
3707 ind[0] = x + y*(nx+1);
3708 ind[1] = x + 1 +y*(nx+1);
3709 ind[2] = x + 1 + (y+1)*(nx+1);
3710 ind[3] = x + (y+1)*(nx+1);
3716 for (
int i = 0; i < nx; i++)
3722 for (
int j = 0; j < ny; j++)
3748 const int NVert = (nx+1) * (ny+1) * (nz+1);
3749 const int NElem = nx * ny * nz * 24;
3750 const int NBdrElem = 2*(nx*ny+nx*nz+ny*nz)*4;
3752 InitMesh(3, 3, NVert, NElem, NBdrElem);
3757 for (
real_t z = 0; z <= nz; z++)
3759 coord[2] = ( z / nz) * sz;
3760 for (
real_t y = 0; y <= ny; y++)
3762 coord[1] = (y / ny) * sy;
3763 for (
real_t x = 0; x <= nx; x++)
3765 coord[0] = (x / nx) * sx;
3771 std::map<std::array<int, 4>,
int> hex_face_verts;
3772 auto VertexIndex = [nx, ny](
int xc,
int yc,
int zc)
3774 return xc + (yc + zc*(ny+1))*(nx+1);
3778 for (
int z = 0; z < nz; z++)
3780 for (
int y = 0; y < ny; y++)
3782 for (
int x = 0; x < nx; x++)
3785 ind[0] = VertexIndex(x , y , z );
3786 ind[1] = VertexIndex(x+1, y , z );
3787 ind[2] = VertexIndex(x+1, y+1, z );
3788 ind[3] = VertexIndex(x , y+1, z );
3789 ind[4] = VertexIndex(x , y , z+1);
3790 ind[5] = VertexIndex(x+1, y , z+1);
3791 ind[6] = VertexIndex(x+1, y+1, z+1);
3792 ind[7] = VertexIndex( x, y+1, z+1);
3798 hex_face_verts.clear();
3807 std::map<std::array<int, 3>,
int> tet_face_count;
3809 std::map<std::array<int, 3>,
int> face_count_map;
3814 return std::array<int, 3>{v[0], v[1], v[2]};
3823 for (
int j = 0; j < el_faces.
Size(); j++)
3826 auto t = get3array(vertidxs);
3827 auto it = tet_face_count.find(
t);
3828 if (it == tet_face_count.end())
3830 tet_face_count.insert({
t, 1});
3831 face_count_map.insert({
t, el_faces[j]});
3840 for (
const auto &edge : tet_face_count)
3842 if (edge.second == 1)
3844 int facenum = (face_count_map.find(edge.first))->second;
3851 ofstream test_stream(
"debug.mesh");
3853 test_stream.close();
3862 bool generate_edges,
bool sfc_ordering)
3886 for (j = 0; j < ny+1; j++)
3888 cy = ((
real_t) j / ny) * sy;
3889 for (i = 0; i < nx+1; i++)
3891 cx = ((
real_t) i / nx) * sx;
3903 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
3905 for (k = 0; k < nx*ny; k++)
3909 ind[0] = i + j*(nx+1);
3910 ind[1] = i + 1 +j*(nx+1);
3911 ind[2] = i + 1 + (j+1)*(nx+1);
3912 ind[3] = i + (j+1)*(nx+1);
3919 for (j = 0; j < ny; j++)
3921 for (i = 0; i < nx; i++)
3923 ind[0] = i + j*(nx+1);
3924 ind[1] = i + 1 +j*(nx+1);
3925 ind[2] = i + 1 + (j+1)*(nx+1);
3926 ind[3] = i + (j+1)*(nx+1);
3935 for (i = 0; i < nx; i++)
3941 for (j = 0; j < ny; j++)
3963 for (j = 0; j < ny+1; j++)
3965 cy = ((
real_t) j / ny) * sy;
3966 for (i = 0; i < nx+1; i++)
3968 cx = ((
real_t) i / nx) * sx;
3977 for (j = 0; j < ny; j++)
3979 for (i = 0; i < nx; i++)
3981 ind[0] = i + j*(nx+1);
3982 ind[1] = i + 1 + (j+1)*(nx+1);
3983 ind[2] = i + (j+1)*(nx+1);
3986 ind[1] = i + 1 + j*(nx+1);
3987 ind[2] = i + 1 + (j+1)*(nx+1);
3995 for (i = 0; i < nx; i++)
4001 for (j = 0; j < ny; j++)
4011 MFEM_ABORT(
"Unsupported element type.");
4017 if (generate_edges == 1)
4055 for (j = 0; j < n+1; j++)
4061 for (j = 0; j < n; j++)
4088 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4140 for (
int i = 0; i <
faces.Size(); i++)
4182 if (
dynamic_cast<const ParMesh*
>(&mesh))
4194 if (mesh.
Nodes && copy_nodes)
4226 int refine,
bool fix_orientation)
4230 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
4231 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
4248 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
4258 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
4294 ref_factors = ref_factor;
4307Mesh::Mesh(
const std::string &filename,
int generate_edges,
int refine,
4308 bool fix_orientation)
4309 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4318 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
4322 Load(imesh, generate_edges, refine, fix_orientation);
4327 bool fix_orientation)
4328 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4331 Load(input, generate_edges, refine, fix_orientation);
4340 "Not enough vertices in external array : "
4341 "len_vertex_data = "<< len_vertex_data <<
", "
4346 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
4351 memcpy(vertex_data,
vertices.GetData(),
4360 int *element_attributes,
int num_elements,
4362 int *boundary_attributes,
int num_boundary_elements,
4364 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4366 if (space_dimension == -1)
4372 num_boundary_elements);
4375 int boundary_index_stride = num_boundary_elements > 0 ?
4379 vertices.MakeRef(
reinterpret_cast<Vertex*
>(vertices_), num_vertices);
4382 for (
int i = 0; i < num_elements; i++)
4385 elements[i]->SetVertices(element_indices + i * element_index_stride);
4386 elements[i]->SetAttribute(element_attributes[i]);
4390 for (
int i = 0; i < num_boundary_elements; i++)
4393 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
4394 boundary[i]->SetAttribute(boundary_attributes[i]);
4410#ifdef MFEM_USE_MEMALLOC
4419 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
4432 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
4435 for (
int i = 0; i < nv; i++)
4448 for (
int j = 0; j < nv; j++)
4520 MFEM_ABORT(
"invalid element type: " << type);
4527 std::string parse_tag)
4529 int curved = 0, read_gf = 1;
4530 bool finalize_topo =
true;
4534 MFEM_ABORT(
"Input stream is not open");
4541 getline(input, mesh_type);
4545 int mfem_version = 0;
4546 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
4547 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
4548 else if (mesh_type ==
"MFEM mesh v1.3") { mfem_version = 13; }
4552 int mfem_nc_version = 0;
4553 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
4554 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
4562 if (mfem_version >= 12 && parse_tag.empty())
4564 parse_tag =
"mfem_mesh_end";
4568 else if (mfem_nc_version)
4570 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4577 MFEM_VERIFY(mfem_nc_version >= 10,
4578 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
4579 "used to load a parallel nonconforming mesh, sorry.");
4582 input, mfem_nc_version, curved, is_nc);
4587 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
4600 else if (mesh_type ==
"linemesh")
4604 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
4606 if (mesh_type ==
"curved_areamesh2")
4612 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
4616 else if (mesh_type ==
"TrueGrid")
4620 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
4622 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
4624 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4625 "Unsupported VTK format");
4626 ReadVTKMesh(input, curved, read_gf, finalize_topo);
4628 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
4632 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
4636 else if (mesh_type ==
"MFEM NURBS mesh v1.1")
4640 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
4645 else if (mesh_type ==
"$MeshFormat")
4650 ((mesh_type.size() > 2 &&
4651 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
4652 (mesh_type.size() > 3 &&
4653 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
4658#ifdef MFEM_USE_NETCDF
4661 MFEM_ABORT(
"NetCDF support requires configuration with"
4662 " MFEM_USE_NETCDF=YES");
4668 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
4669 " Use mfem::named_ifgzstream for input.");
4675 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
4701 bool generate_bdr =
false;
4706 if (curved && read_gf)
4720 if (mfem_version >= 12)
4726 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
4727 getline(input, line);
4733 if (line ==
"mfem_mesh_end") {
break; }
4735 while (line != parse_tag);
4737 else if (mfem_nc_version >= 10)
4742 MFEM_VERIFY(ident ==
"mfem_mesh_end",
4743 "invalid mesh: end of file tag not found");
4750 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4752 int i, j, ie, ib, iv, *v, nv;
4783 for (i = 0; i < num_pieces; i++)
4790 for (i = 0; i < num_pieces; i++)
4796 for (j = 0; j < m->
GetNE(); j++)
4801 for (j = 0; j < m->
GetNBE(); j++)
4806 for (
int k = 0; k < nv; k++)
4808 v[k] = lvert_vert[v[k]];
4813 for (j = 0; j < m->
GetNV(); j++)
4825 for (i = 0; i < num_pieces; i++)
4836 for (i = 0; i < num_pieces; i++)
4840 for (j = 0; j < m->
GetNE(); j++)
4845 for (
int k = 0; k < nv; k++)
4852 for (j = 0; j < m->
GetNBE(); j++)
4857 for (
int k = 0; k < nv; k++)
4864 for (j = 0; j < m->
GetNV(); j++)
4878 for (i = 0; i < num_pieces; i++)
4880 gf_array[i] = mesh_array[i]->
GetNodes();
4893 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
4896 ref_factors = ref_factor;
4907 int orig_ne = orig_mesh.
GetNE();
4908 MFEM_VERIFY(ref_factors.
Size() == orig_ne,
4909 "Number of refinement factors must equal number of elements")
4910 MFEM_VERIFY(orig_ne == 0 ||
4911 ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
4914 "Invalid refinement type. Must use closed basis type.");
4916 int min_ref = orig_ne > 0 ? ref_factors.
Min() : 1;
4917 int max_ref = orig_ne > 0 ? ref_factors.
Max() : 1;
4919 bool var_order = (min_ref != max_ref);
4932 for (
int i = 0; i < orig_ne; i++)
4948 for (
int el = 0; el < orig_ne; el++)
4960 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
4961 for (
int i = 0; i < phys_pts.
Width(); i++)
4970 for (
int k = 0; k < nvert; k++)
4973 v[k] = rdofs[c2h_map[cid]];
4986 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
4997 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
5003 for (
int k = 0; k < nvert; k++)
5006 v[k] = rdofs[c2h_map[cid]];
5028 for (
int iel = 0; iel < orig_ne; iel++)
5037 const int *node_map = NULL;
5040 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
5041 const int *vertex_map = vertex_fec.
GetDofMap(geom);
5042 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
5043 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
5046 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
5049 int iv = vertex_map[iv_lex];
5051 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
5053 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
5056 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
5067 using GeomRef = std::pair<Geometry::Type, int>;
5068 std::map<GeomRef, int> point_matrices_offsets;
5070 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5074 GeomRef id(geom, ref_factors[el_coarse]);
5075 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
5081 point_matrices_offsets[id] = n_point_matrices[geom];
5082 n_point_matrices[geom] += nref_el;
5089 int nmatrices = n_point_matrices[geom];
5096 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
5099 int ref = ref_factors[el_coarse];
5100 int offset = point_matrices_offsets[GeomRef(geom, ref)];
5106 for (
int k = 0; k < nvert; k++)
5138 "Mesh::MakeSimplicial requires a properly oriented input mesh");
5140 "Mesh::MakeSimplicial does not support non-conforming meshes.")
5147 Mesh copy(orig_mesh);
5152 int nv = orig_mesh.
GetNV();
5153 int ne = orig_mesh.
GetNE();
5154 int nbe = orig_mesh.
GetNBE();
5167 int new_ne = 0, new_nbe = 0;
5168 for (
int i=0; i<ne; ++i)
5172 for (
int i=0; i<nbe; ++i)
5181 for (
int i=0; i<nv; ++i)
5191 if (vglobal == NULL)
5194 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
5195 vglobal = vglobal_id.
GetData();
5198 constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
5199 constexpr int quad_ntris = 2, prism_ntets = 3;
5200 static const int quad_trimap[2][nv_tri*quad_ntris] =
5212 static const int prism_rot[nv_prism*nv_prism] =
5221 static const int prism_f[nv_quad] = {1, 2, 5, 4};
5222 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
5236 static const int hex_rot[nv_hex*nv_hex] =
5238 0, 1, 2, 3, 4, 5, 6, 7,
5239 1, 0, 4, 5, 2, 3, 7, 6,
5240 2, 1, 5, 6, 3, 0, 4, 7,
5241 3, 0, 1, 2, 7, 4, 5, 6,
5242 4, 0, 3, 7, 5, 1, 2, 6,
5243 5, 1, 0, 4, 6, 2, 3, 7,
5244 6, 2, 1, 5, 7, 3, 0, 4,
5245 7, 3, 2, 6, 4, 0, 1, 5
5247 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
5248 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
5249 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
5250 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
5251 static const int hex_tetmap0[nv_tet*5] =
5258 static const int hex_tetmap1[nv_tet*6] =
5265 static const int hex_tetmap2[nv_tet*6] =
5272 static const int hex_tetmap3[nv_tet*6] =
5279 static const int *hex_tetmaps[4] =
5281 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
5284 auto find_min = [](
const int*
a,
int n) {
return std::min_element(
a,
a+n)-
a; };
5286 for (
int i=0; i<ne; ++i)
5288 const int *v = orig_mesh.
elements[i]->GetVertices();
5292 if (num_subdivisions[orig_geom] == 1)
5304 for (
int itri=0; itri<quad_ntris; ++itri)
5309 for (
int iv=0; iv<nv_tri; ++iv)
5311 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
5319 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
5322 int irot = find_min(vg, nv_prism);
5323 for (
int iv=0; iv<nv_prism; ++iv)
5325 int jv = prism_rot[iv + irot*nv_prism];
5330 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
5331 int j = find_min(q, nv_quad);
5332 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
5333 for (
int itet=0; itet<prism_ntets; ++itet)
5338 for (
int iv=0; iv<nv_tet; ++iv)
5340 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
5348 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
5352 int irot = find_min(vg, nv_hex);
5353 for (
int iv=0; iv<nv_hex; ++iv)
5355 int jv = hex_rot[iv + irot*nv_hex];
5365 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
5366 j = find_min(q, nv_quad);
5367 if (j == 0 || j == 2) { bitmask += 4; }
5369 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
5370 j = find_min(q, nv_quad);
5371 if (j == 1 || j == 3) { bitmask += 2; }
5373 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
5374 j = find_min(q, nv_quad);
5375 if (j == 0 || j == 2) { bitmask += 1; }
5378 int nrot = num_rot[bitmask];
5379 for (
int k=0; k<nrot; ++k)
5393 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
5394 int ntets = (ndiags == 0) ? 5 : 6;
5395 const int *tetmap = hex_tetmaps[ndiags];
5396 for (
int itet=0; itet<ntets; ++itet)
5401 for (
int iv=0; iv<nv_tet; ++iv)
5403 v2[iv] = vg[tetmap[itet + iv*ntets]];
5412 for (
int i=0; i<nbe; ++i)
5414 const int *v = orig_mesh.
boundary[i]->GetVertices();
5417 if (num_subdivisions[orig_geom] == 1)
5427 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
5429 int iv_min = find_min(vg, nv_quad);
5430 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
5431 for (
int itri=0; itri<quad_ntris; ++itri)
5436 for (
int iv=0; iv<nv_tri; ++iv)
5438 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
5445 MFEM_ABORT(
"Unreachable");
5459 Mesh periodic_mesh(orig_mesh,
true);
5465 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
5470 for (
int j = 0; j < nv; j++)
5476 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
5481 for (
int j = 0; j < nv; j++)
5488 return periodic_mesh;
5492 const std::vector<Vector> &translations,
real_t tol)
const
5496 Vector coord(sdim), at(sdim), dx(sdim);
5497 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
5498 xMax = xMin = xDiff = 0.0;
5501 unordered_set<int> bdr_v;
5502 for (
int be = 0; be <
GetNBE(); be++)
5507 for (
int i = 0; i < dofs.
Size(); i++)
5509 bdr_v.insert(dofs[i]);
5512 for (
int j = 0; j < sdim; j++)
5514 xMax[j] = max(xMax[j], coord[j]);
5515 xMin[j] = min(xMin[j], coord[j]);
5519 add(xMax, -1.0, xMin, xDiff);
5528 unordered_map<int, int> replica2primary;
5530 unordered_map<int, unordered_set<int>> primary2replicas;
5533 std::unique_ptr<KDTreeBase<int,real_t>> kdtree;
5534 if (sdim == 1) { kdtree.reset(
new KDTree1D); }
5535 else if (sdim == 2) { kdtree.reset(
new KDTree2D); }
5536 else if (sdim == 3) { kdtree.reset(
new KDTree3D); }
5537 else { MFEM_ABORT(
"Invalid space dimension."); }
5541 for (
const int v : bdr_v)
5543 primary2replicas[v];
5551 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
5553 if (r ==
p) {
return; }
5554 primary2replicas[
p].insert(r);
5555 replica2primary[r] =
p;
5556 for (
const int s : primary2replicas[r])
5558 primary2replicas[
p].insert(
s);
5559 replica2primary[
s] =
p;
5561 primary2replicas.erase(r);
5564 for (
unsigned int i = 0; i < translations.size(); i++)
5566 for (
int vi : bdr_v)
5569 add(coord, translations[i], at);
5571 const int vj = kdtree->FindClosestPoint(at.GetData());
5573 add(at, -1.0, coord, dx);
5575 if (dx.
Norml2() > dia*tol) {
continue; }
5580 const bool pi = primary2replicas.find(vi) != primary2replicas.end();
5581 const bool pj = primary2replicas.find(vj) != primary2replicas.end();
5587 make_replica(vj, vi);
5592 const int owner_of_vj = replica2primary[vj];
5594 make_replica(vi, owner_of_vj);
5600 const int owner_of_vi = replica2primary[vi];
5601 make_replica(vj, owner_of_vi);
5608 const int owner_of_vi = replica2primary[vi];
5609 const int owner_of_vj = replica2primary[vj];
5610 make_replica(owner_of_vj, owner_of_vi);
5615 std::vector<int> v2v(
GetNV());
5616 for (
size_t i = 0; i < v2v.size(); i++)
5620 for (
const auto &r2p : replica2primary)
5622 v2v[r2p.first] = r2p.second;
5629 MFEM_VERIFY(
NURBSext,
"Mesh::RefineNURBSFromFile: Not a NURBS mesh!");
5630 mfem::out<<
"Refining NURBS from refinement file: "<<ref_file<<endl;
5633 ifstream input(ref_file);
5640 mfem::out<<
"Knot vectors in ref_file: "<<nkv<<endl;
5642 MFEM_ABORT(
"Refine file does not have the correct number of knot vectors");
5647 for (
int kv = 0; kv < nkv; kv++)
5649 knotVec[kv] =
new Vector();
5650 knotVec[kv]->
Load(input);
5658 for (
int kv = 0; kv < nkv; kv++)
5668 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5673 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5690 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5695 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5712 mfem_error(
"Mesh::KnotRemove : Not a NURBS mesh!");
5717 mfem_error(
"Mesh::KnotRemove : KnotVector array size mismatch!");
5740 "Refinement factors must be defined for each dimension");
5742 MFEM_VERIFY(
NURBSext,
"NURBSUniformRefinement is only for NURBS meshes");
5752 cf1 = (cf1 &&
f == 1);
5769 for (
int i=0; i<cf.
Size(); ++i) { cf[i] *= rf[i]; }
5783 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
5807 for (
int i = 0; i <
elements.Size(); i++)
5817 for (
int i = 0; i <
boundary.Size(); i++)
5835 for (
int i = 0; i < vd; i++)
5902 input >> edge_to_knot[j] >> v[0] >> v[1];
5905 edge_to_knot[j] = -1 - edge_to_knot[j];
5925 if (edge_to_knot.
Size() == 0)
5929 constexpr int notset = -9999999;
5930 edge_to_knot = notset;
5942 edge0[0] = 0; edge1[0] = 2;
5943 edge0[1] = 1; edge1[1] = 3;
5951 edge0[0] = 0; edge1[0] = 2;
5952 edge0[1] = 0; edge1[1] = 4;
5953 edge0[2] = 0; edge1[2] = 6;
5955 edge0[3] = 1; edge1[3] = 3;
5956 edge0[4] = 1; edge1[4] = 5;
5957 edge0[5] = 1; edge1[5] = 7;
5959 edge0[6] = 8; edge1[6] = 9;
5960 edge0[7] = 8; edge1[7] = 10;
5961 edge0[8] = 8; edge1[8] = 11;
5969 int e0, e1, v0, v1, df;
5975 const int *v =
elements[
p]->GetVertices();
5976 for (j = 0; j < edges.
Size(); j++)
5979 const int *e =
elements[
p]->GetEdgeVertices(j);
5992 for (j = 0; j < edge1.
Size(); j++)
5994 e0 = edges[edge0[j]];
5995 e1 = edges[edge1[j]];
5996 v0 = edge_to_knot[e0];
5997 v1 = edge_to_knot[e1];
5998 df = flip*oedge[edge0[j]]*oedge[edge1[j]];
6001 if ((v0 == notset) && (v1 == notset))
6003 edge_to_knot[e0] = knot;
6004 edge_to_knot[e1] = knot;
6010 else if ((v0 != notset) && (v1 == notset))
6012 edge_to_knot[e1] = (df >= 0 ? -v0-1 : v0);
6014 else if ((v0 == notset) && (v1 != notset))
6016 edge_to_knot[e0] = (df >= 0 ? -v1-1 : v1);
6036 for (j = 0; j < edge1.
Size(); j++)
6038 e0 = edges[edge0[j]];
6039 e1 = edges[edge1[j]];
6040 v0 = edge_to_knot[e0];
6041 v1 = edge_to_knot[e1];
6042 v0 = ( v0 >= 0 ? v0 : -v0-1);
6043 v1 = ( v1 >= 0 ? v1 : -v1-1);
6049 edge_to_knot[e1] = (oedge[edge1[j]] >= 0 ? v0 : -v0-1);
6053 edge_to_knot[e0] = (oedge[edge0[j]] >= 0 ? v1 : -v1-1);
6061 while (corrections > 0 && passes <
GetNE() + 1);
6064 if (corrections > 0)
6066 mfem::err<<
"Edge_to_knot mapping potentially incorrect"<<endl;
6068 mfem::err<<
" corrections = "<<corrections<<endl;
6078 k = edge_to_knot[j];
6079 cnt[(k >= 0 ? k : -k-1)]++;
6083 for (j = 0; j < cnt.
Size(); j++)
6085 cnt[j] = (cnt[j] > 0 ? k++ : -1);
6090 k = edge_to_knot[j];
6091 edge_to_knot[j] = (k >= 0 ? cnt[k]:-cnt[-k-1]-1);
6095 mfem::out<<
"Generated edge to knot mapping:"<<endl;
6099 k = edge_to_knot[j];
6108 mfem::out<<(k >= 0 ? k:-k-1)<<
" "<< v[0] <<
" "<<v[1]<<endl;
6112 if (corrections > 0 ) {
mfem_error(
"Mesh::LoadPatchTopo");}
6118 if (
p.Size() >= v.
Size())
6120 for (
int d = 0; d < v.
Size(); d++)
6128 for (d = 0; d <
p.Size(); d++)
6132 for ( ; d < v.
Size(); d++)
6175 const bool warn =
true;
6178 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
6182 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
6183 " H1-continuous mesh!\n "
6184 "If this is the desired behavior, you can silence"
6185 " this warning by converting\n "
6186 "the NURBS mesh to high-order mesh in advance by"
6187 " calling the method\n "
6188 "Mesh::SetCurvature().");
6219 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
6238 MFEM_ASSERT(nodes != NULL,
"");
6254 case 1:
return GetNV();
6290#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
6291static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
6296 int i, j, k, wo = 0, fo = 0;
6305 int *vi =
elements[i]->GetVertices();
6308 for (j = 0; j < 3; j++)
6312 for (j = 0; j < 2; j++)
6313 for (k = 0; k < 2; k++)
6315 J(j, k) = v[j+1][k] - v[0][k];
6336 MFEM_ABORT(
"Invalid 2D element type \""
6353 int *vi =
elements[i]->GetVertices();
6359 for (j = 0; j < 4; j++)
6363 for (j = 0; j < 3; j++)
6364 for (k = 0; k < 3; k++)
6366 J(j, k) = v[j+1][k] - v[0][k];
6425 MFEM_ABORT(
"Invalid 3D element type \""
6431#if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
6434 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
6435 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
6439 MFEM_CONTRACT_VAR(fo);
6452 if (test[0] == base[0])
6453 if (test[1] == base[1])
6461 else if (test[0] == base[1])
6462 if (test[1] == base[0])
6471 if (test[1] == base[0])
6482 for (
int j = 0; j < 3; j++)
6483 if (test[aor[j]] != base[j])
6485 mfem::err <<
"Mesh::GetTriOrientation(...)" << endl;
6487 for (
int k = 0; k < 3; k++)
6492 for (
int k = 0; k < 3; k++)
6513 const int oo[6][6] =
6523 int ori_a_c = oo[ori_a_b][ori_b_c];
6529 const int inv_ori[6] = {0, 1, 4, 3, 2, 5};
6530 return inv_ori[ori];
6537 for (i = 0; i < 4; i++)
6538 if (test[i] == base[0])
6545 if (test[(i+1)%4] == base[1])
6554 for (
int j = 0; j < 4; j++)
6555 if (test[aor[j]] != base[j])
6557 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
6559 for (
int k = 0; k < 4; k++)
6564 for (
int k = 0; k < 4; k++)
6573 if (test[(i+1)%4] == base[1])
6590 const int oo[8][8] =
6592 {0, 1, 2, 3, 4, 5, 6, 7},
6593 {1, 0, 3, 2, 5, 4, 7, 6},
6594 {2, 7, 4, 1, 6, 3, 0, 5},
6595 {3, 6, 5, 0, 7, 2, 1, 4},
6596 {4, 5, 6, 7, 0, 1, 2, 3},
6597 {5, 4, 7, 6, 1, 0, 3, 2},
6598 {6, 3, 0, 5, 2, 7, 4, 1},
6599 {7, 2, 1, 4, 3, 6, 5, 0}
6602 int ori_a_c = oo[ori_a_b][ori_b_c];
6608 const int inv_ori[8] = {0, 1, 6, 3, 4, 5, 2, 7};
6609 return inv_ori[ori];
6620 if (test[0] == base[0])
6621 if (test[1] == base[1])
6622 if (test[2] == base[2])
6630 else if (test[2] == base[1])
6631 if (test[3] == base[2])
6640 if (test[1] == base[2])
6648 else if (test[1] == base[0])
6649 if (test[2] == base[1])
6650 if (test[0] == base[2])
6658 else if (test[3] == base[1])
6659 if (test[2] == base[2])
6668 if (test[3] == base[2])
6676 else if (test[2] == base[0])
6677 if (test[3] == base[1])
6678 if (test[0] == base[2])
6686 else if (test[0] == base[1])
6687 if (test[1] == base[2])
6696 if (test[3] == base[2])
6705 if (test[0] == base[1])
6706 if (test[2] == base[2])
6714 else if (test[1] == base[1])
6715 if (test[0] == base[2])
6724 if (test[1] == base[2])
6735 for (
int j = 0; j < 4; j++)
6736 if (test[aor[j]] != base[j])
6761 int *bv =
boundary[i]->GetVertices();
6781 if (
faces_info[fi].Elem2No >= 0) {
continue; }
6784 int *bv =
boundary[i]->GetVertices();
6786 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
6787 const int *fv =
faces[fi]->GetVertices();
6804 MFEM_ABORT(
"Invalid 2D boundary element type \""
6805 << bdr_type <<
"\"");
6810 if (orientation % 2 == 0) {
continue; }
6812 if (!fix_it) {
continue; }
6848 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
6866 MFEM_ASSERT(o >= 0 && o < 2,
"Invalid orientation for Geometry::SEGMENT!");
6878 MFEM_ASSERT(o >= 0 && o < 6,
"Invalid orientation for Geometry::TRIANGLE!");
6891 fip.
x = 1.0 - ip.
x - ip.
y;
6896 fip.
x = 1.0 - ip.
x - ip.
y;
6902 fip.
y = 1.0 - ip.
x - ip.
y;
6907 fip.
y = 1.0 - ip.
x - ip.
y;
6912 MFEM_ASSERT(o >= 0 && o < 8,
"Invalid orientation for Geometry::SQUARE!");
6956 MFEM_ABORT(
"Unsupported face geometry for TransformBdrElementToFace!");
6963 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
6974 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
6993 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
6994 "is not generated.");
6997 const int *v =
elements[i]->GetVertices();
6998 const int ne =
elements[i]->GetNEdges();
7000 for (
int j = 0; j < ne; j++)
7002 const int *e =
elements[i]->GetEdgeVertices(j);
7003 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7014 const int *v =
boundary[i]->GetVertices();
7015 cor[0] = (v[0] < v[1]) ? (1) : (-1);
7028 const int *v =
boundary[i]->GetVertices();
7029 const int ne =
boundary[i]->GetNEdges();
7031 for (
int j = 0; j < ne; j++)
7033 const int *e =
boundary[i]->GetEdgeVertices(j);
7034 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7046 const int *v =
faces[i]->GetVertices();
7047 o[0] = (v[0] < v[1]) ? (1) : (-1);
7059 const int *v =
faces[i]->GetVertices();
7060 const int ne =
faces[i]->GetNEdges();
7062 for (
int j = 0; j < ne; j++)
7064 const int *e =
faces[i]->GetEdgeVertices(j);
7065 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
7093 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
7144 for (j = 0; j < nv; j++)
7156 for (j = 0; j < nv; j++)
7203 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
7207 int n = el_faces.
Size();
7210 for (
int j = 0; j < n; j++)
7214 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
7218 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
7219 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
7236 for (
auto f : elem_faces)
7257 const int *bv =
boundary[i]->GetVertices();
7267 default: MFEM_ABORT(
"invalid geometry");
7276 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
7279 const int *bv =
boundary[bdr_el]->GetVertices();
7287 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
7294 int bdr_el,
int &el,
int &info)
const
7299 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
7302 const int *bv =
boundary[bdr_el]->GetVertices();
7310 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
7337 for (j = 0; j < nv; j++)
7339 pointmat(k, j) =
vertices[v[j]](k);
7354 for (j = 0; j < nv; j++)
7356 pointmat(k, j) =
vertices[v[j]](k);
7368 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
7371 return sqrt(length);
7379 for (
int i = 0; i < elem_array.
Size(); i++)
7384 for (
int i = 0; i < elem_array.
Size(); i++)
7386 const int *v = elem_array[i]->GetVertices();
7387 const int ne = elem_array[i]->GetNEdges();
7388 for (
int j = 0; j < ne; j++)
7390 const int *e = elem_array[i]->GetEdgeVertices(j);
7404 v_to_v.
Push(v[0], v[1]);
7411 const int *v =
elements[i]->GetVertices();
7412 const int ne =
elements[i]->GetNEdges();
7413 for (
int j = 0; j < ne; j++)
7415 const int *e =
elements[i]->GetEdgeVertices(j);
7416 v_to_v.
Push(v[e[0]], v[e[1]]);
7424 int i, NumberOfEdges;
7440 const int *v =
boundary[i]->GetVertices();
7454 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
7458 return NumberOfEdges;
7517 if (
faces[gf] == NULL)
7549 if (
faces[gf] == NULL)
7559 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
7560 "Interior edge found between 2D elements "
7562 <<
" and " << el <<
".");
7563 int *v =
faces[gf]->GetVertices();
7565 if (v[1] == v0 && v[0] == v1)
7569 else if (v[0] == v0 && v[1] == v1)
7579 MFEM_ABORT(
"internal error");
7585 int v0,
int v1,
int v2)
7587 if (
faces[gf] == NULL)
7597 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
7598 "Interior triangular face found connecting elements "
7600 <<
" and " << el <<
".");
7601 int orientation, vv[3] = { v0, v1, v2 };
7608 faces_info[gf].Elem2Inf = 64 * lf + orientation;
7613 int v0,
int v1,
int v2,
int v3)
7625 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
7626 "Interior quadrilateral face found connecting elements "
7628 <<
" and " << el <<
".");
7629 int vv[4] = { v0, v1, v2, v3 };
7649 faces.SetSize(nfaces);
7651 for (
int i = 0; i < nfaces; ++i)
7670 const int ne =
elements[i]->GetNEdges();
7671 for (
int j = 0; j < ne; j++)
7673 const int *e =
elements[i]->GetEdgeVertices(j);
7684 for (
int j = 0; j < 4; j++)
7688 v[fv[0]], v[fv[1]], v[fv[2]]);
7694 for (
int j = 0; j < 2; j++)
7698 v[fv[0]], v[fv[1]], v[fv[2]]);
7700 for (
int j = 2; j < 5; j++)
7704 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7710 for (
int j = 0; j < 1; j++)
7714 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7716 for (
int j = 1; j < 5; j++)
7720 v[fv[0]], v[fv[1]], v[fv[2]]);
7726 for (
int j = 0; j < 6; j++)
7730 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7735 MFEM_ABORT(
"Unexpected type of Element.");
7743 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
7754 nc_faces_info.Reserve(list.masters.Size() + list.slaves.Size());
7761 if (master.index >= nfaces) {
continue; }
7767 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
7768 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
7774 if (slave.index < 0 ||
7775 slave.index >= nfaces ||
7776 slave.master >= nfaces)
7795 list.point_matrices[slave.geom][slave.matrix]));
7804 const int *v =
elements[i]->GetVertices();
7809 for (
int j = 0; j < 4; j++)
7812 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
7818 for (
int j = 0; j < 1; j++)
7821 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7823 for (
int j = 1; j < 5; j++)
7826 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
7832 for (
int j = 0; j < 2; j++)
7835 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
7837 for (
int j = 2; j < 5; j++)
7840 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7848 for (
int j = 0; j < 6; j++)
7851 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7880 for (
int j = 0; j < 4; j++)
7884 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
7890 for (
int j = 0; j < 2; j++)
7894 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
7896 for (
int j = 2; j < 5; j++)
7900 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
7906 for (
int j = 0; j < 1; j++)
7910 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
7912 for (
int j = 1; j < 5; j++)
7916 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
7924 for (
int j = 0; j < 6; j++)
7928 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
7933 MFEM_ABORT(
"Unexpected type of Element.");
7947 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
7952 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
7956 MFEM_ABORT(
"Unexpected type of boundary Element.");
7970void Rotate3(
int &
a,
int &
b,
int &c)
8002 Table *old_elem_vert = NULL;
8013 int *v =
elements[i]->GetVertices();
8015 Rotate3(v[0], v[1], v[2]);
8018 Rotate3(v[1], v[2], v[3]);
8031 int *v =
boundary[i]->GetVertices();
8033 Rotate3(v[0], v[1], v[2]);
8049 delete old_elem_vert;
8065 if (
p[i] < pmin[i]) { pmin[i] =
p[i]; }
8066 if (
p[i] > pmax[i]) { pmax[i] =
p[i]; }
8080 for (
int i =
spaceDim-1; i >= 0; i--)
8082 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
8083 if (idx < 0) { idx = 0; }
8084 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
8085 part = part * nxyz[i] + idx;
8087 partitioning[el] = part;
8090 return partitioning;
8100#ifdef MFEM_USE_METIS
8102 int print_messages = 1;
8105 int init_flag, fin_flag;
8106 MPI_Initialized(&init_flag);
8107 MPI_Finalized(&fin_flag);
8108 if (init_flag && !fin_flag)
8112 if (rank != 0) { print_messages = 0; }
8116 int i, *partitioning;
8126 partitioning[i] = 0;
8133 partitioning[i] = i;
8139#ifndef MFEM_USE_METIS_5
8151 bool freedata =
false;
8153 idx_t *mpartitioning;
8156 if (
sizeof(
idx_t) ==
sizeof(int))
8160 mpartitioning = (
idx_t*) partitioning;
8169 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
8170 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
8171 mpartitioning =
new idx_t[n];
8174#ifndef MFEM_USE_METIS_5
8177 METIS_SetDefaultOptions(options);
8178 options[METIS_OPTION_CONTIG] = 1;
8186 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
8191 if (part_method >= 0 && part_method <= 2)
8193 for (i = 0; i < n; i++)
8199 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
8205 if (part_method == 0 || part_method == 3)
8207#ifndef MFEM_USE_METIS_5
8236 " error in METIS_PartGraphRecursive!");
8243 if (part_method == 1 || part_method == 4)
8245#ifndef MFEM_USE_METIS_5
8274 " error in METIS_PartGraphKway!");
8281 if (part_method == 2 || part_method == 5)
8283#ifndef MFEM_USE_METIS_5
8296 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
8313 " error in METIS_PartGraphKway!");
8321 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
8325 nparts = (int) mparts;
8326 if (mpartitioning != (
idx_t*)partitioning)
8330 partitioning[k] = mpartitioning[k];
8337 delete[] mpartitioning;
8353 auto count_partition_elements = [&]()
8355 for (i = 0; i < nparts; i++)
8363 psize[partitioning[i]].one++;
8367 for (i = 0; i < nparts; i++)
8369 if (psize[i].one == 0) { empty_parts++; }
8373 count_partition_elements();
8381 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
8382 << empty_parts <<
" empty parts!"
8383 <<
" Applying a simple fix ..." << endl;
8388 for (i = nparts-1; i > nparts-1-empty_parts; i--)
8395 for (i = nparts-1; i > nparts-1-empty_parts; i--)
8397 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
8403 partitioning[j] = psize[nparts-1-i].two;
8410 count_partition_elements();
8414 return partitioning;
8418 mfem_error(
"Mesh::GeneratePartitioning(...): "
8419 "MFEM was compiled without Metis.");
8433 int num_elem, *i_elem_elem, *j_elem_elem;
8435 num_elem = elem_elem.
Size();
8436 i_elem_elem = elem_elem.
GetI();
8437 j_elem_elem = elem_elem.
GetJ();
8442 int stack_p, stack_top_p, elem;
8446 for (i = 0; i < num_elem; i++)
8448 if (partitioning[i] > num_part)
8450 num_part = partitioning[i];
8457 for (i = 0; i < num_part; i++)
8464 for (elem = 0; elem < num_elem; elem++)
8466 if (component[elem] >= 0)
8471 component[elem] = num_comp[partitioning[elem]]++;
8473 elem_stack[stack_top_p++] = elem;
8475 for ( ; stack_p < stack_top_p; stack_p++)
8477 i = elem_stack[stack_p];
8478 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
8481 if (partitioning[k] == partitioning[i])
8483 if (component[k] < 0)
8485 component[k] = component[i];
8486 elem_stack[stack_top_p++] = k;
8488 else if (component[k] != component[i])
8500 int i, n_empty, n_mcomp;
8508 n_empty = n_mcomp = 0;
8509 for (i = 0; i < num_comp.
Size(); i++)
8510 if (num_comp[i] == 0)
8514 else if (num_comp[i] > 1)
8521 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
8522 <<
"The following subdomains are empty :\n";
8523 for (i = 0; i < num_comp.
Size(); i++)
8524 if (num_comp[i] == 0)
8532 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
8533 <<
"The following subdomains are NOT connected :\n";
8534 for (i = 0; i < num_comp.
Size(); i++)
8535 if (num_comp[i] > 1)
8541 if (n_empty == 0 && n_mcomp == 0)
8542 mfem::out <<
"Mesh::CheckPartitioning(...) : "
8543 "All subdomains are connected." << endl;
8567 c(0) =
a[0]*
a[3]-
a[1]*
a[2];
8568 c(1) =
a[0]*
b[3]-
a[1]*
b[2]+
b[0]*
a[3]-
b[1]*
a[2];
8569 c(2) =
b[0]*
b[3]-
b[1]*
b[2];
8590 c(0) = (
a[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
8591 a[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
8592 a[2] * (
a[3] *
a[7] -
a[4] *
a[6]));
8594 c(1) = (
b[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
8595 b[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
8596 b[2] * (
a[3] *
a[7] -
a[4] *
a[6]) +
8598 a[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
8599 a[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
8600 a[2] * (
b[3] *
a[7] -
b[4] *
a[6]) +
8602 a[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
8603 a[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
8604 a[2] * (
a[3] *
b[7] -
a[4] *
b[6]));
8606 c(2) = (
a[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
8607 a[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
8608 a[2] * (
b[3] *
b[7] -
b[4] *
b[6]) +
8610 b[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
8611 b[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
8612 b[2] * (
a[3] *
b[7] -
a[4] *
b[6]) +
8614 b[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
8615 b[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
8616 b[2] * (
b[3] *
a[7] -
b[4] *
a[6]));
8618 c(3) = (
b[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
8619 b[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
8620 b[2] * (
b[3] *
b[7] -
b[4] *
b[6]));
8666 real_t a = z(2),
b = z(1), c = z(0);
8674 x(0) = x(1) = -0.5 *
b /
a;
8679 x(0) = -(x(1) = fabs(0.5 * sqrt(D) /
a));
8687 t = -0.5 * (
b + sqrt(D));
8691 t = -0.5 * (
b - sqrt(D));
8705 real_t a = z(2)/z(3),
b = z(1)/z(3), c = z(0)/z(3);
8709 real_t R = (2 *
a *
a *
a - 9 *
a *
b + 27 * c) / 54;
8717 x(0) = x(1) = x(2) = -
a / 3;
8725 x(0) = -2 * sqrtQ -
a / 3;
8726 x(1) = x(2) = sqrtQ -
a / 3;
8730 x(0) = x(1) = - sqrtQ -
a / 3;
8731 x(2) = 2 * sqrtQ -
a / 3;
8738 real_t theta = acos(R / sqrt(Q3));
8741 x0 = A * cos(theta / 3) -
a / 3;
8742 x1 = A * cos((theta + 2.0 * M_PI) / 3) -
a / 3;
8743 x2 = A * cos((theta - 2.0 * M_PI) / 3) -
a / 3;
8768 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
8772 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
8774 x(0) = A + Q / A -
a / 3;
8783 const real_t factor,
const int Dim)
8786 c(0) = c0 * (1.0 - pow(factor, -Dim));
8788 for (
int j = 0; j < nr; j++)
8800 c(0) = c0 * (1.0 - pow(factor, Dim));
8802 for (
int j = 0; j < nr; j++)
8821 const real_t factor = 2.0;
8836 for (
int k = 0; k < nv; k++)
8839 V(j, k) = displacements(v[k]+j*nvs);
8869 for (
int j = 0; j < nv; j++)
8895 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
8898 vertices[i](j) += displacements(j*nv+i);
8906 for (
int i = 0; i < nv; i++)
8909 vert_coord(j*nv+i) =
vertices[i](j);
8915 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
8918 vertices[i](j) = vert_coord(j*nv+i);
8948 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
8965 (*Nodes) += displacements;
8977 node_coord = (*Nodes);
8989 (*Nodes) = node_coord;
9052 for (j = 1; j < n; j++)
9089 int quad_counter = 0;
9111 const int attr =
elements[i]->GetAttribute();
9112 int *v =
elements[i]->GetVertices();
9118 for (
int ei = 0; ei < 3; ei++)
9120 for (
int k = 0; k < 2; k++)
9128 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
9130 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
9132 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
9134 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
9138 const int qe = quad_counter;
9142 for (
int ei = 0; ei < 4; ei++)
9144 for (
int k = 0; k < 2; k++)
9152 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
9154 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
9156 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
9158 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
9162 MFEM_ABORT(
"unknown element type: " << el_type);
9172 const int attr =
boundary[i]->GetAttribute();
9173 int *v =
boundary[i]->GetVertices();
9182 static const real_t A = 0.0, B = 0.5, C = 1.0;
9183 static real_t tri_children[2*3*4] =
9190 static real_t quad_children[2*4*4] =
9204 for (
int i = 0; i <
elements.Size(); i++)
9225 if (!
Nodes || update_nodes)
9255 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
9258 int NumOfQuadFaces = 0;
9264 for (
int i = 0; i <
faces.Size(); i++)
9268 f2qf[i] = NumOfQuadFaces;
9279 int hex_counter = 0;
9282 for (
int i = 0; i <
elements.Size(); i++)
9291 int pyr_counter = 0;
9294 for (
int i = 0; i <
elements.Size(); i++)
9312 DSTable *v_to_v_ptr = v_to_v_p;
9328 std::sort(row_start, J_v2v.
end());
9331 for (
int i = 0; i < J_v2v.
Size(); i++)
9333 e2v[J_v2v[i].two] = i;
9346 it.SetIndex(e2v[it.Index()]);
9356 const int oelem = oface + NumOfQuadFaces;
9369 const int attr =
elements[i]->GetAttribute();
9370 int *v =
elements[i]->GetVertices();
9377 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
9385 for (
int ei = 0; ei < 6; ei++)
9387 for (
int k = 0; k < 2; k++)
9397 const int rt_algo = 1;
9410 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
9411 sqr(J(1,0)-J(1,1)-J(1,2)) +
9412 sqr(J(2,0)-J(2,1)-J(2,2));
9415 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
9416 sqr(J(1,1)-J(1,0)-J(1,2)) +
9417 sqr(J(2,1)-J(2,0)-J(2,2));
9418 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
9420 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
9421 sqr(J(1,2)-J(1,0)-J(1,1)) +
9422 sqr(J(2,2)-J(2,0)-J(2,1));
9423 if (len_sqr < min_len) { rt = 2; }
9428 real_t Em_data[18], Js_data[9], Jp_data[9];
9433 for (
int s = 0;
s < 3;
s++)
9435 for (
int t = 0;
t < 3;
t++)
9437 Em(
t,
s) = 0.5*J(
t,
s);
9440 for (
int t = 0;
t < 3;
t++)
9442 Em(
t,3) = 0.5*(J(
t,0)+J(
t,1));
9443 Em(
t,4) = 0.5*(J(
t,0)+J(
t,2));
9444 Em(
t,5) = 0.5*(J(
t,1)+J(
t,2));
9448 for (
int t = 0;
t < 3;
t++)
9450 Js(
t,0) = Em(
t,5)-Em(
t,0);
9451 Js(
t,1) = Em(
t,1)-Em(
t,0);
9452 Js(
t,2) = Em(
t,2)-Em(
t,0);
9456 for (
int t = 0;
t < 3;
t++)
9458 Js(
t,0) = Em(
t,5)-Em(
t,0);
9459 Js(
t,1) = Em(
t,2)-Em(
t,0);
9460 Js(
t,2) = Em(
t,4)-Em(
t,0);
9464 kappa_min = std::max(ar1, ar2);
9468 for (
int t = 0;
t < 3;
t++)
9470 Js(
t,0) = Em(
t,0)-Em(
t,1);
9471 Js(
t,1) = Em(
t,4)-Em(
t,1);
9472 Js(
t,2) = Em(
t,2)-Em(
t,1);
9476 for (
int t = 0;
t < 3;
t++)
9478 Js(
t,0) = Em(
t,2)-Em(
t,1);
9479 Js(
t,1) = Em(
t,4)-Em(
t,1);
9480 Js(
t,2) = Em(
t,5)-Em(
t,1);
9484 kappa = std::max(ar1, ar2);
9485 if (
kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
9488 for (
int t = 0;
t < 3;
t++)
9490 Js(
t,0) = Em(
t,0)-Em(
t,2);
9491 Js(
t,1) = Em(
t,1)-Em(
t,2);
9492 Js(
t,2) = Em(
t,3)-Em(
t,2);
9496 for (
int t = 0;
t < 3;
t++)
9498 Js(
t,0) = Em(
t,1)-Em(
t,2);
9499 Js(
t,1) = Em(
t,5)-Em(
t,2);
9500 Js(
t,2) = Em(
t,3)-Em(
t,2);
9504 kappa = std::max(ar1, ar2);
9505 if (
kappa < kappa_min) { rt = 2; }
9508 static const int mv_all[3][4][4] =
9510 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
9511 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
9512 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
9514 const int (&mv)[4][4] = mv_all[rt];
9516#ifndef MFEM_USE_MEMALLOC
9518 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
9520 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
9522 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
9524 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
9526 for (
int k = 0; k < 4; k++)
9528 new_elements[j+4+k] =
9529 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
9530 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
9534 new_elements[j+0] = tet =
TetMemory.Alloc();
9535 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
9537 new_elements[j+1] = tet =
TetMemory.Alloc();
9538 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
9540 new_elements[j+2] = tet =
TetMemory.Alloc();
9541 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
9543 new_elements[j+3] = tet =
TetMemory.Alloc();
9544 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
9546 for (
int k = 0; k < 4; k++)
9548 new_elements[j+4+k] = tet =
TetMemory.Alloc();
9549 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
9550 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
9553 for (
int k = 0; k < 4; k++)
9558 for (
int k = 0; k < 4; k++)
9572 for (
int fi = 2; fi < 5; fi++)
9574 for (
int k = 0; k < 4; k++)
9581 for (
int ei = 0; ei < 9; ei++)
9583 for (
int k = 0; k < 2; k++)
9590 const int qf2 = f2qf[
f[2]];
9591 const int qf3 = f2qf[
f[3]];
9592 const int qf4 = f2qf[
f[4]];
9595 new Wedge(v[0], oedge+e[0], oedge+e[2],
9596 oedge+e[6], oface+qf2, oface+qf4, attr);
9599 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
9600 oface+qf3, oface+qf4, oface+qf2, attr);
9603 new Wedge(oedge+e[0], v[1], oedge+e[1],
9604 oface+qf2, oedge+e[7], oface+qf3, attr);
9607 new Wedge(oedge+e[2], oedge+e[1], v[2],
9608 oface+qf4, oface+qf3, oedge+e[8], attr);
9611 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
9612 v[3], oedge+e[3], oedge+e[5], attr);
9615 new Wedge(oface+qf3, oface+qf4, oface+qf2,
9616 oedge+e[4], oedge+e[5], oedge+e[3], attr);
9619 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
9620 oedge+e[3], v[4], oedge+e[4], attr);
9623 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
9624 oedge+e[5], oedge+e[4], v[5], attr);
9633 for (
int fi = 0; fi < 1; fi++)
9635 for (
int k = 0; k < 4; k++)
9642 for (
int ei = 0; ei < 8; ei++)
9644 for (
int k = 0; k < 2; k++)
9651 const int qf0 = f2qf[
f[0]];
9654 new Pyramid(v[0], oedge+e[0], oface+qf0,
9655 oedge+e[3], oedge+e[4], attr);
9658 new Pyramid(oedge+e[0], v[1], oedge+e[1],
9659 oface+qf0, oedge+e[5], attr);
9662 new Pyramid(oface+qf0, oedge+e[1], v[2],
9663 oedge+e[2], oedge+e[6], attr);
9666 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
9667 v[3], oedge+e[7], attr);
9670 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
9671 oedge+e[7], v[4], attr);
9674 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
9675 oedge+e[4], oface+qf0, attr);
9677#ifndef MFEM_USE_MEMALLOC
9679 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
9683 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
9687 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
9691 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
9695 new_elements[j++] = tet =
TetMemory.Alloc();
9696 tet->
Init(oedge+e[0], oedge+e[4], oedge+e[5],
9699 new_elements[j++] = tet =
TetMemory.Alloc();
9700 tet->
Init(oedge+e[1], oedge+e[5], oedge+e[6],
9703 new_elements[j++] = tet =
TetMemory.Alloc();
9704 tet->
Init(oedge+e[2], oedge+e[6], oedge+e[7],
9707 new_elements[j++] = tet =
TetMemory.Alloc();
9708 tet->
Init(oedge+e[3], oedge+e[7], oedge+e[4],
9717 const int he = hex_counter;
9722 if (f2qf.
Size() == 0)
9728 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[
f[k]]; }
9734 for (
int fi = 0; fi < 6; fi++)
9736 for (
int k = 0; k < 4; k++)
9743 for (
int ei = 0; ei < 12; ei++)
9745 for (
int k = 0; k < 2; k++)
9753 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
9754 oedge+e[3], oedge+e[8], oface+qf[1],
9755 oelem+he, oface+qf[4], attr);
9758 oface+qf[0], oface+qf[1], oedge+e[9],
9759 oface+qf[2], oelem+he, attr);
9761 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
9762 oedge+e[2], oelem+he, oface+qf[2],
9763 oedge+e[10], oface+qf[3], attr);
9765 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
9766 v[3], oface+qf[4], oelem+he,
9767 oface+qf[3], oedge+e[11], attr);
9769 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
9770 oface+qf[4], v[4], oedge+e[4],
9771 oface+qf[5], oedge+e[7], attr);
9773 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
9774 oelem+he, oedge+e[4], v[5],
9775 oedge+e[5], oface+qf[5], attr);
9777 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
9778 oface+qf[3], oface+qf[5], oedge+e[5],
9779 v[6], oedge+e[6], attr);
9781 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
9782 oedge+e[11], oedge+e[7], oface+qf[5],
9783 oedge+e[6], v[7], attr);
9788 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
9800 const int attr =
boundary[i]->GetAttribute();
9801 int *v =
boundary[i]->GetVertices();
9808 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
9815 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
9817 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
9819 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
9821 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
9829 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
9831 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
9833 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
9835 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
9839 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
9845 static const real_t A = 0.0, B = 0.5, C = 1.0, D = -1.0;
9846 static real_t tet_children[3*4*16] =
9848 A,A,A, B,A,A, A,B,A, A,A,B,
9849 B,A,A, C,A,A, B,B,A, B,A,B,
9850 A,B,A, B,B,A, A,C,A, A,B,B,
9851 A,A,B, B,A,B, A,B,B, A,A,C,
9856 B,A,A, A,B,B, A,B,A, A,A,B,
9857 B,A,A, A,B,B, A,A,B, B,A,B,
9858 B,A,A, A,B,B, B,A,B, B,B,A,
9859 B,A,A, A,B,B, B,B,A, A,B,A,
9861 A,B,A, B,A,A, B,A,B, A,A,B,
9862 A,B,A, A,A,B, B,A,B, A,B,B,
9863 A,B,A, A,B,B, B,A,B, B,B,A,
9864 A,B,A, B,B,A, B,A,B, B,A,A,
9866 A,A,B, B,A,A, A,B,A, B,B,A,
9867 A,A,B, A,B,A, A,B,B, B,B,A,
9868 A,A,B, A,B,B, B,A,B, B,B,A,
9869 A,A,B, B,A,B, B,A,A, B,B,A
9871 static real_t pyr_children[3*5*10] =
9873 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
9874 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
9875 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
9876 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
9877 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
9878 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
9879 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
9880 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
9881 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
9882 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
9884 static real_t pri_children[3*6*8] =
9886 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
9887 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
9888 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
9889 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
9890 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
9891 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
9892 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
9893 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
9895 static real_t hex_children[3*8*8] =
9897 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
9898 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
9899 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
9900 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
9901 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
9902 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
9903 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
9904 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
9916 for (
int i = 0; i <
elements.Size(); i++)
9947 int i, j, ind, nedges;
9954 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
9968 for (j = 0; j < marked_el.
Size(); j++)
9973 int new_v = cnv + j, new_e = cne + j;
9982 static real_t seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
9984 UseExternalData(seg_children, 1, 2, 3);
9997 int *edge1 =
new int[nedges];
9998 int *edge2 =
new int[nedges];
9999 int *middle =
new int[nedges];
10001 for (i = 0; i < nedges; i++)
10003 edge1[i] = edge2[i] = middle[i] = -1;
10009 for (j = 1; j < v.
Size(); j++)
10011 ind = v_to_v(v[j-1], v[j]);
10012 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
10014 ind = v_to_v(v[0], v[v.
Size()-1]);
10015 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
10019 for (i = 0; i < marked_el.
Size(); i++)
10025 int need_refinement;
10028 need_refinement = 0;
10029 for (i = 0; i < nedges; i++)
10031 if (middle[i] != -1 && edge1[i] != -1)
10033 need_refinement = 1;
10038 while (need_refinement == 1);
10041 int v1[2], v2[2], bisect, temp;
10043 for (i = 0; i < temp; i++)
10046 bisect = v_to_v(v[0], v[1]);
10047 if (middle[bisect] != -1)
10051 v1[0] = v[0]; v1[1] = middle[bisect];
10052 v2[0] = middle[bisect]; v2[1] = v[1];
10058 mfem_error(
"Only bisection of segment is implemented"
10082 MFEM_VERIFY(
GetNE() == 0 ||
10084 "tetrahedral mesh is not marked for refinement:"
10085 " call Finalize(true)");
10092 for (i = 0; i < marked_el.
Size(); i++)
10098 for (i = 0; i < marked_el.
Size(); i++)
10107 for (i = 0; i < marked_el.
Size(); i++)
10124 int need_refinement;
10129 need_refinement = 0;
10137 if (
elements[i]->NeedRefinement(v_to_v))
10139 need_refinement = 1;
10144 while (need_refinement == 1);
10151 need_refinement = 0;
10153 if (
boundary[i]->NeedRefinement(v_to_v))
10155 need_refinement = 1;
10159 while (need_refinement == 1);
10190 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
10191 "not supported. Project the NURBS to Nodes first.");
10201 if (!refinements.
Size())
10222 Swap(*mesh2,
false);
10234 const int *fine,
int nfine,
int op)
10236 real_t error = (op == 3) ? std::pow(elem_error[fine[0]],
10237 2.0) : elem_error[fine[0]];
10239 for (
int i = 1; i < nfine; i++)
10241 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
10243 real_t err_fine = elem_error[fine[i]];
10246 case 0: error = std::min(error, err_fine);
break;
10247 case 1: error += err_fine;
break;
10248 case 2: error = std::max(error, err_fine);
break;
10249 case 3: error += std::pow(err_fine, 2.0);
break;
10250 default: MFEM_ABORT(
"Invalid operation.");
10253 return (op == 3) ? std::sqrt(error) : error;
10257 real_t threshold,
int nc_limit,
int op)
10259 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
10260 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
10261 "Project the NURBS to Nodes first.");
10274 for (
int i = 0; i < dt.
Size(); i++)
10276 if (nc_limit > 0 && !level_ok[i]) {
continue; }
10281 if (error < threshold) { derefs.
Append(i); }
10284 if (!derefs.
Size()) {
return false; }
10291 Swap(*mesh2,
false);
10305 int nc_limit,
int op)
10315 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
10322 int nc_limit,
int op)
10325 for (
int i = 0; i < tmp.
Size(); i++)
10327 tmp[i] = elem_error(i);
10370 : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
10416#ifdef MFEM_USE_MEMALLOC
10444 for (
int i = 0; i < elem_array.
Size(); i++)
10446 if (elem_array[i]->GetGeometryType() == geom)
10451 elem_vtx.
SetSize(nv*num_elems);
10455 for (
int i = 0; i < elem_array.
Size(); i++)
10461 elem_vtx.
Append(loc_vtx);
10469 for (
int i = 0; i < nelem; i++) { list[i] = i; }
10485 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
10497 default: MFEM_ABORT(
"internal error");
10511 bool noInitialCoarsening =
true;
10512 for (
auto f : initialCoarsening)
10514 noInitialCoarsening = (noInitialCoarsening &&
f == 1);
10517 if (noInitialCoarsening)
10537 bool divisible =
true;
10538 for (
int i=0; i<rf.
Size(); ++i)
10541 divisible = divisible && cf * rf[i] == initialCoarsening[i];
10544 MFEM_VERIFY(divisible,
"Invalid coarsening");
10559 int nonconforming,
int nc_limit)
10569 else if (nonconforming < 0)
10590 for (
int i = 0; i < refinements.
Size(); i++)
10592 el_to_refine[i] = refinements[i].index;
10596 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
10597 if (rt == 1 || rt == 2 || rt == 4)
10601 else if (rt == 3 || rt == 5 || rt == 6)
10619 for (
int i = 0; i < el_to_refine.
Size(); i++)
10621 refinements[i] =
Refinement(el_to_refine[i]);
10628 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
10629 "Please project the NURBS to Nodes first, with SetCurvature().");
10632 MFEM_VERIFY(
ncmesh != NULL ||
dynamic_cast<const ParMesh*
>(
this) == NULL,
10633 "Sorry, converting a conforming ParMesh to an NC mesh is "
10641 (simplices_nonconforming && (
meshgen & 0x1)) )
10654 for (
int i = 0; i <
GetNE(); i++)
10661 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
10673 for (
int i = 0; i <
GetNE(); i++)
10676 bool refine =
false;
10677 for (
int j = 0; j < v.
Size(); j++)
10680 for (
int l = 0; l <
spaceDim; l++)
10685 if (dist <= eps*eps) { refine =
true;
break; }
10696 int nonconforming,
int nc_limit)
10698 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
10700 for (
int i = 0; i <
GetNE(); i++)
10702 if (elem_error[i] > threshold)
10716 int nonconforming,
int nc_limit)
10719 elem_error.
Size());
10720 return RefineByError(tmp, threshold, nonconforming, nc_limit);
10725 int *edge1,
int *edge2,
int *middle)
10728 int v[2][4], v_new, bisect,
t;
10740 bisect = v_to_v(vert[0], vert[1]);
10741 MFEM_ASSERT(bisect >= 0,
"");
10743 if (middle[bisect] == -1)
10746 for (
int d = 0; d <
spaceDim; d++)
10754 if (edge1[bisect] == i)
10756 edge1[bisect] = edge2[bisect];
10759 middle[bisect] = v_new;
10763 v_new = middle[bisect];
10766 edge1[bisect] = -1;
10771 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
10772 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
10793 bisect = v_to_v(v[1][0], v[1][1]);
10794 MFEM_ASSERT(bisect >= 0,
"");
10796 if (edge1[bisect] == i)
10800 else if (edge2[bisect] == i)
10809 MFEM_ABORT(
"Bisection for now works only for triangles.");
10816 int v[2][4], v_new, bisect,
t;
10826 "TETRAHEDRON element is not marked for refinement.");
10831 bisect = v_to_v.
FindId(vert[0], vert[1]);
10835 for (
int j = 0; j < 3; j++)
10848 int type, old_redges[2], flag;
10851 int new_type, new_redges[2][2];
10854 new_redges[0][0] = 2;
10855 new_redges[0][1] = 1;
10856 new_redges[1][0] = 2;
10857 new_redges[1][1] = 1;
10858 int tr1 = -1, tr2 = -1;
10859 switch (old_redges[0])
10862 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
10867 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
10871 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
10874 switch (old_redges[1])
10877 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
10882 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
10886 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
10893#ifdef MFEM_USE_MEMALLOC
10929 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
10936 int v[2][3], v_new, bisect,
t;
10947 bisect = v_to_v.
FindId(vert[0], vert[1]);
10948 MFEM_ASSERT(bisect >= 0,
"");
10950 MFEM_ASSERT(v_new != -1,
"");
10954 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
10955 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
10965 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
10971 int *edge1,
int *edge2,
int *middle)
10974 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
10983 bisect[0] = v_to_v(v[0],v[1]);
10984 bisect[1] = v_to_v(v[1],v[2]);
10985 bisect[2] = v_to_v(v[0],v[2]);
10986 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
10988 for (j = 0; j < 3; j++)
10990 if (middle[bisect[j]] == -1)
10993 for (
int d = 0; d <
spaceDim; d++)
11001 if (edge1[bisect[j]] == i)
11003 edge1[bisect[j]] = edge2[bisect[j]];
11006 middle[bisect[j]] = v_new[j];
11010 v_new[j] = middle[bisect[j]];
11013 edge1[bisect[j]] = -1;
11019 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
11020 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
11021 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
11022 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
11056 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
11092 for (
int i = 0; i < elem_geoms.
Size(); i++)
11100 std::map<unsigned, int> mat_no;
11104 for (
int j = 0; j <
elements.Size(); j++)
11107 unsigned code =
elements[j]->GetTransform();
11110 int &matrix = mat_no[code];
11111 if (!matrix) { matrix = mat_no.size(); }
11121 std::map<unsigned, int>::iterator it;
11122 for (it = mat_no.begin(); it != mat_no.end(); ++it)
11136 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
11137 " geom = " << geom);
11147 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
11156 os <<
"areamesh2\n\n";
11160 os <<
"curved_areamesh2\n\n";
11169 os <<
boundary[i]->GetAttribute();
11170 for (j = 0; j < v.
Size(); j++)
11172 os <<
' ' << v[j] + 1;
11184 for (j = 0; j < v.
Size(); j++)
11186 os <<
' ' << v[j] + 1;
11198 for (j = 1; j <
Dim; j++)
11215 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
11223 os <<
"NETGEN_Neutral_Format\n";
11228 for (j = 0; j <
Dim; j++)
11241 os <<
elements[i]->GetAttribute();
11242 for (j = 0; j < nv; j++)
11244 os <<
' ' << ind[j]+1;
11255 os <<
boundary[i]->GetAttribute();
11256 for (j = 0; j < nv; j++)
11258 os <<
' ' << ind[j]+1;
11270 <<
" 0 0 0 0 0 0 0\n"
11271 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
11273 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
11274 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11278 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
11284 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
11285 for (j = 0; j < nv; j++)
11287 os <<
' ' << ind[j]+1;
11296 os <<
boundary[i]->GetAttribute();
11297 for (j = 0; j < nv; j++)
11299 os <<
' ' << ind[j]+1;
11301 os <<
" 1.0 1.0 1.0 1.0\n";
11310 const std::string &comments)
const
11344 os <<
"\n# mesh curvature GridFunction";
11349 os <<
"\nmfem_mesh_end" << endl;
11356 os << (!set_names && section_delimiter.empty()
11357 ?
"MFEM mesh v1.0\n" :
11358 (!set_names ?
"MFEM mesh v1.2\n" :
"MFEM mesh v1.3\n"));
11360 if (set_names && section_delimiter.empty())
11362 section_delimiter =
"mfem_mesh_end";
11366 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
11369 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
11374 "# TETRAHEDRON = 4\n"
11380 os <<
"\ndimension\n" <<
Dim;
11390 os <<
"\nattribute_sets\n";
11402 os <<
"\nbdr_attribute_sets\n";
11427 if (!section_delimiter.empty())
11430 << section_delimiter << endl;
11435 const int version,
const std::string &comments)
const
11437 MFEM_VERIFY(version == 10 || version == 11,
"Invalid NURBS mesh version");
11442 os <<
"MFEM NURBS mesh v" << int(version / 10) <<
"." << version % 10 <<
"\n";
11445 if (!comments.empty()) { os <<
'\n' << comments <<
'\n'; }
11448 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
11454 os <<
"\ndimension\n" <<
Dim
11471 int ki = e_to_k[i];
11476 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
11483 ofstream ofs(fname);
11484 ofs.precision(precision);
11488#ifdef MFEM_USE_ADIOS2
11498 "# vtk DataFile Version 3.0\n"
11499 "Generated by MFEM\n"
11501 "DATASET UNSTRUCTURED_GRID\n";
11514 for ( ; j < 3; j++)
11525 for (
int i = 0; i <
Nodes->
FESpace()->GetNDofs(); i++)
11530 os << (*Nodes)(vdofs[0]);
11534 os <<
' ' << (*Nodes)(vdofs[j]);
11536 for ( ; j < 3; j++)
11550 size +=
elements[i]->GetNVertices() + 1;
11555 const int *v =
elements[i]->GetVertices();
11556 const int nv =
elements[i]->GetNVertices();
11560 for (
int j = 0; j < nv; j++)
11562 os <<
' ' << v[perm ? perm[j] : j];
11575 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
11576 "Point meshes should have a single dof per element");
11577 size += dofs.
Size() + 1;
11582 if (!strcmp(fec_name,
"Linear") ||
11583 !strcmp(fec_name,
"H1_0D_P1") ||
11584 !strcmp(fec_name,
"H1_1D_P1") ||
11585 !strcmp(fec_name,
"H1_2D_P1") ||
11586 !strcmp(fec_name,
"H1_3D_P1"))
11590 else if (!strcmp(fec_name,
"Quadratic") ||
11591 !strcmp(fec_name,
"H1_1D_P2") ||
11592 !strcmp(fec_name,
"H1_2D_P2") ||
11593 !strcmp(fec_name,
"H1_3D_P2"))
11599 mfem::err <<
"Mesh::PrintVTK : can not save '"
11600 << fec_name <<
"' elements!" << endl;
11609 for (
int j = 0; j < dofs.
Size(); j++)
11611 os <<
' ' << dofs[j];
11614 else if (order == 2)
11616 const int *vtk_mfem;
11617 switch (
elements[i]->GetGeometryType())
11631 for (
int j = 0; j < dofs.
Size(); j++)
11633 os <<
' ' << dofs[vtk_mfem[j]];
11643 int vtk_cell_type = 5;
11647 os << vtk_cell_type <<
'\n';
11652 <<
"SCALARS material int\n"
11653 <<
"LOOKUP_TABLE default\n";
11656 os <<
elements[i]->GetAttribute() <<
'\n';
11663 bool high_order_output,
11664 int compression_level,
11667 int ref = (high_order_output &&
Nodes)
11670 fname = fname +
".vtu";
11671 std::fstream os(fname.c_str(),std::ios::out);
11672 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
11673 if (compression_level != 0)
11675 os <<
" compressor=\"vtkZLibDataCompressor\"";
11678 os <<
"<UnstructuredGrid>\n";
11679 PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
11680 os <<
"</Piece>\n";
11681 os <<
"</UnstructuredGrid>\n";
11682 os <<
"</VTKFile>" << std::endl;
11689 bool high_order_output,
11690 int compression_level)
11692 PrintVTU(fname, format, high_order_output, compression_level,
true);
11696 bool high_order_output,
int compression_level,
11704 std::vector<char> buf;
11706 auto get_geom = [&](
int i)
11714 int np = 0, nc_ref = 0;
11715 for (
int i = 0; i < ne; i++)
11724 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
11725 << (high_order_output ? ne : nc_ref) <<
"\">\n";
11728 os <<
"<Points>\n";
11729 os <<
"<DataArray type=\"" << type_str
11730 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
11731 for (
int i = 0; i < ne; i++)
11744 for (
int j = 0; j < pmat.
Width(); j++)
11770 os <<
"</DataArray>" << std::endl;
11771 os <<
"</Points>" << std::endl;
11773 os <<
"<Cells>" << std::endl;
11774 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
11775 << fmt_str <<
"\">" << std::endl;
11777 std::vector<int> offset;
11780 if (high_order_output)
11783 for (
int iel = 0; iel < ne; iel++)
11787 int nnodes = local_connectivity.
Size();
11788 for (
int i=0; i<nnodes; ++i)
11795 offset.push_back(np);
11801 for (
int i = 0; i < ne; i++)
11807 for (
int j = 0; j < RG.
Size(); )
11810 offset.push_back(coff);
11812 for (
int k = 0; k < nv; k++, j++)
11826 os <<
"</DataArray>" << std::endl;
11828 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
11829 << fmt_str <<
"\">" << std::endl;
11831 for (
size_t ii=0; ii<offset.size(); ii++)
11839 os <<
"</DataArray>" << std::endl;
11840 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
11841 << fmt_str <<
"\">" << std::endl;
11843 const int *vtk_geom_map =
11845 for (
int i = 0; i < ne; i++)
11848 uint8_t vtk_cell_type = 5;
11850 vtk_cell_type = vtk_geom_map[geom];
11852 if (high_order_output)
11861 for (
int j = 0; j < RG.
Size(); j += nv)
11871 os <<
"</DataArray>" << std::endl;
11872 os <<
"</Cells>" << std::endl;
11874 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
11875 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
11876 << fmt_str <<
"\">" << std::endl;
11877 for (
int i = 0; i < ne; i++)
11880 if (high_order_output)
11899 os <<
"</DataArray>" << std::endl;
11900 os <<
"</CellData>" << std::endl;
11911 "# vtk DataFile Version 3.0\n"
11912 "Generated by MFEM\n"
11914 "DATASET UNSTRUCTURED_GRID\n";
11919 os <<
"FIELD FieldData 1\n"
11929 np = nc = size = 0;
11930 for (
int i = 0; i <
GetNE(); i++)
11939 os <<
"POINTS " << np <<
" double\n";
11941 for (
int i = 0; i <
GetNE(); i++)
11948 for (
int j = 0; j < pmat.
Width(); j++)
11950 os << pmat(0, j) <<
' ';
11953 os << pmat(1, j) <<
' ';
11965 os << 0.0 <<
' ' << 0.0;
11972 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
11974 for (
int i = 0; i <
GetNE(); i++)
11981 for (
int j = 0; j < RG.
Size(); )
11984 for (
int k = 0; k < nv; k++, j++)
11986 os <<
' ' << np + RG[j];
11992 os <<
"CELL_TYPES " << nc <<
'\n';
11993 for (
int i = 0; i <
GetNE(); i++)
12001 for (
int j = 0; j < RG.
Size(); j += nv)
12003 os << vtk_cell_type <<
'\n';
12007 os <<
"CELL_DATA " << nc <<
'\n'
12008 <<
"SCALARS material int\n"
12009 <<
"LOOKUP_TABLE default\n";
12010 for (
int i = 0; i <
GetNE(); i++)
12018 os << attr <<
'\n';
12025 srand((
unsigned)time(0));
12027 int el0 = (int)floor(
a *
GetNE());
12029 os <<
"SCALARS element_coloring int\n"
12030 <<
"LOOKUP_TABLE default\n";
12031 for (
int i = 0; i <
GetNE(); i++)
12038 os << coloring[i] + 1 <<
'\n';
12044 os <<
"POINT_DATA " << np <<
'\n' << flush;
12049 int delete_el_to_el = (
el_to_el) ? (0) : (1);
12051 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
12054 const int *i_el_el = el_el.
GetI();
12055 const int *j_el_el = el_el.
GetJ();
12060 stack_p = stack_top_p = 0;
12061 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
12063 if (colors[el] != -2)
12069 el_stack[stack_top_p++] = el;
12071 for ( ; stack_p < stack_top_p; stack_p++)
12073 int i = el_stack[stack_p];
12074 int num_nb = i_el_el[i+1] - i_el_el[i];
12075 if (max_num_col < num_nb + 1)
12077 max_num_col = num_nb + 1;
12079 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
12081 int k = j_el_el[j];
12082 if (colors[k] == -2)
12085 el_stack[stack_top_p++] = k;
12093 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
12095 int i = el_stack[stack_p], col;
12097 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
12099 col = colors[j_el_el[j]];
12102 col_marker[col] = 1;
12106 for (col = 0; col < max_num_col; col++)
12107 if (col_marker[col] == 0)
12115 if (delete_el_to_el)
12123 int elem_attr)
const
12125 if (
Dim != 3 &&
Dim != 2) {
return; }
12127 int i, j, k, l, nv, nbe, *v;
12129 os <<
"MFEM mesh v1.0\n";
12133 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12138 "# TETRAHEDRON = 4\n"
12143 os <<
"\ndimension\n" <<
Dim
12147 os << int((elem_attr) ? partitioning[i]+1 :
elements[i]->GetAttribute())
12148 <<
' ' <<
elements[i]->GetGeometryType();
12151 for (j = 0; j < nv; j++)
12163 l = partitioning[l];
12178 os <<
"\nboundary\n" << nbe <<
'\n';
12184 l = partitioning[l];
12187 nv =
faces[i]->GetNVertices();
12188 v =
faces[i]->GetVertices();
12189 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
12190 for (j = 0; j < nv; j++)
12197 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
12198 for (j = nv-1; j >= 0; j--)
12209 nv =
faces[i]->GetNVertices();
12210 v =
faces[i]->GetVertices();
12211 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
12212 for (j = 0; j < nv; j++)
12243 int interior_faces)
12245 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
12246 if (
Dim != 3 &&
Dim != 2) {
return; }
12255 int nv =
elements[i]->GetNVertices();
12256 const int *ind =
elements[i]->GetVertices();
12257 for (
int j = 0; j < nv; j++)
12267 voff[i] = vcount[i-1] + voff[i-1];
12273 vown[i] =
new int[vcount[i]];
12285 int nv =
elements[i]->GetNVertices();
12286 const int *ind =
elements[i]->GetVertices();
12287 for (
int j = 0; j < nv; j++)
12290 vown[ind[j]][vcount[ind[j]]] = i;
12296 vcount[i] = voff[i+1] - voff[i];
12300 for (
int i = 0; i < edge_el.
Size(); i++)
12302 const int *el = edge_el.
GetRow(i);
12305 int k = partitioning[el[0]];
12306 int l = partitioning[el[1]];
12307 if (interior_faces || k != l)
12319 os <<
"areamesh2\n\n" << nbe <<
'\n';
12321 for (
int i = 0; i < edge_el.
Size(); i++)
12323 const int *el = edge_el.
GetRow(i);
12326 int k = partitioning[el[0]];
12327 int l = partitioning[el[1]];
12328 if (interior_faces || k != l)
12333 for (
int j = 0; j < 2; j++)
12334 for (
int s = 0;
s < vcount[ev[j]];
s++)
12335 if (vown[ev[j]][
s] == el[0])
12337 os <<
' ' << voff[ev[j]]+
s+1;
12341 for (
int j = 1; j >= 0; j--)
12342 for (
int s = 0;
s < vcount[ev[j]];
s++)
12343 if (vown[ev[j]][
s] == el[1])
12345 os <<
' ' << voff[ev[j]]+
s+1;
12352 int k = partitioning[el[0]];
12356 for (
int j = 0; j < 2; j++)
12357 for (
int s = 0;
s < vcount[ev[j]];
s++)
12358 if (vown[ev[j]][
s] == el[0])
12360 os <<
' ' << voff[ev[j]]+
s+1;
12370 int nv =
elements[i]->GetNVertices();
12371 const int *ind =
elements[i]->GetVertices();
12372 os << partitioning[i]+1 <<
' ';
12374 for (
int j = 0; j < nv; j++)
12376 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
12377 vown[ind[j]][vcount[ind[j]]] = i;
12384 vcount[i] = voff[i+1] - voff[i];
12390 for (
int k = 0; k < vcount[i]; k++)
12392 for (
int j = 0; j <
Dim; j++)
12402 os <<
"NETGEN_Neutral_Format\n";
12406 for (
int k = 0; k < vcount[i]; k++)
12408 for (
int j = 0; j <
Dim; j++)
12419 int nv =
elements[i]->GetNVertices();
12420 const int *ind =
elements[i]->GetVertices();
12421 os << partitioning[i]+1;
12422 for (
int j = 0; j < nv; j++)
12424 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
12425 vown[ind[j]][vcount[ind[j]]] = i;
12432 vcount[i] = voff[i+1] - voff[i];
12442 int k = partitioning[
faces_info[i].Elem1No];
12443 l = partitioning[l];
12444 if (interior_faces || k != l)
12461 int k = partitioning[
faces_info[i].Elem1No];
12462 l = partitioning[l];
12463 if (interior_faces || k != l)
12465 int nv =
faces[i]->GetNVertices();
12466 const int *ind =
faces[i]->GetVertices();
12468 for (
int j = 0; j < nv; j++)
12469 for (
int s = 0;
s < vcount[ind[j]];
s++)
12472 os <<
' ' << voff[ind[j]]+
s+1;
12476 for (
int j = nv-1; j >= 0; j--)
12477 for (
int s = 0;
s < vcount[ind[j]];
s++)
12480 os <<
' ' << voff[ind[j]]+
s+1;
12487 int k = partitioning[
faces_info[i].Elem1No];
12488 int nv =
faces[i]->GetNVertices();
12489 const int *ind =
faces[i]->GetVertices();
12491 for (
int j = 0; j < nv; j++)
12492 for (
int s = 0;
s < vcount[ind[j]];
s++)
12495 os <<
' ' << voff[ind[j]]+
s+1;
12511 int k = partitioning[
faces_info[i].Elem1No];
12512 l = partitioning[l];
12513 if (interior_faces || k != l)
12526 <<
" 0 0 0 0 0 0 0\n"
12527 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
12528 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
12529 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
12530 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
12533 for (
int k = 0; k < vcount[i]; k++)
12534 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
12539 int nv =
elements[i]->GetNVertices();
12540 const int *ind =
elements[i]->GetVertices();
12541 os << i+1 <<
' ' << partitioning[i]+1;
12542 for (
int j = 0; j < nv; j++)
12544 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
12545 vown[ind[j]][vcount[ind[j]]] = i;
12552 vcount[i] = voff[i+1] - voff[i];
12561 int k = partitioning[
faces_info[i].Elem1No];
12562 l = partitioning[l];
12563 if (interior_faces || k != l)
12565 int nv =
faces[i]->GetNVertices();
12566 const int *ind =
faces[i]->GetVertices();
12568 for (
int j = 0; j < nv; j++)
12569 for (
int s = 0;
s < vcount[ind[j]];
s++)
12572 os <<
' ' << voff[ind[j]]+
s+1;
12574 os <<
" 1.0 1.0 1.0 1.0\n";
12576 for (
int j = nv-1; j >= 0; j--)
12577 for (
int s = 0;
s < vcount[ind[j]];
s++)
12580 os <<
' ' << voff[ind[j]]+
s+1;
12582 os <<
" 1.0 1.0 1.0 1.0\n";
12587 int k = partitioning[
faces_info[i].Elem1No];
12588 int nv =
faces[i]->GetNVertices();
12589 const int *ind =
faces[i]->GetVertices();
12591 for (
int j = 0; j < nv; j++)
12592 for (
int s = 0;
s < vcount[ind[j]];
s++)
12595 os <<
' ' << voff[ind[j]]+
s+1;
12597 os <<
" 1.0 1.0 1.0 1.0\n";
12621 " NURBS mesh is not supported!");
12625 os <<
"MFEM mesh v1.0\n";
12629 "\n#\n# MFEM Geometry Types (see fem/geom.hpp):\n#\n"
12634 "# TETRAHEDRON = 4\n"
12639 os <<
"\ndimension\n" <<
Dim
12647 const int *
const i_AF_f = Aface_face.
GetI();
12648 const int *
const j_AF_f = Aface_face.
GetJ();
12650 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
12651 for (
const int * iface = j_AF_f + i_AF_f[iAF];
12652 iface < j_AF_f + i_AF_f[iAF+1];
12655 os << iAF+1 <<
' ';
12688 int *nbea =
new int[na];
12695 for (i = 0; i < na; i++)
12707 for (k = 0; k < vert.
Size(); k++)
12719 for (k = 0; k < vert.
Size(); k++)
12720 if (vn[vert[k]] == 1)
12725 cg[bea*
spaceDim+j] += pointmat(j,k);
12736 for (k = 0; k < vert.
Size(); k++)
12741 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
12758 int *nbea =
new int[na];
12765 for (i = 0; i < na; i++)
12777 for (k = 0; k < vert.
Size(); k++)
12789 for (k = 0; k < vert.
Size(); k++)
12790 if (vn[vert[k]] == 1)
12795 cg[bea*
spaceDim+j] += pointmat(j,k);
12806 for (k = 0; k < vert.
Size(); k++)
12811 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
12827 for (
int i = 0; i <
vertices.Size(); i++)
12829 for (
int j = 0; j <
spaceDim; j++)
12850 "incompatible vector dimensions");
12858 for (
int d = 0; d <
spaceDim; d++)
12878 for (
int i = 0; i <
GetNE(); i++)
12883 for (
int j = 0; j < nv; j++)
12888 for (
int i = 0; i <
GetNBE(); i++)
12893 for (
int j = 0; j < nv; j++)
12899 for (
int i = 0; i < v2v.
Size(); i++)
12904 v2v[i] = num_vert++;
12908 if (num_vert == v2v.
Size()) {
return; }
12910 Vector nodes_by_element;
12915 for (
int i = 0; i <
GetNE(); i++)
12922 for (
int i = 0; i <
GetNE(); i++)
12931 for (
int i = 0; i <
GetNE(); i++)
12936 for (
int j = 0; j < nv; j++)
12941 for (
int i = 0; i <
GetNBE(); i++)
12946 for (
int j = 0; j < nv; j++)
12970 for (
int i = 0; i <
GetNE(); i++)
12983 int num_bdr_elem = 0;
12984 int new_bel_to_edge_nnz = 0;
12985 for (
int i = 0; i <
GetNBE(); i++)
13001 if (num_bdr_elem ==
GetNBE()) {
return; }
13005 Table *new_bel_to_edge = NULL;
13007 new_be_to_face.
Reserve(num_bdr_elem);
13010 new_bel_to_edge =
new Table;
13011 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
13013 for (
int i = 0; i <
GetNBE(); i++)
13018 int row = new_be_to_face.
Size();
13024 int *new_e = new_bel_to_edge->
GetRow(row);
13025 for (
int j = 0; j < ne; j++)
13029 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
13046 for (
int i = 0; i < attribs.
Size(); i++)
13058#ifdef MFEM_USE_MEMALLOC
13085 const int npts = point_mat.
Width();
13086 if (!npts) {
return 0; }
13087 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
13091 if (!
GetNE()) {
return 0; }
13100 min_dist = std::numeric_limits<real_t>::max();
13104 for (
int i = 0; i <
GetNE(); i++)
13108 for (
int k = 0; k < npts; k++)
13111 if (dist < min_dist(k))
13113 min_dist(k) = dist;
13122 for (
int k = 0; k < npts; k++)
13126 int res = inv_tr->
Transform(pt, ips[k]);
13129 elem_ids[k] = e_idx[k];
13133 if (pts_found != npts)
13137 for (
int k = 0; k < npts; k++)
13139 if (elem_ids[k] != -1) {
continue; }
13143 for (
int v = 0; v < elvertices.
Size(); v++)
13145 int vv = elvertices[v];
13147 const int* els = vtoel->
GetRow(vv);
13148 for (
int e = 0; e < ne; e++)
13150 if (els[e] == e_idx[k]) {
continue; }
13152 int res = inv_tr->
Transform(pt, ips[k]);
13155 elem_ids[k] = els[e];
13167 for (
int e = 0; e < neigh.
Size(); e++)
13173 int res = inv_tr->
Transform(pt, ips[k]);
13186 if (inv_trans == NULL) {
delete inv_tr; }
13188 if (warn && pts_found != npts)
13190 MFEM_WARNING((npts-pts_found) <<
" points were not found");
13205 MFEM_VERIFY(
Dim == 2 ||
Dim == 3,
"Only 2D/3D meshes supported right now.");
13206 MFEM_VERIFY(
Dim ==
spaceDim,
"Surface meshes not currently supported.");
13223 skew(0) = std::atan2(J.
Det(), col1 * col2);
13226 ori(0) = std::atan2(J(1,0), J(0,0));
13233 Vector col1, col2, col3;
13244 col1unit *= 1.0/len1;
13245 col2unit *= 1.0/len2;
13246 col3unit *= 1.0/len3;
13252 aspr(0) = len1/std::sqrt(len2*len3),
13253 aspr(1) = len2/std::sqrt(len1*len3);
13256 aspr(2) = std::sqrt(len1/(len2*len3)),
13257 aspr(3) = std::sqrt(len2/(len1*len3));
13260 Vector crosscol12, crosscol13;
13261 col1.
cross3D(col2, crosscol12);
13262 col1.
cross3D(col3, crosscol13);
13263 skew(0) = std::acos(col1unit*col2unit);
13264 skew(1) = std::acos(col1unit*col3unit);
13265 skew(2) = std::atan(len1*volume/(crosscol12*crosscol13));
13271 for (
int d=0; d<
Dim; d++) { rot(d, 0) = col1unit(d); }
13275 rot1 *= col1unit*col2unit;
13277 col1unit.
cross3D(col2unit, rot1);
13279 for (
int d=0; d <
Dim; d++) { rot(d, 1) = rot2(d); }
13282 for (
int d=0; d <
Dim; d++) { rot(d, 2) = rot1(d); }
13283 real_t delta = sqrt(pow(rot(2,1)-rot(1,2), 2.0) +
13284 pow(rot(0,2)-rot(2,0), 2.0) +
13285 pow(rot(1,0)-rot(0,1), 2.0));
13290 for (
int d = 0; d <
Dim; d++) { Iden(d, d) = 1.0; };
13296 MFEM_ABORT(
"Invalid rotation matrix. Contact TMOP Developers.");
13301 ori(0) = (1./
delta)*(rot(2,1)-rot(1,2));
13302 ori(1) = (1./
delta)*(rot(0,2)-rot(2,0));
13303 ori(2) = (1./
delta)*(rot(1,0)-rot(0,1));
13304 ori(3) = std::acos(0.5*(rot.
Trace()-1.0));
13313 entity_to_vertex(entity_to_vertex_)
13315 int geom_offset = 0;
13329 while (geom_offsets[geom+1] <= bytype_entity_id) { geom++; }
13333 const int geom_elem_id = bytype_entity_id - geom_offsets[geom];
13335 return { geom, nv, v };
13340 os <<
"MFEM mesh v1.2\n";
13344 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
13349 "# TETRAHEDRON = 4\n"
13356 os <<
"\ndimension\n" <<
dim;
13362 "invalid MeshPart state");
13365 "invalid MeshPart state");
13366 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
13368 const int bytype_elem_id = have_element_map ?
13373 for (
int i = 0; i < ent.
num_verts; i++)
13375 os <<
' ' << ent.
verts[i];
13385 "invalid MeshPart state");
13388 "invalid MeshPart state");
13391 const int bytype_bdr_id = have_boundary_map ?
13396 for (
int i = 0; i < ent.
num_verts; i++)
13398 os <<
' ' << ent.
verts[i];
13408 os << sdim <<
'\n';
13412 for (
int d = 1; d < sdim; d++)
13425 os <<
"\nmfem_serial_mesh_end\n";
13429 os <<
"\ncommunication_groups\n";
13430 os <<
"number_of_groups " << num_groups <<
"\n\n";
13432 os <<
"# number of entities in each group, followed by ranks in group\n";
13433 for (
int group_id = 0; group_id < num_groups; ++group_id)
13438 for (
int group_member_index = 0; group_member_index < group_size;
13439 ++group_member_index)
13441 os <<
' ' << group_ptr[group_member_index];
13452 MFEM_VERIFY(g2v.
RowSize(0) == 0,
"internal erroor");
13456 MFEM_VERIFY(g2ev.
RowSize(0) == 0,
"internal erroor");
13461 MFEM_VERIFY(g2tv.
RowSize(0) == 0,
"internal erroor");
13462 MFEM_VERIFY(g2qv.
RowSize(0) == 0,
"internal erroor");
13463 const int total_shared_faces =
13465 os <<
"total_shared_faces " << total_shared_faces <<
'\n';
13467 os <<
"\n# group 0 has no shared entities\n";
13468 for (
int gr = 1; gr < num_groups; gr++)
13471 const int nv = g2v.
RowSize(gr);
13472 const int *sv = g2v.
GetRow(gr);
13473 os <<
"\n# group " << gr <<
"\nshared_vertices " << nv <<
'\n';
13474 for (
int i = 0; i < nv; i++)
13476 os << sv[i] <<
'\n';
13481 const int ne = g2ev.
RowSize(gr)/2;
13482 const int *se = g2ev.
GetRow(gr);
13483 os <<
"\nshared_edges " << ne <<
'\n';
13484 for (
int i = 0; i < ne; i++)
13486 const int *v = se + 2*i;
13487 os << v[0] <<
' ' << v[1] <<
'\n';
13492 const int nt = g2tv.
RowSize(gr)/3;
13493 const int *st = g2tv.
GetRow(gr);
13494 const int nq = g2qv.
RowSize(gr)/4;
13496 os <<
"\nshared_faces " << nt+nq <<
'\n';
13497 for (
int i = 0; i < nt; i++)
13500 const int *v = st + 3*i;
13501 for (
int j = 0; j < 3; j++) { os <<
' ' << v[j]; }
13504 for (
int i = 0; i < nq; i++)
13507 const int *v =
sq + 4*i;
13508 for (
int j = 0; j < 4; j++) { os <<
' ' << v[j]; }
13515 os <<
"\nmfem_mesh_end" << endl;
13532 "invalid MeshPart state");
13535 "invalid MeshPart state");
13537 for (
int nat_elem_id = 0; nat_elem_id <
num_elements; nat_elem_id++)
13539 const int bytype_elem_id = have_element_map ?
13550 static_cast<Tetrahedron*
>(el)->SetRefinementFlag(ref_flag);
13552 mesh->AddElement(el);
13560 "invalid MeshPart state");
13563 "invalid MeshPart state");
13566 const int bytype_bdr_id = have_boundary_map ?
13572 mesh->AddBdrElement(bdr);
13579 MFEM_ASSERT(!
nodes,
"invalid MeshPart state");
13580 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
13588 for (
int vert_id = 0; vert_id <
num_vertices; vert_id++)
13590 mesh->AddVertex(0., 0., 0.);
13595 mesh->FinalizeTopology(
false);
13603 int *partitioning_,
13634 for (
int i = 0; i < boundary_to_part.
Size(); i++)
13636 int face, o, el1, el2;
13639 boundary_to_part[i] =
13645 for (
int i = 0; i < boundary_to_part.
Size(); i++)
13654 for (
int i = 0; i < boundary_to_part.
Size(); i++)
13669 delete vert_element;
13676 MFEM_VERIFY(0 <= part_id && part_id < num_parts,
13677 "invalid part_id = " << part_id
13678 <<
", num_parts = " << num_parts);
13711 mesh_part.
nodes.reset(
nullptr);
13713 mesh_part.
mesh.reset(
nullptr);
13721 int geom_marker = 0, num_geom = 0;
13722 for (
int i = 0; i < num_elems; i++)
13728 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
13730 "overflow in 'entity_to_vertex[geom]', geom: "
13755 if ((geom_marker & (1 << geom)) == 0)
13757 geom_marker |= (1 << geom);
13772 offsets[g] = offset;
13776 for (
int i = 0; i < num_elems; i++)
13787 geom_marker = 0; num_geom = 0;
13788 for (
int i = 0; i < num_bdr_elems; i++)
13794 MFEM_VERIFY(numeric_limits<int>::max() - nv >=
13796 "overflow in 'entity_to_vertex[geom]', geom: "
13800 if ((geom_marker & (1 << geom)) == 0)
13802 geom_marker |= (1 << geom);
13813 offsets[g] = offset;
13817 for (
int i = 0; i < num_bdr_elems; i++)
13828 std::unordered_set<int> vertex_set;
13829 for (
int i = 0; i < num_elems; i++)
13835 vertex_set.insert(v, v + nv);
13837 vertex_loc_to_glob.
SetSize(vertex_set.size());
13838 std::copy(vertex_set.begin(), vertex_set.end(),
13839 vertex_loc_to_glob.
begin());
13841 vertex_loc_to_glob.
Sort();
13851 for (
int i = 0; i < vert_array.
Size(); i++)
13853 const int glob_id = vert_array[i];
13854 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
13855 MFEM_ASSERT(loc_id >= 0,
"internal error: global vertex id not found");
13856 vert_array[i] = loc_id;
13863 MFEM_VERIFY(numeric_limits<int>::max()/sdim >= vertex_loc_to_glob.
Size(),
13864 "overflow in 'vertex_coordinates', num_vertices = "
13865 << vertex_loc_to_glob.
Size() <<
", sdim = " << sdim);
13867 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
13870 for (
int d = 0; d < sdim; d++)
13887 mesh_part.
mesh->NewNodes(*mesh_part.
nodes,
false);
13909 std::unordered_set<int> face_set;
13912 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
13914 const int glob_elem_id = elem_list[loc_elem_id];
13915 const int nfaces = elem_to_face.
RowSize(glob_elem_id);
13916 const int *faces = elem_to_face.
GetRow(glob_elem_id);
13917 face_set.insert(faces, faces + nfaces);
13921 for (
int glob_face_id : face_set)
13925 if (el[1] < 0) {
continue; }
13928 MFEM_ASSERT(el[0] == part_id || el[1] == part_id,
"internal error");
13929 if (el[0] != part_id || el[1] != part_id)
13932 const int group_id = groups.
Insert(group);
13936 shared_faces.
Sort();
13947 std::unordered_set<int> edge_set;
13950 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
13952 const int glob_elem_id = elem_list[loc_elem_id];
13953 const int nedges = elem_to_edge.
RowSize(glob_elem_id);
13954 const int *edges = elem_to_edge.
GetRow(glob_elem_id);
13955 edge_set.insert(edges, edges + nedges);
13959 for (
int glob_edge_id : edge_set)
13965 for (
int j = 0; j < nelem; j++)
13971 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
13972 if (group.
Size() > 1)
13974 const int group_id = groups.
Insert(group);
13978 shared_edges.
Sort();
13989 for (
int i = 0; i < vertex_loc_to_glob.
Size(); i++)
13992 const int glob_vertex_id = vertex_loc_to_glob[i];
13997 for (
int j = 0; j < nelem; j++)
14003 MFEM_ASSERT(gr.
FindSorted(part_id) >= 0,
"internal error");
14004 if (group.
Size() > 1)
14006 const int group_id = groups.
Insert(group);
14013 const int num_groups = groups.
Size();
14019 Table &group__shared_vertex_to_vertex =
14021 group__shared_vertex_to_vertex.
MakeI(num_groups);
14022 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
14024 const int group_id = shared_verts[sv].two;
14027 group__shared_vertex_to_vertex.
MakeJ();
14028 for (
int sv = 0; sv < shared_verts.
Size(); sv++)
14030 const int glob_vertex_id = shared_verts[sv].one;
14031 const int group_id = shared_verts[sv].two;
14032 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(glob_vertex_id);
14033 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
14034 group__shared_vertex_to_vertex.
AddConnection(group_id, loc_vertex_id);
14036 group__shared_vertex_to_vertex.
ShiftUpI();
14041 Table &group__shared_edge_to_vertex =
14043 group__shared_edge_to_vertex.
MakeI(num_groups);
14044 for (
int se = 0; se < shared_edges.
Size(); se++)
14046 const int group_id = shared_edges[se].two;
14049 group__shared_edge_to_vertex.
MakeJ();
14051 for (
int se = 0; se < shared_edges.
Size(); se++)
14053 const int glob_edge_id = shared_edges[se].one;
14054 const int group_id = shared_edges[se].two;
14055 const int *v = edge_to_vertex.
GetRow(glob_edge_id);
14056 for (
int i = 0; i < 2; i++)
14058 const int loc_vertex_id = vertex_loc_to_glob.
FindSorted(v[i]);
14059 MFEM_ASSERT(loc_vertex_id >= 0,
"internal error");
14060 group__shared_edge_to_vertex.
AddConnection(group_id, loc_vertex_id);
14063 group__shared_edge_to_vertex.
ShiftUpI();
14070 Table &group__shared_tria_to_vertex =
14072 Table &group__shared_quad_to_vertex =
14075 group__shared_tria_to_vertex.
MakeI(num_groups);
14076 group__shared_quad_to_vertex.
MakeI(num_groups);
14077 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
14079 const int glob_face_id = shared_faces[sf].one;
14080 const int group_id = shared_faces[sf].two;
14085 group__shared_tria_to_vertex.
MakeJ();
14086 group__shared_quad_to_vertex.
MakeJ();
14087 for (
int sf = 0; sf < shared_faces.
Size(); sf++)
14089 const int glob_face_id = shared_faces[sf].one;
14090 const int group_id = shared_faces[sf].two;
14126 for (
int i = 0; i < vertex_ids.
Size(); i++)
14128 const int glob_id = vertex_ids[i];
14129 const int loc_id = vertex_loc_to_glob.
FindSorted(glob_id);
14130 MFEM_ASSERT(loc_id >= 0,
"internal error");
14131 vertex_ids[i] = loc_id;
14134 AddConnections(group_id, vertex_ids, vertex_ids.
Size());
14136 group__shared_tria_to_vertex.
ShiftUpI();
14137 group__shared_quad_to_vertex.
ShiftUpI();
14141std::unique_ptr<FiniteElementSpace>
14149 return std::unique_ptr<FiniteElementSpace>(
14151 global_fespace.
FEColl(),
14156std::unique_ptr<GridFunction>
14161 std::unique_ptr<GridFunction> local_gf(
new GridFunction(&local_fespace));
14169 for (
int loc_elem_id = 0; loc_elem_id < num_elems; loc_elem_id++)
14171 const int glob_elem_id = elem_list[loc_elem_id];
14174 if (glob_dt) { glob_dt->InvTransformPrimal(loc_vals); }
14177 local_gf->SetSubVector(lvdofs, loc_vals);
14191 "mixed meshes are not supported!");
14192 MFEM_ASSERT(
mesh->
GetNodes(),
"meshes without nodes are not supported!");
14205 Compute(nodes, d_mt);
14208void GeometricFactors::Compute(
const GridFunction &nodes,
14215 const int vdim = fespace->
GetVDim();
14216 const int NE = fespace->
GetNE();
14217 const int ND = fe->
GetDof();
14220 unsigned eval_flags = 0;
14222 Device::GetDeviceMemoryType();
14245 qi->DisableTensorProducts(!use_tensor_products);
14253 Vector Enodes(vdim*ND*NE, my_d_mt);
14254 elem_restr->Mult(nodes, Enodes);
14255 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
14259 qi->Mult(nodes, eval_flags,
X,
J,
detJ);
14275 const int vdim = fespace->
GetVDim();
14289 face_restr->
Mult(*nodes, Fnodes);
14291 unsigned eval_flags = 0;
14300 J.
SetSize(vdim*vdim*NQ*NF, my_d_mt);
14338 V(1) = s * ((ip.
y + layer) / n);
14343 V(2) = s * ((ip.
z + layer) / n);
14352 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
14356 int nvy = (closed) ? (ny) : (ny + 1);
14357 int nvt = mesh->
GetNV() * nvy;
14366 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
14371 for (
int i = 0; i < mesh->
GetNV(); i++)
14374 for (
int j = 0; j < nvy; j++)
14376 vc[1] = sy * (
real_t(j) / ny);
14382 for (
int i = 0; i < mesh->
GetNE(); i++)
14387 for (
int j = 0; j < ny; j++)
14390 qv[0] = vert[0] * nvy + j;
14391 qv[1] = vert[1] * nvy + j;
14392 qv[2] = vert[1] * nvy + (j + 1) % nvy;
14393 qv[3] = vert[0] * nvy + (j + 1) % nvy;
14399 for (
int i = 0; i < mesh->
GetNBE(); i++)
14404 for (
int j = 0; j < ny; j++)
14407 sv[0] = vert[0] * nvy + j;
14408 sv[1] = vert[0] * nvy + (j + 1) % nvy;
14424 for (
int i = 0; i < mesh->
GetNE(); i++)
14430 sv[0] = vert[0] * nvy;
14431 sv[1] = vert[1] * nvy;
14435 sv[0] = vert[1] * nvy + ny;
14436 sv[1] = vert[0] * nvy + ny;
14452 string cname = name;
14453 if (cname ==
"Linear")
14457 else if (cname ==
"Quadratic")
14461 else if (cname ==
"Cubic")
14465 else if (!strncmp(name,
"H1_", 3))
14469 else if (!strncmp(name,
"L2_T", 4))
14473 else if (!strncmp(name,
"L2_", 3))
14480 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
14492 for (
int i = 0; i < mesh->
GetNE(); i++)
14495 for (
int j = ny-1; j >= 0; j--)
14512 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
14517 int nvt = mesh->
GetNV() * nvz;
14522 bool wdgMesh =
false;
14523 bool hexMesh =
false;
14527 for (
int i = 0; i < mesh->
GetNV(); i++)
14531 for (
int j = 0; j < nvz; j++)
14533 vc[2] = sz * (
real_t(j) / nz);
14539 for (
int i = 0; i < mesh->
GetNE(); i++)
14549 for (
int j = 0; j < nz; j++)
14552 pv[0] = vert[0] * nvz + j;
14553 pv[1] = vert[1] * nvz + j;
14554 pv[2] = vert[2] * nvz + j;
14555 pv[3] = vert[0] * nvz + (j + 1) % nvz;
14556 pv[4] = vert[1] * nvz + (j + 1) % nvz;
14557 pv[5] = vert[2] * nvz + (j + 1) % nvz;
14564 for (
int j = 0; j < nz; j++)
14567 hv[0] = vert[0] * nvz + j;
14568 hv[1] = vert[1] * nvz + j;
14569 hv[2] = vert[2] * nvz + j;
14570 hv[3] = vert[3] * nvz + j;
14571 hv[4] = vert[0] * nvz + (j + 1) % nvz;
14572 hv[5] = vert[1] * nvz + (j + 1) % nvz;
14573 hv[6] = vert[2] * nvz + (j + 1) % nvz;
14574 hv[7] = vert[3] * nvz + (j + 1) % nvz;
14576 mesh3d->
AddHex(hv, attr);
14580 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
14581 << geom <<
"\'" << endl;
14587 for (
int i = 0; i < mesh->
GetNBE(); i++)
14592 for (
int j = 0; j < nz; j++)
14595 qv[0] = vert[0] * nvz + j;
14596 qv[1] = vert[1] * nvz + j;
14597 qv[2] = vert[1] * nvz + (j + 1) % nvz;
14598 qv[3] = vert[0] * nvz + (j + 1) % nvz;
14607 for (
int i = 0; i < mesh->
GetNE(); i++)
14618 tv[0] = vert[0] * nvz;
14619 tv[1] = vert[2] * nvz;
14620 tv[2] = vert[1] * nvz;
14624 tv[0] = vert[0] * nvz + nz;
14625 tv[1] = vert[1] * nvz + nz;
14626 tv[2] = vert[2] * nvz + nz;
14634 qv[0] = vert[0] * nvz;
14635 qv[1] = vert[3] * nvz;
14636 qv[2] = vert[2] * nvz;
14637 qv[3] = vert[1] * nvz;
14641 qv[0] = vert[0] * nvz + nz;
14642 qv[1] = vert[1] * nvz + nz;
14643 qv[2] = vert[2] * nvz + nz;
14644 qv[3] = vert[3] * nvz + nz;
14650 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
14651 << geom <<
"\'" << endl;
14657 if ( hexMesh && wdgMesh )
14661 else if ( hexMesh )
14665 else if ( wdgMesh )
14678 string cname = name;
14679 if (cname ==
"Linear")
14683 else if (cname ==
"Quadratic")
14687 else if (cname ==
"Cubic")
14691 else if (!strncmp(name,
"H1_", 3))
14695 else if (!strncmp(name,
"L2_T", 4))
14699 else if (!strncmp(name,
"L2_", 3))
14706 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
14718 for (
int i = 0; i < mesh->
GetNE(); i++)
14721 for (
int j = nz-1; j >= 0; j--)
14742 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
14743 <<
" 0 0 " << i <<
" -1 0\n";
14750 real_t mid[3] = {0, 0, 0};
14751 for (
int j = 0; j < 2; j++)
14753 for (
int k = 0; k <
spaceDim; k++)
14759 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
14760 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
Node::Index insert_node(Float length=1)
Arc::Index insert_arc(Node::Index i, Node::Index j, Float w=1, Float b=1)
uint rank(Node::Index i) const
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int FindSorted(const T &el) const
Do bisection search for 'el' in a sorted array; return -1 if not found.
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
void DeleteAll()
Delete the whole array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
T * GetData()
Returns the data.
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
T * end()
STL-like end. Returns pointer after the last element of the array.
T * begin()
STL-like begin. Returns pointer to the first element of the array.
T & Last()
Return the last element in the array.
bool SetsExist() const
Return true if any named sets are currently defined.
void Print(std::ostream &out=mfem::out, int width=-1) const
Print the contents of the container to an output stream.
void Copy(AttributeSets ©) const
Create a copy of the internal data to the provided copy.
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Piecewise-(bi)cubic continuous finite elements.
int NumberOfEntries() const
Data type dense matrix using column-major storage.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
const real_t * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
real_t * Data() const
Returns the matrix data array.
real_t * GetData() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
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 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.
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 NURBSExtension * GetNURBSext() const
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.
NURBSExtension * StealNURBSext()
int GetNE() const
Returns number of elements in the mesh.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
void SetRelaxedHpConformity(bool relaxed=true)
int GetNFDofs() const
Number of all scalar face-interior dofs.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
const FiniteElementCollection * FEColl() const
Mesh * GetMesh() const
Returns the mesh.
int GetMaxElementOrder() const
Return the maximum polynomial order.
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 vector dimension.
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.
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 and fes.
FiniteElementSpace * FESpace()
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Arbitrary order H1-conforming (continuous) finite elements.
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
int GetId(int p1, int p2)
Get id of item whose parents are p1, p2... Create it if it doesn't exist.
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn't exist.
Data type hexahedron element.
void Recreate(const int n, const int *p)
Create an integer set from C-array 'p' of 'n' integers. Overwrites any existing set data.
int Size()
Return the size of the set.
Class for integration point with weight.
void Get(real_t *p, const int dim) const
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Arbitrary order "L2-conforming" discontinuous finite elements.
Piecewise-(bi/tri)linear continuous finite elements.
int Size()
Return the number of integer sets in the list.
int Insert(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)
Write the list of sets into table 't'.
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_, int *partitioning_=NULL, int part_method=1)
Construct a MeshPartitioner.
std::unique_ptr< GridFunction > ExtractGridFunction(const MeshPart &mesh_part, const GridFunction &global_gf, FiniteElementSpace &local_fespace) const
Construct a local version of the given GridFunction, global_gf, corresponding to the given mesh_part....
void ExtractPart(int part_id, MeshPart &mesh_part) const
Construct a MeshPart corresponding to the given part_id.
List of mesh geometries stored as Array<Geometry::Type>.
int CheckElementOrientation(bool fix_it=true)
Check (and optionally attempt to fix) the orientation of the elements.
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
void GetEdgeOrdering(const DSTable &v_to_v, Array< int > &order)
void GetLocalFaceTransformation(int face_type, int elem_type, IsoparametricTransformation &Transf, int info) const
A helper method that constructs a transformation from the reference space of a face to the reference ...
void NURBSCoarsening(int cf=2, real_t tol=1.0e-12)
void SetVerticesFromNodes(const GridFunction *nodes)
Helper to set vertex coordinates given a high-order curvature function.
void ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
int GetPatchBdrAttribute(int i) const
Return the attribute of patch boundary element i, for a NURBS mesh.
int GetElementToEdgeTable(Table &)
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
void SetVertices(const Vector &vert_coord)
Element * NewElement(int geom)
Operation GetLastOperation() const
Return type of last modification of the mesh.
IsoparametricTransformation Transformation2
int GetNEdges() const
Return the number of edges.
void GetBdrElementFace(int i, int *f, int *o) const
void InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
Begin construction of a mesh.
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
IsoparametricTransformation EdgeTransformation
static FiniteElement * GetTransformationFEforElementType(Element::Type)
Return FiniteElement for reference element of the specified type.
int AddBdrQuad(int v1, int v2, int v3, int v4, int attr=1)
int * CartesianPartitioning(int nxyz[])
Array< int > FindFaceNeighbors(const int elem) const
Returns the sorted, unique indices of elements sharing a face with element elem, including elem.
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
Element::Type GetElementType(int i) const
Returns the type of element i.
void GetLocalSegToQuadTransformation(IsoparametricTransformation &loc, int i) const
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
void BdrBisection(int i, const HashTable< Hashed2 > &)
Bisect a boundary triangle: boundary element with index i is bisected.
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
const Table & ElementToEdgeTable() const
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void ReadNetgen3DMesh(std::istream &input)
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Geometry::Type GetElementGeometry(int i) const
Geometry::Type GetBdrElementGeometry(int i) const
int AddTri(const int *vi, int attr=1)
Adds a triangle to the mesh given by 3 vertices vi.
static Mesh MakeCartesian1D(int n, real_t sx=1.0)
Creates 1D mesh, divided into n equal intervals.
int GetAttribute(int i) const
Return the attribute of element i.
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally,...
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
void UniformRefinement3D_base(Array< int > *f2qf=NULL, DSTable *v_to_v_p=NULL, bool update_nodes=true)
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
long nodes_sequence
Counter for geometric factor invalidation.
virtual void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
IsoparametricTransformation FaceTransformation
Array< NCFaceInfo > nc_faces_info
void MakeRefined_(Mesh &orig_mesh, const Array< int > &ref_factors, int ref_type)
Internal function used in Mesh::MakeRefined.
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
Adds a wedge to the mesh given by 6 vertices v1 through v6.
Array< int > GetFaceToBdrElMap() const
static Mesh MakeCartesian2DWith4TrisPerQuad(int nx, int ny, real_t sx=1.0, real_t sy=1.0)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*4 triangles.
void ReadInlineMesh(std::istream &input, bool generate_edges=false)
void SetPatchAttribute(int i, int attr)
Set the attribute of patch i, for a NURBS mesh.
void FinalizeTetMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a tetrahedral Mesh.
real_t GetLength(int i, int j) const
Return the length of the segment from node i to node j.
const FiniteElementSpace * GetNodalFESpace() const
void AddBdrQuadAsTriangles(const int *vi, int attr=1)
int AddPyramid(int v1, int v2, int v3, int v4, int v5, int attr=1)
Adds a pyramid to the mesh given by 5 vertices v1 through v5.
void Loader(std::istream &input, int generate_edges=0, std::string parse_tag="")
const Table & ElementToElementTable()
void ScaleElements(real_t sf)
void GenerateNCFaceInfo()
void ReadLineMesh(std::istream &input)
void ApplyLocalSlaveTransformation(FaceElementTransformations &FT, const FaceInfo &fi, bool is_ghost) const
real_t AggregateError(const Array< real_t > &elem_error, const int *fine, int nfine, int op)
Derefinement helper.
void CheckPartitioning(int *partitioning_)
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
void GetLocalPtToSegTransformation(IsoparametricTransformation &, int i) const
Used in GetFaceElementTransformations (...)
bool Nonconforming() const
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
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.
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
void GetVertices(Vector &vert_coord) const
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
virtual int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type, does not count master nonconforming face...
void Make1D(int n, real_t sx=1.0)
void RefineNURBSFromFile(std::string ref_file)
const Element * GetElement(int i) const
Return pointer to the i'th element object.
int AddBdrPoint(int v, int attr=1)
void FinalizeWedgeMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a wedge Mesh.
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
void PrintTopo(std::ostream &os, const Array< int > &e_to_k, const int version, const std::string &comment="") const
Write the beginning of a NURBS mesh to os, specifying the NURBS patch topology. Optional file comment...
static Mesh MakeSimplicial(const Mesh &orig_mesh)
void SetPatchBdrAttribute(int i, int attr)
Set the attribute of patch boundary element i, for a NURBS mesh.
int GetNFaces() const
Return the number of faces in a 3D mesh.
int AddVertexAtMeanCenter(const int *vi, const int nverts, int dim=3)
static int GetTetOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
real_t GetGeckoElementOrdering(Array< int > &ordering, int iterations=4, int window=4, int period=2, int seed=0, bool verbose=false, real_t time_limit=0)
static int EncodeFaceInfo(int local_face_index, int orientation)
Given local_face_index and orientation, return the corresponding encoded "face info int".
bool FaceIsTrueInterior(int FaceNo) const
const CoarseFineTransformations & GetRefinementTransforms() const
void Make2D5QuadsFromQuad(int nx, int ny, real_t sx, real_t sy)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*5 quadrilaterals.
void GetLocalQuadToWdgTransformation(IsoparametricTransformation &loc, int i) const
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.
void SetAttribute(int i, int attr)
Set the attribute of element i.
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
static Mesh MakeCartesian2DWith5QuadsPerQuad(int nx, int ny, real_t sx=1.0, real_t sy=1.0)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*5 quadrilaterals.
const FaceGeometricFactors * GetFaceGeometricFactors(const IntegrationRule &ir, const int flags, FaceType type, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors for the faces corresponding to the given integration rule.
void PrintWithPartitioning(int *partitioning, std::ostream &os, int elem_attr=0) const
Prints the mesh with boundary elements given by the boundary of the subdomains, so that the boundary ...
void Clear()
Clear the contents of the Mesh.
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
void ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo, const std::string &xml_prefix="")
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
This function is not public anymore. Use GeneralRefinement instead.
int GetNE() const
Returns number of elements.
void Make3D(int nx, int ny, int nz, Element::Type type, real_t sx, real_t sy, real_t sz, bool sfc_ordering)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
virtual void Save(const std::string &fname, int precision=16) const
int AddTet(int v1, int v2, int v3, int v4, int attr=1)
Adds a tetrahedron to the mesh given by 4 vertices v1 through v4.
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const
int Dimension() const
Dimension of the reference space used within the elements.
void RandomRefinement(real_t prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
void CheckDisplacements(const Vector &displacements, real_t &tmax)
friend class NURBSExtension
void AddTriangleFaceElement(int lf, int gf, int el, int v0, int v1, int v2)
void AddHexAs24TetsWithPoints(int *vi, std::map< std::array< int, 4 >, int > &hex_face_verts, int attr=1)
Adds 24 tetrahedrons to the mesh by splitting a hexahedron.
void GetNode(int i, real_t *coord) const
void ReorderElements(const Array< int > &ordering, bool reorder_vertices=true)
void GreenRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
void UpdateNodes()
Update the nodes of a curved mesh after the topological part of a Mesh::Operation,...
void PrintElementsWithPartitioning(int *partitioning, std::ostream &os, int interior_faces=0)
Mesh & operator=(Mesh &&mesh)
Move assignment operator.
static int InvertQuadOrientation(int ori)
Array< FaceGeometricFactors * > face_geom_factors
void GetLocalTriToPyrTransformation(IsoparametricTransformation &loc, int i) const
static Mesh MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
Create a refined (by any factor) version of orig_mesh.
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Element * ReadElementWithoutAttr(std::istream &input)
int AddBdrSegment(int v1, int v2, int attr=1)
bool DerefineByError(Array< real_t > &elem_error, real_t threshold, int nc_limit=0, int op=1)
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
int AddElement(Element *elem)
static int DecodeFaceInfoLocalIndex(int info)
Given a "face info int", return the local face index.
FaceInformation GetFaceInformation(int f) const
int GetNumFacesWithGhost() const
Return the number of faces (3D), edges (2D) or vertices (1D) including ghost faces.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
void RefineAtVertex(const Vertex &vert, real_t eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
static int InvertTriOrientation(int ori)
STable3D * GetElementToFaceTable(int ret_ftbl=0)
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
void Make3D24TetsFromHex(int nx, int ny, int nz, real_t sx, real_t sy, real_t sz)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz*24 tetrahedrons.
virtual bool NonconformingDerefinement(Array< real_t > &elem_error, real_t threshold, int nc_limit=0, int op=1)
NC version of GeneralDerefinement.
void AddVertexParents(int i, int p1, int p2)
Mark vertex i as nonconforming, with parent vertices p1 and p2.
MFEM_DEPRECATED void GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
Deprecated.
int GetPatchAttribute(int i) const
Return the attribute of patch i, for a NURBS mesh.
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
void Printer(std::ostream &os=mfem::out, std::string section_delimiter="", const std::string &comments="") const
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
ElementTransformation * GetBdrElementTransformation(int i)
Returns a pointer to the transformation defining the i-th boundary element.
IsoparametricTransformation Transformation
void Make2D4TrisFromQuad(int nx, int ny, real_t sx, real_t sy)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny*4 triangles.
void GetLocalTriToWdgTransformation(IsoparametricTransformation &loc, int i) const
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
int SpaceDimension() const
Dimension of the physical space containing the mesh.
void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf, bool spacing=false)
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0, bool sfc_ordering=true)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
void GetCharacteristics(real_t &h_min, real_t &h_max, real_t &kappa_min, real_t &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
void ReadCubit(const std::string &filename, int &curved, int &read_gf)
Load a mesh from a Genesis file.
void SetNode(int i, const real_t *coord)
void GetNodes(Vector &node_coord) const
AttributeSets attribute_sets
Named sets of element attributes.
virtual void UniformRefinement3D()
Refine a mixed 3D mesh uniformly.
static int ComposeTriOrientations(int ori_a_b, int ori_b_c)
void PrintBdrVTU(std::string fname, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0)
void AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
int AddBdrTriangle(int v1, int v2, int v3, int attr=1)
void GetGeometricParametersFromJacobian(const DenseMatrix &J, real_t &volume, Vector &aspr, Vector &skew, Vector &ori) const
Computes geometric parameters associated with a Jacobian matrix in 2D/3D. These parameters are (1) Ar...
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
static int DecodeFaceInfoOrientation(int info)
Given a "face info int", return the face orientation.
void GetHilbertElementOrdering(Array< int > &ordering)
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
void AddQuadAs5QuadsWithPoints(int *vi, int attr=1)
Adds 5 quadrilaterals to the mesh by splitting a quadrilateral given by 4 vertices vi.
void ReadVTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo)
void PrintVTU(std::ostream &os, int ref=1, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0, bool bdr_elements=false)
virtual void UniformRefinement2D()
Refine a mixed 2D mesh uniformly.
static Mesh MakePeriodic(const Mesh &orig_mesh, const std::vector< int > &v2v)
Create a periodic mesh by identifying vertices of orig_mesh.
Element::Type GetFaceElementType(int Face) const
int CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i....
void AverageVertices(const int *indexes, int n, int result)
Averages the vertices with given indexes and saves the result in vertices[result].
real_t GetElementVolume(int i)
static Mesh LoadFromFile(const std::string &filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
static const int vtk_quadratic_hex[27]
void Swap(Mesh &other, bool non_geometry)
Array< Triple< int, int, int > > tmp_vertex_parents
virtual void GenerateBoundaryElements()
static void GetElementArrayEdgeTable(const Array< Element * > &elem_array, const DSTable &v_to_v, Table &el_to_edge)
std::vector< int > CreatePeriodicVertexMapping(const std::vector< Vector > &translations, real_t tol=1e-8) const
Creates a mapping v2v from the vertex indices of the mesh such that coincident vertices under the giv...
void Bisection(int i, const DSTable &, int *, int *, int *)
Bisect a triangle: element with index i is bisected.
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
void UniformRefinement2D_base(bool update_nodes=true)
IsoparametricTransformation BdrTransformation
void GetLocalSegToTriTransformation(IsoparametricTransformation &loc, int i) const
void DegreeElevate(int rel_degree, int degree=16)
int FindCoarseElement(int i)
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
void GetElementCenter(int i, Vector ¢er)
void AddHexAsPyramids(const int *vi, int attr=1)
Adds 6 pyramids to the mesh by splitting a hexahedron given by 8 vertices vi.
static void PrintElementsByGeometry(int dim, const Array< int > &num_elems_by_geom, std::ostream &os)
Auxiliary method used by PrintCharacteristics().
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Adds a hexahedron to the mesh given by 8 vertices v1 through v8.
static IntegrationPoint TransformBdrElementToFace(Geometry::Type geom, int o, const IntegrationPoint &ip)
For the vertex (1D), edge (2D), or face (3D) of a boundary element with the orientation o,...
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
int GetNBE() const
Returns number of boundary elements.
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
void AddQuadFaceElement(int lf, int gf, int el, int v0, int v1, int v2, int v3)
void ReadMFEMMesh(std::istream &input, int version, int &curved)
void KnotRemove(Array< Vector * > &kv)
For NURBS meshes, remove the knots in kv, for each direction.
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr) const
Builds the transformation defining the i-th edge element in EdTr. EdTr must be allocated in advance a...
void GetElementData(const Array< Element * > &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
void PrintSurfaces(const Table &Aface_face, std::ostream &os) const
Print set of disjoint surfaces:
bool IsSlaveFace(const FaceInfo &fi) const
void FreeElement(Element *E)
Array< Element * > boundary
virtual void MarkTetMeshForRefinement(const DSTable &v_to_v)
void GetPointMatrix(int i, DenseMatrix &pointmat) const
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
void GetLocalTriToTetTransformation(IsoparametricTransformation &loc, int i) const
virtual void PrintXG(std::ostream &os=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
NCMesh * ncmesh
Optional nonconforming mesh extension.
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
virtual void NURBSUniformRefinement(int rf=2, real_t tol=1.0e-12)
Refine NURBS mesh, with an optional refinement factor, generally anisotropic.
void DebugDump(std::ostream &os) const
Output an NCMesh-compatible debug dump.
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
bool RefineByError(const Array< real_t > &elem_error, real_t threshold, int nonconforming=-1, int nc_limit=0)
static Mesh MakeCartesian3DWith24TetsPerHex(int nx, int ny, int nz, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz*24 tetrahedrons.
STable3D * GetFacesTable()
Table * GetFaceEdgeTable() const
void EnsureNCMesh(bool simplices_nonconforming=false)
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
virtual MFEM_DEPRECATED void ReorientTetMesh()
void PrintVTK(std::ostream &os)
void MoveNodes(const Vector &displacements)
Array< GeometricFactors * > geom_factors
Optional geometric factors.
void SetMeshGen()
Determine the mesh generator bitmask meshgen, see MeshGenerator().
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
void Make2D(int nx, int ny, Element::Type type, real_t sx, real_t sy, bool generate_edges, bool sfc_ordering)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
void AddHexAsTets(const int *vi, int attr=1)
Adds 6 tetrahedrons to the mesh by splitting a hexahedron given by 8 vertices vi.
FaceElementTransformations FaceElemTr
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
void GetLocalQuadToPyrTransformation(IsoparametricTransformation &loc, int i) const
void AddHexAsWedges(const int *vi, int attr=1)
Adds 2 wedges to the mesh by splitting a hexahedron given by 8 vertices vi.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Element * ReadElement(std::istream &input)
virtual bool HasBoundaryElements() const
Checks if the mesh has boundary elements.
void ScaleSubdomains(real_t sf)
void GetVertexToVertexTable(DSTable &) const
void GetLocalQuadToHexTransformation(IsoparametricTransformation &loc, int i) const
void Transform(void(*f)(const Vector &, Vector &))
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Geometry::Type GetElementBaseGeometry(int i) const
@ LocalSlaveNonconforming
@ SharedSlaveNonconforming
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones.
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
void ChangeVertexDataOwnership(real_t *vertices, int len_vertices, bool zerocopy=false)
Set the internal Vertex array to point to the given vertices array without assuming ownership of the ...
Table * GetVertexToElementTable()
void InitRefinementTransforms()
Table * GetEdgeVertexTable() const
int * GeneratePartitioning(int nparts, int part_method=1)
Table * GetFaceToElementTable() const
void MarkTriMeshForRefinement()
void ReadNetgen2DMesh(std::istream &input, int &curved)
Array< Element * > elements
Array< int > attributes
A list of all unique element attributes used by the Mesh.
void AddPointFaceElement(int lf, int gf, int el)
Used in GenerateFaces()
void RemoveInternalBoundaries()
virtual void NonconformingRefinement(const Array< Refinement > &refinements, int nc_limit=0)
This function is not public anymore. Use GeneralRefinement instead.
void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
void MoveVertices(const Vector &displacements)
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
void DeleteGeometricFactors()
Destroy all GeometricFactors stored by the Mesh.
const Table & ElementToFaceTable() const
void AddQuadAs4TrisWithPoints(int *vi, int attr=1)
Adds 4 triangles to the mesh by splitting a quadrilateral given by 4 vertices vi.
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
void KnotInsert(Array< KnotVector * > &kv)
For NURBS meshes, insert the new knots in kv, for each direction.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
void OnMeshUpdated(Mesh *mesh)
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
const CoarseFineTransformations & GetRefinementTransforms() const
int Dimension() const
Return the dimension of the NCMesh.
BlockArray< Element > elements
static void GridSfcOrdering3D(int width, int height, int depth, Array< int > &coords)
Array< int > leaf_elements
finest elements, in Mesh ordering (+ ghosts)
virtual void LimitNCLevel(int max_nc_level)
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
virtual void Derefine(const Array< int > &derefs)
Array< real_t > coordinates
bool IsGhost(const Element &el) const
Return true if the Element el is a ghost element.
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
void Print(std::ostream &out, const std::string &comments="") const
int spaceDim
dimensions of the elements and the vertex coordinates
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
const Table & GetDerefinementTable()
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)
void GetCoarseningFactors(Array< int > &f) const
void SetPatchAttribute(int i, int attr)
const Array< int > & GetPatchElements(int 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
void SetPatchBdrAttribute(int i, int attr)
void SetCoordsFromPatches(Vector &Nodes)
void Coarsen(int cf=2, real_t tol=1.0e-12)
int GetPatchAttribute(int i) const
void GetElementTopo(Array< Element * > &elements) const
const Array< int > & GetPatchBdrElements(int patch)
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
void GetBdrElementTopo(Array< Element * > &boundary) const
void KnotRemove(Array< Vector * > &kv, real_t tol=1.0e-12)
void GetElementLocalToGlobal(Array< int > &lelem_elem)
void KnotInsert(Array< KnotVector * > &kv)
void ConvertToPatches(const Vector &Nodes)
void SetKnotsFromPatches()
int GetPatchBdrAttribute(int i) const
void DegreeElevate(int rel_degree, int degree=16)
Class for standard nodal finite elements.
Class used to extrude the nodes of a mesh.
void SetLayer(const int l)
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
NodeExtrudeCoefficient(const int dim, const int n_, const real_t s_)
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.
Data type Pyramid element.
Piecewise-(bi)quadratic continuous finite elements.
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
@ VALUES
Evaluate the values at quadrature points.
@ DERIVATIVES
Evaluate the derivatives at quadrature points.
@ DETERMINANTS
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
Data type quadrilateral element.
Symmetric 3D Table stored as an array of rows each of which has a stack of column,...
int Push(int r, int c, int f)
Check to see if this entry is in the table and add it to the table if it is not there....
int NumberOfElements()
Return the number of elements added to the table.
int Push4(int r, int c, int f, int t)
Check to see if this entry is in the table and add it to the table if it is not there....
Data type line segment element.
void GetVertices(Array< int > &v) const override
Get the indices defining the vertices.
void Start()
Start the stopwatch. The elapsed time is not cleared.
double UserTime()
Return the number of user seconds elapsed since the stopwatch was started.
void AddConnections(int r, const int *c, int nc)
void SetSize(int dim, int connections_per_row)
Set the size and the number of connections for the table.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
void AddConnection(int r, int c)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
int Size() const
Returns the number of TYPE I elements.
int Size_of_connections() const
void AddColumnsInRow(int r, int ncol)
void AddAColumnInRow(int r)
void SetDims(int rows, int nnz)
Data type tetrahedron element.
void Init(int ind1, int ind2, int ind3, int ind4, int attr=1, int ref_flag=0)
Initialize the vertex indices and the attribute of a Tetrahedron.
void PushTransform(int tr) override
Add 'tr' to the current chain of coarse-fine transformations.
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag) const
int GetRefinementFlag() const
void SetVertices(const Array< int > &v) override
Set the indices defining the vertices.
void ResetTransform(int tr) override
Set current coarse-fine transformation number.
unsigned GetTransform() const override
Return current coarse-fine transformation.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void GetMarkedFace(const int face, int *fv) const
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
void GetVertices(Array< int > &v) const override
Get the indices defining the vertices.
Data type triangle element.
void SetVertices(const Array< int > &v) override
Set the indices defining the vertices.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void PushTransform(int tr) override
Add 'tr' to the current chain of coarse-fine transformations.
void ResetTransform(int tr) override
Set current coarse-fine transformation number.
unsigned GetTransform() const override
Return current coarse-fine transformation.
void GetVertices(Array< int > &v) const override
Get the indices defining the vertices.
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
A general vector function coefficient.
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
real_t Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
void NewDataAndSize(real_t *d, int s)
Set the Vector data and size, deleting the old data, if owned.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
void cross3D(const Vector &vin, Vector &vout) const
real_t DistanceTo(const real_t *p) const
Compute the Euclidean distance to another vector.
void Print(const Mesh &mesh, const adios2stream::mode print_mode=mode::sync)
const std::string filename
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
int index(int i, int j, int nx, int ny)
void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void METIS_PartGraphVKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
Linear1DFiniteElement SegmentFE
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
PointFiniteElement PointFE
TriLinear3DFiniteElement HexahedronFE
void Swap(Array< T > &, Array< T > &)
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
GeometryRefiner GlobGeometryRefiner
int FindRoots(const Vector &z, Vector &x)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
real_t rand_real()
Generate a random real_t number in the interval [0,1) using rand().
void add(const Vector &v1, const Vector &v2, Vector &v)
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
void ShiftRight(int &a, int &b, int &c)
MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
Mesh * Extrude1D(Mesh *mesh, const int ny, const real_t sy, const bool closed)
Extrude a 1D mesh.
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
void WriteBase64WithSizeAndClear(std::ostream &os, std::vector< char > &buf, int compression_level)
Encode in base 64 (and potentially compress) the given data, write it to the output stream (with a he...
Mesh * Extrude2D(Mesh *mesh, const int nz, const real_t sz)
Extrude a 2D mesh.
VTKFormat
Data array format for VTK and VTU files.
@ ASCII
Data arrays will be written in ASCII format.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
void FindTMax(Vector &c, Vector &x, real_t &tmax, const real_t factor, const int Dim)
BiLinear2DFiniteElement QuadrilateralFE
void WriteBinaryOrASCII(std::ostream &os, std::vector< char > &buf, const T &val, const char *suffix, VTKFormat format)
Write either ASCII data to the stream or binary data to the buffer depending on the given format.
@ byNODES
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
class LinearPyramidFiniteElement PyramidFE
const T & AsConst(const T &a)
Utility function similar to std::as_const in c++17.
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.
@ LEXICOGRAPHIC
Lexicographic ordering for tensor-product FiniteElements.
@ 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.