15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
17 #include "../general/binaryio.hpp"
18 #include "../general/text.hpp"
19 #include "../general/device.hpp"
20 #include "../general/tic_toc.hpp"
21 #include "../general/gecko.hpp"
22 #include "../fem/quadinterpolator.hpp"
37 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
42 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
47 int*,
int*,
int*,
int*,
int*,
idxtype*);
49 int*,
int*,
int*,
int*,
int*,
idxtype*);
51 int*,
int*,
int*,
int*,
int*,
idxtype*);
68 void Mesh::GetElementCenter(
int i,
Vector ¢er)
71 int geom = GetElementBaseGeometry(i);
86 return pow(fabs(J.
Weight()), 1./Dim);
98 double Mesh::GetElementSize(
int i,
int type)
100 return GetElementSize(GetElementTransformation(i), type);
103 double Mesh::GetElementSize(
int i,
const Vector &dir)
107 GetElementJacobian(i, J);
109 return sqrt((d_hat * d_hat) / (dir * dir));
112 double Mesh::GetElementVolume(
int i)
134 for (
int d = 0; d < spaceDim; d++)
143 for (
int i = 0; i < NumOfVertices; i++)
145 coord = GetVertex(i);
146 for (
int d = 0; d < spaceDim; d++)
148 if (coord[d] < min(d)) { min(d) = coord[d]; }
149 if (coord[d] > max(d)) { max(d) = coord[d]; }
155 const bool use_boundary =
false;
156 int ne = use_boundary ? GetNBE() : GetNE();
164 for (
int i = 0; i < ne; i++)
168 GetBdrElementFace(i, &fn, &fo);
170 Tr = GetFaceElementTransformations(fn, 5);
177 T = GetElementTransformation(i);
181 for (
int j = 0; j < pointmat.
Width(); j++)
183 for (
int d = 0; d < pointmat.
Height(); d++)
185 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
186 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
193 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
194 double &kappa_min,
double &kappa_max,
202 sdim = SpaceDimension();
204 if (Vh) { Vh->
SetSize(NumOfElements); }
205 if (Vk) { Vk->
SetSize(NumOfElements); }
208 h_max = kappa_max = -h_min;
209 if (dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
211 for (i = 0; i < NumOfElements; i++)
213 GetElementJacobian(i, J);
214 h = pow(fabs(J.
Weight()), 1.0/
double(dim));
215 kappa = (dim == sdim) ?
217 if (Vh) { (*Vh)(i) = h; }
218 if (Vk) { (*Vk)(i) = kappa; }
220 if (h < h_min) { h_min = h; }
221 if (h > h_max) { h_max = h; }
222 if (kappa < kappa_min) { kappa_min =
kappa; }
223 if (kappa > kappa_max) { kappa_max =
kappa; }
228 void Mesh::PrintElementsByGeometry(
int dim,
232 for (
int g = Geometry::DimStart[dim], first = 1;
233 g < Geometry::DimStart[dim+1]; g++)
235 if (!num_elems_by_geom[g]) {
continue; }
236 if (!first) { os <<
" + "; }
238 os << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
242 void Mesh::PrintCharacteristics(
Vector *Vh,
Vector *Vk, std::ostream &os)
244 double h_min, h_max, kappa_min, kappa_max;
246 os <<
"Mesh Characteristics:";
248 this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
250 Array<int> num_elems_by_geom(Geometry::NumGeom);
251 num_elems_by_geom = 0;
252 for (
int i = 0; i < GetNE(); i++)
254 num_elems_by_geom[GetElementBaseGeometry(i)]++;
258 <<
"Dimension : " << Dimension() <<
'\n'
259 <<
"Space dimension : " << SpaceDimension();
263 <<
"Number of vertices : " << GetNV() <<
'\n'
264 <<
"Number of elements : " << GetNE() <<
'\n'
265 <<
"Number of bdr elem : " << GetNBE() <<
'\n';
270 <<
"Number of vertices : " << GetNV() <<
'\n'
271 <<
"Number of elements : " << GetNE() <<
'\n'
272 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
273 <<
"h_min : " << h_min <<
'\n'
274 <<
"h_max : " << h_max <<
'\n';
279 <<
"Number of vertices : " << GetNV() <<
'\n'
280 <<
"Number of edges : " << GetNEdges() <<
'\n'
281 <<
"Number of elements : " << GetNE() <<
" -- ";
282 PrintElementsByGeometry(2, num_elems_by_geom, os);
284 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
285 <<
"Euler Number : " << EulerNumber2D() <<
'\n'
286 <<
"h_min : " << h_min <<
'\n'
287 <<
"h_max : " << h_max <<
'\n'
288 <<
"kappa_min : " << kappa_min <<
'\n'
289 <<
"kappa_max : " << kappa_max <<
'\n';
293 Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
294 num_bdr_elems_by_geom = 0;
295 for (
int i = 0; i < GetNBE(); i++)
297 num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
299 Array<int> num_faces_by_geom(Geometry::NumGeom);
300 num_faces_by_geom = 0;
301 for (
int i = 0; i < GetNFaces(); i++)
303 num_faces_by_geom[GetFaceBaseGeometry(i)]++;
307 <<
"Number of vertices : " << GetNV() <<
'\n'
308 <<
"Number of edges : " << GetNEdges() <<
'\n'
309 <<
"Number of faces : " << GetNFaces() <<
" -- ";
310 PrintElementsByGeometry(Dim-1, num_faces_by_geom, os);
312 <<
"Number of elements : " << GetNE() <<
" -- ";
313 PrintElementsByGeometry(Dim, num_elems_by_geom, os);
315 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
316 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, os);
318 <<
"Euler Number : " << EulerNumber() <<
'\n'
319 <<
"h_min : " << h_min <<
'\n'
320 <<
"h_max : " << h_max <<
'\n'
321 <<
"kappa_min : " << kappa_min <<
'\n'
322 <<
"kappa_max : " << kappa_max <<
'\n';
324 os <<
'\n' << std::flush;
331 case Element::POINT :
return &
PointFE;
332 case Element::SEGMENT :
return &
SegmentFE;
337 case Element::WEDGE :
return &
WedgeFE;
338 case Element::PYRAMID :
return &
PyramidFE;
340 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
343 MFEM_ABORT(
"Unknown element type");
352 ElTr->
ElementType = ElementTransformation::ELEMENT;
358 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
364 Nodes->FESpace()->GetElementVDofs(i, vdofs);
367 int n = vdofs.
Size()/spaceDim;
369 for (
int k = 0; k < spaceDim; k++)
371 for (
int j = 0; j < n; j++)
373 pm(k,j) = nodes(vdofs[n*k+j]);
376 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
380 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
385 ElTr->
ElementType = ElementTransformation::ELEMENT;
392 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
393 int nv = elements[i]->GetNVertices();
394 const int *v = elements[i]->GetVertices();
395 int n = vertices.Size();
397 for (
int k = 0; k < spaceDim; k++)
399 for (
int j = 0; j < nv; j++)
401 pm(k, j) = nodes(k*n+v[j]);
404 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
408 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
410 Nodes->FESpace()->GetElementVDofs(i, vdofs);
411 int n = vdofs.
Size()/spaceDim;
413 for (
int k = 0; k < spaceDim; k++)
415 for (
int j = 0; j < n; j++)
417 pm(k,j) = nodes(vdofs[n*k+j]);
420 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
426 GetElementTransformation(i, &Transformation);
428 return &Transformation;
433 GetBdrElementTransformation(i, &BdrTransformation);
434 return &BdrTransformation;
441 ElTr->
ElementType = ElementTransformation::BDR_ELEMENT;
447 GetBdrPointMatrix(i, pm);
448 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
458 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
459 int n = vdofs.
Size()/spaceDim;
461 for (
int k = 0; k < spaceDim; k++)
463 for (
int j = 0; j < n; j++)
465 pm(k,j) = nodes(vdofs[n*k+j]);
472 int elem_id, face_info;
473 GetBdrElementAdjacentElement2(i, elem_id, face_info);
475 GetLocalFaceTransformation(GetBdrElementType(i),
476 GetElementType(elem_id),
477 FaceElemTr.Loc1.Transf, face_info);
482 Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
483 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
484 "Mesh requires nodal Finite Element.");
486 FaceElemTr.Loc1.Transf.ElementNo = elem_id;
487 FaceElemTr.Loc1.Transf.mesh =
this;
488 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
489 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
490 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
492 ElTr->
SetFE(face_el);
499 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
507 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
508 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
510 for (
int i = 0; i < spaceDim; i++)
512 for (
int j = 0; j < nv; j++)
514 pm(i, j) = vertices[v[j]](i);
517 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
521 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
527 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
528 int n = vdofs.
Size()/spaceDim;
530 for (
int i = 0; i < spaceDim; i++)
532 for (
int j = 0; j < n; j++)
534 pm(i, j) = nodes(vdofs[n*i+j]);
541 FaceInfo &face_info = faces_info[FaceNo];
546 GetLocalFaceTransformation(face_type,
547 GetElementType(face_info.
Elem1No),
548 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
551 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
553 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
554 "Mesh requires nodal Finite Element.");
557 FaceElemTr.Loc1.Transf.ElementNo = face_info.
Elem1No;
558 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
559 FaceElemTr.Loc1.Transf.mesh =
this;
560 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
561 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
570 GetFaceTransformation(FaceNo, &FaceTransformation);
571 return &FaceTransformation;
578 GetFaceTransformation(EdgeNo, EdTr);
583 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
595 GetEdgeVertices(EdgeNo, v);
598 for (
int i = 0; i < spaceDim; i++)
600 for (
int j = 0; j < nv; j++)
602 pm(i, j) = vertices[v[j]](i);
605 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
609 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
615 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
616 int n = vdofs.
Size()/spaceDim;
618 for (
int i = 0; i < spaceDim; i++)
620 for (
int j = 0; j < n; j++)
622 pm(i, j) = nodes(vdofs[n*i+j]);
625 EdTr->
SetFE(edge_el);
629 MFEM_ABORT(
"Not implemented.");
636 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
637 return &EdgeTransformation;
641 void Mesh::GetLocalPtToSegTransformation(
656 void Mesh::GetLocalSegToTriTransformation(
665 tv = tri_t::Edges[i/64];
666 so = seg_t::Orient[i%64];
669 for (
int j = 0; j < 2; j++)
671 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
672 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
676 void Mesh::GetLocalSegToQuadTransformation(
685 qv = quad_t::Edges[i/64];
686 so = seg_t::Orient[i%64];
689 for (
int j = 0; j < 2; j++)
691 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
692 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
696 void Mesh::GetLocalTriToTetTransformation(
704 const int *tv = tet_t::FaceVert[i/64];
707 const int *to = tri_t::Orient[i%64];
711 for (
int j = 0; j < 3; j++)
714 locpm(0, j) = vert.
x;
715 locpm(1, j) = vert.
y;
716 locpm(2, j) = vert.
z;
720 void Mesh::GetLocalTriToWdgTransformation(
728 MFEM_VERIFY(i < 128,
"Local face index " << i/64
729 <<
" is not a triangular face of a wedge.");
730 const int *pv = pri_t::FaceVert[i/64];
733 const int *to = tri_t::Orient[i%64];
737 for (
int j = 0; j < 3; j++)
740 locpm(0, j) = vert.
x;
741 locpm(1, j) = vert.
y;
742 locpm(2, j) = vert.
z;
746 void Mesh::GetLocalTriToPyrTransformation(
753 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
754 <<
" is not a triangular face of a pyramid.");
755 const int *pv = pyr_t::FaceVert[i/64];
758 const int *to = tri_t::Orient[i%64];
762 for (
int j = 0; j < 3; j++)
765 locpm(0, j) = vert.
x;
766 locpm(1, j) = vert.
y;
767 locpm(2, j) = vert.
z;
771 void Mesh::GetLocalQuadToHexTransformation(
779 const int *hv = hex_t::FaceVert[i/64];
781 const int *qo = quad_t::Orient[i%64];
784 for (
int j = 0; j < 4; j++)
787 locpm(0, j) = vert.
x;
788 locpm(1, j) = vert.
y;
789 locpm(2, j) = vert.
z;
793 void Mesh::GetLocalQuadToWdgTransformation(
801 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
802 <<
" is not a quadrilateral face of a wedge.");
803 const int *pv = pri_t::FaceVert[i/64];
805 const int *qo = quad_t::Orient[i%64];
808 for (
int j = 0; j < 4; j++)
811 locpm(0, j) = vert.
x;
812 locpm(1, j) = vert.
y;
813 locpm(2, j) = vert.
z;
817 void Mesh::GetLocalQuadToPyrTransformation(
824 MFEM_VERIFY(i < 64,
"Local face index " << i/64
825 <<
" is not a quadrilateral face of a pyramid.");
826 const int *pv = pyr_t::FaceVert[i/64];
828 const int *qo = quad_t::Orient[i%64];
831 for (
int j = 0; j < 4; j++)
834 locpm(0, j) = vert.
x;
835 locpm(1, j) = vert.
y;
836 locpm(2, j) = vert.
z;
844 for (
int i = 0; i < geom_factors.Size(); i++)
856 geom_factors.Append(gf);
864 for (
int i = 0; i < face_geom_factors.Size(); i++)
878 face_geom_factors.Append(gf);
882 void Mesh::DeleteGeometricFactors()
884 for (
int i = 0; i < geom_factors.Size(); i++)
886 delete geom_factors[i];
888 geom_factors.SetSize(0);
889 for (
int i = 0; i < face_geom_factors.Size(); i++)
891 delete face_geom_factors[i];
893 face_geom_factors.SetSize(0);
896 void Mesh::GetLocalFaceTransformation(
902 GetLocalPtToSegTransformation(Transf, info);
905 case Element::SEGMENT:
906 if (elem_type == Element::TRIANGLE)
908 GetLocalSegToTriTransformation(Transf, info);
912 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
913 GetLocalSegToQuadTransformation(Transf, info);
917 case Element::TRIANGLE:
918 if (elem_type == Element::TETRAHEDRON)
920 GetLocalTriToTetTransformation(Transf, info);
922 else if (elem_type == Element::WEDGE)
924 GetLocalTriToWdgTransformation(Transf, info);
926 else if (elem_type == Element::PYRAMID)
928 GetLocalTriToPyrTransformation(Transf, info);
932 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
933 "face type " << face_type
934 <<
" and element type " << elem_type <<
"\n");
938 case Element::QUADRILATERAL:
939 if (elem_type == Element::HEXAHEDRON)
941 GetLocalQuadToHexTransformation(Transf, info);
943 else if (elem_type == Element::WEDGE)
945 GetLocalQuadToWdgTransformation(Transf, info);
947 else if (elem_type == Element::PYRAMID)
949 GetLocalQuadToPyrTransformation(Transf, info);
953 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
954 "face type " << face_type
955 <<
" and element type " << elem_type <<
"\n");
964 FaceInfo &face_info = faces_info[FaceNo];
967 FaceElemTr.SetConfigurationMask(cmask);
968 FaceElemTr.Elem1 = NULL;
969 FaceElemTr.Elem2 = NULL;
973 if (mask & FaceElementTransformations::HAVE_ELEM1)
975 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
976 FaceElemTr.Elem1 = &Transformation;
983 FaceElemTr.Elem2No = face_info.
Elem2No;
984 if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
985 FaceElemTr.Elem2No >= 0)
988 if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
989 { MFEM_ABORT(
"NURBS mesh not supported!"); }
991 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
992 FaceElemTr.Elem2 = &Transformation2;
997 if (mask & FaceElementTransformations::HAVE_FACE)
999 GetFaceTransformation(FaceNo, &FaceElemTr);
1004 FaceElemTr.SetGeometryType(GetFaceGeometryType(FaceNo));
1008 int face_type = GetFaceElementType(FaceNo);
1009 if (mask & FaceElementTransformations::HAVE_LOC1)
1011 int elem_type = GetElementType(face_info.
Elem1No);
1012 GetLocalFaceTransformation(face_type, elem_type,
1013 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
1016 if ((mask & FaceElementTransformations::HAVE_LOC2) &&
1017 FaceElemTr.Elem2No >= 0)
1019 int elem_type = GetElementType(face_info.
Elem2No);
1020 GetLocalFaceTransformation(face_type, elem_type,
1021 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
1024 if (Nonconforming() && IsSlaveFace(face_info))
1026 ApplyLocalSlaveTransformation(FaceElemTr, face_info,
false);
1031 FaceElemTr.SetConfigurationMask(cmask);
1037 double dist = FaceElemTr.CheckConsistency();
1040 mfem::out <<
"\nInternal error: face id = " << FaceNo
1041 <<
", dist = " << dist <<
'\n';
1042 FaceElemTr.CheckConsistency(1);
1043 MFEM_ABORT(
"internal error");
1053 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
1059 #ifdef MFEM_THREAD_SAFE
1064 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1065 MFEM_ASSERT(nc_faces_info[fi.
NCFace].Slave,
"internal error");
1076 std::swap(composition(0,0), composition(0,1));
1077 std::swap(composition(1,0), composition(1,1));
1098 int fn = GetBdrFace(BdrElemNo);
1101 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
1105 tr = GetFaceElementTransformations(fn, 21);
1106 tr->
Attribute = boundary[BdrElemNo]->GetAttribute();
1108 tr->
ElementType = ElementTransformation::BDR_FACE;
1113 int Mesh::GetBdrFace(
int BdrElemNo)
const
1118 fn = be_to_face[BdrElemNo];
1122 fn = be_to_edge[BdrElemNo];
1126 fn = boundary[BdrElemNo]->GetVertices()[0];
1137 GetFaceElements(f, &e1, &e2);
1138 GetFaceInfos(f, &inf1, &inf2, &ncface);
1148 if (f < GetNumFaces())
1154 face.
tag = FaceInfoTag::LocalConforming;
1155 face.
topology = FaceTopology::Conforming;
1164 face.
tag = FaceInfoTag::LocalSlaveNonconforming;
1165 face.
topology = FaceTopology::Nonconforming;
1170 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1181 face.
tag = FaceInfoTag::Boundary;
1182 face.
topology = FaceTopology::Boundary;
1191 face.
tag = FaceInfoTag::SharedConforming;
1192 face.
topology = FaceTopology::Conforming;
1204 face.
tag = FaceInfoTag::MasterNonconforming;
1205 face.
topology = FaceTopology::Nonconforming;
1214 face.
tag = FaceInfoTag::SharedSlaveNonconforming;
1215 face.
topology = FaceTopology::Nonconforming;
1230 face.
tag = FaceInfoTag::GhostMaster;
1240 face.
tag = FaceInfoTag::GhostSlave;
1241 face.
topology = FaceTopology::Nonconforming;
1258 case FaceInfoTag::LocalConforming:
1259 res.
Elem1No = element[0].index;
1260 res.Elem2No = element[1].index;
1261 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1262 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1263 res.NCFace = ncface;
1265 case FaceInfoTag::LocalSlaveNonconforming:
1266 res.Elem1No = element[0].index;
1267 res.Elem2No = element[1].index;
1268 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1269 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1270 res.NCFace = ncface;
1272 case FaceInfoTag::Boundary:
1273 res.Elem1No = element[0].index;
1274 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1276 case FaceInfoTag::SharedConforming:
1277 res.Elem1No = element[0].index;
1278 res.Elem2No = -1 - element[1].index;
1279 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1280 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1282 case FaceInfoTag::MasterNonconforming:
1283 res.Elem1No = element[0].index;
1284 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1286 case FaceInfoTag::SharedSlaveNonconforming:
1287 res.Elem1No = element[0].index;
1288 res.Elem2No = -1 - element[1].index;
1289 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1290 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1292 case FaceInfoTag::GhostMaster:
1294 case FaceInfoTag::GhostSlave:
1295 res.Elem1No = element[0].index;
1296 res.Elem2No = -1 - element[1].index;
1297 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1298 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1306 os <<
"face topology=";
1309 case Mesh::FaceTopology::Boundary:
1312 case Mesh::FaceTopology::Conforming:
1315 case Mesh::FaceTopology::Nonconforming:
1316 os <<
"Non-conforming";
1318 case Mesh::FaceTopology::NA:
1322 os <<
"element[0].location=";
1325 case Mesh::ElementLocation::Local:
1328 case Mesh::ElementLocation::FaceNbr:
1331 case Mesh::ElementLocation::NA:
1336 os <<
"element[1].location=";
1339 case Mesh::ElementLocation::Local:
1342 case Mesh::ElementLocation::FaceNbr:
1345 case Mesh::ElementLocation::NA:
1350 os <<
"element[0].conformity=";
1353 case Mesh::ElementConformity::Coincident:
1356 case Mesh::ElementConformity::Superset:
1359 case Mesh::ElementConformity::Subset:
1362 case Mesh::ElementConformity::NA:
1367 os <<
"element[1].conformity=";
1370 case Mesh::ElementConformity::Coincident:
1373 case Mesh::ElementConformity::Superset:
1376 case Mesh::ElementConformity::Subset:
1379 case Mesh::ElementConformity::NA:
1384 os <<
"element[0].index=" << info.
element[0].
index << std::endl
1385 <<
"element[1].index=" << info.
element[1].
index << std::endl
1390 <<
"ncface=" << info.
ncface << std::endl;
1394 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const
1396 *Elem1 = faces_info[Face].Elem1No;
1397 *Elem2 = faces_info[Face].Elem2No;
1400 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const
1402 *Inf1 = faces_info[Face].Elem1Inf;
1403 *Inf2 = faces_info[Face].Elem2Inf;
1406 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2,
int *NCFace)
const
1408 *Inf1 = faces_info[Face].Elem1Inf;
1409 *Inf2 = faces_info[Face].Elem2Inf;
1410 *NCFace = faces_info[Face].NCFace;
1417 case 1:
return Geometry::POINT;
1418 case 2:
return Geometry::SEGMENT;
1420 if (Face < NumOfFaces)
1422 return faces[Face]->GetGeometryType();
1425 const int nc_face_id = faces_info[Face].NCFace;
1426 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1427 return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
1429 return Geometry::INVALID;
1434 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1441 for (
int i = 0; i < NumOfBdrElements; i++)
1443 face_to_be[GetBdrElementEdgeIndex(i)] = i;
1453 NumOfElements = NumOfBdrElements = 0;
1454 NumOfEdges = NumOfFaces = 0;
1455 nbInteriorFaces = -1;
1456 nbBoundaryFaces = -1;
1457 meshgen = mesh_geoms = 0;
1463 last_operation = Mesh::NONE;
1466 void Mesh::InitTables()
1469 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1472 void Mesh::SetEmpty()
1478 void Mesh::DestroyTables()
1483 DeleteGeometricFactors();
1494 void Mesh::DestroyPointers()
1496 if (own_nodes) {
delete Nodes; }
1502 for (
int i = 0; i < NumOfElements; i++)
1504 FreeElement(elements[i]);
1507 for (
int i = 0; i < NumOfBdrElements; i++)
1509 FreeElement(boundary[i]);
1512 for (
int i = 0; i < faces.Size(); i++)
1514 FreeElement(faces[i]);
1520 void Mesh::Destroy()
1524 elements.DeleteAll();
1525 vertices.DeleteAll();
1526 boundary.DeleteAll();
1528 faces_info.DeleteAll();
1529 nc_faces_info.DeleteAll();
1530 be_to_edge.DeleteAll();
1531 be_to_face.DeleteAll();
1539 CoarseFineTr.Clear();
1541 #ifdef MFEM_USE_MEMALLOC
1545 attributes.DeleteAll();
1546 bdr_attributes.DeleteAll();
1549 void Mesh::ResetLazyData()
1551 delete el_to_el; el_to_el = NULL;
1552 delete face_edge; face_edge = NULL;
1553 delete edge_vertex; edge_vertex = NULL;
1554 DeleteGeometricFactors();
1555 nbInteriorFaces = -1;
1556 nbBoundaryFaces = -1;
1559 void Mesh::SetAttributes()
1564 for (
int i = 0; i < attribs.
Size(); i++)
1566 attribs[i] = GetBdrAttribute(i);
1570 attribs.
Copy(bdr_attributes);
1571 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1573 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1577 for (
int i = 0; i < attribs.
Size(); i++)
1579 attribs[i] = GetAttribute(i);
1583 attribs.
Copy(attributes);
1584 if (attributes.Size() > 0 && attributes[0] <= 0)
1586 MFEM_WARNING(
"Non-positive attributes in the domain!");
1590 void Mesh::InitMesh(
int Dim_,
int spaceDim_,
int NVert,
int NElem,
int NBdrElem)
1595 spaceDim = spaceDim_;
1598 vertices.SetSize(NVert);
1601 elements.SetSize(NElem);
1603 NumOfBdrElements = 0;
1604 boundary.SetSize(NBdrElem);
1607 template<
typename T>
1608 static void CheckEnlarge(
Array<T> &array,
int size)
1610 if (size >= array.
Size()) { array.
SetSize(size + 1); }
1613 int Mesh::AddVertex(
double x,
double y,
double z)
1615 CheckEnlarge(vertices, NumOfVertices);
1616 double *v = vertices[NumOfVertices]();
1620 return NumOfVertices++;
1623 int Mesh::AddVertex(
const double *coords)
1625 CheckEnlarge(vertices, NumOfVertices);
1626 vertices[NumOfVertices].SetCoords(spaceDim, coords);
1627 return NumOfVertices++;
1630 void Mesh::AddVertexParents(
int i,
int p1,
int p2)
1636 if (i < vertices.Size())
1638 double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1639 for (
int j = 0; j < 3; j++)
1641 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1646 int Mesh::AddSegment(
int v1,
int v2,
int attr)
1648 CheckEnlarge(elements, NumOfElements);
1649 elements[NumOfElements] =
new Segment(v1, v2, attr);
1650 return NumOfElements++;
1653 int Mesh::AddSegment(
const int *vi,
int attr)
1655 CheckEnlarge(elements, NumOfElements);
1656 elements[NumOfElements] =
new Segment(vi, attr);
1657 return NumOfElements++;
1660 int Mesh::AddTriangle(
int v1,
int v2,
int v3,
int attr)
1662 CheckEnlarge(elements, NumOfElements);
1663 elements[NumOfElements] =
new Triangle(v1, v2, v3, attr);
1664 return NumOfElements++;
1667 int Mesh::AddTriangle(
const int *vi,
int attr)
1669 CheckEnlarge(elements, NumOfElements);
1670 elements[NumOfElements] =
new Triangle(vi, attr);
1671 return NumOfElements++;
1674 int Mesh::AddQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1676 CheckEnlarge(elements, NumOfElements);
1677 elements[NumOfElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1678 return NumOfElements++;
1681 int Mesh::AddQuad(
const int *vi,
int attr)
1683 CheckEnlarge(elements, NumOfElements);
1685 return NumOfElements++;
1688 int Mesh::AddTet(
int v1,
int v2,
int v3,
int v4,
int attr)
1690 int vi[4] = {v1, v2, v3, v4};
1691 return AddTet(vi, attr);
1694 int Mesh::AddTet(
const int *vi,
int attr)
1696 CheckEnlarge(elements, NumOfElements);
1697 #ifdef MFEM_USE_MEMALLOC
1699 tet = TetMemory.Alloc();
1702 elements[NumOfElements] = tet;
1704 elements[NumOfElements] =
new Tetrahedron(vi, attr);
1706 return NumOfElements++;
1709 int Mesh::AddWedge(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int attr)
1711 CheckEnlarge(elements, NumOfElements);
1712 elements[NumOfElements] =
new Wedge(v1, v2, v3, v4, v5, v6, attr);
1713 return NumOfElements++;
1716 int Mesh::AddWedge(
const int *vi,
int attr)
1718 CheckEnlarge(elements, NumOfElements);
1719 elements[NumOfElements] =
new Wedge(vi, attr);
1720 return NumOfElements++;
1723 int Mesh::AddPyramid(
int v1,
int v2,
int v3,
int v4,
int v5,
int attr)
1725 CheckEnlarge(elements, NumOfElements);
1726 elements[NumOfElements] =
new Pyramid(v1, v2, v3, v4, v5, attr);
1727 return NumOfElements++;
1730 int Mesh::AddPyramid(
const int *vi,
int attr)
1732 CheckEnlarge(elements, NumOfElements);
1733 elements[NumOfElements] =
new Pyramid(vi, attr);
1734 return NumOfElements++;
1737 int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1740 CheckEnlarge(elements, NumOfElements);
1741 elements[NumOfElements] =
1742 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1743 return NumOfElements++;
1746 int Mesh::AddHex(
const int *vi,
int attr)
1748 CheckEnlarge(elements, NumOfElements);
1749 elements[NumOfElements] =
new Hexahedron(vi, attr);
1750 return NumOfElements++;
1753 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1755 static const int hex_to_tet[6][4] =
1757 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1758 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1762 for (
int i = 0; i < 6; i++)
1764 for (
int j = 0; j < 4; j++)
1766 ti[j] = vi[hex_to_tet[i][j]];
1772 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1774 static const int hex_to_wdg[2][6] =
1776 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1780 for (
int i = 0; i < 2; i++)
1782 for (
int j = 0; j < 6; j++)
1784 ti[j] = vi[hex_to_wdg[i][j]];
1790 void Mesh::AddHexAsPyramids(
const int *vi,
int attr)
1792 static const int hex_to_pyr[6][5] =
1794 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
1795 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
1799 for (
int i = 0; i < 6; i++)
1801 for (
int j = 0; j < 5; j++)
1803 ti[j] = vi[hex_to_pyr[i][j]];
1805 AddPyramid(ti, attr);
1811 CheckEnlarge(elements, NumOfElements);
1812 elements[NumOfElements] = elem;
1813 return NumOfElements++;
1818 CheckEnlarge(boundary, NumOfBdrElements);
1819 boundary[NumOfBdrElements] = elem;
1820 return NumOfBdrElements++;
1823 int Mesh::AddBdrSegment(
int v1,
int v2,
int attr)
1825 CheckEnlarge(boundary, NumOfBdrElements);
1826 boundary[NumOfBdrElements] =
new Segment(v1, v2, attr);
1827 return NumOfBdrElements++;
1830 int Mesh::AddBdrSegment(
const int *vi,
int attr)
1832 CheckEnlarge(boundary, NumOfBdrElements);
1833 boundary[NumOfBdrElements] =
new Segment(vi, attr);
1834 return NumOfBdrElements++;
1837 int Mesh::AddBdrTriangle(
int v1,
int v2,
int v3,
int attr)
1839 CheckEnlarge(boundary, NumOfBdrElements);
1840 boundary[NumOfBdrElements] =
new Triangle(v1, v2, v3, attr);
1841 return NumOfBdrElements++;
1844 int Mesh::AddBdrTriangle(
const int *vi,
int attr)
1846 CheckEnlarge(boundary, NumOfBdrElements);
1847 boundary[NumOfBdrElements] =
new Triangle(vi, attr);
1848 return NumOfBdrElements++;
1851 int Mesh::AddBdrQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1853 CheckEnlarge(boundary, NumOfBdrElements);
1854 boundary[NumOfBdrElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1855 return NumOfBdrElements++;
1858 int Mesh::AddBdrQuad(
const int *vi,
int attr)
1860 CheckEnlarge(boundary, NumOfBdrElements);
1862 return NumOfBdrElements++;
1865 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1867 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1870 for (
int i = 0; i < 2; i++)
1872 for (
int j = 0; j < 3; j++)
1874 ti[j] = vi[quad_to_tri[i][j]];
1876 AddBdrTriangle(ti, attr);
1880 int Mesh::AddBdrPoint(
int v,
int attr)
1882 CheckEnlarge(boundary, NumOfBdrElements);
1883 boundary[NumOfBdrElements] =
new Point(&v, attr);
1884 return NumOfBdrElements++;
1887 void Mesh::GenerateBoundaryElements()
1890 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1894 for (i = 0; i < boundary.Size(); i++)
1896 FreeElement(boundary[i]);
1906 NumOfBdrElements = 0;
1907 for (i = 0; i < faces_info.Size(); i++)
1909 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1912 boundary.SetSize(NumOfBdrElements);
1913 be2face.
SetSize(NumOfBdrElements);
1914 for (j = i = 0; i < faces_info.Size(); i++)
1916 if (faces_info[i].Elem2No < 0)
1918 boundary[j] = faces[i]->Duplicate(
this);
1925 void Mesh::FinalizeCheck()
1927 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1928 vertices.Size() == 0,
1929 "incorrect number of vertices: preallocated: " << vertices.Size()
1930 <<
", actually added: " << NumOfVertices);
1931 MFEM_VERIFY(elements.Size() == NumOfElements,
1932 "incorrect number of elements: preallocated: " << elements.Size()
1933 <<
", actually added: " << NumOfElements);
1934 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1935 "incorrect number of boundary elements: preallocated: "
1936 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1939 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1942 CheckElementOrientation(fix_orientation);
1946 MarkTriMeshForRefinement();
1951 el_to_edge =
new Table;
1952 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1954 CheckBdrElementOrientation();
1968 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1969 bool fix_orientation)
1972 if (fix_orientation)
1974 CheckElementOrientation(fix_orientation);
1979 el_to_edge =
new Table;
1980 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1982 CheckBdrElementOrientation();
2002 GeckoProgress(
double limit) : limit(limit) { sw.Start(); }
2003 virtual bool quit()
const {
return limit > 0 && sw.UserTime() > limit; }
2006 class GeckoVerboseProgress :
public GeckoProgress
2012 GeckoVerboseProgress(
double limit) : GeckoProgress(limit) {}
2014 virtual void beginorder(
const Graph* graph,
Float cost)
const
2015 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2016 virtual void endorder(
const Graph* graph,
Float cost)
const
2017 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2019 virtual void beginiter(
const Graph* graph,
2022 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
2023 << window << std::flush;
2025 virtual void enditer(
const Graph* graph,
Float mincost,
Float cost)
const
2026 {
mfem::out <<
", cost = " << cost << endl; }
2031 int iterations,
int window,
2032 int period,
int seed,
bool verbose,
2038 GeckoProgress progress(time_limit);
2039 GeckoVerboseProgress vprogress(time_limit);
2042 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2049 const Table &my_el_to_el = ElementToElementTable();
2050 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2052 const int *neighid = my_el_to_el.
GetRow(elemid);
2053 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2055 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2060 graph.
order(&functional, iterations, window, period, seed,
2061 verbose ? &vprogress : &progress);
2067 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2070 return graph.
cost();
2081 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
2082 : coord(coord), dir(dir), points(points), mid(mid) {}
2084 bool operator()(
int i)
const
2086 return (points[3*i + coord] < mid) != dir;
2090 static void HilbertSort2D(
int coord1,
2093 const Array<double> &points,
int *beg,
int *end,
2094 double xmin,
double ymin,
double xmax,
double ymax)
2096 if (end - beg <= 1) {
return; }
2098 double xmid = (xmin + xmax)*0.5;
2099 double ymid = (ymin + ymax)*0.5;
2101 int coord2 = (coord1 + 1) % 2;
2104 int *p0 = beg, *p4 = end;
2105 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2106 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2107 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2111 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2112 ymin, xmin, ymid, xmid);
2114 if (p1 != p0 || p2 != p4)
2116 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2117 xmin, ymid, xmid, ymax);
2119 if (p2 != p0 || p3 != p4)
2121 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2122 xmid, ymid, xmax, ymax);
2126 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2127 ymid, xmax, ymin, xmid);
2131 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2132 const Array<double> &points,
int *beg,
int *end,
2133 double xmin,
double ymin,
double zmin,
2134 double xmax,
double ymax,
double zmax)
2136 if (end - beg <= 1) {
return; }
2138 double xmid = (xmin + xmax)*0.5;
2139 double ymid = (ymin + ymax)*0.5;
2140 double zmid = (zmin + zmax)*0.5;
2142 int coord2 = (coord1 + 1) % 3;
2143 int coord3 = (coord1 + 2) % 3;
2146 int *p0 = beg, *p8 = end;
2147 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2148 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2149 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2150 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2151 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2152 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2153 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2157 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2158 zmin, xmin, ymin, zmid, xmid, ymid);
2160 if (p1 != p0 || p2 != p8)
2162 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2163 ymin, zmid, xmin, ymid, zmax, xmid);
2165 if (p2 != p0 || p3 != p8)
2167 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2168 ymid, zmid, xmin, ymax, zmax, xmid);
2170 if (p3 != p0 || p4 != p8)
2172 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2173 xmin, ymax, zmid, xmid, ymid, zmin);
2175 if (p4 != p0 || p5 != p8)
2177 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2178 xmid, ymax, zmid, xmax, ymid, zmin);
2180 if (p5 != p0 || p6 != p8)
2182 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2183 ymax, zmid, xmax, ymid, zmax, xmid);
2185 if (p6 != p0 || p7 != p8)
2187 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2188 ymid, zmid, xmax, ymin, zmax, xmid);
2192 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2193 zmid, xmax, ymin, zmin, xmid, ymid);
2199 MFEM_VERIFY(spaceDim <= 3,
"");
2202 GetBoundingBox(min, max);
2207 if (spaceDim < 3) { points = 0.0; }
2210 for (
int i = 0; i < GetNE(); i++)
2212 GetElementCenter(i, center);
2213 for (
int j = 0; j < spaceDim; j++)
2215 points[3*i + j] = center(j);
2222 indices.
Sort([&](
int a,
int b)
2223 {
return points[3*
a] < points[3*
b]; });
2225 else if (spaceDim == 2)
2228 HilbertSort2D(0,
false,
false,
2229 points, indices.
begin(), indices.
end(),
2230 min(0), min(1), max(0), max(1));
2235 HilbertSort3D(0,
false,
false,
false,
2236 points, indices.
begin(), indices.
end(),
2237 min(0), min(1), min(2), max(0), max(1), max(2));
2242 for (
int i = 0; i < GetNE(); i++)
2244 ordering[indices[i]] = i;
2249 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
2253 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2258 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
2262 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
2292 old_elem_node_vals.SetSize(GetNE());
2293 nodes_fes = Nodes->FESpace();
2296 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2298 nodes_fes->GetElementVDofs(old_elid, old_dofs);
2300 old_elem_node_vals[old_elid] =
new Vector(vals);
2306 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2308 int new_elid = ordering[old_elid];
2309 new_elements[new_elid] = elements[old_elid];
2314 if (reorder_vertices)
2319 vertex_ordering = -1;
2321 int new_vertex_ind = 0;
2322 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2324 int *elem_vert = elements[new_elid]->GetVertices();
2325 int nv = elements[new_elid]->GetNVertices();
2326 for (
int vi = 0; vi < nv; ++vi)
2328 int old_vertex_ind = elem_vert[vi];
2329 if (vertex_ordering[old_vertex_ind] == -1)
2331 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2332 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
2342 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2344 int *elem_vert = elements[new_elid]->GetVertices();
2345 int nv = elements[new_elid]->GetNVertices();
2346 for (
int vi = 0; vi < nv; ++vi)
2348 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2353 for (
int belid = 0; belid < GetNBE(); ++belid)
2355 int *be_vert = boundary[belid]->GetVertices();
2356 int nv = boundary[belid]->GetNVertices();
2357 for (
int vi = 0; vi < nv; ++vi)
2359 be_vert[vi] = vertex_ordering[be_vert[vi]];
2370 el_to_edge =
new Table;
2371 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2376 GetElementToFaceTable();
2386 last_operation = Mesh::NONE;
2387 nodes_fes->Update(
false);
2390 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2392 int new_elid = ordering[old_elid];
2393 nodes_fes->GetElementVDofs(new_elid, new_dofs);
2394 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
2395 delete old_elem_node_vals[old_elid];
2401 void Mesh::MarkForRefinement()
2407 MarkTriMeshForRefinement();
2411 DSTable v_to_v(NumOfVertices);
2412 GetVertexToVertexTable(v_to_v);
2413 MarkTetMeshForRefinement(v_to_v);
2418 void Mesh::MarkTriMeshForRefinement()
2423 for (
int i = 0; i < NumOfElements; i++)
2425 if (elements[i]->GetType() == Element::TRIANGLE)
2427 GetPointMatrix(i, pmat);
2428 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
2439 for (
int i = 0; i < NumOfVertices; i++)
2444 length_idx[j].one = GetLength(i, it.Column());
2445 length_idx[j].two = j;
2452 for (
int i = 0; i < NumOfEdges; i++)
2454 order[length_idx[i].two] = i;
2458 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
2463 GetEdgeOrdering(v_to_v, order);
2465 for (
int i = 0; i < NumOfElements; i++)
2467 if (elements[i]->GetType() == Element::TETRAHEDRON)
2469 elements[i]->MarkEdge(v_to_v, order);
2472 for (
int i = 0; i < NumOfBdrElements; i++)
2474 if (boundary[i]->GetType() == Element::TRIANGLE)
2476 boundary[i]->MarkEdge(v_to_v, order);
2483 if (*old_v_to_v && *old_elem_vert)
2490 if (*old_v_to_v == NULL)
2492 bool need_v_to_v =
false;
2494 for (
int i = 0; i < GetNEdges(); i++)
2500 if (dofs.
Size() > 0)
2508 *old_v_to_v =
new DSTable(NumOfVertices);
2509 GetVertexToVertexTable(*(*old_v_to_v));
2512 if (*old_elem_vert == NULL)
2514 bool need_elem_vert =
false;
2516 for (
int i = 0; i < GetNE(); i++)
2522 if (dofs.
Size() > 1)
2524 need_elem_vert =
true;
2530 *old_elem_vert =
new Table;
2531 (*old_elem_vert)->
MakeI(GetNE());
2532 for (
int i = 0; i < GetNE(); i++)
2534 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2536 (*old_elem_vert)->MakeJ();
2537 for (
int i = 0; i < GetNE(); i++)
2539 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2540 elements[i]->GetNVertices());
2542 (*old_elem_vert)->ShiftUpI();
2555 const int num_edge_dofs = old_dofs.
Size();
2558 const Vector onodes = *Nodes;
2562 int offset = NumOfVertices * old_dofs.
Size();
2566 if (num_edge_dofs > 0)
2568 DSTable new_v_to_v(NumOfVertices);
2569 GetVertexToVertexTable(new_v_to_v);
2571 for (
int i = 0; i < NumOfVertices; i++)
2575 const int old_i = (*old_v_to_v)(i, it.Column());
2576 const int new_i = it.Index();
2577 if (new_i == old_i) {
continue; }
2579 old_dofs.
SetSize(num_edge_dofs);
2580 new_dofs.
SetSize(num_edge_dofs);
2581 for (
int j = 0; j < num_edge_dofs; j++)
2583 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2584 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2588 for (
int j = 0; j < old_dofs.
Size(); j++)
2590 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2594 offset += NumOfEdges * num_edge_dofs;
2602 Table old_face_vertex;
2603 old_face_vertex.
MakeI(NumOfFaces);
2604 for (
int i = 0; i < NumOfFaces; i++)
2608 old_face_vertex.
MakeJ();
2609 for (
int i = 0; i < NumOfFaces; i++)
2611 faces[i]->GetNVertices());
2615 STable3D *faces_tbl = GetElementToFaceTable(1);
2621 for (
int i = 0; i < NumOfFaces; i++)
2623 const int *old_v = old_face_vertex.
GetRow(i);
2625 switch (old_face_vertex.
RowSize(i))
2628 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2632 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2636 new_fdofs[new_i+1] = old_dofs.
Size();
2641 for (
int i = 0; i < NumOfFaces; i++)
2643 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2646 switch (old_face_vertex.
RowSize(i))
2649 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2650 new_v = faces[new_i]->GetVertices();
2651 new_or = GetTriOrientation(old_v, new_v);
2656 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2657 new_v = faces[new_i]->GetVertices();
2658 new_or = GetQuadOrientation(old_v, new_v);
2665 for (
int j = 0; j < old_dofs.
Size(); j++)
2668 const int old_j = dof_ord[j];
2669 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2673 for (
int j = 0; j < old_dofs.
Size(); j++)
2675 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2697 for (
int i = 0; i < GetNE(); i++)
2699 const int *old_v = old_elem_vert->
GetRow(i);
2700 const int *new_v = elements[i]->GetVertices();
2706 case Geometry::SEGMENT:
2707 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2709 case Geometry::TRIANGLE:
2710 new_or = GetTriOrientation(old_v, new_v);
2712 case Geometry::SQUARE:
2713 new_or = GetQuadOrientation(old_v, new_v);
2715 case Geometry::TETRAHEDRON:
2716 new_or = GetTetOrientation(old_v, new_v);
2720 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2721 <<
" FE collection) are not supported yet!");
2725 MFEM_VERIFY(dof_ord != NULL,
2726 "FE collection '" << fec->
Name()
2727 <<
"' does not define reordering for "
2728 << Geometry::Name[geom] <<
" elements!");
2731 for (
int j = 0; j < new_dofs.
Size(); j++)
2734 const int old_j = dof_ord[j];
2735 new_dofs[old_j] = offset + j;
2737 offset += new_dofs.
Size();
2740 for (
int j = 0; j < old_dofs.
Size(); j++)
2742 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2754 GetElementToFaceTable();
2757 CheckBdrElementOrientation();
2762 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2767 CheckBdrElementOrientation();
2772 last_operation = Mesh::NONE;
2777 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2780 CheckElementOrientation(fix_orientation);
2782 if (NumOfBdrElements == 0)
2784 GetElementToFaceTable();
2786 GenerateBoundaryElements();
2791 DSTable v_to_v(NumOfVertices);
2792 GetVertexToVertexTable(v_to_v);
2793 MarkTetMeshForRefinement(v_to_v);
2796 GetElementToFaceTable();
2799 CheckBdrElementOrientation();
2801 if (generate_edges == 1)
2803 el_to_edge =
new Table;
2804 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2818 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2819 bool fix_orientation)
2822 CheckElementOrientation(fix_orientation);
2824 if (NumOfBdrElements == 0)
2826 GetElementToFaceTable();
2828 GenerateBoundaryElements();
2831 GetElementToFaceTable();
2834 CheckBdrElementOrientation();
2836 if (generate_edges == 1)
2838 el_to_edge =
new Table;
2839 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2853 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2856 CheckElementOrientation(fix_orientation);
2858 GetElementToFaceTable();
2861 if (NumOfBdrElements == 0)
2863 GenerateBoundaryElements();
2866 CheckBdrElementOrientation();
2870 el_to_edge =
new Table;
2871 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2883 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2887 Finalize(refine, fix_orientation);
2890 void Mesh::FinalizeTopology(
bool generate_bdr)
2902 bool generate_edges =
true;
2904 if (spaceDim == 0) { spaceDim = Dim; }
2905 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2909 if (tmp_vertex_parents.Size())
2911 MFEM_VERIFY(ncmesh == NULL,
"");
2912 ncmesh =
new NCMesh(
this);
2916 InitFromNCMesh(*ncmesh);
2917 ncmesh->OnMeshUpdated(
this);
2918 GenerateNCFaceInfo();
2922 tmp_vertex_parents.DeleteAll();
2932 GetElementToFaceTable();
2934 if (NumOfBdrElements == 0 && generate_bdr)
2936 GenerateBoundaryElements();
2937 GetElementToFaceTable();
2946 if (Dim > 1 && generate_edges)
2949 if (!el_to_edge) { el_to_edge =
new Table; }
2950 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2954 if (NumOfBdrElements == 0 && generate_bdr)
2956 GenerateBoundaryElements();
2973 ncmesh->OnMeshUpdated(
this);
2976 GenerateNCFaceInfo();
2983 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2985 if (NURBSext || ncmesh)
2987 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
2988 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
2997 const bool check_orientation =
true;
2998 const bool curved = (Nodes != NULL);
2999 const bool may_change_topology =
3000 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
3001 ( check_orientation && fix_orientation &&
3002 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
3005 Table *old_elem_vert = NULL;
3007 if (curved && may_change_topology)
3009 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3012 if (check_orientation)
3015 CheckElementOrientation(fix_orientation);
3019 MarkForRefinement();
3022 if (may_change_topology)
3026 DoNodeReorder(old_v_to_v, old_elem_vert);
3027 delete old_elem_vert;
3040 CheckBdrElementOrientation();
3045 if (Dim >= 2 && Dim == spaceDim)
3047 const int num_faces = GetNumFaces();
3048 for (
int i = 0; i < num_faces; i++)
3050 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
3051 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
3052 " Interior face with incompatible orientations.");
3059 double sx,
double sy,
double sz,
bool sfc_ordering)
3063 int NVert, NElem, NBdrElem;
3065 NVert = (nx+1) * (ny+1) * (nz+1);
3066 NElem = nx * ny * nz;
3067 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3068 if (type == Element::TETRAHEDRON)
3073 else if (type == Element::WEDGE)
3076 NBdrElem += 2*nx*ny;
3078 else if (type == Element::PYRAMID)
3081 NVert += nx * ny * nz;
3084 InitMesh(3, 3, NVert, NElem, NBdrElem);
3090 for (z = 0; z <= nz; z++)
3092 coord[2] = ((double) z / nz) * sz;
3093 for (y = 0; y <= ny; y++)
3095 coord[1] = ((double) y / ny) * sy;
3096 for (x = 0; x <= nx; x++)
3098 coord[0] = ((double) x / nx) * sx;
3103 if (type == Element::PYRAMID)
3105 for (z = 0; z < nz; z++)
3107 coord[2] = (((double) z + 0.5) / nz) * sz;
3108 for (y = 0; y < ny; y++)
3110 coord[1] = (((double) y + 0.5 ) / ny) * sy;
3111 for (x = 0; x < nx; x++)
3113 coord[0] = (((double) x + 0.5 ) / nx) * sx;
3120 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3121 #define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3124 if (sfc_ordering && type == Element::HEXAHEDRON)
3127 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
3128 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3130 for (
int k = 0; k < nx*ny*nz; k++)
3137 ind[0] = VTX(x , y , z );
3138 ind[1] = VTX(x+1, y , z );
3139 ind[2] = VTX(x+1, y+1, z );
3140 ind[3] = VTX(x , y+1, z );
3141 ind[4] = VTX(x , y , z+1);
3142 ind[5] = VTX(x+1, y , z+1);
3143 ind[6] = VTX(x+1, y+1, z+1);
3144 ind[7] = VTX(x , y+1, z+1);
3152 for (z = 0; z < nz; z++)
3154 for (y = 0; y < ny; y++)
3156 for (x = 0; x < nx; x++)
3159 ind[0] = VTX(x , y , z );
3160 ind[1] = VTX(x+1, y , z );
3161 ind[2] = VTX(x+1, y+1, z );
3162 ind[3] = VTX(x , y+1, z );
3163 ind[4] = VTX(x , y , z+1);
3164 ind[5] = VTX(x+1, y , z+1);
3165 ind[6] = VTX(x+1, y+1, z+1);
3166 ind[7] = VTX( x, y+1, z+1);
3168 if (type == Element::TETRAHEDRON)
3170 AddHexAsTets(ind, 1);
3172 else if (type == Element::WEDGE)
3174 AddHexAsWedges(ind, 1);
3176 else if (type == Element::PYRAMID)
3178 ind[8] = VTXP(x, y, z);
3179 AddHexAsPyramids(ind, 1);
3192 for (y = 0; y < ny; y++)
3194 for (x = 0; x < nx; x++)
3197 ind[0] = VTX(x , y , 0);
3198 ind[1] = VTX(x , y+1, 0);
3199 ind[2] = VTX(x+1, y+1, 0);
3200 ind[3] = VTX(x+1, y , 0);
3202 if (type == Element::TETRAHEDRON)
3204 AddBdrQuadAsTriangles(ind, 1);
3206 else if (type == Element::WEDGE)
3208 AddBdrQuadAsTriangles(ind, 1);
3217 for (y = 0; y < ny; y++)
3219 for (x = 0; x < nx; x++)
3222 ind[0] = VTX(x , y , nz);
3223 ind[1] = VTX(x+1, y , nz);
3224 ind[2] = VTX(x+1, y+1, nz);
3225 ind[3] = VTX(x , y+1, nz);
3227 if (type == Element::TETRAHEDRON)
3229 AddBdrQuadAsTriangles(ind, 6);
3231 else if (type == Element::WEDGE)
3233 AddBdrQuadAsTriangles(ind, 6);
3242 for (z = 0; z < nz; z++)
3244 for (y = 0; y < ny; y++)
3247 ind[0] = VTX(0 , y , z );
3248 ind[1] = VTX(0 , y , z+1);
3249 ind[2] = VTX(0 , y+1, z+1);
3250 ind[3] = VTX(0 , y+1, z );
3252 if (type == Element::TETRAHEDRON)
3254 AddBdrQuadAsTriangles(ind, 5);
3263 for (z = 0; z < nz; z++)
3265 for (y = 0; y < ny; y++)
3268 ind[0] = VTX(nx, y , z );
3269 ind[1] = VTX(nx, y+1, z );
3270 ind[2] = VTX(nx, y+1, z+1);
3271 ind[3] = VTX(nx, y , z+1);
3273 if (type == Element::TETRAHEDRON)
3275 AddBdrQuadAsTriangles(ind, 3);
3284 for (x = 0; x < nx; x++)
3286 for (z = 0; z < nz; z++)
3289 ind[0] = VTX(x , 0, z );
3290 ind[1] = VTX(x+1, 0, z );
3291 ind[2] = VTX(x+1, 0, z+1);
3292 ind[3] = VTX(x , 0, z+1);
3294 if (type == Element::TETRAHEDRON)
3296 AddBdrQuadAsTriangles(ind, 2);
3305 for (x = 0; x < nx; x++)
3307 for (z = 0; z < nz; z++)
3310 ind[0] = VTX(x , ny, z );
3311 ind[1] = VTX(x , ny, z+1);
3312 ind[2] = VTX(x+1, ny, z+1);
3313 ind[3] = VTX(x+1, ny, z );
3315 if (type == Element::TETRAHEDRON)
3317 AddBdrQuadAsTriangles(ind, 4);
3329 ofstream test_stream(
"debug.mesh");
3331 test_stream.close();
3340 double sx,
double sy,
3341 bool generate_edges,
bool sfc_ordering)
3350 if (type == Element::QUADRILATERAL)
3352 NumOfVertices = (nx+1) * (ny+1);
3353 NumOfElements = nx * ny;
3354 NumOfBdrElements = 2 * nx + 2 * ny;
3356 vertices.SetSize(NumOfVertices);
3357 elements.SetSize(NumOfElements);
3358 boundary.SetSize(NumOfBdrElements);
3365 for (j = 0; j < ny+1; j++)
3367 cy = ((double) j / ny) * sy;
3368 for (i = 0; i < nx+1; i++)
3370 cx = ((double) i / nx) * sx;
3371 vertices[k](0) = cx;
3372 vertices[k](1) = cy;
3381 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
3382 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
3384 for (k = 0; k < nx*ny; k++)
3388 ind[0] = i + j*(nx+1);
3389 ind[1] = i + 1 +j*(nx+1);
3390 ind[2] = i + 1 + (j+1)*(nx+1);
3391 ind[3] = i + (j+1)*(nx+1);
3398 for (j = 0; j < ny; j++)
3400 for (i = 0; i < nx; i++)
3402 ind[0] = i + j*(nx+1);
3403 ind[1] = i + 1 +j*(nx+1);
3404 ind[2] = i + 1 + (j+1)*(nx+1);
3405 ind[3] = i + (j+1)*(nx+1);
3414 for (i = 0; i < nx; i++)
3416 boundary[i] =
new Segment(i, i+1, 1);
3417 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3420 for (j = 0; j < ny; j++)
3422 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3423 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3427 else if (type == Element::TRIANGLE)
3429 NumOfVertices = (nx+1) * (ny+1);
3430 NumOfElements = 2 * nx * ny;
3431 NumOfBdrElements = 2 * nx + 2 * ny;
3433 vertices.SetSize(NumOfVertices);
3434 elements.SetSize(NumOfElements);
3435 boundary.SetSize(NumOfBdrElements);
3442 for (j = 0; j < ny+1; j++)
3444 cy = ((double) j / ny) * sy;
3445 for (i = 0; i < nx+1; i++)
3447 cx = ((double) i / nx) * sx;
3448 vertices[k](0) = cx;
3449 vertices[k](1) = cy;
3456 for (j = 0; j < ny; j++)
3458 for (i = 0; i < nx; i++)
3460 ind[0] = i + j*(nx+1);
3461 ind[1] = i + 1 + (j+1)*(nx+1);
3462 ind[2] = i + (j+1)*(nx+1);
3465 ind[1] = i + 1 + j*(nx+1);
3466 ind[2] = i + 1 + (j+1)*(nx+1);
3474 for (i = 0; i < nx; i++)
3476 boundary[i] =
new Segment(i, i+1, 1);
3477 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3480 for (j = 0; j < ny; j++)
3482 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3483 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3490 MFEM_ABORT(
"Unsupported element type.");
3494 CheckElementOrientation();
3496 if (generate_edges == 1)
3498 el_to_edge =
new Table;
3499 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3501 CheckBdrElementOrientation();
3510 attributes.Append(1);
3511 bdr_attributes.Append(1); bdr_attributes.Append(2);
3512 bdr_attributes.Append(3); bdr_attributes.Append(4);
3517 void Mesh::Make1D(
int n,
double sx)
3526 NumOfVertices = n + 1;
3528 NumOfBdrElements = 2;
3529 vertices.SetSize(NumOfVertices);
3530 elements.SetSize(NumOfElements);
3531 boundary.SetSize(NumOfBdrElements);
3534 for (j = 0; j < n+1; j++)
3536 vertices[j](0) = ((
double) j / n) * sx;
3540 for (j = 0; j < n; j++)
3542 elements[j] =
new Segment(j, j+1, 1);
3547 boundary[0] =
new Point(ind, 1);
3549 boundary[1] =
new Point(ind, 2);
3557 attributes.Append(1);
3558 bdr_attributes.Append(1); bdr_attributes.Append(2);
3561 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
3579 last_operation = Mesh::NONE;
3582 elements.SetSize(NumOfElements);
3583 for (
int i = 0; i < NumOfElements; i++)
3585 elements[i] = mesh.
elements[i]->Duplicate(
this);
3592 boundary.SetSize(NumOfBdrElements);
3593 for (
int i = 0; i < NumOfBdrElements; i++)
3595 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
3614 faces.SetSize(mesh.
faces.Size());
3615 for (
int i = 0; i < faces.Size(); i++)
3618 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
3652 if (dynamic_cast<const ParMesh*>(&mesh))
3664 if (mesh.
Nodes && copy_nodes)
3669 FiniteElementCollection::New(fec->
Name());
3673 Nodes->MakeOwner(fec_copy);
3674 *Nodes = *mesh.
Nodes;
3696 bool fix_orientation)
3700 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
3701 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
3715 double sx,
double sy,
bool sfc_ordering)
3718 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
3725 double sx,
double sy,
double sz,
bool sfc_ordering)
3728 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
3737 ref_factors = ref_factor;
3751 bool fix_orientation)
3760 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
3764 Load(imesh, generate_edges, refine, fix_orientation);
3769 bool fix_orientation)
3772 Load(input, generate_edges, refine, fix_orientation);
3781 "Not enough vertices in external array : "
3782 "len_vertex_data = "<< len_vertex_data <<
", "
3785 if (vertex_data == (
double *)(
vertices.GetData()))
3787 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
3792 memcpy(vertex_data,
vertices.GetData(),
3801 int *element_attributes,
int num_elements,
3803 int *boundary_attributes,
int num_boundary_elements,
3806 if (space_dimension == -1)
3811 InitMesh(dimension, space_dimension, 0, num_elements,
3812 num_boundary_elements);
3815 int boundary_index_stride = num_boundary_elements > 0 ?
3819 vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
3822 for (
int i = 0; i < num_elements; i++)
3825 elements[i]->SetVertices(element_indices + i * element_index_stride);
3826 elements[i]->SetAttribute(element_attributes[i]);
3830 for (
int i = 0; i < num_boundary_elements; i++)
3833 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3834 boundary[i]->SetAttribute(boundary_attributes[i]);
3850 #ifdef MFEM_USE_MEMALLOC
3859 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3872 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3875 for (
int i = 0; i < nv; i++)
3888 for (
int j = 0; j < nv; j++)
3960 MFEM_ABORT(
"invalid element type: " << type);
3967 std::string parse_tag)
3969 int curved = 0, read_gf = 1;
3970 bool finalize_topo =
true;
3974 MFEM_ABORT(
"Input stream is not open");
3981 getline(input, mesh_type);
3985 int mfem_version = 0;
3986 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
3987 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
3991 int mfem_nc_version = 0;
3992 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
3993 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
4001 if (mfem_version == 12 && parse_tag.empty())
4003 parse_tag =
"mfem_mesh_end";
4007 else if (mfem_nc_version)
4009 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4016 MFEM_VERIFY(mfem_nc_version >= 10,
4017 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
4018 "used to load a parallel nonconforming mesh, sorry.");
4021 input, mfem_nc_version, curved, is_nc);
4026 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
4039 else if (mesh_type ==
"linemesh")
4043 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
4045 if (mesh_type ==
"curved_areamesh2")
4051 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
4055 else if (mesh_type ==
"TrueGrid")
4059 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
4061 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
4063 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4064 "Unsupported VTK format");
4065 ReadVTKMesh(input, curved, read_gf, finalize_topo);
4067 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
4071 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
4075 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
4080 else if (mesh_type ==
"$MeshFormat")
4085 ((mesh_type.size() > 2 &&
4086 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
4087 (mesh_type.size() > 3 &&
4088 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
4093 #ifdef MFEM_USE_NETCDF
4096 MFEM_ABORT(
"NetCDF support requires configuration with"
4097 " MFEM_USE_NETCDF=YES");
4103 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
4104 " Use mfem::named_ifgzstream for input.");
4110 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
4136 bool generate_bdr =
false;
4141 if (curved && read_gf)
4155 if (mfem_version == 12)
4161 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
4162 getline(input, line);
4168 if (line ==
"mfem_mesh_end") {
break; }
4170 while (line != parse_tag);
4172 else if (mfem_nc_version >= 10)
4177 MFEM_VERIFY(ident ==
"mfem_mesh_end",
4178 "invalid mesh: end of file tag not found");
4186 int i, j, ie, ib, iv, *v, nv;
4217 for (i = 0; i < num_pieces; i++)
4224 for (i = 0; i < num_pieces; i++)
4230 for (j = 0; j < m->
GetNE(); j++)
4235 for (j = 0; j < m->
GetNBE(); j++)
4240 for (
int k = 0; k < nv; k++)
4242 v[k] = lvert_vert[v[k]];
4247 for (j = 0; j < m->
GetNV(); j++)
4259 for (i = 0; i < num_pieces; i++)
4270 for (i = 0; i < num_pieces; i++)
4274 for (j = 0; j < m->
GetNE(); j++)
4279 for (
int k = 0; k < nv; k++)
4286 for (j = 0; j < m->
GetNBE(); j++)
4291 for (
int k = 0; k < nv; k++)
4298 for (j = 0; j < m->
GetNV(); j++)
4312 for (i = 0; i < num_pieces; i++)
4314 gf_array[i] = mesh_array[i]->
GetNodes();
4329 ref_factors = ref_factor;
4340 int orig_ne = orig_mesh.
GetNE();
4341 MFEM_VERIFY(ref_factors.
Size() == orig_ne && orig_ne > 0,
4342 "Number of refinement factors must equal number of elements")
4343 MFEM_VERIFY(ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
4346 "Invalid refinement type. Must use closed basis type.");
4348 int min_ref = ref_factors.
Min();
4349 int max_ref = ref_factors.
Max();
4351 bool var_order = (min_ref != max_ref);
4364 for (
int i = 0; i < orig_ne; i++)
4382 for (
int el = 0; el < orig_ne; el++)
4394 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
4395 for (
int i = 0; i < phys_pts.
Width(); i++)
4404 for (
int k = 0; k < nvert; k++)
4407 v[k] = rdofs[c2h_map[cid]];
4420 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
4431 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
4437 for (
int k = 0; k < nvert; k++)
4440 v[k] = rdofs[c2h_map[cid]];
4462 for (
int iel = 0; iel < orig_ne; iel++)
4471 const int *node_map = NULL;
4474 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
4475 const int *vertex_map = vertex_fec.
GetDofMap(geom);
4476 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
4477 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
4480 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
4483 int iv = vertex_map[iv_lex];
4485 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
4487 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
4490 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
4501 using GeomRef = std::pair<Geometry::Type, int>;
4502 std::map<GeomRef, int> point_matrices_offsets;
4504 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4508 GeomRef id(geom, ref_factors[el_coarse]);
4509 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
4515 point_matrices_offsets[id] = n_point_matrices[geom];
4516 n_point_matrices[geom] += nref_el;
4523 int nmatrices = n_point_matrices[geom];
4530 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4533 int ref = ref_factors[el_coarse];
4534 int offset = point_matrices_offsets[GeomRef(geom, ref)];
4540 for (
int k = 0; k < nvert; k++)
4572 "Mesh::MakeSimplicial requires a properly oriented input mesh");
4574 "Mesh::MakeSimplicial does not support non-conforming meshes.")
4581 Mesh copy(orig_mesh);
4586 int nv = orig_mesh.
GetNV();
4587 int ne = orig_mesh.
GetNE();
4588 int nbe = orig_mesh.
GetNBE();
4601 int new_ne = 0, new_nbe = 0;
4602 for (
int i=0; i<ne; ++i)
4606 for (
int i=0; i<nbe; ++i)
4615 for (
int i=0; i<nv; ++i)
4625 if (vglobal == NULL)
4628 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
4629 vglobal = vglobal_id.
GetData();
4632 constexpr
int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
4633 constexpr
int quad_ntris = 2, prism_ntets = 3;
4634 static const int quad_trimap[2][nv_tri*quad_ntris] =
4646 static const int prism_rot[nv_prism*nv_prism] =
4655 static const int prism_f[nv_quad] = {1, 2, 5, 4};
4656 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
4670 static const int hex_rot[nv_hex*nv_hex] =
4672 0, 1, 2, 3, 4, 5, 6, 7,
4673 1, 0, 4, 5, 2, 3, 7, 6,
4674 2, 1, 5, 6, 3, 0, 4, 7,
4675 3, 0, 1, 2, 7, 4, 5, 6,
4676 4, 0, 3, 7, 5, 1, 2, 6,
4677 5, 1, 0, 4, 6, 2, 3, 7,
4678 6, 2, 1, 5, 7, 3, 0, 4,
4679 7, 3, 2, 6, 4, 0, 1, 5
4681 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
4682 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
4683 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
4684 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
4685 static const int hex_tetmap0[nv_tet*5] =
4692 static const int hex_tetmap1[nv_tet*6] =
4699 static const int hex_tetmap2[nv_tet*6] =
4706 static const int hex_tetmap3[nv_tet*6] =
4713 static const int *hex_tetmaps[4] =
4715 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
4718 auto find_min = [](
const int*
a,
int n) {
return std::min_element(a,a+n)-
a; };
4720 for (
int i=0; i<ne; ++i)
4722 const int *v = orig_mesh.
elements[i]->GetVertices();
4726 if (num_subdivisions[orig_geom] == 1)
4738 for (
int itri=0; itri<quad_ntris; ++itri)
4743 for (
int iv=0; iv<nv_tri; ++iv)
4745 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
4753 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
4756 int irot = find_min(vg, nv_prism);
4757 for (
int iv=0; iv<nv_prism; ++iv)
4759 int jv = prism_rot[iv + irot*nv_prism];
4764 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
4765 int j = find_min(q, nv_quad);
4766 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
4767 for (
int itet=0; itet<prism_ntets; ++itet)
4772 for (
int iv=0; iv<nv_tet; ++iv)
4774 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
4782 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
4786 int irot = find_min(vg, nv_hex);
4787 for (
int iv=0; iv<nv_hex; ++iv)
4789 int jv = hex_rot[iv + irot*nv_hex];
4799 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
4800 j = find_min(q, nv_quad);
4801 if (j == 0 || j == 2) { bitmask += 4; }
4803 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
4804 j = find_min(q, nv_quad);
4805 if (j == 1 || j == 3) { bitmask += 2; }
4807 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
4808 j = find_min(q, nv_quad);
4809 if (j == 0 || j == 2) { bitmask += 1; }
4812 int nrot = num_rot[bitmask];
4813 for (
int k=0; k<nrot; ++k)
4827 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
4828 int ntets = (ndiags == 0) ? 5 : 6;
4829 const int *tetmap = hex_tetmaps[ndiags];
4830 for (
int itet=0; itet<ntets; ++itet)
4835 for (
int iv=0; iv<nv_tet; ++iv)
4837 v2[iv] = vg[tetmap[itet + iv*ntets]];
4846 for (
int i=0; i<nbe; ++i)
4848 const int *v = orig_mesh.
boundary[i]->GetVertices();
4851 if (num_subdivisions[orig_geom] == 1)
4861 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
4863 int iv_min = find_min(vg, nv_quad);
4864 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
4865 for (
int itri=0; itri<quad_ntris; ++itri)
4870 for (
int iv=0; iv<nv_tri; ++iv)
4872 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
4879 MFEM_ABORT(
"Unreachable");
4893 Mesh periodic_mesh(orig_mesh,
true);
4899 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
4904 for (
int j = 0; j < nv; j++)
4910 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
4915 for (
int j = 0; j < nv; j++)
4922 return periodic_mesh;
4926 const std::vector<Vector> &translations,
double tol)
const
4930 Vector coord(sdim), at(sdim), dx(sdim);
4931 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
4932 xMax = xMin = xDiff = 0.0;
4936 for (
int be = 0; be <
GetNBE(); be++)
4941 for (
int i = 0; i < dofs.
Size(); i++)
4943 bdr_v.insert(dofs[i]);
4946 for (
int j = 0; j < sdim; j++)
4948 xMax[j] = max(xMax[j], coord[j]);
4949 xMin[j] = min(xMin[j], coord[j]);
4953 add(xMax, -1.0, xMin, xDiff);
4954 double dia = xDiff.
Norml2();
4962 map<int, int> replica2primary;
4964 map<int, set<int>> primary2replicas;
4968 for (
int v : bdr_v) { primary2replicas[v]; }
4972 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
4974 if (r ==
p) {
return; }
4975 primary2replicas[
p].insert(r);
4976 replica2primary[r] =
p;
4977 for (
int s : primary2replicas[r])
4979 primary2replicas[
p].insert(
s);
4980 replica2primary[
s] =
p;
4982 primary2replicas.erase(r);
4985 for (
unsigned int i = 0; i < translations.size(); i++)
4987 for (
int vi : bdr_v)
4990 add(coord, translations[i], at);
4992 for (
int vj : bdr_v)
4995 add(at, -1.0, coord, dx);
4996 if (dx.
Norml2() > dia*tol) {
continue; }
5001 bool pi = primary2replicas.find(vi) != primary2replicas.end();
5002 bool pj = primary2replicas.find(vj) != primary2replicas.end();
5008 make_replica(vj, vi);
5013 int owner_of_vj = replica2primary[vj];
5015 make_replica(vi, owner_of_vj);
5021 int owner_of_vi = replica2primary[vi];
5022 make_replica(vj, owner_of_vi);
5029 int owner_of_vi = replica2primary[vi];
5030 int owner_of_vj = replica2primary[vj];
5031 make_replica(owner_of_vj, owner_of_vi);
5038 std::vector<int> v2v(
GetNV());
5039 for (
size_t i = 0; i < v2v.size(); i++)
5043 for (
auto &&r2p : replica2primary)
5045 v2v[r2p.first] = r2p.second;
5054 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5059 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5076 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5081 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5111 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
5135 for (
int i = 0; i <
elements.Size(); i++)
5145 for (
int i = 0; i <
boundary.Size(); i++)
5162 for (
int i = 0; i < vd; i++)
5216 boundary.SetSize(NumOfBdrElements);
5227 edge_to_knot.
SetSize(NumOfEdges);
5231 input >> edge_to_knot[j] >> v[0] >> v[1];
5234 edge_to_knot[j] = -1 - edge_to_knot[j];
5252 for (
int d = 0; d < v.
Size(); d++)
5260 for (d = 0; d < p.
Size(); d++)
5264 for ( ; d < v.
Size(); d++)
5296 if (dynamic_cast<const H1_FECollection*>(fec)
5297 || dynamic_cast<const L2_FECollection*>(fec))
5306 #ifndef MFEM_USE_MPI
5307 const bool warn =
true;
5310 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
5314 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
5315 " H1-continuous mesh!\n "
5316 "If this is the desired behavior, you can silence"
5317 " this warning by converting\n "
5318 "the NURBS mesh to high-order mesh in advance by"
5319 " calling the method\n "
5320 "Mesh::SetCurvature().");
5345 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
5364 MFEM_ASSERT(nodes != NULL,
"");
5380 case 1:
return GetNV();
5416 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5417 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
5422 int i, j, k, wo = 0, fo = 0;
5431 int *vi =
elements[i]->GetVertices();
5434 for (j = 0; j < 3; j++)
5438 for (j = 0; j < 2; j++)
5439 for (k = 0; k < 2; k++)
5441 J(j, k) = v[j+1][k] - v[0][k];
5462 MFEM_ABORT(
"Invalid 2D element type \""
5479 int *vi =
elements[i]->GetVertices();
5485 for (j = 0; j < 4; j++)
5489 for (j = 0; j < 3; j++)
5490 for (k = 0; k < 3; k++)
5492 J(j, k) = v[j+1][k] - v[0][k];
5551 MFEM_ABORT(
"Invalid 3D element type \""
5557 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5560 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
5561 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
5576 if (test[0] == base[0])
5577 if (test[1] == base[1])
5585 else if (test[0] == base[1])
5586 if (test[1] == base[0])
5595 if (test[1] == base[0])
5606 for (
int j = 0; j < 3; j++)
5607 if (test[aor[j]] != base[j])
5620 for (i = 0; i < 4; i++)
5621 if (test[i] == base[0])
5628 if (test[(i+1)%4] == base[1])
5637 for (
int j = 0; j < 4; j++)
5638 if (test[aor[j]] != base[j])
5640 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
5642 for (
int k = 0; k < 4; k++)
5647 for (
int k = 0; k < 4; k++)
5656 if (test[(i+1)%4] == base[1])
5672 if (test[0] == base[0])
5673 if (test[1] == base[1])
5674 if (test[2] == base[2])
5682 else if (test[2] == base[1])
5683 if (test[3] == base[2])
5692 if (test[1] == base[2])
5700 else if (test[1] == base[0])
5701 if (test[2] == base[1])
5702 if (test[0] == base[2])
5710 else if (test[3] == base[1])
5711 if (test[2] == base[2])
5720 if (test[3] == base[2])
5728 else if (test[2] == base[0])
5729 if (test[3] == base[1])
5730 if (test[0] == base[2])
5738 else if (test[0] == base[1])
5739 if (test[1] == base[2])
5748 if (test[3] == base[2])
5757 if (test[0] == base[1])
5758 if (test[2] == base[2])
5766 else if (test[1] == base[1])
5767 if (test[0] == base[2])
5776 if (test[1] == base[2])
5787 for (
int j = 0; j < 4; j++)
5788 if (test[aor[j]] != base[j])
5813 int *bv =
boundary[i]->GetVertices();
5819 mfem::Swap<int>(bv[0], bv[1]);
5833 if (
faces_info[fi].Elem2No >= 0) {
continue; }
5836 int *bv =
boundary[i]->GetVertices();
5838 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
5839 const int *fv =
faces[fi]->GetVertices();
5856 MFEM_ABORT(
"Invalid 2D boundary element type \""
5857 << bdr_type <<
"\"");
5862 if (orientation % 2 == 0) {
continue; }
5864 if (!fix_it) {
continue; }
5872 mfem::Swap<int>(bv[0], bv[1]);
5876 mfem::Swap<int>(be[1], be[2]);
5882 mfem::Swap<int>(bv[0], bv[2]);
5886 mfem::Swap<int>(be[0], be[1]);
5887 mfem::Swap<int>(be[2], be[3]);
5900 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
5910 MFEM_ASSERT(0 <= dim && dim <=
Dim,
"invalid dim: " << dim);
5921 MFEM_ASSERT(0 <= dim && dim <=
Dim,
"invalid dim: " << dim);
5940 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
5941 "is not generated.");
5944 const int *v =
elements[i]->GetVertices();
5945 const int ne =
elements[i]->GetNEdges();
5947 for (
int j = 0; j < ne; j++)
5949 const int *e =
elements[i]->GetEdgeVertices(j);
5950 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5961 const int *v =
boundary[i]->GetVertices();
5962 cor[0] = (v[0] < v[1]) ? (1) : (-1);
5975 const int *v =
boundary[i]->GetVertices();
5976 const int ne =
boundary[i]->GetNEdges();
5978 for (
int j = 0; j < ne; j++)
5980 const int *e =
boundary[i]->GetEdgeVertices(j);
5981 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5993 const int *v =
faces[i]->GetVertices();
5994 o[0] = (v[0] < v[1]) ? (1) : (-1);
6006 const int *v =
faces[i]->GetVertices();
6007 const int ne =
faces[i]->GetNEdges();
6009 for (
int j = 0; j < ne; j++)
6011 const int *e =
faces[i]->GetEdgeVertices(j);
6012 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6040 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
6091 for (j = 0; j < nv; j++)
6103 for (j = 0; j < nv; j++)
6150 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
6154 int n = el_faces.
Size();
6157 for (
int j = 0; j < n; j++)
6161 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
6165 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
6166 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
6190 MFEM_ABORT(
"invalid geometry");
6198 case 1:
return boundary[i]->GetVertices()[0];
6201 default: MFEM_ABORT(
"invalid dimension!");
6211 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6214 const int *bv =
boundary[bdr_el]->GetVertices();
6222 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6233 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6236 const int *bv =
boundary[bdr_el]->GetVertices();
6244 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6271 for (j = 0; j < nv; j++)
6273 pointmat(k, j) =
vertices[v[j]](k);
6288 for (j = 0; j < nv; j++)
6290 pointmat(k, j) =
vertices[v[j]](k);
6302 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
6305 return sqrt(length);
6313 for (
int i = 0; i < elem_array.
Size(); i++)
6318 for (
int i = 0; i < elem_array.
Size(); i++)
6320 const int *v = elem_array[i]->GetVertices();
6321 const int ne = elem_array[i]->GetNEdges();
6322 for (
int j = 0; j < ne; j++)
6324 const int *e = elem_array[i]->GetEdgeVertices(j);
6338 v_to_v.
Push(v[0], v[1]);
6345 const int *v =
elements[i]->GetVertices();
6346 const int ne =
elements[i]->GetNEdges();
6347 for (
int j = 0; j < ne; j++)
6349 const int *e =
elements[i]->GetEdgeVertices(j);
6350 v_to_v.
Push(v[e[0]], v[e[1]]);
6358 int i, NumberOfEdges;
6374 const int *v =
boundary[i]->GetVertices();
6375 be_to_f[i] = v_to_v(v[0], v[1]);
6388 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
6392 return NumberOfEdges;
6483 if (
faces[gf] == NULL)
6493 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6494 "Interior edge found between 2D elements "
6496 <<
" and " << el <<
".");
6497 int *v =
faces[gf]->GetVertices();
6499 if ( v[1] == v0 && v[0] == v1 )
6503 else if ( v[0] == v0 && v[1] == v1 )
6513 MFEM_ABORT(
"internal error");
6519 int v0,
int v1,
int v2)
6521 if (
faces[gf] == NULL)
6531 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6532 "Interior triangular face found connecting elements "
6534 <<
" and " << el <<
".");
6535 int orientation, vv[3] = { v0, v1, v2 };
6542 faces_info[gf].Elem2Inf = 64 * lf + orientation;
6547 int v0,
int v1,
int v2,
int v3)
6559 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6560 "Interior quadrilateral face found connecting elements "
6562 <<
" and " << el <<
".");
6563 int vv[4] = { v0, v1, v2, v3 };
6577 for (i = 0; i <
faces.Size(); i++)
6583 faces.SetSize(nfaces);
6585 for (i = 0; i < nfaces; i++)
6593 const int *v =
elements[i]->GetVertices();
6603 const int ne =
elements[i]->GetNEdges();
6604 for (
int j = 0; j < ne; j++)
6606 const int *e =
elements[i]->GetEdgeVertices(j);
6617 for (
int j = 0; j < 4; j++)
6621 v[fv[0]], v[fv[1]], v[fv[2]]);
6627 for (
int j = 0; j < 2; j++)
6631 v[fv[0]], v[fv[1]], v[fv[2]]);
6633 for (
int j = 2; j < 5; j++)
6637 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6643 for (
int j = 0; j < 1; j++)
6647 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6649 for (
int j = 1; j < 5; j++)
6653 v[fv[0]], v[fv[1]], v[fv[2]]);
6659 for (
int j = 0; j < 6; j++)
6663 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6668 MFEM_ABORT(
"Unexpected type of Element.");
6676 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
6692 for (
int i = 0; i < list.
masters.Size(); i++)
6695 if (master.
index >= nfaces) {
continue; }
6701 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
6702 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
6706 for (
int i = 0; i < list.
slaves.Size(); i++)
6710 if (slave.
index < 0 ||
6711 slave.
index >= nfaces ||
6740 const int *v =
elements[i]->GetVertices();
6745 for (
int j = 0; j < 4; j++)
6748 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6754 for (
int j = 0; j < 1; j++)
6757 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6759 for (
int j = 1; j < 5; j++)
6762 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6768 for (
int j = 0; j < 2; j++)
6771 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6773 for (
int j = 2; j < 5; j++)
6776 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6784 for (
int j = 0; j < 6; j++)
6787 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6792 MFEM_ABORT(
"Unexpected type of Element.");
6816 for (
int j = 0; j < 4; j++)
6820 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6826 for (
int j = 0; j < 2; j++)
6830 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6832 for (
int j = 2; j < 5; j++)
6836 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6842 for (
int j = 0; j < 1; j++)
6846 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6848 for (
int j = 1; j < 5; j++)
6852 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6860 for (
int j = 0; j < 6; j++)
6864 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6869 MFEM_ABORT(
"Unexpected type of Element.");
6882 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
6887 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
6891 MFEM_ABORT(
"Unexpected type of boundary Element.");
6905 void Rotate3(
int &
a,
int &
b,
int &c)
6937 Table *old_elem_vert = NULL;
6948 int *v =
elements[i]->GetVertices();
6950 Rotate3(v[0], v[1], v[2]);
6953 Rotate3(v[1], v[2], v[3]);
6966 int *v =
boundary[i]->GetVertices();
6968 Rotate3(v[0], v[1], v[2]);
6984 delete old_elem_vert;
7000 if (p[i] < pmin[i]) { pmin[i] = p[i]; }
7001 if (p[i] > pmax[i]) { pmax[i] = p[i]; }
7015 for (
int i =
spaceDim-1; i >= 0; i--)
7017 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
7018 if (idx < 0) { idx = 0; }
7019 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
7020 part = part * nxyz[i] + idx;
7022 partitioning[el] = part;
7025 return partitioning;
7035 #ifdef MFEM_USE_METIS
7037 int print_messages = 1;
7040 int init_flag, fin_flag;
7041 MPI_Initialized(&init_flag);
7042 MPI_Finalized(&fin_flag);
7043 if (init_flag && !fin_flag)
7047 if (rank != 0) { print_messages = 0; }
7051 int i, *partitioning;
7061 partitioning[i] = 0;
7068 partitioning[i] = i;
7074 #ifndef MFEM_USE_METIS_5
7086 bool freedata =
false;
7088 idx_t *mpartitioning;
7091 if (
sizeof(
idx_t) ==
sizeof(int))
7095 mpartitioning = (
idx_t*) partitioning;
7104 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
7105 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
7106 mpartitioning =
new idx_t[n];
7109 #ifndef MFEM_USE_METIS_5
7112 METIS_SetDefaultOptions(options);
7113 options[METIS_OPTION_CONTIG] = 1;
7121 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
7126 if (part_method >= 0 && part_method <= 2)
7128 for (i = 0; i < n; i++)
7134 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
7140 if (part_method == 0 || part_method == 3)
7142 #ifndef MFEM_USE_METIS_5
7171 " error in METIS_PartGraphRecursive!");
7178 if (part_method == 1 || part_method == 4)
7180 #ifndef MFEM_USE_METIS_5
7209 " error in METIS_PartGraphKway!");
7216 if (part_method == 2 || part_method == 5)
7218 #ifndef MFEM_USE_METIS_5
7231 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
7248 " error in METIS_PartGraphKway!");
7256 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
7260 nparts = (int) mparts;
7261 if (mpartitioning != (
idx_t*)partitioning)
7265 partitioning[k] = mpartitioning[k];
7272 delete[] mpartitioning;
7288 auto count_partition_elements = [&]()
7290 for (i = 0; i < nparts; i++)
7298 psize[partitioning[i]].one++;
7302 for (i = 0; i < nparts; i++)
7304 if (psize[i].one == 0) { empty_parts++; }
7308 count_partition_elements();
7316 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
7317 << empty_parts <<
" empty parts!"
7318 <<
" Applying a simple fix ..." << endl;
7321 SortPairs<int,int>(psize, nparts);
7323 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7330 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7332 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
7338 partitioning[j] = psize[nparts-1-i].two;
7345 count_partition_elements();
7349 return partitioning;
7353 mfem_error(
"Mesh::GeneratePartitioning(...): "
7354 "MFEM was compiled without Metis.");
7368 int num_elem, *i_elem_elem, *j_elem_elem;
7370 num_elem = elem_elem.
Size();
7371 i_elem_elem = elem_elem.
GetI();
7372 j_elem_elem = elem_elem.
GetJ();
7377 int stack_p, stack_top_p, elem;
7381 for (i = 0; i < num_elem; i++)
7383 if (partitioning[i] > num_part)
7385 num_part = partitioning[i];
7392 for (i = 0; i < num_part; i++)
7399 for (elem = 0; elem < num_elem; elem++)
7401 if (component[elem] >= 0)
7406 component[elem] = num_comp[partitioning[elem]]++;
7408 elem_stack[stack_top_p++] = elem;
7410 for ( ; stack_p < stack_top_p; stack_p++)
7412 i = elem_stack[stack_p];
7413 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
7416 if (partitioning[k] == partitioning[i])
7418 if (component[k] < 0)
7420 component[k] = component[i];
7421 elem_stack[stack_top_p++] = k;
7423 else if (component[k] != component[i])
7435 int i, n_empty, n_mcomp;
7443 n_empty = n_mcomp = 0;
7444 for (i = 0; i < num_comp.
Size(); i++)
7445 if (num_comp[i] == 0)
7449 else if (num_comp[i] > 1)
7456 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
7457 <<
"The following subdomains are empty :\n";
7458 for (i = 0; i < num_comp.
Size(); i++)
7459 if (num_comp[i] == 0)
7467 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
7468 <<
"The following subdomains are NOT connected :\n";
7469 for (i = 0; i < num_comp.
Size(); i++)
7470 if (num_comp[i] > 1)
7476 if (n_empty == 0 && n_mcomp == 0)
7477 mfem::out <<
"Mesh::CheckPartitioning(...) : "
7478 "All subdomains are connected." << endl;
7492 const double *a = A.
Data();
7493 const double *b = B.
Data();
7502 c(0) = a[0]*a[3]-a[1]*a[2];
7503 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
7504 c(2) = b[0]*b[3]-b[1]*b[2];
7525 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
7526 a[1] * (a[5] * a[6] - a[3] * a[8]) +
7527 a[2] * (a[3] * a[7] - a[4] * a[6]));
7529 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
7530 b[1] * (a[5] * a[6] - a[3] * a[8]) +
7531 b[2] * (a[3] * a[7] - a[4] * a[6]) +
7533 a[0] * (b[4] * a[8] - b[5] * a[7]) +
7534 a[1] * (b[5] * a[6] - b[3] * a[8]) +
7535 a[2] * (b[3] * a[7] - b[4] * a[6]) +
7537 a[0] * (a[4] * b[8] - a[5] * b[7]) +
7538 a[1] * (a[5] * b[6] - a[3] * b[8]) +
7539 a[2] * (a[3] * b[7] - a[4] * b[6]));
7541 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
7542 a[1] * (b[5] * b[6] - b[3] * b[8]) +
7543 a[2] * (b[3] * b[7] - b[4] * b[6]) +
7545 b[0] * (a[4] * b[8] - a[5] * b[7]) +
7546 b[1] * (a[5] * b[6] - a[3] * b[8]) +
7547 b[2] * (a[3] * b[7] - a[4] * b[6]) +
7549 b[0] * (b[4] * a[8] - b[5] * a[7]) +
7550 b[1] * (b[5] * a[6] - b[3] * a[8]) +
7551 b[2] * (b[3] * a[7] - b[4] * a[6]));
7553 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
7554 b[1] * (b[5] * b[6] - b[3] * b[8]) +
7555 b[2] * (b[3] * b[7] - b[4] * b[6]));
7601 double a = z(2), b = z(1), c = z(0);
7602 double D = b*b-4*a*c;
7609 x(0) = x(1) = -0.5 * b /
a;
7614 x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
7622 t = -0.5 * (b + sqrt(D));
7626 t = -0.5 * (b - sqrt(D));
7632 Swap<double>(x(0), x(1));
7640 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
7643 double Q = (a * a - 3 *
b) / 9;
7644 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
7645 double Q3 = Q * Q * Q;
7652 x(0) = x(1) = x(2) = - a / 3;
7656 double sqrtQ = sqrt(Q);
7660 x(0) = -2 * sqrtQ - a / 3;
7661 x(1) = x(2) = sqrtQ - a / 3;
7665 x(0) = x(1) = - sqrtQ - a / 3;
7666 x(2) = 2 * sqrtQ - a / 3;
7673 double theta = acos(R / sqrt(Q3));
7674 double A = -2 * sqrt(Q);
7676 x0 = A * cos(theta / 3) - a / 3;
7677 x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
7678 x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
7683 Swap<double>(x0, x1);
7687 Swap<double>(x1, x2);
7690 Swap<double>(x0, x1);
7703 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
7707 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
7709 x(0) = A + Q / A - a / 3;
7718 const double factor,
const int Dim)
7720 const double c0 = c(0);
7721 c(0) = c0 * (1.0 - pow(factor, -Dim));
7723 for (
int j = 0; j < nr; j++)
7735 c(0) = c0 * (1.0 - pow(factor, Dim));
7737 for (
int j = 0; j < nr; j++)
7756 const double factor = 2.0;
7771 for (
int k = 0; k < nv; k++)
7774 V(j, k) = displacements(v[k]+j*nvs);
7804 for (
int j = 0; j < nv; j++)
7830 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7833 vertices[i](j) += displacements(j*nv+i);
7841 for (
int i = 0; i < nv; i++)
7844 vert_coord(j*nv+i) =
vertices[i](j);
7850 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7853 vertices[i](j) = vert_coord(j*nv+i);
7864 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
7883 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
7900 (*Nodes) += displacements;
7912 node_coord = (*Nodes);
7924 (*Nodes) = node_coord;
7959 mfem::Swap<GridFunction*>(
Nodes, nodes);
7979 for (j = 1; j < n; j++)
8016 int quad_counter = 0;
8031 vertices.SetSize(oelem + quad_counter);
8032 new_elements.
SetSize(4 * NumOfElements);
8038 const int attr =
elements[i]->GetAttribute();
8039 int *v =
elements[i]->GetVertices();
8045 for (
int ei = 0; ei < 3; ei++)
8047 for (
int k = 0; k < 2; k++)
8055 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8057 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8059 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8061 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8065 const int qe = quad_counter;
8069 for (
int ei = 0; ei < 4; ei++)
8071 for (
int k = 0; k < 2; k++)
8079 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
8081 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
8083 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
8085 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
8089 MFEM_ABORT(
"unknown element type: " << el_type);
8099 const int attr =
boundary[i]->GetAttribute();
8100 int *v =
boundary[i]->GetVertices();
8109 static const double A = 0.0, B = 0.5, C = 1.0;
8110 static double tri_children[2*3*4] =
8117 static double quad_children[2*4*4] =
8131 for (
int i = 0; i <
elements.Size(); i++)
8152 if (!
Nodes || update_nodes)
8160 static inline double sqr(
const double &x)
8182 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
8185 int NumOfQuadFaces = 0;
8191 for (
int i = 0; i <
faces.Size(); i++)
8195 f2qf[i] = NumOfQuadFaces;
8202 NumOfQuadFaces =
faces.Size();
8206 int hex_counter = 0;
8209 for (
int i = 0; i <
elements.Size(); i++)
8218 int pyr_counter = 0;
8221 for (
int i = 0; i <
elements.Size(); i++)
8239 DSTable *v_to_v_ptr = v_to_v_p;
8255 std::sort(row_start, J_v2v.
end());
8258 for (
int i = 0; i < J_v2v.
Size(); i++)
8260 e2v[J_v2v[i].two] = i;
8273 it.SetIndex(e2v[it.Index()]);
8283 const int oelem = oface + NumOfQuadFaces;
8288 vertices.SetSize(oelem + hex_counter);
8296 const int attr =
elements[i]->GetAttribute();
8297 int *v =
elements[i]->GetVertices();
8304 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8312 for (
int ei = 0; ei < 6; ei++)
8314 for (
int k = 0; k < 2; k++)
8324 const int rt_algo = 1;
8335 double len_sqr, min_len;
8337 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
8338 sqr(J(1,0)-J(1,1)-J(1,2)) +
8339 sqr(J(2,0)-J(2,1)-J(2,2));
8342 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
8343 sqr(J(1,1)-J(1,0)-J(1,2)) +
8344 sqr(J(2,1)-J(2,0)-J(2,2));
8345 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
8347 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
8348 sqr(J(1,2)-J(1,0)-J(1,1)) +
8349 sqr(J(2,2)-J(2,0)-J(2,1));
8350 if (len_sqr < min_len) { rt = 2; }
8355 double Em_data[18], Js_data[9], Jp_data[9];
8358 double ar1, ar2,
kappa, kappa_min;
8360 for (
int s = 0;
s < 3;
s++)
8362 for (
int t = 0;
t < 3;
t++)
8364 Em(
t,
s) = 0.5*J(
t,
s);
8367 for (
int t = 0;
t < 3;
t++)
8369 Em(
t,3) = 0.5*(J(
t,0)+J(
t,1));
8370 Em(
t,4) = 0.5*(J(
t,0)+J(
t,2));
8371 Em(
t,5) = 0.5*(J(
t,1)+J(
t,2));
8375 for (
int t = 0;
t < 3;
t++)
8377 Js(
t,0) = Em(
t,5)-Em(
t,0);
8378 Js(
t,1) = Em(
t,1)-Em(
t,0);
8379 Js(
t,2) = Em(
t,2)-Em(
t,0);
8383 for (
int t = 0;
t < 3;
t++)
8385 Js(
t,0) = Em(
t,5)-Em(
t,0);
8386 Js(
t,1) = Em(
t,2)-Em(
t,0);
8387 Js(
t,2) = Em(
t,4)-Em(
t,0);
8391 kappa_min = std::max(ar1, ar2);
8395 for (
int t = 0;
t < 3;
t++)
8397 Js(
t,0) = Em(
t,0)-Em(
t,1);
8398 Js(
t,1) = Em(
t,4)-Em(
t,1);
8399 Js(
t,2) = Em(
t,2)-Em(
t,1);
8403 for (
int t = 0;
t < 3;
t++)
8405 Js(
t,0) = Em(
t,2)-Em(
t,1);
8406 Js(
t,1) = Em(
t,4)-Em(
t,1);
8407 Js(
t,2) = Em(
t,5)-Em(
t,1);
8411 kappa = std::max(ar1, ar2);
8412 if (kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
8415 for (
int t = 0;
t < 3;
t++)
8417 Js(
t,0) = Em(
t,0)-Em(
t,2);
8418 Js(
t,1) = Em(
t,1)-Em(
t,2);
8419 Js(
t,2) = Em(
t,3)-Em(
t,2);
8423 for (
int t = 0;
t < 3;
t++)
8425 Js(
t,0) = Em(
t,1)-Em(
t,2);
8426 Js(
t,1) = Em(
t,5)-Em(
t,2);
8427 Js(
t,2) = Em(
t,3)-Em(
t,2);
8431 kappa = std::max(ar1, ar2);
8432 if (kappa < kappa_min) { rt = 2; }
8435 static const int mv_all[3][4][4] =
8437 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
8438 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
8439 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
8441 const int (&mv)[4][4] = mv_all[rt];
8443 #ifndef MFEM_USE_MEMALLOC
8445 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8447 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8449 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8451 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8453 for (
int k = 0; k < 4; k++)
8455 new_elements[j+4+k] =
8456 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8457 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8461 new_elements[j+0] = tet =
TetMemory.Alloc();
8462 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8464 new_elements[j+1] = tet =
TetMemory.Alloc();
8465 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8467 new_elements[j+2] = tet =
TetMemory.Alloc();
8468 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8470 new_elements[j+3] = tet =
TetMemory.Alloc();
8471 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8473 for (
int k = 0; k < 4; k++)
8475 new_elements[j+4+k] = tet =
TetMemory.Alloc();
8476 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8477 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8480 for (
int k = 0; k < 4; k++)
8485 for (
int k = 0; k < 4; k++)
8499 for (
int fi = 2; fi < 5; fi++)
8501 for (
int k = 0; k < 4; k++)
8508 for (
int ei = 0; ei < 9; ei++)
8510 for (
int k = 0; k < 2; k++)
8517 const int qf2 = f2qf[f[2]];
8518 const int qf3 = f2qf[f[3]];
8519 const int qf4 = f2qf[f[4]];
8522 new Wedge(v[0], oedge+e[0], oedge+e[2],
8523 oedge+e[6], oface+qf2, oface+qf4, attr);
8526 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
8527 oface+qf3, oface+qf4, oface+qf2, attr);
8530 new Wedge(oedge+e[0], v[1], oedge+e[1],
8531 oface+qf2, oedge+e[7], oface+qf3, attr);
8534 new Wedge(oedge+e[2], oedge+e[1], v[2],
8535 oface+qf4, oface+qf3, oedge+e[8], attr);
8538 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
8539 v[3], oedge+e[3], oedge+e[5], attr);
8542 new Wedge(oface+qf3, oface+qf4, oface+qf2,
8543 oedge+e[4], oedge+e[5], oedge+e[3], attr);
8546 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
8547 oedge+e[3], v[4], oedge+e[4], attr);
8550 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
8551 oedge+e[5], oedge+e[4], v[5], attr);
8560 for (
int fi = 0; fi < 1; fi++)
8562 for (
int k = 0; k < 4; k++)
8569 for (
int ei = 0; ei < 8; ei++)
8571 for (
int k = 0; k < 2; k++)
8578 const int qf0 = f2qf[f[0]];
8581 new Pyramid(v[0], oedge+e[0], oface+qf0,
8582 oedge+e[3], oedge+e[4], attr);
8585 new Pyramid(oedge+e[0], v[1], oedge+e[1],
8586 oface+qf0, oedge+e[5], attr);
8589 new Pyramid(oface+qf0, oedge+e[1], v[2],
8590 oedge+e[2], oedge+e[6], attr);
8593 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
8594 v[3], oedge+e[7], attr);
8597 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
8598 oedge+e[7], v[4], attr);
8601 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
8602 oedge+e[4], oface+qf0, attr);
8605 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
8609 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
8613 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
8617 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
8625 const int he = hex_counter;
8630 if (f2qf.
Size() == 0)
8636 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[f[k]]; }
8642 for (
int fi = 0; fi < 6; fi++)
8644 for (
int k = 0; k < 4; k++)
8651 for (
int ei = 0; ei < 12; ei++)
8653 for (
int k = 0; k < 2; k++)
8661 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
8662 oedge+e[3], oedge+e[8], oface+qf[1],
8663 oelem+he, oface+qf[4], attr);
8666 oface+qf[0], oface+qf[1], oedge+e[9],
8667 oface+qf[2], oelem+he, attr);
8669 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
8670 oedge+e[2], oelem+he, oface+qf[2],
8671 oedge+e[10], oface+qf[3], attr);
8673 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
8674 v[3], oface+qf[4], oelem+he,
8675 oface+qf[3], oedge+e[11], attr);
8677 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
8678 oface+qf[4], v[4], oedge+e[4],
8679 oface+qf[5], oedge+e[7], attr);
8681 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
8682 oelem+he, oedge+e[4], v[5],
8683 oedge+e[5], oface+qf[5], attr);
8685 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
8686 oface+qf[3], oface+qf[5], oedge+e[5],
8687 v[6], oedge+e[6], attr);
8689 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
8690 oedge+e[11], oedge+e[7], oface+qf[5],
8691 oedge+e[6], v[7], attr);
8696 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
8708 const int attr =
boundary[i]->GetAttribute();
8709 int *v =
boundary[i]->GetVertices();
8716 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8723 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8725 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8727 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8729 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8737 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
8739 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
8741 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
8743 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
8747 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
8753 static const double A = 0.0, B = 0.5, C = 1.0, D = -1.0;
8754 static double tet_children[3*4*16] =
8756 A,A,A, B,A,A, A,B,A, A,A,B,
8757 B,A,A, C,A,A, B,B,A, B,A,B,
8758 A,B,A, B,B,A, A,C,A, A,B,B,
8759 A,A,B, B,A,B, A,B,B, A,A,C,
8764 B,A,A, A,B,B, A,B,A, A,A,B,
8765 B,A,A, A,B,B, A,A,B, B,A,B,
8766 B,A,A, A,B,B, B,A,B, B,B,A,
8767 B,A,A, A,B,B, B,B,A, A,B,A,
8769 A,B,A, B,A,A, B,A,B, A,A,B,
8770 A,B,A, A,A,B, B,A,B, A,B,B,
8771 A,B,A, A,B,B, B,A,B, B,B,A,
8772 A,B,A, B,B,A, B,A,B, B,A,A,
8774 A,A,B, B,A,A, A,B,A, B,B,A,
8775 A,A,B, A,B,A, A,B,B, B,B,A,
8776 A,A,B, A,B,B, B,A,B, B,B,A,
8777 A,A,B, B,A,B, B,A,A, B,B,A
8779 static double pyr_children[3*5*10] =
8781 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
8782 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
8783 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
8784 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
8785 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
8786 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
8787 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
8788 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
8789 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
8790 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
8792 static double pri_children[3*6*8] =
8794 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
8795 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
8796 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
8797 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
8798 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
8799 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
8800 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
8801 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
8803 static double hex_children[3*8*8] =
8805 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
8806 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
8807 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
8808 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
8809 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
8810 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
8811 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
8812 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
8824 for (
int i = 0; i <
elements.Size(); i++)
8835 NumOfElements = 8 * NumOfElements + 2 * pyr_counter;
8855 int i, j, ind, nedges;
8862 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
8876 for (j = 0; j < marked_el.
Size(); j++)
8881 int new_v = cnv + j, new_e = cne + j;
8890 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
8892 UseExternalData(seg_children, 1, 2, 3);
8905 int *edge1 =
new int[nedges];
8906 int *edge2 =
new int[nedges];
8907 int *middle =
new int[nedges];
8909 for (i = 0; i < nedges; i++)
8911 edge1[i] = edge2[i] = middle[i] = -1;
8917 for (j = 1; j < v.
Size(); j++)
8919 ind = v_to_v(v[j-1], v[j]);
8920 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8922 ind = v_to_v(v[0], v[v.
Size()-1]);
8923 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8927 for (i = 0; i < marked_el.
Size(); i++)
8933 int need_refinement;
8936 need_refinement = 0;
8937 for (i = 0; i < nedges; i++)
8939 if (middle[i] != -1 && edge1[i] != -1)
8941 need_refinement = 1;
8946 while (need_refinement == 1);
8949 int v1[2], v2[2], bisect, temp;
8951 for (i = 0; i < temp; i++)
8954 bisect = v_to_v(v[0], v[1]);
8955 if (middle[bisect] != -1)
8959 v1[0] = v[0]; v1[1] = middle[bisect];
8960 v2[0] = middle[bisect]; v2[1] = v[1];
8966 mfem_error(
"Only bisection of segment is implemented"
8990 MFEM_VERIFY(
GetNE() == 0 ||
8992 "tetrahedral mesh is not marked for refinement:"
8993 " call Finalize(true)");
9000 for (i = 0; i < marked_el.
Size(); i++)
9006 for (i = 0; i < marked_el.
Size(); i++)
9015 for (i = 0; i < marked_el.
Size(); i++)
9032 int need_refinement;
9037 need_refinement = 0;
9045 if (
elements[i]->NeedRefinement(v_to_v))
9047 need_refinement = 1;
9052 while (need_refinement == 1);
9059 need_refinement = 0;
9061 if (
boundary[i]->NeedRefinement(v_to_v))
9063 need_refinement = 1;
9067 while (need_refinement == 1);
9098 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
9099 "not supported. Project the NURBS to Nodes first.");
9109 if (!refinements.
Size())
9130 Swap(*mesh2,
false);
9146 const int *fine,
int nfine,
int op)
9149 for (
int i = 0; i < nfine; i++)
9151 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
9153 double err_fine = elem_error[fine[i]];
9156 case 0: error = std::min(error, err_fine);
break;
9157 case 1: error += err_fine;
break;
9158 case 2: error = std::max(error, err_fine);
break;
9165 double threshold,
int nc_limit,
int op)
9167 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
9168 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
9169 "Project the NURBS to Nodes first.");
9182 for (
int i = 0; i < dt.
Size(); i++)
9184 if (nc_limit > 0 && !level_ok[i]) {
continue; }
9189 if (error < threshold) { derefs.
Append(i); }
9192 if (!derefs.
Size()) {
return false; }
9199 Swap(*mesh2,
false);
9213 int nc_limit,
int op)
9223 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
9230 int nc_limit,
int op)
9233 for (
int i = 0; i < tmp.Size(); i++)
9235 tmp[i] = elem_error(i);
9320 #ifdef MFEM_USE_MEMALLOC
9347 for (
int i = 0; i < elem_array.
Size(); i++)
9349 if (elem_array[i]->GetGeometryType() == geom)
9354 elem_vtx.
SetSize(nv*num_elems);
9358 for (
int i = 0; i < elem_array.
Size(); i++)
9364 elem_vtx.
Append(loc_vtx);
9372 for (
int i = 0; i < nelem; i++) { list[i] = i; }
9388 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
9400 default: MFEM_ABORT(
"internal error");
9406 int nonconforming,
int nc_limit)
9416 else if (nonconforming < 0)
9437 for (
int i = 0; i < refinements.
Size(); i++)
9439 el_to_refine[i] = refinements[i].index;
9443 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
9444 if (rt == 1 || rt == 2 || rt == 4)
9448 else if (rt == 3 || rt == 5 || rt == 6)
9466 for (
int i = 0; i < el_to_refine.
Size(); i++)
9468 refinements[i] =
Refinement(el_to_refine[i]);
9475 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
9476 "Please project the NURBS to Nodes first, with SetCurvature().");
9479 MFEM_VERIFY(
ncmesh != NULL || dynamic_cast<const ParMesh*>(
this) == NULL,
9480 "Sorry, converting a conforming ParMesh to an NC mesh is "
9488 (simplices_nonconforming && (
meshgen & 0x1)) )
9501 for (
int i = 0; i <
GetNE(); i++)
9503 if ((
double) rand() / RAND_MAX <
prob)
9508 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
9520 for (
int i = 0; i <
GetNE(); i++)
9523 bool refine =
false;
9524 for (
int j = 0; j < v.
Size(); j++)
9529 double d = vert(l) -
vertices[v[j]](l);
9532 if (dist <= eps*eps) { refine =
true;
break; }
9543 int nonconforming,
int nc_limit)
9545 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
9547 for (
int i = 0; i <
GetNE(); i++)
9549 if (elem_error[i] > threshold)
9563 int nonconforming,
int nc_limit)
9567 return RefineByError(tmp, threshold, nonconforming, nc_limit);
9572 int *edge1,
int *edge2,
int *middle)
9575 int v[2][4], v_new, bisect,
t;
9587 bisect = v_to_v(vert[0], vert[1]);
9588 MFEM_ASSERT(bisect >= 0,
"");
9590 if (middle[bisect] == -1)
9601 if (edge1[bisect] == i)
9603 edge1[bisect] = edge2[bisect];
9606 middle[bisect] = v_new;
9610 v_new = middle[bisect];
9618 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9619 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9640 bisect = v_to_v(v[1][0], v[1][1]);
9641 MFEM_ASSERT(bisect >= 0,
"");
9643 if (edge1[bisect] == i)
9647 else if (edge2[bisect] == i)
9656 MFEM_ABORT(
"Bisection for now works only for triangles.");
9663 int v[2][4], v_new, bisect,
t;
9670 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
9674 "TETRAHEDRON element is not marked for refinement.");
9679 bisect = v_to_v.
FindId(vert[0], vert[1]);
9683 for (j = 0; j < 3; j++)
9700 new_redges[0][0] = 2;
9701 new_redges[0][1] = 1;
9702 new_redges[1][0] = 2;
9703 new_redges[1][1] = 1;
9704 int tr1 = -1, tr2 = -1;
9705 switch (old_redges[0])
9708 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
9713 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
9717 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
9720 switch (old_redges[1])
9723 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
9728 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
9732 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
9739 #ifdef MFEM_USE_MEMALLOC
9775 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
9782 int v[2][3], v_new, bisect,
t;
9793 bisect = v_to_v.
FindId(vert[0], vert[1]);
9794 MFEM_ASSERT(bisect >= 0,
"");
9796 MFEM_ASSERT(v_new != -1,
"");
9800 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9801 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9811 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
9817 int *edge1,
int *edge2,
int *middle)
9820 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
9829 bisect[0] = v_to_v(v[0],v[1]);
9830 bisect[1] = v_to_v(v[1],v[2]);
9831 bisect[2] = v_to_v(v[0],v[2]);
9832 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
9834 for (j = 0; j < 3; j++)
9836 if (middle[bisect[j]] == -1)
9847 if (edge1[bisect[j]] == i)
9849 edge1[bisect[j]] = edge2[bisect[j]];
9852 middle[bisect[j]] = v_new[j];
9856 v_new[j] = middle[bisect[j]];
9859 edge1[bisect[j]] = -1;
9865 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
9866 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
9867 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
9868 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
9883 tri2->ResetTransform(code);
9884 tri3->ResetTransform(code);
9888 tri2->PushTransform(1);
9889 tri3->PushTransform(2);
9902 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
9938 for (
int i = 0; i < elem_geoms.
Size(); i++)
9946 std::map<unsigned, int> mat_no;
9950 for (
int j = 0; j <
elements.Size(); j++)
9953 unsigned code =
elements[j]->GetTransform();
9956 int &matrix = mat_no[code];
9957 if (!matrix) { matrix = mat_no.size(); }
9967 std::map<unsigned, int>::iterator it;
9968 for (it = mat_no.begin(); it != mat_no.end(); ++it)
9982 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
9983 " geom = " << geom);
9993 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
10002 os <<
"areamesh2\n\n";
10006 os <<
"curved_areamesh2\n\n";
10015 os <<
boundary[i]->GetAttribute();
10016 for (j = 0; j < v.
Size(); j++)
10018 os <<
' ' << v[j] + 1;
10030 for (j = 0; j < v.
Size(); j++)
10032 os <<
' ' << v[j] + 1;
10044 for (j = 1; j <
Dim; j++)
10061 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
10069 os <<
"NETGEN_Neutral_Format\n";
10074 for (j = 0; j <
Dim; j++)
10087 os <<
elements[i]->GetAttribute();
10088 for (j = 0; j < nv; j++)
10090 os <<
' ' << ind[j]+1;
10101 os <<
boundary[i]->GetAttribute();
10102 for (j = 0; j < nv; j++)
10104 os <<
' ' << ind[j]+1;
10116 <<
" 0 0 0 0 0 0 0\n"
10117 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
10119 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
10120 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
10124 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
10130 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
10131 for (j = 0; j < nv; j++)
10133 os <<
' ' << ind[j]+1;
10142 os <<
boundary[i]->GetAttribute();
10143 for (j = 0; j < nv; j++)
10145 os <<
' ' << ind[j]+1;
10147 os <<
" 1.0 1.0 1.0 1.0\n";
10180 os <<
"\n# mesh curvature GridFunction";
10185 os <<
"\nmfem_mesh_end" << endl;
10190 os << (section_delimiter.empty()
10191 ?
"MFEM mesh v1.0\n" :
"MFEM mesh v1.2\n");
10195 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10200 "# TETRAHEDRON = 4\n"
10206 os <<
"\ndimension\n" <<
Dim;
10241 if (!section_delimiter.empty())
10243 os << section_delimiter << endl;
10252 os <<
"MFEM NURBS mesh v1.0\n";
10256 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10262 os <<
"\ndimension\n" <<
Dim
10279 int ki = e_to_k[i];
10284 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
10291 ofstream ofs(fname);
10292 ofs.precision(precision);
10296 #ifdef MFEM_USE_ADIOS2
10306 "# vtk DataFile Version 3.0\n"
10307 "Generated by MFEM\n"
10309 "DATASET UNSTRUCTURED_GRID\n";
10322 for ( ; j < 3; j++)
10338 os << (*Nodes)(vdofs[0]);
10342 os <<
' ' << (*Nodes)(vdofs[j]);
10344 for ( ; j < 3; j++)
10358 size +=
elements[i]->GetNVertices() + 1;
10360 os <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
10363 const int *v =
elements[i]->GetVertices();
10364 const int nv =
elements[i]->GetNVertices();
10368 for (
int j = 0; j < nv; j++)
10370 os <<
' ' << v[perm ? perm[j] : j];
10383 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
10384 "Point meshes should have a single dof per element");
10385 size += dofs.
Size() + 1;
10387 os <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
10390 if (!strcmp(fec_name,
"Linear") ||
10391 !strcmp(fec_name,
"H1_0D_P1") ||
10392 !strcmp(fec_name,
"H1_1D_P1") ||
10393 !strcmp(fec_name,
"H1_2D_P1") ||
10394 !strcmp(fec_name,
"H1_3D_P1"))
10398 else if (!strcmp(fec_name,
"Quadratic") ||
10399 !strcmp(fec_name,
"H1_1D_P2") ||
10400 !strcmp(fec_name,
"H1_2D_P2") ||
10401 !strcmp(fec_name,
"H1_3D_P2"))
10407 mfem::err <<
"Mesh::PrintVTK : can not save '"
10408 << fec_name <<
"' elements!" << endl;
10417 for (
int j = 0; j < dofs.
Size(); j++)
10419 os <<
' ' << dofs[j];
10422 else if (order == 2)
10424 const int *vtk_mfem;
10425 switch (
elements[i]->GetGeometryType())
10439 for (
int j = 0; j < dofs.
Size(); j++)
10441 os <<
' ' << dofs[vtk_mfem[j]];
10451 int vtk_cell_type = 5;
10455 os << vtk_cell_type <<
'\n';
10459 os <<
"CELL_DATA " << NumOfElements <<
'\n'
10460 <<
"SCALARS material int\n"
10461 <<
"LOOKUP_TABLE default\n";
10464 os <<
elements[i]->GetAttribute() <<
'\n';
10471 bool high_order_output,
10472 int compression_level,
10475 int ref = (high_order_output &&
Nodes)
10478 fname = fname +
".vtu";
10480 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
10481 if (compression_level != 0)
10483 os <<
" compressor=\"vtkZLibDataCompressor\"";
10486 os <<
"<UnstructuredGrid>\n";
10487 PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
10488 os <<
"</Piece>\n";
10489 os <<
"</UnstructuredGrid>\n";
10490 os <<
"</VTKFile>" << std::endl;
10497 bool high_order_output,
10498 int compression_level)
10500 PrintVTU(fname, format, high_order_output, compression_level,
true);
10504 bool high_order_output,
int compression_level,
10512 std::vector<char> buf;
10514 auto get_geom = [&](
int i)
10522 int np = 0, nc_ref = 0;
10523 for (
int i = 0; i < ne; i++)
10532 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
10533 << (high_order_output ? ne : nc_ref) <<
"\">\n";
10536 os <<
"<Points>\n";
10537 os <<
"<DataArray type=\"" << type_str
10538 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
10539 for (
int i = 0; i < ne; i++)
10552 for (
int j = 0; j < pmat.
Width(); j++)
10578 os <<
"</DataArray>" << std::endl;
10579 os <<
"</Points>" << std::endl;
10581 os <<
"<Cells>" << std::endl;
10582 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
10583 << fmt_str <<
"\">" << std::endl;
10585 std::vector<int> offset;
10588 if (high_order_output)
10591 for (
int iel = 0; iel < ne; iel++)
10595 int nnodes = local_connectivity.
Size();
10596 for (
int i=0; i<nnodes; ++i)
10603 offset.push_back(np);
10609 for (
int i = 0; i < ne; i++)
10615 for (
int j = 0; j < RG.
Size(); )
10618 offset.push_back(coff);
10620 for (
int k = 0; k < nv; k++, j++)
10634 os <<
"</DataArray>" << std::endl;
10636 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
10637 << fmt_str <<
"\">" << std::endl;
10639 for (
size_t ii=0; ii<offset.size(); ii++)
10647 os <<
"</DataArray>" << std::endl;
10648 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
10649 << fmt_str <<
"\">" << std::endl;
10651 const int *vtk_geom_map =
10653 for (
int i = 0; i < ne; i++)
10656 uint8_t vtk_cell_type = 5;
10658 vtk_cell_type = vtk_geom_map[geom];
10660 if (high_order_output)
10669 for (
int j = 0; j < RG.
Size(); j += nv)
10679 os <<
"</DataArray>" << std::endl;
10680 os <<
"</Cells>" << std::endl;
10682 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
10683 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
10684 << fmt_str <<
"\">" << std::endl;
10685 for (
int i = 0; i < ne; i++)
10688 if (high_order_output)
10707 os <<
"</DataArray>" << std::endl;
10708 os <<
"</CellData>" << std::endl;
10719 "# vtk DataFile Version 3.0\n"
10720 "Generated by MFEM\n"
10722 "DATASET UNSTRUCTURED_GRID\n";
10727 os <<
"FIELD FieldData 1\n"
10737 np = nc = size = 0;
10738 for (
int i = 0; i <
GetNE(); i++)
10747 os <<
"POINTS " << np <<
" double\n";
10749 for (
int i = 0; i <
GetNE(); i++)
10756 for (
int j = 0; j < pmat.
Width(); j++)
10758 os << pmat(0, j) <<
' ';
10761 os << pmat(1, j) <<
' ';
10773 os << 0.0 <<
' ' << 0.0;
10780 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
10782 for (
int i = 0; i <
GetNE(); i++)
10789 for (
int j = 0; j < RG.
Size(); )
10792 for (
int k = 0; k < nv; k++, j++)
10794 os <<
' ' << np + RG[j];
10800 os <<
"CELL_TYPES " << nc <<
'\n';
10801 for (
int i = 0; i <
GetNE(); i++)
10809 for (
int j = 0; j < RG.
Size(); j += nv)
10811 os << vtk_cell_type <<
'\n';
10815 os <<
"CELL_DATA " << nc <<
'\n'
10816 <<
"SCALARS material int\n"
10817 <<
"LOOKUP_TABLE default\n";
10818 for (
int i = 0; i <
GetNE(); i++)
10826 os << attr <<
'\n';
10833 srand((
unsigned)time(0));
10834 double a = double(rand()) / (double(RAND_MAX) + 1.);
10835 int el0 = (int)floor(a *
GetNE());
10837 os <<
"SCALARS element_coloring int\n"
10838 <<
"LOOKUP_TABLE default\n";
10839 for (
int i = 0; i <
GetNE(); i++)
10846 os << coloring[i] + 1 <<
'\n';
10852 os <<
"POINT_DATA " << np <<
'\n' << flush;
10857 int delete_el_to_el = (
el_to_el) ? (0) : (1);
10859 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
10862 const int *i_el_el = el_el.
GetI();
10863 const int *j_el_el = el_el.
GetJ();
10868 stack_p = stack_top_p = 0;
10869 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
10871 if (colors[el] != -2)
10877 el_stack[stack_top_p++] = el;
10879 for ( ; stack_p < stack_top_p; stack_p++)
10881 int i = el_stack[stack_p];
10882 int num_nb = i_el_el[i+1] - i_el_el[i];
10883 if (max_num_col < num_nb + 1)
10885 max_num_col = num_nb + 1;
10887 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10889 int k = j_el_el[j];
10890 if (colors[k] == -2)
10893 el_stack[stack_top_p++] = k;
10901 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
10903 int i = el_stack[stack_p], col;
10905 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10907 col = colors[j_el_el[j]];
10910 col_marker[col] = 1;
10914 for (col = 0; col < max_num_col; col++)
10915 if (col_marker[col] == 0)
10923 if (delete_el_to_el)
10931 int elem_attr)
const
10933 if (
Dim != 3 &&
Dim != 2) {
return; }
10935 int i, j, k, l, nv, nbe, *v;
10937 os <<
"MFEM mesh v1.0\n";
10941 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10946 "# TETRAHEDRON = 4\n"
10951 os <<
"\ndimension\n" <<
Dim
10956 <<
' ' <<
elements[i]->GetGeometryType();
10959 for (j = 0; j < nv; j++)
10971 l = partitioning[l];
10986 os <<
"\nboundary\n" << nbe <<
'\n';
10992 l = partitioning[l];
10995 nv =
faces[i]->GetNVertices();
10996 v =
faces[i]->GetVertices();
10997 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
10998 for (j = 0; j < nv; j++)
11005 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
11006 for (j = nv-1; j >= 0; j--)
11017 nv =
faces[i]->GetNVertices();
11018 v =
faces[i]->GetVertices();
11019 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
11020 for (j = 0; j < nv; j++)
11051 int interior_faces)
11053 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
11054 if (
Dim != 3 &&
Dim != 2) {
return; }
11063 int nv =
elements[i]->GetNVertices();
11064 const int *ind =
elements[i]->GetVertices();
11065 for (
int j = 0; j < nv; j++)
11071 int *voff =
new int[NumOfVertices+1];
11075 voff[i] = vcount[i-1] + voff[i-1];
11081 vown[i] =
new int[vcount[i]];
11093 int nv =
elements[i]->GetNVertices();
11094 const int *ind =
elements[i]->GetVertices();
11095 for (
int j = 0; j < nv; j++)
11098 vown[ind[j]][vcount[ind[j]]] = i;
11104 vcount[i] = voff[i+1] - voff[i];
11108 for (
int i = 0; i < edge_el.
Size(); i++)
11110 const int *el = edge_el.
GetRow(i);
11113 int k = partitioning[el[0]];
11114 int l = partitioning[el[1]];
11115 if (interior_faces || k != l)
11127 os <<
"areamesh2\n\n" << nbe <<
'\n';
11129 for (
int i = 0; i < edge_el.
Size(); i++)
11131 const int *el = edge_el.
GetRow(i);
11134 int k = partitioning[el[0]];
11135 int l = partitioning[el[1]];
11136 if (interior_faces || k != l)
11141 for (
int j = 0; j < 2; j++)
11142 for (
int s = 0;
s < vcount[ev[j]];
s++)
11143 if (vown[ev[j]][
s] == el[0])
11145 os <<
' ' << voff[ev[j]]+
s+1;
11149 for (
int j = 1; j >= 0; j--)
11150 for (
int s = 0;
s < vcount[ev[j]];
s++)
11151 if (vown[ev[j]][
s] == el[1])
11153 os <<
' ' << voff[ev[j]]+
s+1;
11160 int k = partitioning[el[0]];
11164 for (
int j = 0; j < 2; j++)
11165 for (
int s = 0;
s < vcount[ev[j]];
s++)
11166 if (vown[ev[j]][
s] == el[0])
11168 os <<
' ' << voff[ev[j]]+
s+1;
11175 os << NumOfElements <<
'\n';
11178 int nv =
elements[i]->GetNVertices();
11179 const int *ind =
elements[i]->GetVertices();
11180 os << partitioning[i]+1 <<
' ';
11182 for (
int j = 0; j < nv; j++)
11184 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11185 vown[ind[j]][vcount[ind[j]]] = i;
11192 vcount[i] = voff[i+1] - voff[i];
11198 for (
int k = 0; k < vcount[i]; k++)
11200 for (
int j = 0; j <
Dim; j++)
11210 os <<
"NETGEN_Neutral_Format\n";
11214 for (
int k = 0; k < vcount[i]; k++)
11216 for (
int j = 0; j <
Dim; j++)
11224 os << NumOfElements <<
'\n';
11227 int nv =
elements[i]->GetNVertices();
11228 const int *ind =
elements[i]->GetVertices();
11229 os << partitioning[i]+1;
11230 for (
int j = 0; j < nv; j++)
11232 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11233 vown[ind[j]][vcount[ind[j]]] = i;
11240 vcount[i] = voff[i+1] - voff[i];
11250 int k = partitioning[
faces_info[i].Elem1No];
11251 l = partitioning[l];
11252 if (interior_faces || k != l)
11269 int k = partitioning[
faces_info[i].Elem1No];
11270 l = partitioning[l];
11271 if (interior_faces || k != l)
11273 int nv =
faces[i]->GetNVertices();
11274 const int *ind =
faces[i]->GetVertices();
11276 for (
int j = 0; j < nv; j++)
11277 for (
int s = 0;
s < vcount[ind[j]];
s++)
11278 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11280 os <<
' ' << voff[ind[j]]+
s+1;
11284 for (
int j = nv-1; j >= 0; j--)
11285 for (
int s = 0;
s < vcount[ind[j]];
s++)
11286 if (vown[ind[j]][
s] == faces_info[i].Elem2No)
11288 os <<
' ' << voff[ind[j]]+
s+1;
11295 int k = partitioning[
faces_info[i].Elem1No];
11296 int nv =
faces[i]->GetNVertices();
11297 const int *ind =
faces[i]->GetVertices();
11299 for (
int j = 0; j < nv; j++)
11300 for (
int s = 0;
s < vcount[ind[j]];
s++)
11301 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11303 os <<
' ' << voff[ind[j]]+
s+1;
11319 int k = partitioning[
faces_info[i].Elem1No];
11320 l = partitioning[l];
11321 if (interior_faces || k != l)
11334 <<
" 0 0 0 0 0 0 0\n"
11335 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
11336 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
11337 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
11338 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11341 for (
int k = 0; k < vcount[i]; k++)
11342 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
11347 int nv =
elements[i]->GetNVertices();
11348 const int *ind =
elements[i]->GetVertices();
11349 os << i+1 <<
' ' << partitioning[i]+1;
11350 for (
int j = 0; j < nv; j++)
11352 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11353 vown[ind[j]][vcount[ind[j]]] = i;
11360 vcount[i] = voff[i+1] - voff[i];
11369 int k = partitioning[
faces_info[i].Elem1No];
11370 l = partitioning[l];
11371 if (interior_faces || k != l)
11373 int nv =
faces[i]->GetNVertices();
11374 const int *ind =
faces[i]->GetVertices();
11376 for (
int j = 0; j < nv; j++)
11377 for (
int s = 0;
s < vcount[ind[j]];
s++)
11378 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11380 os <<
' ' << voff[ind[j]]+
s+1;
11382 os <<
" 1.0 1.0 1.0 1.0\n";
11384 for (
int j = nv-1; j >= 0; j--)
11385 for (
int s = 0;
s < vcount[ind[j]];
s++)
11386 if (vown[ind[j]][
s] == faces_info[i].Elem2No)
11388 os <<
' ' << voff[ind[j]]+
s+1;
11390 os <<
" 1.0 1.0 1.0 1.0\n";
11395 int k = partitioning[
faces_info[i].Elem1No];
11396 int nv =
faces[i]->GetNVertices();
11397 const int *ind =
faces[i]->GetVertices();
11399 for (
int j = 0; j < nv; j++)
11400 for (
int s = 0;
s < vcount[ind[j]];
s++)
11401 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11403 os <<
' ' << voff[ind[j]]+
s+1;
11405 os <<
" 1.0 1.0 1.0 1.0\n";
11429 " NURBS mesh is not supported!");
11433 os <<
"MFEM mesh v1.0\n";
11437 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
11442 "# TETRAHEDRON = 4\n"
11447 os <<
"\ndimension\n" <<
Dim
11455 const int *
const i_AF_f = Aface_face.
GetI();
11456 const int *
const j_AF_f = Aface_face.
GetJ();
11458 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
11459 for (
const int * iface = j_AF_f + i_AF_f[iAF];
11460 iface < j_AF_f + i_AF_f[iAF+1];
11463 os << iAF+1 <<
' ';
11495 double *cg =
new double[na*
spaceDim];
11496 int *nbea =
new int[na];
11503 for (i = 0; i < na; i++)
11507 cg[i*spaceDim+j] = 0.0;
11515 for (k = 0; k < vert.
Size(); k++)
11527 for (k = 0; k < vert.
Size(); k++)
11528 if (vn[vert[k]] == 1)
11533 cg[bea*spaceDim+j] += pointmat(j,k);
11544 for (k = 0; k < vert.
Size(); k++)
11549 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
11565 double *cg =
new double[na*
spaceDim];
11566 int *nbea =
new int[na];
11573 for (i = 0; i < na; i++)
11577 cg[i*spaceDim+j] = 0.0;
11585 for (k = 0; k < vert.
Size(); k++)
11597 for (k = 0; k < vert.
Size(); k++)
11598 if (vn[vert[k]] == 1)
11603 cg[bea*spaceDim+j] += pointmat(j,k);
11614 for (k = 0; k < vert.
Size(); k++)
11619 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
11635 for (
int i = 0; i <
vertices.Size(); i++)
11637 for (
int j = 0; j <
spaceDim; j++)
11649 xnew.ProjectCoefficient(f_pert);
11657 "incompatible vector dimensions");
11665 for (
int d = 0; d <
spaceDim; d++)
11667 vertices[i](d) = xnew(d + spaceDim*i);
11684 for (
int i = 0; i <
GetNE(); i++)
11689 for (
int j = 0; j < nv; j++)
11694 for (
int i = 0; i <
GetNBE(); i++)
11699 for (
int j = 0; j < nv; j++)
11705 for (
int i = 0; i < v2v.
Size(); i++)
11710 v2v[i] = num_vert++;
11714 if (num_vert == v2v.
Size()) {
return; }
11716 Vector nodes_by_element;
11721 for (
int i = 0; i <
GetNE(); i++)
11728 for (
int i = 0; i <
GetNE(); i++)
11737 for (
int i = 0; i <
GetNE(); i++)
11742 for (
int j = 0; j < nv; j++)
11747 for (
int i = 0; i <
GetNBE(); i++)
11752 for (
int j = 0; j < nv; j++)
11776 for (
int i = 0; i <
GetNE(); i++)
11789 int num_bdr_elem = 0;
11790 int new_bel_to_edge_nnz = 0;
11791 for (
int i = 0; i <
GetNBE(); i++)
11807 if (num_bdr_elem ==
GetNBE()) {
return; }
11811 Table *new_bel_to_edge = NULL;
11815 new_be_to_edge.
Reserve(num_bdr_elem);
11819 new_be_to_face.
Reserve(num_bdr_elem);
11820 new_bel_to_edge =
new Table;
11821 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
11823 for (
int i = 0; i <
GetNBE(); i++)
11834 int row = new_be_to_face.
Size();
11838 int *new_e = new_bel_to_edge->
GetRow(row);
11839 for (
int j = 0; j < ne; j++)
11843 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
11863 for (
int i = 0; i < attribs.
Size(); i++)
11875 #ifdef MFEM_USE_MEMALLOC
11902 const int npts = point_mat.
Width();
11903 if (!npts) {
return 0; }
11904 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
11908 if (!
GetNE()) {
return 0; }
11910 double *data = point_mat.
GetData();
11917 min_dist = std::numeric_limits<double>::max();
11921 for (
int i = 0; i <
GetNE(); i++)
11925 for (
int k = 0; k < npts; k++)
11928 if (dist < min_dist(k))
11930 min_dist(k) = dist;
11939 for (
int k = 0; k < npts; k++)
11943 int res = inv_tr->
Transform(pt, ips[k]);
11946 elem_ids[k] = e_idx[k];
11950 if (pts_found != npts)
11954 for (
int k = 0; k < npts; k++)
11956 if (elem_ids[k] != -1) {
continue; }
11960 for (
int v = 0; v < elvertices.
Size(); v++)
11962 int vv = elvertices[v];
11964 const int* els = vtoel->
GetRow(vv);
11965 for (
int e = 0; e < ne; e++)
11967 if (els[e] == e_idx[k]) {
continue; }
11969 int res = inv_tr->
Transform(pt, ips[k]);
11972 elem_ids[k] = els[e];
11984 for (
int e = 0; e < neigh.
Size(); e++)
11990 int res = inv_tr->
Transform(pt, ips[k]);
12003 if (inv_trans == NULL) {
delete inv_tr; }
12005 if (warn && pts_found != npts)
12007 MFEM_WARNING((npts-pts_found) <<
" points were not found");
12021 "mixed meshes are not supported!");
12022 MFEM_ASSERT(mesh->
GetNodes(),
"meshes without nodes are not supported!");
12035 Compute(nodes, d_mt);
12038 void GeometricFactors::Compute(
const GridFunction &nodes,
12045 const int vdim = fespace->
GetVDim();
12046 const int NE = fespace->
GetNE();
12047 const int ND = fe->
GetDof();
12050 unsigned eval_flags = 0;
12060 J.
SetSize(dim*vdim*NQ*NE, my_d_mt);
12075 qi->DisableTensorProducts(!use_tensor_products);
12083 Vector Enodes(vdim*ND*NE, my_d_mt);
12084 elem_restr->Mult(nodes, Enodes);
12085 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
12089 qi->Mult(nodes, eval_flags,
X,
J,
detJ);
12105 const int vdim = fespace->
GetVDim();
12119 face_restr->
Mult(*nodes, Fnodes);
12121 unsigned eval_flags = 0;
12130 J.
SetSize(vdim*vdim*NQ*NF, my_d_mt);
12168 V(1) = s * ((ip.
y + layer) / n);
12173 V(2) = s * ((ip.
z + layer) / n);
12182 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
12186 int nvy = (closed) ? (ny) : (ny + 1);
12187 int nvt = mesh->
GetNV() * nvy;
12196 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
12201 for (
int i = 0; i < mesh->
GetNV(); i++)
12204 for (
int j = 0; j < nvy; j++)
12206 vc[1] = sy * (double(j) / ny);
12212 for (
int i = 0; i < mesh->
GetNE(); i++)
12217 for (
int j = 0; j < ny; j++)
12220 qv[0] = vert[0] * nvy + j;
12221 qv[1] = vert[1] * nvy + j;
12222 qv[2] = vert[1] * nvy + (j + 1) % nvy;
12223 qv[3] = vert[0] * nvy + (j + 1) % nvy;
12229 for (
int i = 0; i < mesh->
GetNBE(); i++)
12234 for (
int j = 0; j < ny; j++)
12237 sv[0] = vert[0] * nvy + j;
12238 sv[1] = vert[0] * nvy + (j + 1) % nvy;
12242 Swap<int>(sv[0], sv[1]);
12254 for (
int i = 0; i < mesh->
GetNE(); i++)
12260 sv[0] = vert[0] * nvy;
12261 sv[1] = vert[1] * nvy;
12265 sv[0] = vert[1] * nvy + ny;
12266 sv[1] = vert[0] * nvy + ny;
12282 string cname = name;
12283 if (cname ==
"Linear")
12287 else if (cname ==
"Quadratic")
12291 else if (cname ==
"Cubic")
12295 else if (!strncmp(name,
"H1_", 3))
12299 else if (!strncmp(name,
"L2_T", 4))
12303 else if (!strncmp(name,
"L2_", 3))
12310 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
12322 for (
int i = 0; i < mesh->
GetNE(); i++)
12325 for (
int j = ny-1; j >= 0; j--)
12342 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
12347 int nvt = mesh->
GetNV() * nvz;
12352 bool wdgMesh =
false;
12353 bool hexMesh =
false;
12357 for (
int i = 0; i < mesh->
GetNV(); i++)
12361 for (
int j = 0; j < nvz; j++)
12363 vc[2] = sz * (double(j) / nz);
12369 for (
int i = 0; i < mesh->
GetNE(); i++)
12379 for (
int j = 0; j < nz; j++)
12382 pv[0] = vert[0] * nvz + j;
12383 pv[1] = vert[1] * nvz + j;
12384 pv[2] = vert[2] * nvz + j;
12385 pv[3] = vert[0] * nvz + (j + 1) % nvz;
12386 pv[4] = vert[1] * nvz + (j + 1) % nvz;
12387 pv[5] = vert[2] * nvz + (j + 1) % nvz;
12394 for (
int j = 0; j < nz; j++)
12397 hv[0] = vert[0] * nvz + j;
12398 hv[1] = vert[1] * nvz + j;
12399 hv[2] = vert[2] * nvz + j;
12400 hv[3] = vert[3] * nvz + j;
12401 hv[4] = vert[0] * nvz + (j + 1) % nvz;
12402 hv[5] = vert[1] * nvz + (j + 1) % nvz;
12403 hv[6] = vert[2] * nvz + (j + 1) % nvz;
12404 hv[7] = vert[3] * nvz + (j + 1) % nvz;
12406 mesh3d->
AddHex(hv, attr);
12410 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
12411 << geom <<
"\'" << endl;
12417 for (
int i = 0; i < mesh->
GetNBE(); i++)
12422 for (
int j = 0; j < nz; j++)
12425 qv[0] = vert[0] * nvz + j;
12426 qv[1] = vert[1] * nvz + j;
12427 qv[2] = vert[1] * nvz + (j + 1) % nvz;
12428 qv[3] = vert[0] * nvz + (j + 1) % nvz;
12437 for (
int i = 0; i < mesh->
GetNE(); i++)
12448 tv[0] = vert[0] * nvz;
12449 tv[1] = vert[2] * nvz;
12450 tv[2] = vert[1] * nvz;
12454 tv[0] = vert[0] * nvz + nz;
12455 tv[1] = vert[1] * nvz + nz;
12456 tv[2] = vert[2] * nvz + nz;
12464 qv[0] = vert[0] * nvz;
12465 qv[1] = vert[3] * nvz;
12466 qv[2] = vert[2] * nvz;
12467 qv[3] = vert[1] * nvz;
12471 qv[0] = vert[0] * nvz + nz;
12472 qv[1] = vert[1] * nvz + nz;
12473 qv[2] = vert[2] * nvz + nz;
12474 qv[3] = vert[3] * nvz + nz;
12480 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
12481 << geom <<
"\'" << endl;
12487 if ( hexMesh && wdgMesh )
12491 else if ( hexMesh )
12495 else if ( wdgMesh )
12508 string cname = name;
12509 if (cname ==
"Linear")
12513 else if (cname ==
"Quadratic")
12517 else if (cname ==
"Cubic")
12521 else if (!strncmp(name,
"H1_", 3))
12525 else if (!strncmp(name,
"L2_T", 4))
12529 else if (!strncmp(name,
"L2_", 3))
12536 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
12548 for (
int i = 0; i < mesh->
GetNE(); i++)
12551 for (
int j = nz-1; j >= 0; j--)
12572 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
12573 <<
" 0 0 " << i <<
" -1 0\n";
12580 double mid[3] = {0, 0, 0};
12581 for (
int j = 0; j < 2; j++)
12583 for (
int k = 0; k <
spaceDim; k++)
12588 os << NumOfVertices+i <<
" "
12589 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
12590 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
int GetNPoints() const
Returns the number of the points in the integration rule.
static Mesh MakePeriodic(const Mesh &orig_mesh, const std::vector< int > &v2v)
Create a periodic mesh by identifying vertices of orig_mesh.
Abstract class for all finite elements.
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
const IntegrationRule * IntRule
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Geometry::Type GetGeometryType() const
void Loader(std::istream &input, int generate_edges=0, std::string parse_tag="")
Mesh * Make3D(int nsteps, double rstep, double aspect, int order, bool sfc)
Ordering::Type GetOrdering() const
Return the ordering method.
int Size() const
Return the logical size of the array.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
void SetCoordsFromPatches(Vector &Nodes)
void PrintSurfaces(const Table &Aface_face, std::ostream &out) const
Print set of disjoint surfaces:
void Get(double *p, const int dim) const
virtual void Save(const char *fname, int precision=16) const
void GetPointMatrix(int i, DenseMatrix &pointmat) const
static const int vtk_quadratic_hex[27]
int * CartesianPartitioning(int nxyz[])
void METIS_PartGraphVKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
int GetDim() const
Returns the reference space dimension for the finite element.
int GetNDofs() const
Returns number of degrees of freedom.
Arc::Index insert_arc(Node::Index i, Node::Index j, Float w=1, Float b=1)
FaceInformation GetFaceInformation(int f) const
Class for an integration rule - an Array of IntegrationPoint.
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
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.
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
void ScaleElements(double sf)
static const int HighOrderMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to arbitrary-order Lagrange VTK geometries.
Class for grid function - Vector with associated FE space.
int DofToVDof(int dof, int vd, int ndofs=-1) const
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T...
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
T * end()
STL-like end. Returns pointer after the last element of the array.
void ReadVTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo)
void UseExternalData(double *ext_data, int i, int j, int k)
void FreeElement(Element *E)
void Mult(const Vector &x, Vector &y) const override=0
Extract the face degrees of freedom from x into y.
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
int CheckElementOrientation(bool fix_it=true)
Check (and optionally attempt to fix) the orientation of the elements.
void SetVertices(const Vector &vert_coord)
virtual Element * Duplicate(Mesh *m) const =0
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
static const int vtk_quadratic_tet[10]
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
void Make2D(int nx, int ny, Element::Type type, double sx, double sy, bool generate_edges, bool sfc_ordering)
Vector J
Jacobians of the element transformations at all quadrature points.
void AddColumnsInRow(int r, int ncol)
Vector X
Mapped (physical) coordinates of all quadrature points.
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 ...
static Mesh MakeSimplicial(const Mesh &orig_mesh)
Base class for vector Coefficients that optionally depend on time and space.
void UniformRefinement3D_base(Array< int > *f2qf=NULL, DSTable *v_to_v_p=NULL, bool update_nodes=true)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
class LinearPyramidFiniteElement PyramidFE
Array< Element * > boundary
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
int * GeneratePartitioning(int nparts, int part_method=1)
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
CoarseFineTransformations CoarseFineTr
virtual void LimitNCLevel(int max_nc_level)
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
void MoveVertices(const Vector &displacements)
void SetSize(int s)
Resize the vector to size s.
int GetNBE() const
Returns number of boundary elements.
static FiniteElement * GetTransformationFEforElementType(Element::Type)
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'.
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
FaceGeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, FaceType type, MemoryType d_mt=MemoryType::DEFAULT)
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, double sx=1.0, double sy=1.0, double sz=1.0, bool sfc_ordering=true)
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Lists all edges/faces in the nonconforming mesh.
virtual void UniformRefinement2D()
Refine a mixed 2D mesh uniformly.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Element::Type GetElementType(int i) const
Returns the type of element i.
double Norml2() const
Returns the l2 norm of the vector.
void ReadNetgen2DMesh(std::istream &input, int &curved)
void ShiftRight(int &a, int &b, int &c)
void SetDims(int rows, int nnz)
static const int FaceVert[NumFaces][MaxFaceVert]
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
static Mesh MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
Create a refined (by any factor) version of orig_mesh.
unsigned matrix
index into NCList::point_matrices[geom]
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
T * GetData()
Returns the data.
Vector detJ
Determinants of the Jacobians at all quadrature points.
static const int FaceVert[NumFaces][MaxFaceVert]
static const int Edges[NumEdges][2]
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
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. Returns the number assigned to the table entry.
Piecewise-(bi/tri)linear continuous finite elements.
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
void Print(std::ostream &out) const
I/O: Print the mesh in "MFEM NC mesh v1.0" format.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Mesh * Extrude1D(Mesh *mesh, const int ny, const double sy, const bool closed)
Extrude a 1D mesh.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
void Transform(void(*f)(const Vector &, Vector &))
int GetBdrElementEdgeIndex(int i) const
void AverageVertices(const int *indexes, int n, int result)
Averages the vertices with given indexes and saves the result in vertices[result].
int GetNE() const
Returns number of elements.
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
void FinalizeWedgeMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a wedge Mesh.
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
Evaluate the derivatives at quadrature points.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
NodeExtrudeCoefficient(const int dim, const int n_, const double s_)
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
void GetElementLocalToGlobal(Array< int > &lelem_elem)
void SetMeshGen()
Determine the mesh generator bitmask meshgen, see MeshGenerator().
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
Data type quadrilateral element.
void GetVertices(Vector &vert_coord) const
void ReadNetgen3DMesh(std::istream &input)
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Data arrays will be written in ASCII format.
double * GetData() const
Returns the matrix data array.
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
void ReadInlineMesh(std::istream &input, bool generate_edges=false)
void GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i...
double * GetData() const
Return a pointer to the beginning of the Vector data.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
void RemoveInternalBoundaries()
void DebugDump(std::ostream &out) const
Output an NCMesh-compatible debug dump.
Mesh * Extrude2D(Mesh *mesh, const int nz, const double sz)
Extrude a 2D mesh.
Geometry::Type GetElementBaseGeometry(int i) const
void Print(const Mesh &mesh, const adios2stream::mode print_mode=mode::sync)
int Size_of_connections() const
int spaceDim
dimensions of the elements and the vertex coordinates
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
friend class NURBSExtension
Geometry::Type GetBdrElementGeometry(int i) const
void GetVertexToVertexTable(DSTable &) const
void RedRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
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...
void ReadCubit(const char *filename, int &curved, int &read_gf)
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
void add(const Vector &v1, const Vector &v2, Vector &v)
void DeleteAll()
Delete the whole array.
void AddConnections(int r, const int *c, int nc)
int master
master number (in Mesh numbering)
void InitRefinementTransforms()
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
void UniformRefinement2D_base(bool update_nodes=true)
Array< NCFaceInfo > nc_faces_info
Element * ReadElement(std::istream &)
void KnotInsert(Array< KnotVector * > &kv)
void OnMeshUpdated(Mesh *mesh)
std::vector< int > CreatePeriodicVertexMapping(const std::vector< Vector > &translations, double tol=1e-8) const
Creates a mapping v2v from the vertex indices of the mesh such that coincident vertices under the giv...
Element * NewElement(int geom)
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
void SetVerticesFromNodes(const GridFunction *nodes)
Helper to set vertex coordinates given a high-order curvature function.
Structure for storing face geometric factors: coordinates, Jacobians, determinants of the Jacobians...
ElementTransformation * GetBdrElementTransformation(int i)
Returns the transformation defining the i-th boundary element.
virtual void UpdateMeshPointer(Mesh *new_mesh)
int AddBdrTriangle(int v1, int v2, int v3, int attr=1)
Mesh * Make2D(int nsteps, double rstep, double phi, double aspect, int order, bool sfc)
void PrintTopo(std::ostream &out, const Array< int > &e_to_k) const
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void MakeRefined_(Mesh &orig_mesh, const Array< int > ref_factors, int ref_type)
Internal function used in Mesh::MakeRefined.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Piecewise-(bi)cubic continuous finite elements.
int GetNE() const
Returns number of elements in the mesh.
uint rank(Node::Index i) const
PointFiniteElement PointFE
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void GetElementJacobian(int i, DenseMatrix &J)
Data type hexahedron element.
static int GetTetOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void GetVertexDofs(int i, Array< int > &dofs) const
int AddVertex(double x, double y=0.0, double z=0.0)
bool Nonconforming() const
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Vector J
Jacobians of the element transformations at all quadrature points.
Data type Pyramid element.
static const int Edges[NumEdges][2]
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void MoveNodes(const Vector &displacements)
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 ...
double f(const Vector &xvec)
Symmetric 3D Table stored as an array of rows each of which has a stack of column, floor, number nodes. The number of the node is assigned by counting the nodes from zero as they are pushed into the table. Diagonals of any kind are not allowed so the row, column and floor must all be different for each node. Only one node is stored for all 6 symmetric entries that are indexable by unique triplets of row, column, and floor.
const Table & GetDerefinementTable()
void PrintVTK(std::ostream &os)
Native ordering as defined by the FiniteElement.
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void Make3D(int nx, int ny, int nz, Element::Type type, double sx, double sy, double sz, bool sfc_ordering)
void AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
virtual void SetAttributes()
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
void EnsureNCMesh(bool simplices_nonconforming=false)
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
Mesh * GetMesh() const
Returns the mesh.
Vector detJ
Determinants of the Jacobians at all quadrature points.
int FindCoarseElement(int i)
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
virtual void Derefine(const Array< int > &derefs)
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
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.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
static const int Map[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to linear VTK geometries.
void AddQuadFaceElement(int lf, int gf, int el, int v0, int v1, int v2, int v3)
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
static const int NumVerts[NumGeom]
void CheckPartitioning(int *partitioning_)
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void AddConnection(int r, int c)
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
void ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
STable3D * GetElementToFaceTable(int ret_ftbl=0)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Type
Constants for the classes derived from Element.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const
void CheckDisplacements(const Vector &displacements, double &tmax)
This structure stores the low level information necessary to interpret the configuration of elements ...
NURBSExtension * StealNURBSext()
Array< DenseMatrix * > point_matrices[Geometry::NumGeom]
List of unique point matrices for each slave geometry.
VTKFormat
Data array format for VTK and VTU files.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
int AddBdrSegment(int v1, int v2, int attr=1)
int AddElement(Element *elem)
The parameter elem should be allocated using the NewElement() method.
int GetMaxElementOrder() const
Return the maximum polynomial order.
void ReadLineMesh(std::istream &input)
void GetBdrElementTopo(Array< Element * > &boundary) const
virtual void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
GeometryRefiner GlobGeometryRefiner
const CoarseFineTransformations & GetRefinementTransforms()
void SetLayer(const int l)
int GetAttribute() const
Return element's attribute.
int GetElementToEdgeTable(Table &, Array< int > &)
T * begin()
STL-like begin. Returns pointer to the first element of the array.
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...
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
void SetKnotsFromPatches()
void RandomRefinement(double prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Data type triangle element.
const Element * GetElement(int i) const
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
void GetElementTopo(Array< Element * > &elements) const
signed char local
local number within 'element'
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
int GetVDim()
Returns dimension of the vector.
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
int Size() const
Returns the number of TYPE I elements.
int GetVDim() const
Returns vector dimension.
Vector normal
Normals at all quadrature points.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
static const int Edges[NumEdges][2]
FiniteElementSpace * FESpace()
virtual MFEM_DEPRECATED void ReorientTetMesh()
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
static const int QuadraticMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to legacy quadratic VTK geometries/.
void GetColumn(int c, Vector &col) const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
void ReadMFEMMesh(std::istream &input, int version, int &curved)
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
This function is not public anymore. Use GeneralRefinement instead.
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
void FindPartitioningComponents(Table &elem_elem, const Array< int > &partitioning, Array< int > &component, Array< int > &num_comp)
Data type tetrahedron element.
void GetElementInteriorDofs(int i, Array< int > &dofs) const
List of mesh geometries stored as Array<Geometry::Type>.
A general vector function coefficient.
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)
int CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
int SpaceDimension() const
GeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, MemoryType d_mt=MemoryType::DEFAULT)
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
void PrintElementsWithPartitioning(int *partitioning, std::ostream &out, int interior_faces=0)
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
static const int Edges[NumEdges][2]
void AddPointFaceElement(int lf, int gf, int el)
Used in GenerateFaces()
void Make1D(int n, double sx=1.0)
Creates a 1D mesh for the interval [0,sx] divided into n equal intervals.
double * Data() const
Returns the matrix data array.
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]...
void Swap(Array< T > &, Array< T > &)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
virtual void NonconformingRefinement(const Array< Refinement > &refinements, int nc_limit=0)
This function is not public anymore. Use GeneralRefinement instead.
virtual void PushTransform(int tr)
Add 'tr' to the current chain of coarse-fine transformations.
void FinalizeMesh(int refine=0, bool fix_orientation=true)
Finalize the construction of any type of Mesh.
double p(const Vector &x, double t)
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
const NURBSExtension * GetNURBSext() const
int SpaceDimension() const
Return the space dimension of the NCMesh.
Nonconforming edge/face within a bigger edge/face.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
static Mesh LoadFromFile(const char *filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
static void PrintElement(const Element *, std::ostream &)
virtual void NURBSUniformRefinement()
Refine NURBS mesh.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
const IntegrationRule * IntRule
MemoryType
Memory types supported by MFEM.
Vector X
Mapped (physical) coordinates of all quadrature points.
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
static void GetElementArrayEdgeTable(const Array< Element * > &elem_array, const DSTable &v_to_v, Table &el_to_edge)
void AddAColumnInRow(int r)
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
int AddBdrQuad(int v1, int v2, int v3, int v4, int attr=1)
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
int FindRoots(const Vector &z, Vector &x)
bool IsSlaveFace(const FaceInfo &fi) const
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. The entry is addressed by the three smallest values of (r,c,f,t). Returns the number assigned to the table entry.
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, double sx=1.0, double sy=1.0, bool sfc_ordering=true)
void InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
Begin construction of a mesh.
void ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo, const std::string &xml_prefix="")
Helper struct for defining a connectivity table, see Table::MakeFromList.
Geometry::Type GetElementGeometry(int i) const
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
class Linear3DFiniteElement TetrahedronFE
A class to initialize the size of a Tensor.
static Mesh MakeCartesian1D(int n, double sx=1.0)
Array< Element * > elements
const Table & ElementToElementTable()
double GetLength(int i, int j) const
Return the length of the segment from node i to node j.
void SetRelaxedHpConformity(bool relaxed=true)
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Linear2DFiniteElement TriangleFE
static const int FaceVert[NumFaces][MaxFaceVert]
const CoarseFineTransformations & GetRefinementTransforms()
bool RefineByError(const Array< double > &elem_error, double threshold, int nonconforming=-1, int nc_limit=0)
Evaluate the values at quadrature points.
Operation GetLastOperation() const
Return type of last modification of the mesh.
Table * GetFaceToElementTable() const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
NURBSExtension * NURBSext
Optional NURBS mesh extension.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void PrintBdrVTU(std::string fname, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0)
void GetNode(int i, double *coord) const
virtual const char * Name() const
const Table & ElementToFaceTable() const
static const int Edges[NumEdges][2]
Class for integration point with weight.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
Element * ReadElementWithoutAttr(std::istream &)
void Swap(Mesh &other, bool non_geometry)
int GetNumFacesWithGhost() const
Return the number of faces (3D), edges (2D) or vertices (1D) including ghost faces.
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
static const int vtk_quadratic_wedge[18]
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
void DisableTensorProducts(bool disable=true) const
Disable the use of tensor product evaluations, for tensor-product elements, e.g. quads and hexes...
static const int DimStart[MaxDim+2]
const FiniteElementSpace * GetNodalFESpace() const
class LinearWedgeFiniteElement WedgeFE
virtual void Print(std::ostream &os=mfem::out) const
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Array< int > leaf_elements
finest elements, in Mesh ordering (+ ghosts)
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
static const int Orient[NumOrient][NumVert]
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
void DegreeElevate(int rel_degree, int degree=16)
signed char geom
Geometry::Type (faces only) (char to save RAM)
void FindTMax(Vector &c, Vector &x, double &tmax, const double factor, const int Dim)
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int index(int i, int j, int nx, int ny)
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
static const int Edges[NumEdges][2]
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
Lexicographic ordering for tensor-product FiniteElements.
Array< FaceInfo > faces_info
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
virtual int GetNVertices() const =0
int parent
Coarse Element index in the coarse mesh.
void SetAttribute(const int attr)
Set element's attribute.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
static const int Orient[NumOrient][NumVert]
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
STable3D * GetFacesTable()
Piecewise-(bi)quadratic continuous finite elements.
int Dimension() const
Return the dimension of the NCMesh.
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
NCMesh * ncmesh
Optional nonconforming mesh extension.
Geometry::Type GetBdrElementBaseGeometry(int i) const
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
virtual void UniformRefinement3D()
Refine a mixed 3D mesh uniformly.
int AddBdrElement(Element *elem)
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
static void PrintElementWithoutAttr(const Element *, std::ostream &)
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
BlockArray< Element > elements
virtual bool NonconformingDerefinement(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
NC version of GeneralDerefinement.
void Clear()
Clear the contents of the Mesh.
int GetNEdges() const
Return the number of edges.
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
double AggregateError(const Array< double > &elem_error, const int *fine, int nfine, int op)
Derefinement helper.
int GetNFaces() const
Return the number of faces in a 3D mesh.
BiLinear2DFiniteElement QuadrilateralFE
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...
Mesh & operator=(Mesh &&mesh)
Move assignment operstor.
virtual void PrintXG(std::ostream &os=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
Evaluate the derivatives at quadrature points.
void ReadTrueGridMesh(std::istream &input)
const FiniteElementCollection * FEColl() const
void GetNodes(Vector &node_coord) const
virtual int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type, does not count master nonconforming face...
int NumberOfEntries() const
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
void SetNode(int i, const double *coord)
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
virtual void PushTransform(int tr)
Add 'tr' to the current chain of coarse-fine transformations.
void CreateVTKElementConnectivity(Array< int > &con, Geometry::Type geom, int ref)
Create the VTK element connectivity array for a given element geometry and refinement level...
int GetId(int p1, int p2)
Get id of item whose parents are p1, p2... Create it if it doesn't exist.
Arbitrary order H1-conforming (continuous) finite elements.
void XYZ_VectorFunction(const Vector &p, Vector &v)
void GetElementColoring(Array< int > &colors, int el0=0)
void ChangeVertexDataOwnership(double *vertices, int len_vertices, bool zerocopy=false)
Set the internal Vertex array to point to the given vertices array without assuming ownership of the ...
const std::string filename
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
MemAlloc< Tetrahedron, 1024 > TetMemory
Table * GetFaceEdgeTable() const
Returns the face-to-edge Table (3D)
bool IsGhost(const Element &el) const
Return true if the Element el is a ghost element.
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
Node::Index insert_node(Float length=1)
TriLinear3DFiniteElement HexahedronFE
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
void SetNodes(const Vector &node_coord)
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.
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
static const int * VertexPermutation[Geometry::NUM_GEOMETRIES]
Permutation from MFEM's vertex ordering to VTK's vertex ordering.
void ScaleSubdomains(double sf)
void DegreeElevate(int rel_degree, int degree=16)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
Base class for operators that extracts Face degrees of freedom.
void Bisection(int i, const DSTable &, int *, int *, int *)
Bisect a triangle: element with index i is bisected.
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn't exist.
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...
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Rank 3 tensor (array of matrices)
Class for parallel meshes.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Abstract data type element.
int GetAttribute(int i) const
Return the attribute of element i.
void AddTriangleFaceElement(int lf, int gf, int el, int v0, int v1, int v2)
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...
const Table & ElementToEdgeTable() const
Data type line segment element.
Linear1DFiniteElement SegmentFE
Array< GeometricFactors * > geom_factors
Optional geometric factors.
int GetNFDofs() const
Number of all scalar face-interior dofs.
void BdrBisection(int i, const HashTable< Hashed2 > &)
Bisect a boundary triangle: boundary element with index i is bisected.
MPI_Comm GetGlobalMPI_Comm()
Get MFEM's "global" MPI communicator.
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
int NumberOfElements()
Return the number of elements added to the table.
void GenerateNCFaceInfo()
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
virtual Type GetType() const =0
Returns element's type.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Table * GetVertexToElementTable()
The returned Table must be destroyed by the caller.
const Element * GetBdrElement(int i) const
void ConvertToPatches(const Vector &Nodes)
static const int Orient[NumOrient][NumVert]
void GetElementData(const Array< Element * > &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void UpdateNodes()
Update the nodes of a curved mesh after refinement.
Defines the position of a fine element within a coarse element.
void KnotInsert(Array< KnotVector * > &kv)
void Print(std::ostream &out) const
virtual void Refine(const Array< Refinement > &refinements)
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
void GreenRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
Arbitrary order "L2-conforming" discontinuous finite elements.
Evaluate the values at quadrature points.
double f(const Vector &p)
void Printer(std::ostream &out=mfem::out, std::string section_delimiter="") const
static const int FaceVert[NumFaces][MaxFaceVert]
void NodesUpdated()
This function should be called after the mesh node coordinates have changed, e.g. after the mesh has ...
Class used to extrude the nodes of a mesh.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...