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);
179 T->Transform(RefG->
RefPts, pointmat);
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));
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; }
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[GetFaceGeometry(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(GetFaceGeometry(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;
1470 face_to_elem = NULL;
1473 void Mesh::SetEmpty()
1479 void Mesh::DestroyTables()
1484 DeleteGeometricFactors();
1494 delete face_to_elem;
1495 face_to_elem = NULL;
1498 void Mesh::DestroyPointers()
1500 if (own_nodes) {
delete Nodes; }
1506 for (
int i = 0; i < NumOfElements; i++)
1508 FreeElement(elements[i]);
1511 for (
int i = 0; i < NumOfBdrElements; i++)
1513 FreeElement(boundary[i]);
1516 for (
int i = 0; i < faces.Size(); i++)
1518 FreeElement(faces[i]);
1524 void Mesh::Destroy()
1528 elements.DeleteAll();
1529 vertices.DeleteAll();
1530 boundary.DeleteAll();
1532 faces_info.DeleteAll();
1533 nc_faces_info.DeleteAll();
1534 be_to_edge.DeleteAll();
1535 be_to_face.DeleteAll();
1543 CoarseFineTr.Clear();
1545 #ifdef MFEM_USE_MEMALLOC 1549 attributes.DeleteAll();
1550 bdr_attributes.DeleteAll();
1553 void Mesh::ResetLazyData()
1555 delete el_to_el; el_to_el = NULL;
1556 delete face_edge; face_edge = NULL;
1557 delete face_to_elem; face_to_elem = NULL;
1558 delete edge_vertex; edge_vertex = NULL;
1559 DeleteGeometricFactors();
1560 nbInteriorFaces = -1;
1561 nbBoundaryFaces = -1;
1564 void Mesh::SetAttributes()
1569 for (
int i = 0; i < attribs.
Size(); i++)
1571 attribs[i] = GetBdrAttribute(i);
1575 attribs.
Copy(bdr_attributes);
1576 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1578 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1582 for (
int i = 0; i < attribs.
Size(); i++)
1584 attribs[i] = GetAttribute(i);
1588 attribs.
Copy(attributes);
1589 if (attributes.Size() > 0 && attributes[0] <= 0)
1591 MFEM_WARNING(
"Non-positive attributes in the domain!");
1595 void Mesh::InitMesh(
int Dim_,
int spaceDim_,
int NVert,
int NElem,
int NBdrElem)
1600 spaceDim = spaceDim_;
1603 vertices.SetSize(NVert);
1606 elements.SetSize(NElem);
1608 NumOfBdrElements = 0;
1609 boundary.SetSize(NBdrElem);
1612 template<
typename T>
1613 static void CheckEnlarge(
Array<T> &array,
int size)
1615 if (size >= array.
Size()) { array.
SetSize(size + 1); }
1618 int Mesh::AddVertex(
double x,
double y,
double z)
1620 CheckEnlarge(vertices, NumOfVertices);
1621 double *v = vertices[NumOfVertices]();
1625 return NumOfVertices++;
1628 int Mesh::AddVertex(
const double *coords)
1630 CheckEnlarge(vertices, NumOfVertices);
1631 vertices[NumOfVertices].SetCoords(spaceDim, coords);
1632 return NumOfVertices++;
1637 MFEM_ASSERT(coords.
Size() >= spaceDim,
1638 "invalid 'coords' size: " << coords.
Size());
1639 return AddVertex(coords.
GetData());
1642 void Mesh::AddVertexParents(
int i,
int p1,
int p2)
1648 if (i < vertices.Size())
1650 double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1651 for (
int j = 0; j < 3; j++)
1653 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1658 int Mesh::AddSegment(
int v1,
int v2,
int attr)
1660 CheckEnlarge(elements, NumOfElements);
1661 elements[NumOfElements] =
new Segment(v1, v2, attr);
1662 return NumOfElements++;
1665 int Mesh::AddSegment(
const int *vi,
int attr)
1667 CheckEnlarge(elements, NumOfElements);
1668 elements[NumOfElements] =
new Segment(vi, attr);
1669 return NumOfElements++;
1672 int Mesh::AddTriangle(
int v1,
int v2,
int v3,
int attr)
1674 CheckEnlarge(elements, NumOfElements);
1675 elements[NumOfElements] =
new Triangle(v1, v2, v3, attr);
1676 return NumOfElements++;
1679 int Mesh::AddTriangle(
const int *vi,
int attr)
1681 CheckEnlarge(elements, NumOfElements);
1682 elements[NumOfElements] =
new Triangle(vi, attr);
1683 return NumOfElements++;
1686 int Mesh::AddQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1688 CheckEnlarge(elements, NumOfElements);
1689 elements[NumOfElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1690 return NumOfElements++;
1693 int Mesh::AddQuad(
const int *vi,
int attr)
1695 CheckEnlarge(elements, NumOfElements);
1697 return NumOfElements++;
1700 int Mesh::AddTet(
int v1,
int v2,
int v3,
int v4,
int attr)
1702 int vi[4] = {v1, v2, v3, v4};
1703 return AddTet(vi, attr);
1706 int Mesh::AddTet(
const int *vi,
int attr)
1708 CheckEnlarge(elements, NumOfElements);
1709 #ifdef MFEM_USE_MEMALLOC 1711 tet = TetMemory.Alloc();
1714 elements[NumOfElements] = tet;
1716 elements[NumOfElements] =
new Tetrahedron(vi, attr);
1718 return NumOfElements++;
1721 int Mesh::AddWedge(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int attr)
1723 CheckEnlarge(elements, NumOfElements);
1724 elements[NumOfElements] =
new Wedge(v1, v2, v3, v4, v5, v6, attr);
1725 return NumOfElements++;
1728 int Mesh::AddWedge(
const int *vi,
int attr)
1730 CheckEnlarge(elements, NumOfElements);
1731 elements[NumOfElements] =
new Wedge(vi, attr);
1732 return NumOfElements++;
1735 int Mesh::AddPyramid(
int v1,
int v2,
int v3,
int v4,
int v5,
int attr)
1737 CheckEnlarge(elements, NumOfElements);
1738 elements[NumOfElements] =
new Pyramid(v1, v2, v3, v4, v5, attr);
1739 return NumOfElements++;
1742 int Mesh::AddPyramid(
const int *vi,
int attr)
1744 CheckEnlarge(elements, NumOfElements);
1745 elements[NumOfElements] =
new Pyramid(vi, attr);
1746 return NumOfElements++;
1749 int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1752 CheckEnlarge(elements, NumOfElements);
1753 elements[NumOfElements] =
1754 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1755 return NumOfElements++;
1758 int Mesh::AddHex(
const int *vi,
int attr)
1760 CheckEnlarge(elements, NumOfElements);
1761 elements[NumOfElements] =
new Hexahedron(vi, attr);
1762 return NumOfElements++;
1765 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1767 static const int hex_to_tet[6][4] =
1769 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1770 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1774 for (
int i = 0; i < 6; i++)
1776 for (
int j = 0; j < 4; j++)
1778 ti[j] = vi[hex_to_tet[i][j]];
1784 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1786 static const int hex_to_wdg[2][6] =
1788 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1792 for (
int i = 0; i < 2; i++)
1794 for (
int j = 0; j < 6; j++)
1796 ti[j] = vi[hex_to_wdg[i][j]];
1802 void Mesh::AddHexAsPyramids(
const int *vi,
int attr)
1804 static const int hex_to_pyr[6][5] =
1806 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
1807 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
1811 for (
int i = 0; i < 6; i++)
1813 for (
int j = 0; j < 5; j++)
1815 ti[j] = vi[hex_to_pyr[i][j]];
1817 AddPyramid(ti, attr);
1823 CheckEnlarge(elements, NumOfElements);
1824 elements[NumOfElements] = elem;
1825 return NumOfElements++;
1830 CheckEnlarge(boundary, NumOfBdrElements);
1831 boundary[NumOfBdrElements] = elem;
1832 return NumOfBdrElements++;
1835 int Mesh::AddBdrSegment(
int v1,
int v2,
int attr)
1837 CheckEnlarge(boundary, NumOfBdrElements);
1838 boundary[NumOfBdrElements] =
new Segment(v1, v2, attr);
1839 return NumOfBdrElements++;
1842 int Mesh::AddBdrSegment(
const int *vi,
int attr)
1844 CheckEnlarge(boundary, NumOfBdrElements);
1845 boundary[NumOfBdrElements] =
new Segment(vi, attr);
1846 return NumOfBdrElements++;
1849 int Mesh::AddBdrTriangle(
int v1,
int v2,
int v3,
int attr)
1851 CheckEnlarge(boundary, NumOfBdrElements);
1852 boundary[NumOfBdrElements] =
new Triangle(v1, v2, v3, attr);
1853 return NumOfBdrElements++;
1856 int Mesh::AddBdrTriangle(
const int *vi,
int attr)
1858 CheckEnlarge(boundary, NumOfBdrElements);
1859 boundary[NumOfBdrElements] =
new Triangle(vi, attr);
1860 return NumOfBdrElements++;
1863 int Mesh::AddBdrQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1865 CheckEnlarge(boundary, NumOfBdrElements);
1866 boundary[NumOfBdrElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1867 return NumOfBdrElements++;
1870 int Mesh::AddBdrQuad(
const int *vi,
int attr)
1872 CheckEnlarge(boundary, NumOfBdrElements);
1874 return NumOfBdrElements++;
1877 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1879 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1882 for (
int i = 0; i < 2; i++)
1884 for (
int j = 0; j < 3; j++)
1886 ti[j] = vi[quad_to_tri[i][j]];
1888 AddBdrTriangle(ti, attr);
1892 int Mesh::AddBdrPoint(
int v,
int attr)
1894 CheckEnlarge(boundary, NumOfBdrElements);
1895 boundary[NumOfBdrElements] =
new Point(&v, attr);
1896 return NumOfBdrElements++;
1899 void Mesh::GenerateBoundaryElements()
1902 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1906 for (i = 0; i < boundary.Size(); i++)
1908 FreeElement(boundary[i]);
1918 NumOfBdrElements = 0;
1919 for (i = 0; i < faces_info.Size(); i++)
1921 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1924 boundary.SetSize(NumOfBdrElements);
1925 be2face.
SetSize(NumOfBdrElements);
1926 for (j = i = 0; i < faces_info.Size(); i++)
1928 if (faces_info[i].Elem2No < 0)
1930 boundary[j] = faces[i]->Duplicate(
this);
1937 void Mesh::FinalizeCheck()
1939 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1940 vertices.Size() == 0,
1941 "incorrect number of vertices: preallocated: " << vertices.Size()
1942 <<
", actually added: " << NumOfVertices);
1943 MFEM_VERIFY(elements.Size() == NumOfElements,
1944 "incorrect number of elements: preallocated: " << elements.Size()
1945 <<
", actually added: " << NumOfElements);
1946 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1947 "incorrect number of boundary elements: preallocated: " 1948 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1951 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1954 CheckElementOrientation(fix_orientation);
1958 MarkTriMeshForRefinement();
1963 el_to_edge =
new Table;
1964 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1966 CheckBdrElementOrientation();
1980 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1981 bool fix_orientation)
1984 if (fix_orientation)
1986 CheckElementOrientation(fix_orientation);
1991 el_to_edge =
new Table;
1992 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1994 CheckBdrElementOrientation();
2014 GeckoProgress(
double limit) : limit(limit) { sw.
Start(); }
2015 virtual bool quit()
const {
return limit > 0 && sw.
UserTime() > limit; }
2018 class GeckoVerboseProgress :
public GeckoProgress
2024 GeckoVerboseProgress(
double limit) : GeckoProgress(limit) {}
2026 virtual void beginorder(
const Graph* graph,
Float cost)
const 2027 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2028 virtual void endorder(
const Graph* graph,
Float cost)
const 2029 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2031 virtual void beginiter(
const Graph* graph,
2034 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window " 2035 << window << std::flush;
2037 virtual void enditer(
const Graph* graph,
Float mincost,
Float cost)
const 2038 {
mfem::out <<
", cost = " << cost << endl; }
2043 int iterations,
int window,
2044 int period,
int seed,
bool verbose,
2050 GeckoProgress progress(time_limit);
2051 GeckoVerboseProgress vprogress(time_limit);
2054 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2061 const Table &my_el_to_el = ElementToElementTable();
2062 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2064 const int *neighid = my_el_to_el.
GetRow(elemid);
2065 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2067 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2072 graph.
order(&functional, iterations, window, period, seed,
2073 verbose ? &vprogress : &progress);
2079 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2082 return graph.
cost();
2093 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
2094 : coord(coord), dir(dir), points(points), mid(mid) {}
2096 bool operator()(
int i)
const 2098 return (points[3*i + coord] < mid) != dir;
2102 static void HilbertSort2D(
int coord1,
2105 const Array<double> &points,
int *beg,
int *end,
2106 double xmin,
double ymin,
double xmax,
double ymax)
2108 if (end - beg <= 1) {
return; }
2110 double xmid = (xmin + xmax)*0.5;
2111 double ymid = (ymin + ymax)*0.5;
2113 int coord2 = (coord1 + 1) % 2;
2116 int *p0 = beg, *p4 = end;
2117 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2118 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2119 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2123 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2124 ymin, xmin, ymid, xmid);
2126 if (p1 != p0 || p2 != p4)
2128 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2129 xmin, ymid, xmid, ymax);
2131 if (p2 != p0 || p3 != p4)
2133 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2134 xmid, ymid, xmax, ymax);
2138 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2139 ymid, xmax, ymin, xmid);
2143 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2144 const Array<double> &points,
int *beg,
int *end,
2145 double xmin,
double ymin,
double zmin,
2146 double xmax,
double ymax,
double zmax)
2148 if (end - beg <= 1) {
return; }
2150 double xmid = (xmin + xmax)*0.5;
2151 double ymid = (ymin + ymax)*0.5;
2152 double zmid = (zmin + zmax)*0.5;
2154 int coord2 = (coord1 + 1) % 3;
2155 int coord3 = (coord1 + 2) % 3;
2158 int *p0 = beg, *p8 = end;
2159 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2160 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2161 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2162 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2163 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2164 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2165 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2169 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2170 zmin, xmin, ymin, zmid, xmid, ymid);
2172 if (p1 != p0 || p2 != p8)
2174 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2175 ymin, zmid, xmin, ymid, zmax, xmid);
2177 if (p2 != p0 || p3 != p8)
2179 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2180 ymid, zmid, xmin, ymax, zmax, xmid);
2182 if (p3 != p0 || p4 != p8)
2184 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2185 xmin, ymax, zmid, xmid, ymid, zmin);
2187 if (p4 != p0 || p5 != p8)
2189 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2190 xmid, ymax, zmid, xmax, ymid, zmin);
2192 if (p5 != p0 || p6 != p8)
2194 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2195 ymax, zmid, xmax, ymid, zmax, xmid);
2197 if (p6 != p0 || p7 != p8)
2199 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2200 ymid, zmid, xmax, ymin, zmax, xmid);
2204 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2205 zmid, xmax, ymin, zmin, xmid, ymid);
2211 MFEM_VERIFY(spaceDim <= 3,
"");
2214 GetBoundingBox(min, max);
2219 if (spaceDim < 3) { points = 0.0; }
2222 for (
int i = 0; i < GetNE(); i++)
2224 GetElementCenter(i, center);
2225 for (
int j = 0; j < spaceDim; j++)
2227 points[3*i + j] = center(j);
2234 indices.
Sort([&](
int a,
int b)
2235 {
return points[3*
a] < points[3*
b]; });
2237 else if (spaceDim == 2)
2240 HilbertSort2D(0,
false,
false,
2241 points, indices.
begin(), indices.
end(),
2242 min(0), min(1), max(0), max(1));
2247 HilbertSort3D(0,
false,
false,
false,
2248 points, indices.
begin(), indices.
end(),
2249 min(0), min(1), min(2), max(0), max(1), max(2));
2254 for (
int i = 0; i < GetNE(); i++)
2256 ordering[indices[i]] = i;
2261 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
2265 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2270 MFEM_WARNING(
"element reordering of non-conforming meshes is not" 2274 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
2304 old_elem_node_vals.SetSize(GetNE());
2305 nodes_fes = Nodes->FESpace();
2308 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2310 nodes_fes->GetElementVDofs(old_elid, old_dofs);
2312 old_elem_node_vals[old_elid] =
new Vector(vals);
2318 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2320 int new_elid = ordering[old_elid];
2321 new_elements[new_elid] = elements[old_elid];
2326 if (reorder_vertices)
2331 vertex_ordering = -1;
2333 int new_vertex_ind = 0;
2334 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2336 int *elem_vert = elements[new_elid]->GetVertices();
2337 int nv = elements[new_elid]->GetNVertices();
2338 for (
int vi = 0; vi < nv; ++vi)
2340 int old_vertex_ind = elem_vert[vi];
2341 if (vertex_ordering[old_vertex_ind] == -1)
2343 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2344 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
2354 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2356 int *elem_vert = elements[new_elid]->GetVertices();
2357 int nv = elements[new_elid]->GetNVertices();
2358 for (
int vi = 0; vi < nv; ++vi)
2360 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2365 for (
int belid = 0; belid < GetNBE(); ++belid)
2367 int *be_vert = boundary[belid]->GetVertices();
2368 int nv = boundary[belid]->GetNVertices();
2369 for (
int vi = 0; vi < nv; ++vi)
2371 be_vert[vi] = vertex_ordering[be_vert[vi]];
2382 el_to_edge =
new Table;
2383 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2388 GetElementToFaceTable();
2398 last_operation = Mesh::NONE;
2399 nodes_fes->Update(
false);
2402 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2404 int new_elid = ordering[old_elid];
2405 nodes_fes->GetElementVDofs(new_elid, new_dofs);
2406 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
2407 delete old_elem_node_vals[old_elid];
2413 void Mesh::MarkForRefinement()
2419 MarkTriMeshForRefinement();
2423 DSTable v_to_v(NumOfVertices);
2424 GetVertexToVertexTable(v_to_v);
2425 MarkTetMeshForRefinement(v_to_v);
2430 void Mesh::MarkTriMeshForRefinement()
2435 for (
int i = 0; i < NumOfElements; i++)
2437 if (elements[i]->GetType() == Element::TRIANGLE)
2439 GetPointMatrix(i, pmat);
2440 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
2451 for (
int i = 0; i < NumOfVertices; i++)
2456 length_idx[j].one = GetLength(i, it.Column());
2457 length_idx[j].two = j;
2464 for (
int i = 0; i < NumOfEdges; i++)
2466 order[length_idx[i].two] = i;
2470 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
2475 GetEdgeOrdering(v_to_v, order);
2477 for (
int i = 0; i < NumOfElements; i++)
2479 if (elements[i]->GetType() == Element::TETRAHEDRON)
2481 elements[i]->MarkEdge(v_to_v, order);
2484 for (
int i = 0; i < NumOfBdrElements; i++)
2486 if (boundary[i]->GetType() == Element::TRIANGLE)
2488 boundary[i]->MarkEdge(v_to_v, order);
2495 if (*old_v_to_v && *old_elem_vert)
2502 if (*old_v_to_v == NULL)
2504 bool need_v_to_v =
false;
2506 for (
int i = 0; i < GetNEdges(); i++)
2512 if (dofs.
Size() > 0)
2520 *old_v_to_v =
new DSTable(NumOfVertices);
2521 GetVertexToVertexTable(*(*old_v_to_v));
2524 if (*old_elem_vert == NULL)
2526 bool need_elem_vert =
false;
2528 for (
int i = 0; i < GetNE(); i++)
2534 if (dofs.
Size() > 1)
2536 need_elem_vert =
true;
2542 *old_elem_vert =
new Table;
2543 (*old_elem_vert)->
MakeI(GetNE());
2544 for (
int i = 0; i < GetNE(); i++)
2546 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2548 (*old_elem_vert)->MakeJ();
2549 for (
int i = 0; i < GetNE(); i++)
2551 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2552 elements[i]->GetNVertices());
2554 (*old_elem_vert)->ShiftUpI();
2567 const int num_edge_dofs = old_dofs.
Size();
2570 const Vector onodes = *Nodes;
2574 int offset = NumOfVertices * old_dofs.
Size();
2578 if (num_edge_dofs > 0)
2580 DSTable new_v_to_v(NumOfVertices);
2581 GetVertexToVertexTable(new_v_to_v);
2583 for (
int i = 0; i < NumOfVertices; i++)
2587 const int old_i = (*old_v_to_v)(i, it.Column());
2588 const int new_i = it.Index();
2589 if (new_i == old_i) {
continue; }
2591 old_dofs.
SetSize(num_edge_dofs);
2592 new_dofs.
SetSize(num_edge_dofs);
2593 for (
int j = 0; j < num_edge_dofs; j++)
2595 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2596 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2600 for (
int j = 0; j < old_dofs.
Size(); j++)
2602 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2606 offset += NumOfEdges * num_edge_dofs;
2614 Table old_face_vertex;
2615 old_face_vertex.
MakeI(NumOfFaces);
2616 for (
int i = 0; i < NumOfFaces; i++)
2620 old_face_vertex.
MakeJ();
2621 for (
int i = 0; i < NumOfFaces; i++)
2623 faces[i]->GetNVertices());
2627 STable3D *faces_tbl = GetElementToFaceTable(1);
2633 for (
int i = 0; i < NumOfFaces; i++)
2635 const int *old_v = old_face_vertex.
GetRow(i);
2637 switch (old_face_vertex.
RowSize(i))
2640 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2644 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2648 new_fdofs[new_i+1] = old_dofs.
Size();
2653 for (
int i = 0; i < NumOfFaces; i++)
2655 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2658 switch (old_face_vertex.
RowSize(i))
2661 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2662 new_v = faces[new_i]->GetVertices();
2663 new_or = GetTriOrientation(old_v, new_v);
2668 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2669 new_v = faces[new_i]->GetVertices();
2670 new_or = GetQuadOrientation(old_v, new_v);
2677 for (
int j = 0; j < old_dofs.
Size(); j++)
2680 const int old_j = dof_ord[j];
2681 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2685 for (
int j = 0; j < old_dofs.
Size(); j++)
2687 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2709 for (
int i = 0; i < GetNE(); i++)
2711 const int *old_v = old_elem_vert->
GetRow(i);
2712 const int *new_v = elements[i]->GetVertices();
2718 case Geometry::SEGMENT:
2719 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2721 case Geometry::TRIANGLE:
2722 new_or = GetTriOrientation(old_v, new_v);
2724 case Geometry::SQUARE:
2725 new_or = GetQuadOrientation(old_v, new_v);
2727 case Geometry::TETRAHEDRON:
2728 new_or = GetTetOrientation(old_v, new_v);
2732 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2733 <<
" FE collection) are not supported yet!");
2737 MFEM_VERIFY(dof_ord != NULL,
2738 "FE collection '" << fec->
Name()
2739 <<
"' does not define reordering for " 2740 << Geometry::Name[geom] <<
" elements!");
2743 for (
int j = 0; j < new_dofs.
Size(); j++)
2746 const int old_j = dof_ord[j];
2747 new_dofs[old_j] = offset + j;
2749 offset += new_dofs.
Size();
2752 for (
int j = 0; j < old_dofs.
Size(); j++)
2754 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2766 GetElementToFaceTable();
2769 CheckBdrElementOrientation();
2774 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2779 CheckBdrElementOrientation();
2784 last_operation = Mesh::NONE;
2789 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2792 CheckElementOrientation(fix_orientation);
2794 if (NumOfBdrElements == 0)
2796 GetElementToFaceTable();
2798 GenerateBoundaryElements();
2803 DSTable v_to_v(NumOfVertices);
2804 GetVertexToVertexTable(v_to_v);
2805 MarkTetMeshForRefinement(v_to_v);
2808 GetElementToFaceTable();
2811 CheckBdrElementOrientation();
2813 if (generate_edges == 1)
2815 el_to_edge =
new Table;
2816 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2830 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2831 bool fix_orientation)
2834 CheckElementOrientation(fix_orientation);
2836 if (NumOfBdrElements == 0)
2838 GetElementToFaceTable();
2840 GenerateBoundaryElements();
2843 GetElementToFaceTable();
2846 CheckBdrElementOrientation();
2848 if (generate_edges == 1)
2850 el_to_edge =
new Table;
2851 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2865 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2868 CheckElementOrientation(fix_orientation);
2870 GetElementToFaceTable();
2873 if (NumOfBdrElements == 0)
2875 GenerateBoundaryElements();
2878 CheckBdrElementOrientation();
2882 el_to_edge =
new Table;
2883 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2895 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2899 Finalize(refine, fix_orientation);
2902 void Mesh::FinalizeTopology(
bool generate_bdr)
2914 bool generate_edges =
true;
2916 if (spaceDim == 0) { spaceDim = Dim; }
2917 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2921 if (tmp_vertex_parents.Size())
2923 MFEM_VERIFY(ncmesh == NULL,
"");
2924 ncmesh =
new NCMesh(
this);
2928 InitFromNCMesh(*ncmesh);
2929 ncmesh->OnMeshUpdated(
this);
2930 GenerateNCFaceInfo();
2934 tmp_vertex_parents.DeleteAll();
2944 GetElementToFaceTable();
2946 if (NumOfBdrElements == 0 && generate_bdr)
2948 GenerateBoundaryElements();
2949 GetElementToFaceTable();
2958 if (Dim > 1 && generate_edges)
2961 if (!el_to_edge) { el_to_edge =
new Table; }
2962 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2966 if (NumOfBdrElements == 0 && generate_bdr)
2968 GenerateBoundaryElements();
2985 ncmesh->OnMeshUpdated(
this);
2988 GenerateNCFaceInfo();
2995 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2997 if (NURBSext || ncmesh)
2999 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
3000 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
3009 const bool check_orientation =
true;
3010 const bool curved = (Nodes != NULL);
3011 const bool may_change_topology =
3012 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
3013 ( check_orientation && fix_orientation &&
3014 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
3017 Table *old_elem_vert = NULL;
3019 if (curved && may_change_topology)
3021 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3024 if (check_orientation)
3027 CheckElementOrientation(fix_orientation);
3031 MarkForRefinement();
3034 if (may_change_topology)
3038 DoNodeReorder(old_v_to_v, old_elem_vert);
3039 delete old_elem_vert;
3052 CheckBdrElementOrientation();
3057 if (Dim >= 2 && Dim == spaceDim)
3059 const int num_faces = GetNumFaces();
3060 for (
int i = 0; i < num_faces; i++)
3062 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
3063 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology." 3064 " Interior face with incompatible orientations.");
3071 double sx,
double sy,
double sz,
bool sfc_ordering)
3075 int NVert, NElem, NBdrElem;
3077 NVert = (nx+1) * (ny+1) * (nz+1);
3078 NElem = nx * ny * nz;
3079 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3080 if (type == Element::TETRAHEDRON)
3085 else if (type == Element::WEDGE)
3088 NBdrElem += 2*nx*ny;
3090 else if (type == Element::PYRAMID)
3093 NVert += nx * ny * nz;
3096 InitMesh(3, 3, NVert, NElem, NBdrElem);
3102 for (z = 0; z <= nz; z++)
3104 coord[2] = ((double) z / nz) * sz;
3105 for (y = 0; y <= ny; y++)
3107 coord[1] = ((double) y / ny) * sy;
3108 for (x = 0; x <= nx; x++)
3110 coord[0] = ((double) x / nx) * sx;
3115 if (type == Element::PYRAMID)
3117 for (z = 0; z < nz; z++)
3119 coord[2] = (((double) z + 0.5) / nz) * sz;
3120 for (y = 0; y < ny; y++)
3122 coord[1] = (((double) y + 0.5 ) / ny) * sy;
3123 for (x = 0; x < nx; x++)
3125 coord[0] = (((double) x + 0.5 ) / nx) * sx;
3132 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1)) 3133 #define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx) 3136 if (sfc_ordering && type == Element::HEXAHEDRON)
3139 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
3140 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3142 for (
int k = 0; k < nx*ny*nz; k++)
3149 ind[0] = VTX(x , y , z );
3150 ind[1] = VTX(x+1, y , z );
3151 ind[2] = VTX(x+1, y+1, z );
3152 ind[3] = VTX(x , y+1, z );
3153 ind[4] = VTX(x , y , z+1);
3154 ind[5] = VTX(x+1, y , z+1);
3155 ind[6] = VTX(x+1, y+1, z+1);
3156 ind[7] = VTX(x , y+1, z+1);
3164 for (z = 0; z < nz; z++)
3166 for (y = 0; y < ny; y++)
3168 for (x = 0; x < nx; x++)
3171 ind[0] = VTX(x , y , z );
3172 ind[1] = VTX(x+1, y , z );
3173 ind[2] = VTX(x+1, y+1, z );
3174 ind[3] = VTX(x , y+1, z );
3175 ind[4] = VTX(x , y , z+1);
3176 ind[5] = VTX(x+1, y , z+1);
3177 ind[6] = VTX(x+1, y+1, z+1);
3178 ind[7] = VTX( x, y+1, z+1);
3180 if (type == Element::TETRAHEDRON)
3182 AddHexAsTets(ind, 1);
3184 else if (type == Element::WEDGE)
3186 AddHexAsWedges(ind, 1);
3188 else if (type == Element::PYRAMID)
3190 ind[8] = VTXP(x, y, z);
3191 AddHexAsPyramids(ind, 1);
3204 for (y = 0; y < ny; y++)
3206 for (x = 0; x < nx; x++)
3209 ind[0] = VTX(x , y , 0);
3210 ind[1] = VTX(x , y+1, 0);
3211 ind[2] = VTX(x+1, y+1, 0);
3212 ind[3] = VTX(x+1, y , 0);
3214 if (type == Element::TETRAHEDRON)
3216 AddBdrQuadAsTriangles(ind, 1);
3218 else if (type == Element::WEDGE)
3220 AddBdrQuadAsTriangles(ind, 1);
3229 for (y = 0; y < ny; y++)
3231 for (x = 0; x < nx; x++)
3234 ind[0] = VTX(x , y , nz);
3235 ind[1] = VTX(x+1, y , nz);
3236 ind[2] = VTX(x+1, y+1, nz);
3237 ind[3] = VTX(x , y+1, nz);
3239 if (type == Element::TETRAHEDRON)
3241 AddBdrQuadAsTriangles(ind, 6);
3243 else if (type == Element::WEDGE)
3245 AddBdrQuadAsTriangles(ind, 6);
3254 for (z = 0; z < nz; z++)
3256 for (y = 0; y < ny; y++)
3259 ind[0] = VTX(0 , y , z );
3260 ind[1] = VTX(0 , y , z+1);
3261 ind[2] = VTX(0 , y+1, z+1);
3262 ind[3] = VTX(0 , y+1, z );
3264 if (type == Element::TETRAHEDRON)
3266 AddBdrQuadAsTriangles(ind, 5);
3275 for (z = 0; z < nz; z++)
3277 for (y = 0; y < ny; y++)
3280 ind[0] = VTX(nx, y , z );
3281 ind[1] = VTX(nx, y+1, z );
3282 ind[2] = VTX(nx, y+1, z+1);
3283 ind[3] = VTX(nx, y , z+1);
3285 if (type == Element::TETRAHEDRON)
3287 AddBdrQuadAsTriangles(ind, 3);
3296 for (x = 0; x < nx; x++)
3298 for (z = 0; z < nz; z++)
3301 ind[0] = VTX(x , 0, z );
3302 ind[1] = VTX(x+1, 0, z );
3303 ind[2] = VTX(x+1, 0, z+1);
3304 ind[3] = VTX(x , 0, z+1);
3306 if (type == Element::TETRAHEDRON)
3308 AddBdrQuadAsTriangles(ind, 2);
3317 for (x = 0; x < nx; x++)
3319 for (z = 0; z < nz; z++)
3322 ind[0] = VTX(x , ny, z );
3323 ind[1] = VTX(x , ny, z+1);
3324 ind[2] = VTX(x+1, ny, z+1);
3325 ind[3] = VTX(x+1, ny, z );
3327 if (type == Element::TETRAHEDRON)
3329 AddBdrQuadAsTriangles(ind, 4);
3341 ofstream test_stream(
"debug.mesh");
3343 test_stream.close();
3352 double sx,
double sy,
3353 bool generate_edges,
bool sfc_ordering)
3362 if (type == Element::QUADRILATERAL)
3364 NumOfVertices = (nx+1) * (ny+1);
3365 NumOfElements = nx * ny;
3366 NumOfBdrElements = 2 * nx + 2 * ny;
3368 vertices.SetSize(NumOfVertices);
3369 elements.SetSize(NumOfElements);
3370 boundary.SetSize(NumOfBdrElements);
3377 for (j = 0; j < ny+1; j++)
3379 cy = ((double) j / ny) * sy;
3380 for (i = 0; i < nx+1; i++)
3382 cx = ((double) i / nx) * sx;
3383 vertices[k](0) = cx;
3384 vertices[k](1) = cy;
3393 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
3394 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
3396 for (k = 0; k < nx*ny; k++)
3400 ind[0] = i + j*(nx+1);
3401 ind[1] = i + 1 +j*(nx+1);
3402 ind[2] = i + 1 + (j+1)*(nx+1);
3403 ind[3] = i + (j+1)*(nx+1);
3410 for (j = 0; j < ny; j++)
3412 for (i = 0; i < nx; i++)
3414 ind[0] = i + j*(nx+1);
3415 ind[1] = i + 1 +j*(nx+1);
3416 ind[2] = i + 1 + (j+1)*(nx+1);
3417 ind[3] = i + (j+1)*(nx+1);
3426 for (i = 0; i < nx; i++)
3428 boundary[i] =
new Segment(i, i+1, 1);
3429 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3432 for (j = 0; j < ny; j++)
3434 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3435 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3439 else if (type == Element::TRIANGLE)
3441 NumOfVertices = (nx+1) * (ny+1);
3442 NumOfElements = 2 * nx * ny;
3443 NumOfBdrElements = 2 * nx + 2 * ny;
3445 vertices.SetSize(NumOfVertices);
3446 elements.SetSize(NumOfElements);
3447 boundary.SetSize(NumOfBdrElements);
3454 for (j = 0; j < ny+1; j++)
3456 cy = ((double) j / ny) * sy;
3457 for (i = 0; i < nx+1; i++)
3459 cx = ((double) i / nx) * sx;
3460 vertices[k](0) = cx;
3461 vertices[k](1) = cy;
3468 for (j = 0; j < ny; j++)
3470 for (i = 0; i < nx; i++)
3472 ind[0] = i + j*(nx+1);
3473 ind[1] = i + 1 + (j+1)*(nx+1);
3474 ind[2] = i + (j+1)*(nx+1);
3477 ind[1] = i + 1 + j*(nx+1);
3478 ind[2] = i + 1 + (j+1)*(nx+1);
3486 for (i = 0; i < nx; i++)
3488 boundary[i] =
new Segment(i, i+1, 1);
3489 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3492 for (j = 0; j < ny; j++)
3494 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3495 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3502 MFEM_ABORT(
"Unsupported element type.");
3506 CheckElementOrientation();
3508 if (generate_edges == 1)
3510 el_to_edge =
new Table;
3511 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3513 CheckBdrElementOrientation();
3522 attributes.Append(1);
3523 bdr_attributes.Append(1); bdr_attributes.Append(2);
3524 bdr_attributes.Append(3); bdr_attributes.Append(4);
3529 void Mesh::Make1D(
int n,
double sx)
3538 NumOfVertices = n + 1;
3540 NumOfBdrElements = 2;
3541 vertices.SetSize(NumOfVertices);
3542 elements.SetSize(NumOfElements);
3543 boundary.SetSize(NumOfBdrElements);
3546 for (j = 0; j < n+1; j++)
3548 vertices[j](0) = ((
double) j / n) * sx;
3552 for (j = 0; j < n; j++)
3554 elements[j] =
new Segment(j, j+1, 1);
3559 boundary[0] =
new Point(ind, 1);
3561 boundary[1] =
new Point(ind, 2);
3569 attributes.Append(1);
3570 bdr_attributes.Append(1); bdr_attributes.Append(2);
3573 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
3591 last_operation = Mesh::NONE;
3594 elements.SetSize(NumOfElements);
3595 for (
int i = 0; i < NumOfElements; i++)
3597 elements[i] = mesh.
elements[i]->Duplicate(
this);
3604 boundary.SetSize(NumOfBdrElements);
3605 for (
int i = 0; i < NumOfBdrElements; i++)
3607 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
3626 faces.SetSize(mesh.
faces.Size());
3627 for (
int i = 0; i < faces.Size(); i++)
3630 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
3640 face_to_elem = NULL;
3665 if (dynamic_cast<const ParMesh*>(&mesh))
3677 if (mesh.
Nodes && copy_nodes)
3682 FiniteElementCollection::New(fec->
Name());
3686 Nodes->MakeOwner(fec_copy);
3687 *Nodes = *mesh.
Nodes;
3709 bool fix_orientation)
3713 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
3714 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
3728 double sx,
double sy,
bool sfc_ordering)
3731 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
3738 double sx,
double sy,
double sz,
bool sfc_ordering)
3741 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
3750 ref_factors = ref_factor;
3764 bool fix_orientation)
3773 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
3777 Load(imesh, generate_edges, refine, fix_orientation);
3782 bool fix_orientation)
3785 Load(input, generate_edges, refine, fix_orientation);
3794 "Not enough vertices in external array : " 3795 "len_vertex_data = "<< len_vertex_data <<
", " 3798 if (vertex_data == (
double *)(
vertices.GetData()))
3800 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
3805 memcpy(vertex_data,
vertices.GetData(),
3814 int *element_attributes,
int num_elements,
3816 int *boundary_attributes,
int num_boundary_elements,
3819 if (space_dimension == -1)
3825 num_boundary_elements);
3828 int boundary_index_stride = num_boundary_elements > 0 ?
3832 vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
3835 for (
int i = 0; i < num_elements; i++)
3838 elements[i]->SetVertices(element_indices + i * element_index_stride);
3839 elements[i]->SetAttribute(element_attributes[i]);
3843 for (
int i = 0; i < num_boundary_elements; i++)
3846 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3847 boundary[i]->SetAttribute(boundary_attributes[i]);
3863 #ifdef MFEM_USE_MEMALLOC 3872 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3885 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3888 for (
int i = 0; i < nv; i++)
3901 for (
int j = 0; j < nv; j++)
3973 MFEM_ABORT(
"invalid element type: " << type);
3980 std::string parse_tag)
3982 int curved = 0, read_gf = 1;
3983 bool finalize_topo =
true;
3987 MFEM_ABORT(
"Input stream is not open");
3994 getline(input, mesh_type);
3998 int mfem_version = 0;
3999 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
4000 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
4004 int mfem_nc_version = 0;
4005 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
4006 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
4014 if (mfem_version == 12 && parse_tag.empty())
4016 parse_tag =
"mfem_mesh_end";
4020 else if (mfem_nc_version)
4022 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4029 MFEM_VERIFY(mfem_nc_version >= 10,
4030 "Legacy nonconforming format (MFEM mesh v1.1) cannot be " 4031 "used to load a parallel nonconforming mesh, sorry.");
4034 input, mfem_nc_version, curved, is_nc);
4039 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
4052 else if (mesh_type ==
"linemesh")
4056 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
4058 if (mesh_type ==
"curved_areamesh2")
4064 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
4068 else if (mesh_type ==
"TrueGrid")
4072 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
4074 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
4076 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4077 "Unsupported VTK format");
4078 ReadVTKMesh(input, curved, read_gf, finalize_topo);
4080 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
4084 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
4088 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
4093 else if (mesh_type ==
"$MeshFormat")
4098 ((mesh_type.size() > 2 &&
4099 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
4100 (mesh_type.size() > 3 &&
4101 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
4106 #ifdef MFEM_USE_NETCDF 4109 MFEM_ABORT(
"NetCDF support requires configuration with" 4110 " MFEM_USE_NETCDF=YES");
4116 MFEM_ABORT(
"Can not determine Cubit mesh filename!" 4117 " Use mfem::named_ifgzstream for input.");
4123 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
4149 bool generate_bdr =
false;
4154 if (curved && read_gf)
4168 if (mfem_version == 12)
4174 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
4175 getline(input, line);
4181 if (line ==
"mfem_mesh_end") {
break; }
4183 while (line != parse_tag);
4185 else if (mfem_nc_version >= 10)
4190 MFEM_VERIFY(ident ==
"mfem_mesh_end",
4191 "invalid mesh: end of file tag not found");
4199 int i, j, ie, ib, iv, *v, nv;
4230 for (i = 0; i < num_pieces; i++)
4237 for (i = 0; i < num_pieces; i++)
4243 for (j = 0; j < m->
GetNE(); j++)
4248 for (j = 0; j < m->
GetNBE(); j++)
4253 for (
int k = 0; k < nv; k++)
4255 v[k] = lvert_vert[v[k]];
4260 for (j = 0; j < m->
GetNV(); j++)
4272 for (i = 0; i < num_pieces; i++)
4283 for (i = 0; i < num_pieces; i++)
4287 for (j = 0; j < m->
GetNE(); j++)
4292 for (
int k = 0; k < nv; k++)
4299 for (j = 0; j < m->
GetNBE(); j++)
4304 for (
int k = 0; k < nv; k++)
4311 for (j = 0; j < m->
GetNV(); j++)
4325 for (i = 0; i < num_pieces; i++)
4327 gf_array[i] = mesh_array[i]->
GetNodes();
4342 ref_factors = ref_factor;
4353 int orig_ne = orig_mesh.
GetNE();
4354 MFEM_VERIFY(ref_factors.
Size() == orig_ne && orig_ne > 0,
4355 "Number of refinement factors must equal number of elements")
4356 MFEM_VERIFY(ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
4359 "Invalid refinement type. Must use closed basis type.");
4361 int min_ref = ref_factors.
Min();
4362 int max_ref = ref_factors.
Max();
4364 bool var_order = (min_ref != max_ref);
4377 for (
int i = 0; i < orig_ne; i++)
4395 for (
int el = 0; el < orig_ne; el++)
4407 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
4408 for (
int i = 0; i < phys_pts.
Width(); i++)
4417 for (
int k = 0; k < nvert; k++)
4420 v[k] = rdofs[c2h_map[cid]];
4433 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
4444 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
4450 for (
int k = 0; k < nvert; k++)
4453 v[k] = rdofs[c2h_map[cid]];
4475 for (
int iel = 0; iel < orig_ne; iel++)
4484 const int *node_map = NULL;
4487 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
4488 const int *vertex_map = vertex_fec.
GetDofMap(geom);
4489 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
4490 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
4493 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
4496 int iv = vertex_map[iv_lex];
4498 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
4500 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
4503 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
4514 using GeomRef = std::pair<Geometry::Type, int>;
4515 std::map<GeomRef, int> point_matrices_offsets;
4517 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4521 GeomRef id(geom, ref_factors[el_coarse]);
4522 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
4528 point_matrices_offsets[id] = n_point_matrices[geom];
4529 n_point_matrices[geom] += nref_el;
4536 int nmatrices = n_point_matrices[geom];
4543 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4546 int ref = ref_factors[el_coarse];
4547 int offset = point_matrices_offsets[GeomRef(geom, ref)];
4553 for (
int k = 0; k < nvert; k++)
4585 "Mesh::MakeSimplicial requires a properly oriented input mesh");
4587 "Mesh::MakeSimplicial does not support non-conforming meshes.")
4594 Mesh copy(orig_mesh);
4599 int nv = orig_mesh.
GetNV();
4600 int ne = orig_mesh.
GetNE();
4601 int nbe = orig_mesh.
GetNBE();
4614 int new_ne = 0, new_nbe = 0;
4615 for (
int i=0; i<ne; ++i)
4619 for (
int i=0; i<nbe; ++i)
4628 for (
int i=0; i<nv; ++i)
4638 if (vglobal == NULL)
4641 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
4642 vglobal = vglobal_id.
GetData();
4645 constexpr
int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
4646 constexpr
int quad_ntris = 2, prism_ntets = 3;
4647 static const int quad_trimap[2][nv_tri*quad_ntris] =
4659 static const int prism_rot[nv_prism*nv_prism] =
4668 static const int prism_f[nv_quad] = {1, 2, 5, 4};
4669 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
4683 static const int hex_rot[nv_hex*nv_hex] =
4685 0, 1, 2, 3, 4, 5, 6, 7,
4686 1, 0, 4, 5, 2, 3, 7, 6,
4687 2, 1, 5, 6, 3, 0, 4, 7,
4688 3, 0, 1, 2, 7, 4, 5, 6,
4689 4, 0, 3, 7, 5, 1, 2, 6,
4690 5, 1, 0, 4, 6, 2, 3, 7,
4691 6, 2, 1, 5, 7, 3, 0, 4,
4692 7, 3, 2, 6, 4, 0, 1, 5
4694 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
4695 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
4696 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
4697 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
4698 static const int hex_tetmap0[nv_tet*5] =
4705 static const int hex_tetmap1[nv_tet*6] =
4712 static const int hex_tetmap2[nv_tet*6] =
4719 static const int hex_tetmap3[nv_tet*6] =
4726 static const int *hex_tetmaps[4] =
4728 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
4731 auto find_min = [](
const int*
a,
int n) {
return std::min_element(
a,
a+n)-
a; };
4733 for (
int i=0; i<ne; ++i)
4735 const int *v = orig_mesh.
elements[i]->GetVertices();
4739 if (num_subdivisions[orig_geom] == 1)
4751 for (
int itri=0; itri<quad_ntris; ++itri)
4756 for (
int iv=0; iv<nv_tri; ++iv)
4758 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
4766 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
4769 int irot = find_min(vg, nv_prism);
4770 for (
int iv=0; iv<nv_prism; ++iv)
4772 int jv = prism_rot[iv + irot*nv_prism];
4777 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
4778 int j = find_min(q, nv_quad);
4779 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
4780 for (
int itet=0; itet<prism_ntets; ++itet)
4785 for (
int iv=0; iv<nv_tet; ++iv)
4787 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
4795 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
4799 int irot = find_min(vg, nv_hex);
4800 for (
int iv=0; iv<nv_hex; ++iv)
4802 int jv = hex_rot[iv + irot*nv_hex];
4812 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
4813 j = find_min(q, nv_quad);
4814 if (j == 0 || j == 2) { bitmask += 4; }
4816 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
4817 j = find_min(q, nv_quad);
4818 if (j == 1 || j == 3) { bitmask += 2; }
4820 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
4821 j = find_min(q, nv_quad);
4822 if (j == 0 || j == 2) { bitmask += 1; }
4825 int nrot = num_rot[bitmask];
4826 for (
int k=0; k<nrot; ++k)
4840 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
4841 int ntets = (ndiags == 0) ? 5 : 6;
4842 const int *tetmap = hex_tetmaps[ndiags];
4843 for (
int itet=0; itet<ntets; ++itet)
4848 for (
int iv=0; iv<nv_tet; ++iv)
4850 v2[iv] = vg[tetmap[itet + iv*ntets]];
4859 for (
int i=0; i<nbe; ++i)
4861 const int *v = orig_mesh.
boundary[i]->GetVertices();
4864 if (num_subdivisions[orig_geom] == 1)
4874 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
4876 int iv_min = find_min(vg, nv_quad);
4877 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
4878 for (
int itri=0; itri<quad_ntris; ++itri)
4883 for (
int iv=0; iv<nv_tri; ++iv)
4885 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
4892 MFEM_ABORT(
"Unreachable");
4906 Mesh periodic_mesh(orig_mesh,
true);
4912 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
4917 for (
int j = 0; j < nv; j++)
4923 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
4928 for (
int j = 0; j < nv; j++)
4935 return periodic_mesh;
4939 const std::vector<Vector> &translations,
double tol)
const 4943 Vector coord(sdim), at(sdim), dx(sdim);
4944 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
4945 xMax = xMin = xDiff = 0.0;
4949 for (
int be = 0; be <
GetNBE(); be++)
4954 for (
int i = 0; i < dofs.
Size(); i++)
4956 bdr_v.insert(dofs[i]);
4959 for (
int j = 0; j < sdim; j++)
4961 xMax[j] = max(xMax[j], coord[j]);
4962 xMin[j] = min(xMin[j], coord[j]);
4966 add(xMax, -1.0, xMin, xDiff);
4967 double dia = xDiff.
Norml2();
4975 map<int, int> replica2primary;
4977 map<int, set<int>> primary2replicas;
4981 for (
int v : bdr_v) { primary2replicas[v]; }
4985 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
4987 if (r ==
p) {
return; }
4988 primary2replicas[
p].insert(r);
4989 replica2primary[r] =
p;
4990 for (
int s : primary2replicas[r])
4992 primary2replicas[
p].insert(
s);
4993 replica2primary[
s] =
p;
4995 primary2replicas.erase(r);
4998 for (
unsigned int i = 0; i < translations.size(); i++)
5000 for (
int vi : bdr_v)
5003 add(coord, translations[i], at);
5005 for (
int vj : bdr_v)
5008 add(at, -1.0, coord, dx);
5009 if (dx.
Norml2() > dia*tol) {
continue; }
5014 bool pi = primary2replicas.find(vi) != primary2replicas.end();
5015 bool pj = primary2replicas.find(vj) != primary2replicas.end();
5021 make_replica(vj, vi);
5026 int owner_of_vj = replica2primary[vj];
5028 make_replica(vi, owner_of_vj);
5034 int owner_of_vi = replica2primary[vi];
5035 make_replica(vj, owner_of_vi);
5042 int owner_of_vi = replica2primary[vi];
5043 int owner_of_vj = replica2primary[vj];
5044 make_replica(owner_of_vj, owner_of_vi);
5051 std::vector<int> v2v(
GetNV());
5052 for (
size_t i = 0; i < v2v.size(); i++)
5056 for (
auto &&r2p : replica2primary)
5058 v2v[r2p.first] = r2p.second;
5067 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5072 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5089 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5094 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5124 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
5148 for (
int i = 0; i <
elements.Size(); i++)
5158 for (
int i = 0; i <
boundary.Size(); i++)
5175 for (
int i = 0; i < vd; i++)
5244 input >> edge_to_knot[j] >> v[0] >> v[1];
5247 edge_to_knot[j] = -1 - edge_to_knot[j];
5263 if (
p.Size() >= v.
Size())
5265 for (
int d = 0; d < v.
Size(); d++)
5273 for (d = 0; d <
p.Size(); d++)
5277 for ( ; d < v.
Size(); d++)
5309 if (dynamic_cast<const H1_FECollection*>(fec)
5310 || dynamic_cast<const L2_FECollection*>(fec))
5319 #ifndef MFEM_USE_MPI 5320 const bool warn =
true;
5323 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
5327 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
5328 " H1-continuous mesh!\n " 5329 "If this is the desired behavior, you can silence" 5330 " this warning by converting\n " 5331 "the NURBS mesh to high-order mesh in advance by" 5332 " calling the method\n " 5333 "Mesh::SetCurvature().");
5358 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
5377 MFEM_ASSERT(nodes != NULL,
"");
5393 case 1:
return GetNV();
5429 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG)) 5430 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
5435 int i, j, k, wo = 0, fo = 0;
5444 int *vi =
elements[i]->GetVertices();
5447 for (j = 0; j < 3; j++)
5451 for (j = 0; j < 2; j++)
5452 for (k = 0; k < 2; k++)
5454 J(j, k) = v[j+1][k] - v[0][k];
5475 MFEM_ABORT(
"Invalid 2D element type \"" 5492 int *vi =
elements[i]->GetVertices();
5498 for (j = 0; j < 4; j++)
5502 for (j = 0; j < 3; j++)
5503 for (k = 0; k < 3; k++)
5505 J(j, k) = v[j+1][k] - v[0][k];
5564 MFEM_ABORT(
"Invalid 3D element type \"" 5570 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG)) 5573 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / " 5574 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
5589 if (test[0] == base[0])
5590 if (test[1] == base[1])
5598 else if (test[0] == base[1])
5599 if (test[1] == base[0])
5608 if (test[1] == base[0])
5619 for (
int j = 0; j < 3; j++)
5620 if (test[aor[j]] != base[j])
5633 for (i = 0; i < 4; i++)
5634 if (test[i] == base[0])
5641 if (test[(i+1)%4] == base[1])
5650 for (
int j = 0; j < 4; j++)
5651 if (test[aor[j]] != base[j])
5653 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
5655 for (
int k = 0; k < 4; k++)
5660 for (
int k = 0; k < 4; k++)
5669 if (test[(i+1)%4] == base[1])
5685 if (test[0] == base[0])
5686 if (test[1] == base[1])
5687 if (test[2] == base[2])
5695 else if (test[2] == base[1])
5696 if (test[3] == base[2])
5705 if (test[1] == base[2])
5713 else if (test[1] == base[0])
5714 if (test[2] == base[1])
5715 if (test[0] == base[2])
5723 else if (test[3] == base[1])
5724 if (test[2] == base[2])
5733 if (test[3] == base[2])
5741 else if (test[2] == base[0])
5742 if (test[3] == base[1])
5743 if (test[0] == base[2])
5751 else if (test[0] == base[1])
5752 if (test[1] == base[2])
5761 if (test[3] == base[2])
5770 if (test[0] == base[1])
5771 if (test[2] == base[2])
5779 else if (test[1] == base[1])
5780 if (test[0] == base[2])
5789 if (test[1] == base[2])
5800 for (
int j = 0; j < 4; j++)
5801 if (test[aor[j]] != base[j])
5826 int *bv =
boundary[i]->GetVertices();
5832 mfem::Swap<int>(bv[0], bv[1]);
5846 if (
faces_info[fi].Elem2No >= 0) {
continue; }
5849 int *bv =
boundary[i]->GetVertices();
5851 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
5852 const int *fv =
faces[fi]->GetVertices();
5869 MFEM_ABORT(
"Invalid 2D boundary element type \"" 5870 << bdr_type <<
"\"");
5875 if (orientation % 2 == 0) {
continue; }
5877 if (!fix_it) {
continue; }
5885 mfem::Swap<int>(bv[0], bv[1]);
5889 mfem::Swap<int>(be[1], be[2]);
5895 mfem::Swap<int>(bv[0], bv[2]);
5899 mfem::Swap<int>(be[0], be[1]);
5900 mfem::Swap<int>(be[2], be[3]);
5913 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / " 5923 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
5934 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
5953 mfem_error(
"Mesh::GetElementEdges(...) element to edge table " 5954 "is not generated.");
5957 const int *v =
elements[i]->GetVertices();
5958 const int ne =
elements[i]->GetNEdges();
5960 for (
int j = 0; j < ne; j++)
5962 const int *e =
elements[i]->GetEdgeVertices(j);
5963 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5974 const int *v =
boundary[i]->GetVertices();
5975 cor[0] = (v[0] < v[1]) ? (1) : (-1);
5988 const int *v =
boundary[i]->GetVertices();
5989 const int ne =
boundary[i]->GetNEdges();
5991 for (
int j = 0; j < ne; j++)
5993 const int *e =
boundary[i]->GetEdgeVertices(j);
5994 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6006 const int *v =
faces[i]->GetVertices();
6007 o[0] = (v[0] < v[1]) ? (1) : (-1);
6019 const int *v =
faces[i]->GetVertices();
6020 const int ne =
faces[i]->GetNEdges();
6022 for (
int j = 0; j < ne; j++)
6024 const int *e =
faces[i]->GetEdgeVertices(j);
6025 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6053 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
6104 for (j = 0; j < nv; j++)
6116 for (j = 0; j < nv; j++)
6163 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
6167 int n = el_faces.
Size();
6170 for (
int j = 0; j < n; j++)
6174 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
6178 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
6179 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
6196 for (
auto f : elem_faces)
6231 MFEM_ABORT(
"invalid geometry");
6239 case 1:
return boundary[i]->GetVertices()[0];
6242 default: MFEM_ABORT(
"invalid dimension!");
6252 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6255 const int *bv =
boundary[bdr_el]->GetVertices();
6263 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6274 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6277 const int *bv =
boundary[bdr_el]->GetVertices();
6285 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6312 for (j = 0; j < nv; j++)
6314 pointmat(k, j) =
vertices[v[j]](k);
6329 for (j = 0; j < nv; j++)
6331 pointmat(k, j) =
vertices[v[j]](k);
6343 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
6346 return sqrt(length);
6354 for (
int i = 0; i < elem_array.
Size(); i++)
6359 for (
int i = 0; i < elem_array.
Size(); i++)
6361 const int *v = elem_array[i]->GetVertices();
6362 const int ne = elem_array[i]->GetNEdges();
6363 for (
int j = 0; j < ne; j++)
6365 const int *e = elem_array[i]->GetEdgeVertices(j);
6379 v_to_v.
Push(v[0], v[1]);
6386 const int *v =
elements[i]->GetVertices();
6387 const int ne =
elements[i]->GetNEdges();
6388 for (
int j = 0; j < ne; j++)
6390 const int *e =
elements[i]->GetEdgeVertices(j);
6391 v_to_v.
Push(v[e[0]], v[e[1]]);
6399 int i, NumberOfEdges;
6415 const int *v =
boundary[i]->GetVertices();
6416 be_to_f[i] = v_to_v(v[0], v[1]);
6429 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
6433 return NumberOfEdges;
6524 if (
faces[gf] == NULL)
6534 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. " 6535 "Interior edge found between 2D elements " 6537 <<
" and " << el <<
".");
6538 int *v =
faces[gf]->GetVertices();
6540 if ( v[1] == v0 && v[0] == v1 )
6544 else if ( v[0] == v0 && v[1] == v1 )
6554 MFEM_ABORT(
"internal error");
6560 int v0,
int v1,
int v2)
6562 if (
faces[gf] == NULL)
6572 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. " 6573 "Interior triangular face found connecting elements " 6575 <<
" and " << el <<
".");
6576 int orientation, vv[3] = { v0, v1, v2 };
6583 faces_info[gf].Elem2Inf = 64 * lf + orientation;
6588 int v0,
int v1,
int v2,
int v3)
6600 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. " 6601 "Interior quadrilateral face found connecting elements " 6603 <<
" and " << el <<
".");
6604 int vv[4] = { v0, v1, v2, v3 };
6618 for (i = 0; i <
faces.Size(); i++)
6624 faces.SetSize(nfaces);
6626 for (i = 0; i < nfaces; i++)
6634 const int *v =
elements[i]->GetVertices();
6644 const int ne =
elements[i]->GetNEdges();
6645 for (
int j = 0; j < ne; j++)
6647 const int *e =
elements[i]->GetEdgeVertices(j);
6658 for (
int j = 0; j < 4; j++)
6662 v[fv[0]], v[fv[1]], v[fv[2]]);
6668 for (
int j = 0; j < 2; j++)
6672 v[fv[0]], v[fv[1]], v[fv[2]]);
6674 for (
int j = 2; j < 5; j++)
6678 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6684 for (
int j = 0; j < 1; j++)
6688 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6690 for (
int j = 1; j < 5; j++)
6694 v[fv[0]], v[fv[1]], v[fv[2]]);
6700 for (
int j = 0; j < 6; j++)
6704 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6709 MFEM_ABORT(
"Unexpected type of Element.");
6717 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
6733 for (
int i = 0; i < list.
masters.Size(); i++)
6736 if (master.
index >= nfaces) {
continue; }
6742 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
6743 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
6747 for (
int i = 0; i < list.
slaves.Size(); i++)
6751 if (slave.
index < 0 ||
6752 slave.
index >= nfaces ||
6781 const int *v =
elements[i]->GetVertices();
6786 for (
int j = 0; j < 4; j++)
6789 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6795 for (
int j = 0; j < 1; j++)
6798 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6800 for (
int j = 1; j < 5; j++)
6803 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6809 for (
int j = 0; j < 2; j++)
6812 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6814 for (
int j = 2; j < 5; j++)
6817 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6825 for (
int j = 0; j < 6; j++)
6828 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6833 MFEM_ABORT(
"Unexpected type of Element.");
6857 for (
int j = 0; j < 4; j++)
6861 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6867 for (
int j = 0; j < 2; j++)
6871 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6873 for (
int j = 2; j < 5; j++)
6877 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6883 for (
int j = 0; j < 1; j++)
6887 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6889 for (
int j = 1; j < 5; j++)
6893 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6901 for (
int j = 0; j < 6; j++)
6905 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6910 MFEM_ABORT(
"Unexpected type of Element.");
6923 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
6928 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
6932 MFEM_ABORT(
"Unexpected type of boundary Element.");
6946 void Rotate3(
int &
a,
int &
b,
int &c)
6978 Table *old_elem_vert = NULL;
6989 int *v =
elements[i]->GetVertices();
6991 Rotate3(v[0], v[1], v[2]);
6994 Rotate3(v[1], v[2], v[3]);
7007 int *v =
boundary[i]->GetVertices();
7009 Rotate3(v[0], v[1], v[2]);
7025 delete old_elem_vert;
7041 if (
p[i] < pmin[i]) { pmin[i] =
p[i]; }
7042 if (
p[i] > pmax[i]) { pmax[i] =
p[i]; }
7056 for (
int i =
spaceDim-1; i >= 0; i--)
7058 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
7059 if (idx < 0) { idx = 0; }
7060 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
7061 part = part * nxyz[i] + idx;
7063 partitioning[el] = part;
7066 return partitioning;
7076 #ifdef MFEM_USE_METIS 7078 int print_messages = 1;
7081 int init_flag, fin_flag;
7082 MPI_Initialized(&init_flag);
7083 MPI_Finalized(&fin_flag);
7084 if (init_flag && !fin_flag)
7088 if (rank != 0) { print_messages = 0; }
7092 int i, *partitioning;
7102 partitioning[i] = 0;
7109 partitioning[i] = i;
7115 #ifndef MFEM_USE_METIS_5 7127 bool freedata =
false;
7129 idx_t *mpartitioning;
7132 if (
sizeof(
idx_t) ==
sizeof(int))
7136 mpartitioning = (
idx_t*) partitioning;
7145 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
7146 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
7147 mpartitioning =
new idx_t[n];
7150 #ifndef MFEM_USE_METIS_5 7153 METIS_SetDefaultOptions(options);
7154 options[METIS_OPTION_CONTIG] = 1;
7162 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
7167 if (part_method >= 0 && part_method <= 2)
7169 for (i = 0; i < n; i++)
7175 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
7181 if (part_method == 0 || part_method == 3)
7183 #ifndef MFEM_USE_METIS_5 7212 " error in METIS_PartGraphRecursive!");
7219 if (part_method == 1 || part_method == 4)
7221 #ifndef MFEM_USE_METIS_5 7250 " error in METIS_PartGraphKway!");
7257 if (part_method == 2 || part_method == 5)
7259 #ifndef MFEM_USE_METIS_5 7272 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
7289 " error in METIS_PartGraphKway!");
7297 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = " 7301 nparts = (int) mparts;
7302 if (mpartitioning != (
idx_t*)partitioning)
7306 partitioning[k] = mpartitioning[k];
7313 delete[] mpartitioning;
7329 auto count_partition_elements = [&]()
7331 for (i = 0; i < nparts; i++)
7339 psize[partitioning[i]].one++;
7343 for (i = 0; i < nparts; i++)
7345 if (psize[i].one == 0) { empty_parts++; }
7349 count_partition_elements();
7357 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned " 7358 << empty_parts <<
" empty parts!" 7359 <<
" Applying a simple fix ..." << endl;
7362 SortPairs<int,int>(psize, nparts);
7364 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7371 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7373 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
7379 partitioning[j] = psize[nparts-1-i].two;
7386 count_partition_elements();
7390 return partitioning;
7394 mfem_error(
"Mesh::GeneratePartitioning(...): " 7395 "MFEM was compiled without Metis.");
7409 int num_elem, *i_elem_elem, *j_elem_elem;
7411 num_elem = elem_elem.
Size();
7412 i_elem_elem = elem_elem.
GetI();
7413 j_elem_elem = elem_elem.
GetJ();
7418 int stack_p, stack_top_p, elem;
7422 for (i = 0; i < num_elem; i++)
7424 if (partitioning[i] > num_part)
7426 num_part = partitioning[i];
7433 for (i = 0; i < num_part; i++)
7440 for (elem = 0; elem < num_elem; elem++)
7442 if (component[elem] >= 0)
7447 component[elem] = num_comp[partitioning[elem]]++;
7449 elem_stack[stack_top_p++] = elem;
7451 for ( ; stack_p < stack_top_p; stack_p++)
7453 i = elem_stack[stack_p];
7454 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
7457 if (partitioning[k] == partitioning[i])
7459 if (component[k] < 0)
7461 component[k] = component[i];
7462 elem_stack[stack_top_p++] = k;
7464 else if (component[k] != component[i])
7476 int i, n_empty, n_mcomp;
7484 n_empty = n_mcomp = 0;
7485 for (i = 0; i < num_comp.
Size(); i++)
7486 if (num_comp[i] == 0)
7490 else if (num_comp[i] > 1)
7497 mfem::out <<
"Mesh::CheckPartitioning(...) :\n" 7498 <<
"The following subdomains are empty :\n";
7499 for (i = 0; i < num_comp.
Size(); i++)
7500 if (num_comp[i] == 0)
7508 mfem::out <<
"Mesh::CheckPartitioning(...) :\n" 7509 <<
"The following subdomains are NOT connected :\n";
7510 for (i = 0; i < num_comp.
Size(); i++)
7511 if (num_comp[i] > 1)
7517 if (n_empty == 0 && n_mcomp == 0)
7518 mfem::out <<
"Mesh::CheckPartitioning(...) : " 7519 "All subdomains are connected." << endl;
7533 const double *
a = A.
Data();
7534 const double *
b = B.
Data();
7543 c(0) =
a[0]*
a[3]-
a[1]*
a[2];
7544 c(1) =
a[0]*
b[3]-
a[1]*
b[2]+
b[0]*
a[3]-
b[1]*
a[2];
7545 c(2) =
b[0]*
b[3]-
b[1]*
b[2];
7566 c(0) = (
a[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
7567 a[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
7568 a[2] * (
a[3] *
a[7] -
a[4] *
a[6]));
7570 c(1) = (
b[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
7571 b[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
7572 b[2] * (
a[3] *
a[7] -
a[4] *
a[6]) +
7574 a[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
7575 a[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
7576 a[2] * (
b[3] *
a[7] -
b[4] *
a[6]) +
7578 a[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
7579 a[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
7580 a[2] * (
a[3] *
b[7] -
a[4] *
b[6]));
7582 c(2) = (
a[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
7583 a[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
7584 a[2] * (
b[3] *
b[7] -
b[4] *
b[6]) +
7586 b[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
7587 b[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
7588 b[2] * (
a[3] *
b[7] -
a[4] *
b[6]) +
7590 b[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
7591 b[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
7592 b[2] * (
b[3] *
a[7] -
b[4] *
a[6]));
7594 c(3) = (
b[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
7595 b[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
7596 b[2] * (
b[3] *
b[7] -
b[4] *
b[6]));
7642 double a = z(2),
b = z(1), c = z(0);
7643 double D =
b*
b-4*
a*c;
7650 x(0) = x(1) = -0.5 *
b /
a;
7655 x(0) = -(x(1) = fabs(0.5 * sqrt(D) /
a));
7663 t = -0.5 * (
b + sqrt(D));
7667 t = -0.5 * (
b - sqrt(D));
7673 Swap<double>(x(0), x(1));
7681 double a = z(2)/z(3),
b = z(1)/z(3), c = z(0)/z(3);
7684 double Q = (
a *
a - 3 *
b) / 9;
7685 double R = (2 *
a *
a *
a - 9 *
a *
b + 27 * c) / 54;
7686 double Q3 = Q * Q * Q;
7693 x(0) = x(1) = x(2) = -
a / 3;
7697 double sqrtQ = sqrt(Q);
7701 x(0) = -2 * sqrtQ -
a / 3;
7702 x(1) = x(2) = sqrtQ -
a / 3;
7706 x(0) = x(1) = - sqrtQ -
a / 3;
7707 x(2) = 2 * sqrtQ -
a / 3;
7714 double theta = acos(R / sqrt(Q3));
7715 double A = -2 * sqrt(Q);
7717 x0 = A * cos(theta / 3) -
a / 3;
7718 x1 = A * cos((theta + 2.0 * M_PI) / 3) -
a / 3;
7719 x2 = A * cos((theta - 2.0 * M_PI) / 3) -
a / 3;
7724 Swap<double>(x0, x1);
7728 Swap<double>(x1, x2);
7731 Swap<double>(x0, x1);
7744 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
7748 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
7750 x(0) = A + Q / A -
a / 3;
7759 const double factor,
const int Dim)
7761 const double c0 = c(0);
7762 c(0) = c0 * (1.0 - pow(factor, -Dim));
7764 for (
int j = 0; j < nr; j++)
7776 c(0) = c0 * (1.0 - pow(factor, Dim));
7778 for (
int j = 0; j < nr; j++)
7797 const double factor = 2.0;
7812 for (
int k = 0; k < nv; k++)
7815 V(j, k) = displacements(v[k]+j*nvs);
7845 for (
int j = 0; j < nv; j++)
7871 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7874 vertices[i](j) += displacements(j*nv+i);
7882 for (
int i = 0; i < nv; i++)
7885 vert_coord(j*nv+i) =
vertices[i](j);
7891 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7894 vertices[i](j) = vert_coord(j*nv+i);
7905 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
7924 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
7941 (*Nodes) += displacements;
7953 node_coord = (*Nodes);
7965 (*Nodes) = node_coord;
8000 mfem::Swap<GridFunction*>(
Nodes, nodes);
8020 for (j = 1; j < n; j++)
8057 int quad_counter = 0;
8072 vertices.SetSize(oelem + quad_counter);
8079 const int attr =
elements[i]->GetAttribute();
8080 int *v =
elements[i]->GetVertices();
8086 for (
int ei = 0; ei < 3; ei++)
8088 for (
int k = 0; k < 2; k++)
8096 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8098 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8100 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8102 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8106 const int qe = quad_counter;
8110 for (
int ei = 0; ei < 4; ei++)
8112 for (
int k = 0; k < 2; k++)
8120 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
8122 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
8124 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
8126 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
8130 MFEM_ABORT(
"unknown element type: " << el_type);
8140 const int attr =
boundary[i]->GetAttribute();
8141 int *v =
boundary[i]->GetVertices();
8150 static const double A = 0.0, B = 0.5, C = 1.0;
8151 static double tri_children[2*3*4] =
8158 static double quad_children[2*4*4] =
8172 for (
int i = 0; i <
elements.Size(); i++)
8193 if (!
Nodes || update_nodes)
8201 static inline double sqr(
const double &x)
8223 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
8226 int NumOfQuadFaces = 0;
8232 for (
int i = 0; i <
faces.Size(); i++)
8236 f2qf[i] = NumOfQuadFaces;
8243 NumOfQuadFaces =
faces.Size();
8247 int hex_counter = 0;
8250 for (
int i = 0; i <
elements.Size(); i++)
8259 int pyr_counter = 0;
8262 for (
int i = 0; i <
elements.Size(); i++)
8280 DSTable *v_to_v_ptr = v_to_v_p;
8296 std::sort(row_start, J_v2v.
end());
8299 for (
int i = 0; i < J_v2v.
Size(); i++)
8301 e2v[J_v2v[i].two] = i;
8314 it.SetIndex(e2v[it.Index()]);
8324 const int oelem = oface + NumOfQuadFaces;
8329 vertices.SetSize(oelem + hex_counter);
8337 const int attr =
elements[i]->GetAttribute();
8338 int *v =
elements[i]->GetVertices();
8345 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8353 for (
int ei = 0; ei < 6; ei++)
8355 for (
int k = 0; k < 2; k++)
8365 const int rt_algo = 1;
8376 double len_sqr, min_len;
8378 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
8379 sqr(J(1,0)-J(1,1)-J(1,2)) +
8380 sqr(J(2,0)-J(2,1)-J(2,2));
8383 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
8384 sqr(J(1,1)-J(1,0)-J(1,2)) +
8385 sqr(J(2,1)-J(2,0)-J(2,2));
8386 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
8388 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
8389 sqr(J(1,2)-J(1,0)-J(1,1)) +
8390 sqr(J(2,2)-J(2,0)-J(2,1));
8391 if (len_sqr < min_len) { rt = 2; }
8396 double Em_data[18], Js_data[9], Jp_data[9];
8399 double ar1, ar2,
kappa, kappa_min;
8401 for (
int s = 0;
s < 3;
s++)
8403 for (
int t = 0;
t < 3;
t++)
8405 Em(
t,
s) = 0.5*J(
t,
s);
8408 for (
int t = 0;
t < 3;
t++)
8410 Em(
t,3) = 0.5*(J(
t,0)+J(
t,1));
8411 Em(
t,4) = 0.5*(J(
t,0)+J(
t,2));
8412 Em(
t,5) = 0.5*(J(
t,1)+J(
t,2));
8416 for (
int t = 0;
t < 3;
t++)
8418 Js(
t,0) = Em(
t,5)-Em(
t,0);
8419 Js(
t,1) = Em(
t,1)-Em(
t,0);
8420 Js(
t,2) = Em(
t,2)-Em(
t,0);
8424 for (
int t = 0;
t < 3;
t++)
8426 Js(
t,0) = Em(
t,5)-Em(
t,0);
8427 Js(
t,1) = Em(
t,2)-Em(
t,0);
8428 Js(
t,2) = Em(
t,4)-Em(
t,0);
8432 kappa_min = std::max(ar1, ar2);
8436 for (
int t = 0;
t < 3;
t++)
8438 Js(
t,0) = Em(
t,0)-Em(
t,1);
8439 Js(
t,1) = Em(
t,4)-Em(
t,1);
8440 Js(
t,2) = Em(
t,2)-Em(
t,1);
8444 for (
int t = 0;
t < 3;
t++)
8446 Js(
t,0) = Em(
t,2)-Em(
t,1);
8447 Js(
t,1) = Em(
t,4)-Em(
t,1);
8448 Js(
t,2) = Em(
t,5)-Em(
t,1);
8452 kappa = std::max(ar1, ar2);
8453 if (
kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
8456 for (
int t = 0;
t < 3;
t++)
8458 Js(
t,0) = Em(
t,0)-Em(
t,2);
8459 Js(
t,1) = Em(
t,1)-Em(
t,2);
8460 Js(
t,2) = Em(
t,3)-Em(
t,2);
8464 for (
int t = 0;
t < 3;
t++)
8466 Js(
t,0) = Em(
t,1)-Em(
t,2);
8467 Js(
t,1) = Em(
t,5)-Em(
t,2);
8468 Js(
t,2) = Em(
t,3)-Em(
t,2);
8472 kappa = std::max(ar1, ar2);
8473 if (
kappa < kappa_min) { rt = 2; }
8476 static const int mv_all[3][4][4] =
8478 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
8479 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
8480 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
8482 const int (&mv)[4][4] = mv_all[rt];
8484 #ifndef MFEM_USE_MEMALLOC 8486 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8488 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8490 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8492 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8494 for (
int k = 0; k < 4; k++)
8496 new_elements[j+4+k] =
8497 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8498 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8502 new_elements[j+0] = tet =
TetMemory.Alloc();
8503 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8505 new_elements[j+1] = tet =
TetMemory.Alloc();
8506 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8508 new_elements[j+2] = tet =
TetMemory.Alloc();
8509 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8511 new_elements[j+3] = tet =
TetMemory.Alloc();
8512 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8514 for (
int k = 0; k < 4; k++)
8516 new_elements[j+4+k] = tet =
TetMemory.Alloc();
8517 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8518 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8521 for (
int k = 0; k < 4; k++)
8526 for (
int k = 0; k < 4; k++)
8540 for (
int fi = 2; fi < 5; fi++)
8542 for (
int k = 0; k < 4; k++)
8549 for (
int ei = 0; ei < 9; ei++)
8551 for (
int k = 0; k < 2; k++)
8558 const int qf2 = f2qf[
f[2]];
8559 const int qf3 = f2qf[
f[3]];
8560 const int qf4 = f2qf[
f[4]];
8563 new Wedge(v[0], oedge+e[0], oedge+e[2],
8564 oedge+e[6], oface+qf2, oface+qf4, attr);
8567 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
8568 oface+qf3, oface+qf4, oface+qf2, attr);
8571 new Wedge(oedge+e[0], v[1], oedge+e[1],
8572 oface+qf2, oedge+e[7], oface+qf3, attr);
8575 new Wedge(oedge+e[2], oedge+e[1], v[2],
8576 oface+qf4, oface+qf3, oedge+e[8], attr);
8579 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
8580 v[3], oedge+e[3], oedge+e[5], attr);
8583 new Wedge(oface+qf3, oface+qf4, oface+qf2,
8584 oedge+e[4], oedge+e[5], oedge+e[3], attr);
8587 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
8588 oedge+e[3], v[4], oedge+e[4], attr);
8591 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
8592 oedge+e[5], oedge+e[4], v[5], attr);
8601 for (
int fi = 0; fi < 1; fi++)
8603 for (
int k = 0; k < 4; k++)
8610 for (
int ei = 0; ei < 8; ei++)
8612 for (
int k = 0; k < 2; k++)
8619 const int qf0 = f2qf[
f[0]];
8622 new Pyramid(v[0], oedge+e[0], oface+qf0,
8623 oedge+e[3], oedge+e[4], attr);
8626 new Pyramid(oedge+e[0], v[1], oedge+e[1],
8627 oface+qf0, oedge+e[5], attr);
8630 new Pyramid(oface+qf0, oedge+e[1], v[2],
8631 oedge+e[2], oedge+e[6], attr);
8634 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
8635 v[3], oedge+e[7], attr);
8638 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
8639 oedge+e[7], v[4], attr);
8642 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
8643 oedge+e[4], oface+qf0, attr);
8645 #ifndef MFEM_USE_MEMALLOC 8647 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
8651 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
8655 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
8659 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
8663 new_elements[j++] = tet =
TetMemory.Alloc();
8664 tet->
Init(oedge+e[0], oedge+e[4], oedge+e[5],
8667 new_elements[j++] = tet =
TetMemory.Alloc();
8668 tet->
Init(oedge+e[1], oedge+e[5], oedge+e[6],
8671 new_elements[j++] = tet =
TetMemory.Alloc();
8672 tet->
Init(oedge+e[2], oedge+e[6], oedge+e[7],
8675 new_elements[j++] = tet =
TetMemory.Alloc();
8676 tet->
Init(oedge+e[3], oedge+e[7], oedge+e[4],
8685 const int he = hex_counter;
8690 if (f2qf.
Size() == 0)
8696 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[
f[k]]; }
8702 for (
int fi = 0; fi < 6; fi++)
8704 for (
int k = 0; k < 4; k++)
8711 for (
int ei = 0; ei < 12; ei++)
8713 for (
int k = 0; k < 2; k++)
8721 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
8722 oedge+e[3], oedge+e[8], oface+qf[1],
8723 oelem+he, oface+qf[4], attr);
8726 oface+qf[0], oface+qf[1], oedge+e[9],
8727 oface+qf[2], oelem+he, attr);
8729 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
8730 oedge+e[2], oelem+he, oface+qf[2],
8731 oedge+e[10], oface+qf[3], attr);
8733 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
8734 v[3], oface+qf[4], oelem+he,
8735 oface+qf[3], oedge+e[11], attr);
8737 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
8738 oface+qf[4], v[4], oedge+e[4],
8739 oface+qf[5], oedge+e[7], attr);
8741 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
8742 oelem+he, oedge+e[4], v[5],
8743 oedge+e[5], oface+qf[5], attr);
8745 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
8746 oface+qf[3], oface+qf[5], oedge+e[5],
8747 v[6], oedge+e[6], attr);
8749 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
8750 oedge+e[11], oedge+e[7], oface+qf[5],
8751 oedge+e[6], v[7], attr);
8756 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
8768 const int attr =
boundary[i]->GetAttribute();
8769 int *v =
boundary[i]->GetVertices();
8776 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8783 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8785 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8787 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8789 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8797 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
8799 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
8801 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
8803 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
8807 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
8813 static const double A = 0.0, B = 0.5, C = 1.0, D = -1.0;
8814 static double tet_children[3*4*16] =
8816 A,A,A, B,A,A, A,B,A, A,A,B,
8817 B,A,A, C,A,A, B,B,A, B,A,B,
8818 A,B,A, B,B,A, A,C,A, A,B,B,
8819 A,A,B, B,A,B, A,B,B, A,A,C,
8824 B,A,A, A,B,B, A,B,A, A,A,B,
8825 B,A,A, A,B,B, A,A,B, B,A,B,
8826 B,A,A, A,B,B, B,A,B, B,B,A,
8827 B,A,A, A,B,B, B,B,A, A,B,A,
8829 A,B,A, B,A,A, B,A,B, A,A,B,
8830 A,B,A, A,A,B, B,A,B, A,B,B,
8831 A,B,A, A,B,B, B,A,B, B,B,A,
8832 A,B,A, B,B,A, B,A,B, B,A,A,
8834 A,A,B, B,A,A, A,B,A, B,B,A,
8835 A,A,B, A,B,A, A,B,B, B,B,A,
8836 A,A,B, A,B,B, B,A,B, B,B,A,
8837 A,A,B, B,A,B, B,A,A, B,B,A
8839 static double pyr_children[3*5*10] =
8841 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
8842 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
8843 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
8844 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
8845 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
8846 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
8847 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
8848 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
8849 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
8850 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
8852 static double pri_children[3*6*8] =
8854 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
8855 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
8856 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
8857 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
8858 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
8859 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
8860 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
8861 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
8863 static double hex_children[3*8*8] =
8865 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
8866 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
8867 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
8868 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
8869 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
8870 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
8871 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
8872 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
8884 for (
int i = 0; i <
elements.Size(); i++)
8915 int i, j, ind, nedges;
8922 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
8936 for (j = 0; j < marked_el.
Size(); j++)
8941 int new_v = cnv + j, new_e = cne + j;
8950 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
8952 UseExternalData(seg_children, 1, 2, 3);
8965 int *edge1 =
new int[nedges];
8966 int *edge2 =
new int[nedges];
8967 int *middle =
new int[nedges];
8969 for (i = 0; i < nedges; i++)
8971 edge1[i] = edge2[i] = middle[i] = -1;
8977 for (j = 1; j < v.
Size(); j++)
8979 ind = v_to_v(v[j-1], v[j]);
8980 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8982 ind = v_to_v(v[0], v[v.
Size()-1]);
8983 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8987 for (i = 0; i < marked_el.
Size(); i++)
8993 int need_refinement;
8996 need_refinement = 0;
8997 for (i = 0; i < nedges; i++)
8999 if (middle[i] != -1 && edge1[i] != -1)
9001 need_refinement = 1;
9006 while (need_refinement == 1);
9009 int v1[2], v2[2], bisect, temp;
9011 for (i = 0; i < temp; i++)
9014 bisect = v_to_v(v[0], v[1]);
9015 if (middle[bisect] != -1)
9019 v1[0] = v[0]; v1[1] = middle[bisect];
9020 v2[0] = middle[bisect]; v2[1] = v[1];
9026 mfem_error(
"Only bisection of segment is implemented" 9050 MFEM_VERIFY(
GetNE() == 0 ||
9052 "tetrahedral mesh is not marked for refinement:" 9053 " call Finalize(true)");
9060 for (i = 0; i < marked_el.
Size(); i++)
9066 for (i = 0; i < marked_el.
Size(); i++)
9075 for (i = 0; i < marked_el.
Size(); i++)
9092 int need_refinement;
9097 need_refinement = 0;
9105 if (
elements[i]->NeedRefinement(v_to_v))
9107 need_refinement = 1;
9112 while (need_refinement == 1);
9119 need_refinement = 0;
9121 if (
boundary[i]->NeedRefinement(v_to_v))
9123 need_refinement = 1;
9127 while (need_refinement == 1);
9158 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is " 9159 "not supported. Project the NURBS to Nodes first.");
9169 if (!refinements.
Size())
9190 Swap(*mesh2,
false);
9206 const int *fine,
int nfine,
int op)
9208 double error = elem_error[fine[0]];
9210 for (
int i = 1; i < nfine; i++)
9212 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
9214 double err_fine = elem_error[fine[i]];
9217 case 0: error = std::min(error, err_fine);
break;
9218 case 1: error += err_fine;
break;
9219 case 2: error = std::max(error, err_fine);
break;
9226 double threshold,
int nc_limit,
int op)
9228 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
9229 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. " 9230 "Project the NURBS to Nodes first.");
9243 for (
int i = 0; i < dt.
Size(); i++)
9245 if (nc_limit > 0 && !level_ok[i]) {
continue; }
9250 if (error < threshold) { derefs.
Append(i); }
9253 if (!derefs.
Size()) {
return false; }
9260 Swap(*mesh2,
false);
9274 int nc_limit,
int op)
9284 MFEM_ABORT(
"Derefinement is currently supported for non-conforming " 9291 int nc_limit,
int op)
9294 for (
int i = 0; i < tmp.Size(); i++)
9296 tmp[i] = elem_error(i);
9382 #ifdef MFEM_USE_MEMALLOC 9409 for (
int i = 0; i < elem_array.
Size(); i++)
9411 if (elem_array[i]->GetGeometryType() == geom)
9416 elem_vtx.
SetSize(nv*num_elems);
9420 for (
int i = 0; i < elem_array.
Size(); i++)
9426 elem_vtx.
Append(loc_vtx);
9434 for (
int i = 0; i < nelem; i++) { list[i] = i; }
9450 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
9462 default: MFEM_ABORT(
"internal error");
9468 int nonconforming,
int nc_limit)
9478 else if (nonconforming < 0)
9499 for (
int i = 0; i < refinements.
Size(); i++)
9501 el_to_refine[i] = refinements[i].index;
9505 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
9506 if (rt == 1 || rt == 2 || rt == 4)
9510 else if (rt == 3 || rt == 5 || rt == 6)
9528 for (
int i = 0; i < el_to_refine.
Size(); i++)
9530 refinements[i] =
Refinement(el_to_refine[i]);
9537 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. " 9538 "Please project the NURBS to Nodes first, with SetCurvature().");
9541 MFEM_VERIFY(
ncmesh != NULL || dynamic_cast<const ParMesh*>(
this) == NULL,
9542 "Sorry, converting a conforming ParMesh to an NC mesh is " 9550 (simplices_nonconforming && (
meshgen & 0x1)) )
9563 for (
int i = 0; i <
GetNE(); i++)
9565 if ((
double) rand() / RAND_MAX <
prob)
9570 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
9582 for (
int i = 0; i <
GetNE(); i++)
9585 bool refine =
false;
9586 for (
int j = 0; j < v.
Size(); j++)
9591 double d = vert(l) -
vertices[v[j]](l);
9594 if (dist <= eps*eps) { refine =
true;
break; }
9605 int nonconforming,
int nc_limit)
9607 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
9609 for (
int i = 0; i <
GetNE(); i++)
9611 if (elem_error[i] > threshold)
9625 int nonconforming,
int nc_limit)
9629 return RefineByError(tmp, threshold, nonconforming, nc_limit);
9634 int *edge1,
int *edge2,
int *middle)
9637 int v[2][4], v_new, bisect,
t;
9649 bisect = v_to_v(vert[0], vert[1]);
9650 MFEM_ASSERT(bisect >= 0,
"");
9652 if (middle[bisect] == -1)
9663 if (edge1[bisect] == i)
9665 edge1[bisect] = edge2[bisect];
9668 middle[bisect] = v_new;
9672 v_new = middle[bisect];
9680 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9681 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9702 bisect = v_to_v(v[1][0], v[1][1]);
9703 MFEM_ASSERT(bisect >= 0,
"");
9705 if (edge1[bisect] == i)
9709 else if (edge2[bisect] == i)
9718 MFEM_ABORT(
"Bisection for now works only for triangles.");
9725 int v[2][4], v_new, bisect,
t;
9732 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
9736 "TETRAHEDRON element is not marked for refinement.");
9741 bisect = v_to_v.
FindId(vert[0], vert[1]);
9745 for (j = 0; j < 3; j++)
9762 new_redges[0][0] = 2;
9763 new_redges[0][1] = 1;
9764 new_redges[1][0] = 2;
9765 new_redges[1][1] = 1;
9766 int tr1 = -1, tr2 = -1;
9767 switch (old_redges[0])
9770 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
9775 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
9779 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
9782 switch (old_redges[1])
9785 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
9790 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
9794 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
9801 #ifdef MFEM_USE_MEMALLOC 9837 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
9844 int v[2][3], v_new, bisect,
t;
9855 bisect = v_to_v.
FindId(vert[0], vert[1]);
9856 MFEM_ASSERT(bisect >= 0,
"");
9858 MFEM_ASSERT(v_new != -1,
"");
9862 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9863 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9873 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for" 9879 int *edge1,
int *edge2,
int *middle)
9882 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
9891 bisect[0] = v_to_v(v[0],v[1]);
9892 bisect[1] = v_to_v(v[1],v[2]);
9893 bisect[2] = v_to_v(v[0],v[2]);
9894 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
9896 for (j = 0; j < 3; j++)
9898 if (middle[bisect[j]] == -1)
9909 if (edge1[bisect[j]] == i)
9911 edge1[bisect[j]] = edge2[bisect[j]];
9914 middle[bisect[j]] = v_new[j];
9918 v_new[j] = middle[bisect[j]];
9921 edge1[bisect[j]] = -1;
9927 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
9928 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
9929 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
9930 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
9945 tri2->ResetTransform(code);
9946 tri3->ResetTransform(code);
9950 tri2->PushTransform(1);
9951 tri3->PushTransform(2);
9964 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
10000 for (
int i = 0; i < elem_geoms.
Size(); i++)
10008 std::map<unsigned, int> mat_no;
10012 for (
int j = 0; j <
elements.Size(); j++)
10015 unsigned code =
elements[j]->GetTransform();
10018 int &matrix = mat_no[code];
10019 if (!matrix) { matrix = mat_no.size(); }
10029 std::map<unsigned, int>::iterator it;
10030 for (it = mat_no.begin(); it != mat_no.end(); ++it)
10044 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for" 10045 " geom = " << geom);
10055 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
10064 os <<
"areamesh2\n\n";
10068 os <<
"curved_areamesh2\n\n";
10077 os <<
boundary[i]->GetAttribute();
10078 for (j = 0; j < v.
Size(); j++)
10080 os <<
' ' << v[j] + 1;
10092 for (j = 0; j < v.
Size(); j++)
10094 os <<
' ' << v[j] + 1;
10106 for (j = 1; j <
Dim; j++)
10123 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
10131 os <<
"NETGEN_Neutral_Format\n";
10136 for (j = 0; j <
Dim; j++)
10149 os <<
elements[i]->GetAttribute();
10150 for (j = 0; j < nv; j++)
10152 os <<
' ' << ind[j]+1;
10163 os <<
boundary[i]->GetAttribute();
10164 for (j = 0; j < nv; j++)
10166 os <<
' ' << ind[j]+1;
10178 <<
" 0 0 0 0 0 0 0\n" 10179 <<
"0 0 0 1 0 0 0 0 0 0 0\n" 10181 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n" 10182 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
10186 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
10192 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
10193 for (j = 0; j < nv; j++)
10195 os <<
' ' << ind[j]+1;
10204 os <<
boundary[i]->GetAttribute();
10205 for (j = 0; j < nv; j++)
10207 os <<
' ' << ind[j]+1;
10209 os <<
" 1.0 1.0 1.0 1.0\n";
10242 os <<
"\n# mesh curvature GridFunction";
10247 os <<
"\nmfem_mesh_end" << endl;
10252 os << (section_delimiter.empty()
10253 ?
"MFEM mesh v1.0\n" :
"MFEM mesh v1.2\n");
10257 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 10262 "# TETRAHEDRON = 4\n" 10268 os <<
"\ndimension\n" <<
Dim;
10303 if (!section_delimiter.empty())
10305 os << section_delimiter << endl;
10314 os <<
"MFEM NURBS mesh v1.0\n";
10318 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 10324 os <<
"\ndimension\n" <<
Dim 10341 int ki = e_to_k[i];
10346 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
10353 ofstream ofs(fname);
10354 ofs.precision(precision);
10358 #ifdef MFEM_USE_ADIOS2 10368 "# vtk DataFile Version 3.0\n" 10369 "Generated by MFEM\n" 10371 "DATASET UNSTRUCTURED_GRID\n";
10384 for ( ; j < 3; j++)
10400 os << (*Nodes)(vdofs[0]);
10404 os <<
' ' << (*Nodes)(vdofs[j]);
10406 for ( ; j < 3; j++)
10420 size +=
elements[i]->GetNVertices() + 1;
10425 const int *v =
elements[i]->GetVertices();
10426 const int nv =
elements[i]->GetNVertices();
10430 for (
int j = 0; j < nv; j++)
10432 os <<
' ' << v[perm ? perm[j] : j];
10445 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
10446 "Point meshes should have a single dof per element");
10447 size += dofs.
Size() + 1;
10452 if (!strcmp(fec_name,
"Linear") ||
10453 !strcmp(fec_name,
"H1_0D_P1") ||
10454 !strcmp(fec_name,
"H1_1D_P1") ||
10455 !strcmp(fec_name,
"H1_2D_P1") ||
10456 !strcmp(fec_name,
"H1_3D_P1"))
10460 else if (!strcmp(fec_name,
"Quadratic") ||
10461 !strcmp(fec_name,
"H1_1D_P2") ||
10462 !strcmp(fec_name,
"H1_2D_P2") ||
10463 !strcmp(fec_name,
"H1_3D_P2"))
10469 mfem::err <<
"Mesh::PrintVTK : can not save '" 10470 << fec_name <<
"' elements!" << endl;
10479 for (
int j = 0; j < dofs.
Size(); j++)
10481 os <<
' ' << dofs[j];
10484 else if (order == 2)
10486 const int *vtk_mfem;
10487 switch (
elements[i]->GetGeometryType())
10501 for (
int j = 0; j < dofs.
Size(); j++)
10503 os <<
' ' << dofs[vtk_mfem[j]];
10513 int vtk_cell_type = 5;
10517 os << vtk_cell_type <<
'\n';
10522 <<
"SCALARS material int\n" 10523 <<
"LOOKUP_TABLE default\n";
10526 os <<
elements[i]->GetAttribute() <<
'\n';
10533 bool high_order_output,
10534 int compression_level,
10537 int ref = (high_order_output &&
Nodes)
10540 fname = fname +
".vtu";
10542 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
10543 if (compression_level != 0)
10545 os <<
" compressor=\"vtkZLibDataCompressor\"";
10548 os <<
"<UnstructuredGrid>\n";
10549 PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
10550 os <<
"</Piece>\n";
10551 os <<
"</UnstructuredGrid>\n";
10552 os <<
"</VTKFile>" << std::endl;
10559 bool high_order_output,
10560 int compression_level)
10562 PrintVTU(fname, format, high_order_output, compression_level,
true);
10566 bool high_order_output,
int compression_level,
10574 std::vector<char> buf;
10576 auto get_geom = [&](
int i)
10584 int np = 0, nc_ref = 0;
10585 for (
int i = 0; i < ne; i++)
10594 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\"" 10595 << (high_order_output ? ne : nc_ref) <<
"\">\n";
10598 os <<
"<Points>\n";
10599 os <<
"<DataArray type=\"" << type_str
10600 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
10601 for (
int i = 0; i < ne; i++)
10614 for (
int j = 0; j < pmat.
Width(); j++)
10640 os <<
"</DataArray>" << std::endl;
10641 os <<
"</Points>" << std::endl;
10643 os <<
"<Cells>" << std::endl;
10644 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"" 10645 << fmt_str <<
"\">" << std::endl;
10647 std::vector<int> offset;
10650 if (high_order_output)
10653 for (
int iel = 0; iel < ne; iel++)
10657 int nnodes = local_connectivity.
Size();
10658 for (
int i=0; i<nnodes; ++i)
10665 offset.push_back(np);
10671 for (
int i = 0; i < ne; i++)
10677 for (
int j = 0; j < RG.
Size(); )
10680 offset.push_back(coff);
10682 for (
int k = 0; k < nv; k++, j++)
10696 os <<
"</DataArray>" << std::endl;
10698 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\"" 10699 << fmt_str <<
"\">" << std::endl;
10701 for (
size_t ii=0; ii<offset.size(); ii++)
10709 os <<
"</DataArray>" << std::endl;
10710 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\"" 10711 << fmt_str <<
"\">" << std::endl;
10713 const int *vtk_geom_map =
10715 for (
int i = 0; i < ne; i++)
10718 uint8_t vtk_cell_type = 5;
10720 vtk_cell_type = vtk_geom_map[geom];
10722 if (high_order_output)
10731 for (
int j = 0; j < RG.
Size(); j += nv)
10741 os <<
"</DataArray>" << std::endl;
10742 os <<
"</Cells>" << std::endl;
10744 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
10745 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\"" 10746 << fmt_str <<
"\">" << std::endl;
10747 for (
int i = 0; i < ne; i++)
10750 if (high_order_output)
10769 os <<
"</DataArray>" << std::endl;
10770 os <<
"</CellData>" << std::endl;
10781 "# vtk DataFile Version 3.0\n" 10782 "Generated by MFEM\n" 10784 "DATASET UNSTRUCTURED_GRID\n";
10789 os <<
"FIELD FieldData 1\n" 10799 np = nc = size = 0;
10800 for (
int i = 0; i <
GetNE(); i++)
10809 os <<
"POINTS " << np <<
" double\n";
10811 for (
int i = 0; i <
GetNE(); i++)
10818 for (
int j = 0; j < pmat.
Width(); j++)
10820 os << pmat(0, j) <<
' ';
10823 os << pmat(1, j) <<
' ';
10835 os << 0.0 <<
' ' << 0.0;
10842 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
10844 for (
int i = 0; i <
GetNE(); i++)
10851 for (
int j = 0; j < RG.
Size(); )
10854 for (
int k = 0; k < nv; k++, j++)
10856 os <<
' ' << np + RG[j];
10862 os <<
"CELL_TYPES " << nc <<
'\n';
10863 for (
int i = 0; i <
GetNE(); i++)
10871 for (
int j = 0; j < RG.
Size(); j += nv)
10873 os << vtk_cell_type <<
'\n';
10877 os <<
"CELL_DATA " << nc <<
'\n' 10878 <<
"SCALARS material int\n" 10879 <<
"LOOKUP_TABLE default\n";
10880 for (
int i = 0; i <
GetNE(); i++)
10888 os << attr <<
'\n';
10895 srand((
unsigned)time(0));
10896 double a = double(rand()) / (double(RAND_MAX) + 1.);
10897 int el0 = (int)floor(
a *
GetNE());
10899 os <<
"SCALARS element_coloring int\n" 10900 <<
"LOOKUP_TABLE default\n";
10901 for (
int i = 0; i <
GetNE(); i++)
10908 os << coloring[i] + 1 <<
'\n';
10914 os <<
"POINT_DATA " << np <<
'\n' << flush;
10919 int delete_el_to_el = (
el_to_el) ? (0) : (1);
10921 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
10924 const int *i_el_el = el_el.
GetI();
10925 const int *j_el_el = el_el.
GetJ();
10930 stack_p = stack_top_p = 0;
10931 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
10933 if (colors[el] != -2)
10939 el_stack[stack_top_p++] = el;
10941 for ( ; stack_p < stack_top_p; stack_p++)
10943 int i = el_stack[stack_p];
10944 int num_nb = i_el_el[i+1] - i_el_el[i];
10945 if (max_num_col < num_nb + 1)
10947 max_num_col = num_nb + 1;
10949 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10951 int k = j_el_el[j];
10952 if (colors[k] == -2)
10955 el_stack[stack_top_p++] = k;
10963 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
10965 int i = el_stack[stack_p], col;
10967 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10969 col = colors[j_el_el[j]];
10972 col_marker[col] = 1;
10976 for (col = 0; col < max_num_col; col++)
10977 if (col_marker[col] == 0)
10985 if (delete_el_to_el)
10993 int elem_attr)
const 10995 if (
Dim != 3 &&
Dim != 2) {
return; }
10997 int i, j, k, l, nv, nbe, *v;
10999 os <<
"MFEM mesh v1.0\n";
11003 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 11008 "# TETRAHEDRON = 4\n" 11013 os <<
"\ndimension\n" <<
Dim 11018 <<
' ' <<
elements[i]->GetGeometryType();
11021 for (j = 0; j < nv; j++)
11033 l = partitioning[l];
11048 os <<
"\nboundary\n" << nbe <<
'\n';
11054 l = partitioning[l];
11057 nv =
faces[i]->GetNVertices();
11058 v =
faces[i]->GetVertices();
11059 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
11060 for (j = 0; j < nv; j++)
11067 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
11068 for (j = nv-1; j >= 0; j--)
11079 nv =
faces[i]->GetNVertices();
11080 v =
faces[i]->GetVertices();
11081 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
11082 for (j = 0; j < nv; j++)
11113 int interior_faces)
11115 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
11116 if (
Dim != 3 &&
Dim != 2) {
return; }
11125 int nv =
elements[i]->GetNVertices();
11126 const int *ind =
elements[i]->GetVertices();
11127 for (
int j = 0; j < nv; j++)
11137 voff[i] = vcount[i-1] + voff[i-1];
11143 vown[i] =
new int[vcount[i]];
11155 int nv =
elements[i]->GetNVertices();
11156 const int *ind =
elements[i]->GetVertices();
11157 for (
int j = 0; j < nv; j++)
11160 vown[ind[j]][vcount[ind[j]]] = i;
11166 vcount[i] = voff[i+1] - voff[i];
11170 for (
int i = 0; i < edge_el.
Size(); i++)
11172 const int *el = edge_el.
GetRow(i);
11175 int k = partitioning[el[0]];
11176 int l = partitioning[el[1]];
11177 if (interior_faces || k != l)
11189 os <<
"areamesh2\n\n" << nbe <<
'\n';
11191 for (
int i = 0; i < edge_el.
Size(); i++)
11193 const int *el = edge_el.
GetRow(i);
11196 int k = partitioning[el[0]];
11197 int l = partitioning[el[1]];
11198 if (interior_faces || k != l)
11203 for (
int j = 0; j < 2; j++)
11204 for (
int s = 0;
s < vcount[ev[j]];
s++)
11205 if (vown[ev[j]][
s] == el[0])
11207 os <<
' ' << voff[ev[j]]+
s+1;
11211 for (
int j = 1; j >= 0; j--)
11212 for (
int s = 0;
s < vcount[ev[j]];
s++)
11213 if (vown[ev[j]][
s] == el[1])
11215 os <<
' ' << voff[ev[j]]+
s+1;
11222 int k = partitioning[el[0]];
11226 for (
int j = 0; j < 2; j++)
11227 for (
int s = 0;
s < vcount[ev[j]];
s++)
11228 if (vown[ev[j]][
s] == el[0])
11230 os <<
' ' << voff[ev[j]]+
s+1;
11240 int nv =
elements[i]->GetNVertices();
11241 const int *ind =
elements[i]->GetVertices();
11242 os << partitioning[i]+1 <<
' ';
11244 for (
int j = 0; j < nv; j++)
11246 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11247 vown[ind[j]][vcount[ind[j]]] = i;
11254 vcount[i] = voff[i+1] - voff[i];
11260 for (
int k = 0; k < vcount[i]; k++)
11262 for (
int j = 0; j <
Dim; j++)
11272 os <<
"NETGEN_Neutral_Format\n";
11276 for (
int k = 0; k < vcount[i]; k++)
11278 for (
int j = 0; j <
Dim; j++)
11289 int nv =
elements[i]->GetNVertices();
11290 const int *ind =
elements[i]->GetVertices();
11291 os << partitioning[i]+1;
11292 for (
int j = 0; j < nv; j++)
11294 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11295 vown[ind[j]][vcount[ind[j]]] = i;
11302 vcount[i] = voff[i+1] - voff[i];
11312 int k = partitioning[
faces_info[i].Elem1No];
11313 l = partitioning[l];
11314 if (interior_faces || k != l)
11331 int k = partitioning[
faces_info[i].Elem1No];
11332 l = partitioning[l];
11333 if (interior_faces || k != l)
11335 int nv =
faces[i]->GetNVertices();
11336 const int *ind =
faces[i]->GetVertices();
11338 for (
int j = 0; j < nv; j++)
11339 for (
int s = 0;
s < vcount[ind[j]];
s++)
11342 os <<
' ' << voff[ind[j]]+
s+1;
11346 for (
int j = nv-1; j >= 0; j--)
11347 for (
int s = 0;
s < vcount[ind[j]];
s++)
11350 os <<
' ' << voff[ind[j]]+
s+1;
11357 int k = partitioning[
faces_info[i].Elem1No];
11358 int nv =
faces[i]->GetNVertices();
11359 const int *ind =
faces[i]->GetVertices();
11361 for (
int j = 0; j < nv; j++)
11362 for (
int s = 0;
s < vcount[ind[j]];
s++)
11365 os <<
' ' << voff[ind[j]]+
s+1;
11381 int k = partitioning[
faces_info[i].Elem1No];
11382 l = partitioning[l];
11383 if (interior_faces || k != l)
11396 <<
" 0 0 0 0 0 0 0\n" 11397 <<
"0 0 0 1 0 0 0 0 0 0 0\n" 11398 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n" 11399 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n" 11400 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11403 for (
int k = 0; k < vcount[i]; k++)
11404 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' ' 11409 int nv =
elements[i]->GetNVertices();
11410 const int *ind =
elements[i]->GetVertices();
11411 os << i+1 <<
' ' << partitioning[i]+1;
11412 for (
int j = 0; j < nv; j++)
11414 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11415 vown[ind[j]][vcount[ind[j]]] = i;
11422 vcount[i] = voff[i+1] - voff[i];
11431 int k = partitioning[
faces_info[i].Elem1No];
11432 l = partitioning[l];
11433 if (interior_faces || k != l)
11435 int nv =
faces[i]->GetNVertices();
11436 const int *ind =
faces[i]->GetVertices();
11438 for (
int j = 0; j < nv; j++)
11439 for (
int s = 0;
s < vcount[ind[j]];
s++)
11442 os <<
' ' << voff[ind[j]]+
s+1;
11444 os <<
" 1.0 1.0 1.0 1.0\n";
11446 for (
int j = nv-1; j >= 0; j--)
11447 for (
int s = 0;
s < vcount[ind[j]];
s++)
11450 os <<
' ' << voff[ind[j]]+
s+1;
11452 os <<
" 1.0 1.0 1.0 1.0\n";
11457 int k = partitioning[
faces_info[i].Elem1No];
11458 int nv =
faces[i]->GetNVertices();
11459 const int *ind =
faces[i]->GetVertices();
11461 for (
int j = 0; j < nv; j++)
11462 for (
int s = 0;
s < vcount[ind[j]];
s++)
11465 os <<
' ' << voff[ind[j]]+
s+1;
11467 os <<
" 1.0 1.0 1.0 1.0\n";
11491 " NURBS mesh is not supported!");
11495 os <<
"MFEM mesh v1.0\n";
11499 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 11504 "# TETRAHEDRON = 4\n" 11509 os <<
"\ndimension\n" <<
Dim 11517 const int *
const i_AF_f = Aface_face.
GetI();
11518 const int *
const j_AF_f = Aface_face.
GetJ();
11520 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
11521 for (
const int * iface = j_AF_f + i_AF_f[iAF];
11522 iface < j_AF_f + i_AF_f[iAF+1];
11525 os << iAF+1 <<
' ';
11557 double *cg =
new double[na*
spaceDim];
11558 int *nbea =
new int[na];
11565 for (i = 0; i < na; i++)
11577 for (k = 0; k < vert.
Size(); k++)
11589 for (k = 0; k < vert.
Size(); k++)
11590 if (vn[vert[k]] == 1)
11595 cg[bea*
spaceDim+j] += pointmat(j,k);
11606 for (k = 0; k < vert.
Size(); k++)
11611 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
11627 double *cg =
new double[na*
spaceDim];
11628 int *nbea =
new int[na];
11635 for (i = 0; i < na; i++)
11647 for (k = 0; k < vert.
Size(); k++)
11659 for (k = 0; k < vert.
Size(); k++)
11660 if (vn[vert[k]] == 1)
11665 cg[bea*
spaceDim+j] += pointmat(j,k);
11676 for (k = 0; k < vert.
Size(); k++)
11681 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
11697 for (
int i = 0; i <
vertices.Size(); i++)
11699 for (
int j = 0; j <
spaceDim; j++)
11711 xnew.ProjectCoefficient(f_pert);
11720 "incompatible vector dimensions");
11728 for (
int d = 0; d <
spaceDim; d++)
11748 for (
int i = 0; i <
GetNE(); i++)
11753 for (
int j = 0; j < nv; j++)
11758 for (
int i = 0; i <
GetNBE(); i++)
11763 for (
int j = 0; j < nv; j++)
11769 for (
int i = 0; i < v2v.
Size(); i++)
11774 v2v[i] = num_vert++;
11778 if (num_vert == v2v.
Size()) {
return; }
11780 Vector nodes_by_element;
11785 for (
int i = 0; i <
GetNE(); i++)
11792 for (
int i = 0; i <
GetNE(); i++)
11801 for (
int i = 0; i <
GetNE(); i++)
11806 for (
int j = 0; j < nv; j++)
11811 for (
int i = 0; i <
GetNBE(); i++)
11816 for (
int j = 0; j < nv; j++)
11840 for (
int i = 0; i <
GetNE(); i++)
11853 int num_bdr_elem = 0;
11854 int new_bel_to_edge_nnz = 0;
11855 for (
int i = 0; i <
GetNBE(); i++)
11871 if (num_bdr_elem ==
GetNBE()) {
return; }
11875 Table *new_bel_to_edge = NULL;
11879 new_be_to_edge.
Reserve(num_bdr_elem);
11883 new_be_to_face.
Reserve(num_bdr_elem);
11884 new_bel_to_edge =
new Table;
11885 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
11887 for (
int i = 0; i <
GetNBE(); i++)
11898 int row = new_be_to_face.
Size();
11902 int *new_e = new_bel_to_edge->
GetRow(row);
11903 for (
int j = 0; j < ne; j++)
11907 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
11927 for (
int i = 0; i < attribs.
Size(); i++)
11939 #ifdef MFEM_USE_MEMALLOC 11966 const int npts = point_mat.
Width();
11967 if (!npts) {
return 0; }
11968 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
11972 if (!
GetNE()) {
return 0; }
11974 double *data = point_mat.
GetData();
11981 min_dist = std::numeric_limits<double>::max();
11985 for (
int i = 0; i <
GetNE(); i++)
11989 for (
int k = 0; k < npts; k++)
11992 if (dist < min_dist(k))
11994 min_dist(k) = dist;
12003 for (
int k = 0; k < npts; k++)
12007 int res = inv_tr->
Transform(pt, ips[k]);
12010 elem_ids[k] = e_idx[k];
12014 if (pts_found != npts)
12018 for (
int k = 0; k < npts; k++)
12020 if (elem_ids[k] != -1) {
continue; }
12024 for (
int v = 0; v < elvertices.
Size(); v++)
12026 int vv = elvertices[v];
12028 const int* els = vtoel->
GetRow(vv);
12029 for (
int e = 0; e < ne; e++)
12031 if (els[e] == e_idx[k]) {
continue; }
12033 int res = inv_tr->
Transform(pt, ips[k]);
12036 elem_ids[k] = els[e];
12048 for (
int e = 0; e < neigh.
Size(); e++)
12054 int res = inv_tr->
Transform(pt, ips[k]);
12067 if (inv_trans == NULL) {
delete inv_tr; }
12069 if (warn && pts_found != npts)
12071 MFEM_WARNING((npts-pts_found) <<
" points were not found");
12085 "mixed meshes are not supported!");
12086 MFEM_ASSERT(
mesh->
GetNodes(),
"meshes without nodes are not supported!");
12099 Compute(nodes, d_mt);
12102 void GeometricFactors::Compute(
const GridFunction &nodes,
12109 const int vdim = fespace->
GetVDim();
12110 const int NE = fespace->
GetNE();
12111 const int ND = fe->
GetDof();
12114 unsigned eval_flags = 0;
12139 qi->DisableTensorProducts(!use_tensor_products);
12147 Vector Enodes(vdim*ND*NE, my_d_mt);
12148 elem_restr->Mult(nodes, Enodes);
12149 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
12153 qi->Mult(nodes, eval_flags,
X,
J,
detJ);
12169 const int vdim = fespace->
GetVDim();
12183 face_restr->
Mult(*nodes, Fnodes);
12185 unsigned eval_flags = 0;
12194 J.
SetSize(vdim*vdim*NQ*NF, my_d_mt);
12228 T.Transform(ip, tip);
12232 V(1) = s * ((ip.
y + layer) / n);
12237 V(2) = s * ((ip.
z + layer) / n);
12246 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
12250 int nvy = (closed) ? (ny) : (ny + 1);
12251 int nvt = mesh->
GetNV() * nvy;
12260 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
12265 for (
int i = 0; i < mesh->
GetNV(); i++)
12268 for (
int j = 0; j < nvy; j++)
12270 vc[1] = sy * (double(j) / ny);
12276 for (
int i = 0; i < mesh->
GetNE(); i++)
12281 for (
int j = 0; j < ny; j++)
12284 qv[0] = vert[0] * nvy + j;
12285 qv[1] = vert[1] * nvy + j;
12286 qv[2] = vert[1] * nvy + (j + 1) % nvy;
12287 qv[3] = vert[0] * nvy + (j + 1) % nvy;
12293 for (
int i = 0; i < mesh->
GetNBE(); i++)
12298 for (
int j = 0; j < ny; j++)
12301 sv[0] = vert[0] * nvy + j;
12302 sv[1] = vert[0] * nvy + (j + 1) % nvy;
12306 Swap<int>(sv[0], sv[1]);
12318 for (
int i = 0; i < mesh->
GetNE(); i++)
12324 sv[0] = vert[0] * nvy;
12325 sv[1] = vert[1] * nvy;
12329 sv[0] = vert[1] * nvy + ny;
12330 sv[1] = vert[0] * nvy + ny;
12346 string cname = name;
12347 if (cname ==
"Linear")
12351 else if (cname ==
"Quadratic")
12355 else if (cname ==
"Cubic")
12359 else if (!strncmp(name,
"H1_", 3))
12363 else if (!strncmp(name,
"L2_T", 4))
12367 else if (!strncmp(name,
"L2_", 3))
12374 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : " 12386 for (
int i = 0; i < mesh->
GetNE(); i++)
12389 for (
int j = ny-1; j >= 0; j--)
12406 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
12411 int nvt = mesh->
GetNV() * nvz;
12416 bool wdgMesh =
false;
12417 bool hexMesh =
false;
12421 for (
int i = 0; i < mesh->
GetNV(); i++)
12425 for (
int j = 0; j < nvz; j++)
12427 vc[2] = sz * (double(j) / nz);
12433 for (
int i = 0; i < mesh->
GetNE(); i++)
12443 for (
int j = 0; j < nz; j++)
12446 pv[0] = vert[0] * nvz + j;
12447 pv[1] = vert[1] * nvz + j;
12448 pv[2] = vert[2] * nvz + j;
12449 pv[3] = vert[0] * nvz + (j + 1) % nvz;
12450 pv[4] = vert[1] * nvz + (j + 1) % nvz;
12451 pv[5] = vert[2] * nvz + (j + 1) % nvz;
12458 for (
int j = 0; j < nz; j++)
12461 hv[0] = vert[0] * nvz + j;
12462 hv[1] = vert[1] * nvz + j;
12463 hv[2] = vert[2] * nvz + j;
12464 hv[3] = vert[3] * nvz + j;
12465 hv[4] = vert[0] * nvz + (j + 1) % nvz;
12466 hv[5] = vert[1] * nvz + (j + 1) % nvz;
12467 hv[6] = vert[2] * nvz + (j + 1) % nvz;
12468 hv[7] = vert[3] * nvz + (j + 1) % nvz;
12470 mesh3d->
AddHex(hv, attr);
12474 mfem::err <<
"Extrude2D : Invalid 2D element type \'" 12475 << geom <<
"\'" << endl;
12481 for (
int i = 0; i < mesh->
GetNBE(); i++)
12486 for (
int j = 0; j < nz; j++)
12489 qv[0] = vert[0] * nvz + j;
12490 qv[1] = vert[1] * nvz + j;
12491 qv[2] = vert[1] * nvz + (j + 1) % nvz;
12492 qv[3] = vert[0] * nvz + (j + 1) % nvz;
12501 for (
int i = 0; i < mesh->
GetNE(); i++)
12512 tv[0] = vert[0] * nvz;
12513 tv[1] = vert[2] * nvz;
12514 tv[2] = vert[1] * nvz;
12518 tv[0] = vert[0] * nvz + nz;
12519 tv[1] = vert[1] * nvz + nz;
12520 tv[2] = vert[2] * nvz + nz;
12528 qv[0] = vert[0] * nvz;
12529 qv[1] = vert[3] * nvz;
12530 qv[2] = vert[2] * nvz;
12531 qv[3] = vert[1] * nvz;
12535 qv[0] = vert[0] * nvz + nz;
12536 qv[1] = vert[1] * nvz + nz;
12537 qv[2] = vert[2] * nvz + nz;
12538 qv[3] = vert[3] * nvz + nz;
12544 mfem::err <<
"Extrude2D : Invalid 2D element type \'" 12545 << geom <<
"\'" << endl;
12551 if ( hexMesh && wdgMesh )
12555 else if ( hexMesh )
12559 else if ( wdgMesh )
12572 string cname = name;
12573 if (cname ==
"Linear")
12577 else if (cname ==
"Quadratic")
12581 else if (cname ==
"Cubic")
12585 else if (!strncmp(name,
"H1_", 3))
12589 else if (!strncmp(name,
"L2_T", 4))
12593 else if (!strncmp(name,
"L2_", 3))
12600 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : " 12612 for (
int i = 0; i < mesh->
GetNE(); i++)
12615 for (
int j = nz-1; j >= 0; j--)
12636 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
12637 <<
" 0 0 " << i <<
" -1 0\n";
12644 double mid[3] = {0, 0, 0};
12645 for (
int j = 0; j < 2; j++)
12647 for (
int k = 0; k <
spaceDim; k++)
12653 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" " 12654 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
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.
const IntegrationRule * IntRule
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
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)
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
void SetCoordsFromPatches(Vector &Nodes)
void Print(std::ostream &out) const
I/O: Print the mesh in "MFEM NC mesh v1.0" format.
BiLinear2DFiniteElement QuadrilateralFE
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 GetNPoints() const
Returns the number of the points in the integration rule.
Arc::Index insert_arc(Node::Index i, Node::Index j, Float w=1, Float b=1)
Class for an integration rule - an Array of IntegrationPoint.
const FiniteElementSpace * GetNodalFESpace() const
void Init(int ind1, int ind2, int ind3, int ind4, int attr=1, int ref_flag=0)
Initialize the vertex indices and the attribute of a Tetrahedron.
void 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.
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T...
T * end()
STL-like end. Returns pointer after the last element of the array.
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 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.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
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)
Linear1DFiniteElement SegmentFE
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 GetVertices(Array< int > &v) const =0
Returns element's vertices.
class LinearPyramidFiniteElement PyramidFE
Geometry::Type GetElementBaseGeometry(int i) const
int GetBdrElementEdgeIndex(int i) const
const Table & ElementToFaceTable() const
Array< Element * > boundary
int * GeneratePartitioning(int nparts, int part_method=1)
CoarseFineTransformations CoarseFineTr
virtual void LimitNCLevel(int max_nc_level)
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void MoveVertices(const Vector &displacements)
void SetSize(int s)
Resize the vector to size s.
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
static FiniteElement * GetTransformationFEforElementType(Element::Type)
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.
int Dimension() const
Return the dimension of the NCMesh.
Lists all edges/faces in the nonconforming mesh.
virtual void UniformRefinement2D()
Refine a mixed 2D mesh uniformly.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones...
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
void DisableTensorProducts(bool disable=true) const
Disable the use of tensor product evaluations, for tensor-product elements, e.g. quads and hexes...
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
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]
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'.
void GetVertexToVertexTable(DSTable &) const
int GetAttribute() const
Return element's attribute.
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]
T * GetData()
Returns the data.
Vector detJ
Determinants of the Jacobians at all quadrature points.
bool Nonconforming() const
static const int FaceVert[NumFaces][MaxFaceVert]
int Size() const
Returns the size of the vector.
static const int Edges[NumEdges][2]
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.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Piecewise-(bi/tri)linear continuous finite elements.
Data type dense matrix using column-major storage.
const NURBSExtension * GetNURBSext() const
double * Data() const
Returns the matrix data array.
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Mesh * Extrude1D(Mesh *mesh, const int ny, const double sy, const bool closed)
Extrude a 1D mesh.
void Transform(void(*f)(const Vector &, Vector &))
bool IsSlaveFace(const FaceInfo &fi) const
int GetNDofs() const
Returns number of degrees of freedom.
void AverageVertices(const int *indexes, int n, int result)
Averages the vertices with given indexes and saves the result in vertices[result].
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 GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
void GetBdrElementTopo(Array< Element *> &boundary) const
void GetElementData(const Array< Element *> &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
int GetNEdges() const
Return the number of edges.
const Element * GetElement(int i) const
TriLinear3DFiniteElement HexahedronFE
Evaluate the derivatives at quadrature points.
bool IsGhost(const Element &el) const
Return true if the Element el is a ghost element.
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().
void KnotInsert(Array< KnotVector *> &kv)
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
Data type quadrilateral element.
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.
void Print(std::ostream &out) const
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 RemoveInternalBoundaries()
void GetVertexDofs(int i, Array< int > &dofs) const
Mesh * Extrude2D(Mesh *mesh, const int nz, const double sz)
Extrude a 2D mesh.
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
void Print(const Mesh &mesh, const adios2stream::mode print_mode=mode::sync)
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...
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...
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
void GetPointMatrix(int i, DenseMatrix &pointmat) const
int GetNBE() const
Returns number of boundary elements.
friend class NURBSExtension
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)
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.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
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 &)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
void OnMeshUpdated(Mesh *mesh)
Element * NewElement(int geom)
void SetVerticesFromNodes(const GridFunction *nodes)
Helper to set vertex coordinates given a high-order curvature function.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Structure for storing face geometric factors: coordinates, Jacobians, determinants of the Jacobians...
ElementTransformation * GetBdrElementTransformation(int i)
void GetElementTopo(Array< Element *> &elements) const
virtual void UpdateMeshPointer(Mesh *new_mesh)
void DebugDump(std::ostream &out) const
Output an NCMesh-compatible debug dump.
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)
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void MakeRefined_(Mesh &orig_mesh, const Array< int > ref_factors, int ref_type)
Internal function used in Mesh::MakeRefined.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
int GetAttribute(int i) const
Return the attribute of element i.
Piecewise-(bi)cubic continuous finite elements.
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 GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i...
int AddVertex(double x, double y=0.0, double z=0.0)
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
Vector J
Jacobians of the element transformations at all quadrature points.
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
Data type Pyramid element.
static const int Edges[NumEdges][2]
void MoveNodes(const Vector &displacements)
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
GeometryRefiner GlobGeometryRefiner
double * GetData() const
Returns the matrix data array.
double f(const Vector &xvec)
int GetMaxElementOrder() const
Return the maximum polynomial order.
Geometry::Type GetGeometryType() const
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.
void GetVertices(Vector &vert_coord) const
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)
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)
void GetNode(int i, double *coord) const
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
const FiniteElementCollection * FEColl() const
Vector detJ
Determinants of the Jacobians at all quadrature points.
int FindCoarseElement(int i)
int SpaceDimension() const
Return the space dimension of the NCMesh.
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
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)
virtual void PrintXG(std::ostream &os=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
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 KnotInsert(Array< KnotVector *> &kv)
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)
Set the curvature of the mesh nodes using the given polynomial degree.
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Type
Constants for the classes derived from Element.
void Get(double *p, const int dim) 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.
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.
PointFiniteElement PointFE
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
VTKFormat
Data array format for VTK and VTU files.
virtual const char * Name() const
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
int AddBdrSegment(int v1, int v2, int attr=1)
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...
int AddElement(Element *elem)
The parameter elem should be allocated using the NewElement() method.
void ReadLineMesh(std::istream &input)
virtual void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
const CoarseFineTransformations & GetRefinementTransforms()
void SetLayer(const int l)
Table * GetFaceToElementTable() const
Array< int > FindFaceNeighbors(const int elem) const
Returns the sorted, unique indices of elements sharing a face with element elem, including elem...
int GetElementToEdgeTable(Table &, Array< int > &)
T * begin()
STL-like begin. Returns pointer to the first element of the array.
void SetKnotsFromPatches()
void RandomRefinement(double prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
Data type triangle element.
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
signed char local
local number within 'element'
void GetColumn(int c, Vector &col) const
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.
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()
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
int GetNE() const
Returns number of elements in the mesh.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
static const int QuadraticMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to legacy quadratic VTK geometries/.
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
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 ReadMFEMMesh(std::istream &input, int version, int &curved)
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
This function is not public anymore. Use GeneralRefinement instead.
void Start()
Start the stopwatch. The elapsed time is not cleared.
Table * GetFaceEdgeTable() const
Returns the face-to-edge Table (3D)
void FindPartitioningComponents(Table &elem_elem, const Array< int > &partitioning, Array< int > &component, Array< int > &num_comp)
Data type tetrahedron element.
double * GetData() const
Return a pointer to the beginning of the Vector data.
List of mesh geometries stored as Array<Geometry::Type>.
int GetVDim() const
Returns vector dimension.
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)
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
int CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
Geometry::Type GetElementGeometry(int i) const
GeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, MemoryType d_mt=MemoryType::DEFAULT)
void PrintElementsWithPartitioning(int *partitioning, std::ostream &out, int interior_faces=0)
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
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.
Mesh * GetMesh() const
Returns the mesh.
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.
const Table & ElementToEdgeTable() const
double p(const Vector &x, double t)
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
int GetDim() const
Returns the reference space dimension for the finite element.
Nonconforming edge/face within a bigger edge/face.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
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.
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
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.
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)
void PrintTopo(std::ostream &out, const Array< int > &e_to_k) 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.
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
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.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE
A class to initialize the size of a Tensor.
static Mesh MakeCartesian1D(int n, double sx=1.0)
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
Array< Element * > elements
const Table & ElementToElementTable()
void SetRelaxedHpConformity(bool relaxed=true)
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
MFEM_EXPORT Linear2DFiniteElement TriangleFE
int SpaceDimension() const
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)
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
Evaluate the values at quadrature points.
Geometry::Type GetBdrElementGeometry(int i) const
int GetNE() const
Returns number of elements.
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
NURBSExtension * NURBSext
Optional NURBS mesh extension.
void PrintBdrVTU(std::string fname, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0)
static const int Edges[NumEdges][2]
Class for integration point with weight.
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
Element::Type GetElementType(int i) const
Returns the type of element i.
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Element * ReadElementWithoutAttr(std::istream &)
void Swap(Mesh &other, bool non_geometry)
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
static const int vtk_quadratic_wedge[18]
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
static const int DimStart[MaxDim+2]
virtual int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type, does not count master nonconforming face...
int Size() const
Returns the number of TYPE I elements.
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
Ordering::Type GetOrdering() const
Return the ordering method.
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.
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn't exist.
static const int Orient[NumOrient][NumVert]
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
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)
int index(int i, int j, int nx, int ny)
int NumberOfEntries() const
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
static const int Edges[NumEdges][2]
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]
double Norml2() const
Returns the l2 norm of the vector.
STable3D * GetFacesTable()
Piecewise-(bi)quadratic continuous finite elements.
int Size_of_connections() const
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
NCMesh * ncmesh
Optional nonconforming mesh extension.
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)
int Size() const
Return the logical size of the array.
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
FaceInformation GetFaceInformation(int f) const
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.
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
double AggregateError(const Array< double > &elem_error, const int *fine, int nfine, int op)
Derefinement helper.
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...
int DofToVDof(int dof, int vd, int ndofs=-1) const
Mesh & operator=(Mesh &&mesh)
Move assignment operstor.
Evaluate the derivatives at quadrature points.
int GetNFaces() const
Return the number of faces in a 3D mesh.
void ReadTrueGridMesh(std::istream &input)
virtual void Print(std::ostream &os=mfem::out) const
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
static void GetElementArrayEdgeTable(const Array< Element *> &elem_array, const DSTable &v_to_v, Table &el_to_edge)
virtual void Save(const char *fname, int precision=16) 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 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.
void PrintSurfaces(const Table &Aface_face, std::ostream &out) const
Print set of disjoint surfaces:
Arbitrary order H1-conforming (continuous) finite elements.
void XYZ_VectorFunction(const Vector &p, Vector &v)
void GetElementColoring(Array< int > &colors, int el0=0)
void GetNodes(Vector &node_coord) const
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 ...
int GetNumFacesWithGhost() const
Return the number of faces (3D), edges (2D) or vertices (1D) including ghost faces.
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
void Printer(std::ostream &out=mfem::out, std::string section_delimiter="") const
const std::string filename
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
void GetElementInteriorDofs(int i, Array< int > &dofs) const
MemAlloc< Tetrahedron, 1024 > TetMemory
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
Node::Index insert_node(Float length=1)
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
const Element * GetBdrElement(int i) const
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Operation GetLastOperation() const
Return type of last modification of the mesh.
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.
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)
std::ostream & operator<<(std::ostream &os, const Mesh &mesh)
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
void DegreeElevate(int rel_degree, int degree=16)
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.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
uint rank(Node::Index i) const
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Rank 3 tensor (array of matrices)
Class for parallel meshes.
Geometry::Type GetBdrElementBaseGeometry(int i) const
Abstract data type element.
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...
Data type line segment element.
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 DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
int NumberOfElements()
Return the number of elements added to the table.
void GenerateNCFaceInfo()
virtual Type GetType() const =0
Returns element's type.
double GetLength(int i, int j) const
Return the length of the segment from node i to node j.
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.
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
void ConvertToPatches(const Vector &Nodes)
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 ...
static const int Orient[NumOrient][NumVert]
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 the topological part of a Mesh::Operation, such as refinement, has been performed.
Defines the position of a fine element within a coarse element.
void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const
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...
virtual void Refine(const Array< Refinement > &refinements)
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)
static const int FaceVert[NumFaces][MaxFaceVert]
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally, e.g. by modifying the internal nodal GridFunction returned by GetNodes().
Class used to extrude the nodes of a mesh.