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*);
75 void Mesh::GetElementCenter(
int i,
Vector ¢er)
78 int geom = GetElementBaseGeometry(i);
93 return pow(fabs(J.
Weight()), 1./Dim);
105 double Mesh::GetElementSize(
int i,
int type)
107 return GetElementSize(GetElementTransformation(i), type);
110 double Mesh::GetElementSize(
int i,
const Vector &dir)
114 GetElementJacobian(i, J);
116 return sqrt((d_hat * d_hat) / (dir * dir));
119 double Mesh::GetElementVolume(
int i)
141 for (
int d = 0; d < spaceDim; d++)
150 for (
int i = 0; i < NumOfVertices; i++)
152 coord = GetVertex(i);
153 for (
int d = 0; d < spaceDim; d++)
155 if (coord[d] < min(d)) { min(d) = coord[d]; }
156 if (coord[d] > max(d)) { max(d) = coord[d]; }
162 const bool use_boundary =
false;
163 int ne = use_boundary ? GetNBE() : GetNE();
171 for (
int i = 0; i < ne; i++)
175 GetBdrElementFace(i, &fn, &fo);
177 Tr = GetFaceElementTransformations(fn, 5);
184 T = GetElementTransformation(i);
186 T->Transform(RefG->
RefPts, pointmat);
188 for (
int j = 0; j < pointmat.
Width(); j++)
190 for (
int d = 0; d < pointmat.
Height(); d++)
192 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
193 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
200 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
201 double &kappa_min,
double &kappa_max,
209 sdim = SpaceDimension();
211 if (Vh) { Vh->
SetSize(NumOfElements); }
212 if (Vk) { Vk->
SetSize(NumOfElements); }
215 h_max = kappa_max = -h_min;
216 if (
dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
218 for (i = 0; i < NumOfElements; i++)
220 GetElementJacobian(i, J);
221 h = pow(fabs(J.
Weight()), 1.0/
double(
dim));
224 if (Vh) { (*Vh)(i) = h; }
225 if (Vk) { (*Vk)(i) =
kappa; }
227 if (h < h_min) { h_min = h; }
228 if (h > h_max) { h_max = h; }
235 void Mesh::PrintElementsByGeometry(
int dim,
239 for (
int g = Geometry::DimStart[
dim], first = 1;
240 g < Geometry::DimStart[
dim+1]; g++)
242 if (!num_elems_by_geom[g]) {
continue; }
243 if (!first) { os <<
" + "; }
245 os << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
249 void Mesh::PrintCharacteristics(
Vector *Vh,
Vector *Vk, std::ostream &os)
251 double h_min, h_max, kappa_min, kappa_max;
253 os <<
"Mesh Characteristics:";
255 this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
257 Array<int> num_elems_by_geom(Geometry::NumGeom);
258 num_elems_by_geom = 0;
259 for (
int i = 0; i < GetNE(); i++)
261 num_elems_by_geom[GetElementBaseGeometry(i)]++;
265 <<
"Dimension : " << Dimension() <<
'\n' 266 <<
"Space dimension : " << SpaceDimension();
270 <<
"Number of vertices : " << GetNV() <<
'\n' 271 <<
"Number of elements : " << GetNE() <<
'\n' 272 <<
"Number of bdr elem : " << GetNBE() <<
'\n';
277 <<
"Number of vertices : " << GetNV() <<
'\n' 278 <<
"Number of elements : " << GetNE() <<
'\n' 279 <<
"Number of bdr elem : " << GetNBE() <<
'\n' 280 <<
"h_min : " << h_min <<
'\n' 281 <<
"h_max : " << h_max <<
'\n';
286 <<
"Number of vertices : " << GetNV() <<
'\n' 287 <<
"Number of edges : " << GetNEdges() <<
'\n' 288 <<
"Number of elements : " << GetNE() <<
" -- ";
289 PrintElementsByGeometry(2, num_elems_by_geom, os);
291 <<
"Number of bdr elem : " << GetNBE() <<
'\n' 292 <<
"Euler Number : " << EulerNumber2D() <<
'\n' 293 <<
"h_min : " << h_min <<
'\n' 294 <<
"h_max : " << h_max <<
'\n' 295 <<
"kappa_min : " << kappa_min <<
'\n' 296 <<
"kappa_max : " << kappa_max <<
'\n';
300 Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
301 num_bdr_elems_by_geom = 0;
302 for (
int i = 0; i < GetNBE(); i++)
304 num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
306 Array<int> num_faces_by_geom(Geometry::NumGeom);
307 num_faces_by_geom = 0;
308 for (
int i = 0; i < GetNFaces(); i++)
310 num_faces_by_geom[GetFaceGeometry(i)]++;
314 <<
"Number of vertices : " << GetNV() <<
'\n' 315 <<
"Number of edges : " << GetNEdges() <<
'\n' 316 <<
"Number of faces : " << GetNFaces() <<
" -- ";
317 PrintElementsByGeometry(Dim-1, num_faces_by_geom, os);
319 <<
"Number of elements : " << GetNE() <<
" -- ";
320 PrintElementsByGeometry(Dim, num_elems_by_geom, os);
322 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
323 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, os);
325 <<
"Euler Number : " << EulerNumber() <<
'\n' 326 <<
"h_min : " << h_min <<
'\n' 327 <<
"h_max : " << h_max <<
'\n' 328 <<
"kappa_min : " << kappa_min <<
'\n' 329 <<
"kappa_max : " << kappa_max <<
'\n';
331 os <<
'\n' << std::flush;
338 case Element::POINT :
return &
PointFE;
339 case Element::SEGMENT :
return &
SegmentFE;
344 case Element::WEDGE :
return &
WedgeFE;
345 case Element::PYRAMID :
return &
PyramidFE;
347 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
350 MFEM_ABORT(
"Unknown element type");
359 ElTr->
ElementType = ElementTransformation::ELEMENT;
365 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
371 Nodes->FESpace()->GetElementVDofs(i, vdofs);
374 int n = vdofs.
Size()/spaceDim;
376 for (
int k = 0; k < spaceDim; k++)
378 for (
int j = 0; j < n; j++)
380 pm(k,j) = nodes(vdofs[n*k+j]);
383 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
387 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
392 ElTr->
ElementType = ElementTransformation::ELEMENT;
399 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
400 int nv = elements[i]->GetNVertices();
401 const int *v = elements[i]->GetVertices();
402 int n = vertices.Size();
404 for (
int k = 0; k < spaceDim; k++)
406 for (
int j = 0; j < nv; j++)
408 pm(k, j) = nodes(k*n+v[j]);
411 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
415 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
417 Nodes->FESpace()->GetElementVDofs(i, vdofs);
418 int n = vdofs.
Size()/spaceDim;
420 for (
int k = 0; k < spaceDim; k++)
422 for (
int j = 0; j < n; j++)
424 pm(k,j) = nodes(vdofs[n*k+j]);
427 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
433 GetElementTransformation(i, &Transformation);
435 return &Transformation;
440 GetBdrElementTransformation(i, &BdrTransformation);
441 return &BdrTransformation;
448 ElTr->
ElementType = ElementTransformation::BDR_ELEMENT;
454 GetBdrPointMatrix(i, pm);
455 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
465 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
466 int n = vdofs.
Size()/spaceDim;
468 for (
int k = 0; k < spaceDim; k++)
470 for (
int j = 0; j < n; j++)
472 pm(k,j) = nodes(vdofs[n*k+j]);
479 int elem_id, face_info;
480 GetBdrElementAdjacentElement2(i, elem_id, face_info);
482 GetLocalFaceTransformation(GetBdrElementType(i),
483 GetElementType(elem_id),
484 FaceElemTr.Loc1.Transf, face_info);
489 Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
490 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
491 "Mesh requires nodal Finite Element.");
493 FaceElemTr.Loc1.Transf.ElementNo = elem_id;
494 FaceElemTr.Loc1.Transf.mesh =
this;
495 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
496 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
497 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
499 ElTr->
SetFE(face_el);
506 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
514 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
515 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
517 for (
int i = 0; i < spaceDim; i++)
519 for (
int j = 0; j < nv; j++)
521 pm(i, j) = vertices[v[j]](i);
524 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
528 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
534 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
535 int n = vdofs.
Size()/spaceDim;
537 for (
int i = 0; i < spaceDim; i++)
539 for (
int j = 0; j < n; j++)
541 pm(i, j) = nodes(vdofs[n*i+j]);
548 FaceInfo &face_info = faces_info[FaceNo];
553 GetLocalFaceTransformation(face_type,
554 GetElementType(face_info.
Elem1No),
555 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
558 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
560 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
561 "Mesh requires nodal Finite Element.");
564 FaceElemTr.Loc1.Transf.ElementNo = face_info.
Elem1No;
565 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
566 FaceElemTr.Loc1.Transf.mesh =
this;
567 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
568 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
577 GetFaceTransformation(FaceNo, &FaceTransformation);
578 return &FaceTransformation;
585 GetFaceTransformation(EdgeNo, EdTr);
590 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
602 GetEdgeVertices(EdgeNo, v);
605 for (
int i = 0; i < spaceDim; i++)
607 for (
int j = 0; j < nv; j++)
609 pm(i, j) = vertices[v[j]](i);
612 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
616 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
622 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
623 int n = vdofs.
Size()/spaceDim;
625 for (
int i = 0; i < spaceDim; i++)
627 for (
int j = 0; j < n; j++)
629 pm(i, j) = nodes(vdofs[n*i+j]);
632 EdTr->
SetFE(edge_el);
636 MFEM_ABORT(
"Not implemented.");
643 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
644 return &EdgeTransformation;
648 void Mesh::GetLocalPtToSegTransformation(
663 void Mesh::GetLocalSegToTriTransformation(
672 tv = tri_t::Edges[i/64];
673 so = seg_t::Orient[i%64];
676 for (
int j = 0; j < 2; j++)
678 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
679 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
683 void Mesh::GetLocalSegToQuadTransformation(
692 qv = quad_t::Edges[i/64];
693 so = seg_t::Orient[i%64];
696 for (
int j = 0; j < 2; j++)
698 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
699 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
703 void Mesh::GetLocalTriToTetTransformation(
711 const int *tv = tet_t::FaceVert[i/64];
714 const int *to = tri_t::Orient[i%64];
718 for (
int j = 0; j < 3; j++)
721 locpm(0, j) = vert.
x;
722 locpm(1, j) = vert.
y;
723 locpm(2, j) = vert.
z;
727 void Mesh::GetLocalTriToWdgTransformation(
735 MFEM_VERIFY(i < 128,
"Local face index " << i/64
736 <<
" is not a triangular face of a wedge.");
737 const int *pv = pri_t::FaceVert[i/64];
740 const int *to = tri_t::Orient[i%64];
744 for (
int j = 0; j < 3; j++)
747 locpm(0, j) = vert.
x;
748 locpm(1, j) = vert.
y;
749 locpm(2, j) = vert.
z;
753 void Mesh::GetLocalTriToPyrTransformation(
760 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
761 <<
" is not a triangular face of a pyramid.");
762 const int *pv = pyr_t::FaceVert[i/64];
765 const int *to = tri_t::Orient[i%64];
769 for (
int j = 0; j < 3; j++)
772 locpm(0, j) = vert.
x;
773 locpm(1, j) = vert.
y;
774 locpm(2, j) = vert.
z;
778 void Mesh::GetLocalQuadToHexTransformation(
786 const int *hv = hex_t::FaceVert[i/64];
788 const int *qo = quad_t::Orient[i%64];
791 for (
int j = 0; j < 4; j++)
794 locpm(0, j) = vert.
x;
795 locpm(1, j) = vert.
y;
796 locpm(2, j) = vert.
z;
800 void Mesh::GetLocalQuadToWdgTransformation(
808 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
809 <<
" is not a quadrilateral face of a wedge.");
810 const int *pv = pri_t::FaceVert[i/64];
812 const int *qo = quad_t::Orient[i%64];
815 for (
int j = 0; j < 4; j++)
818 locpm(0, j) = vert.
x;
819 locpm(1, j) = vert.
y;
820 locpm(2, j) = vert.
z;
824 void Mesh::GetLocalQuadToPyrTransformation(
831 MFEM_VERIFY(i < 64,
"Local face index " << i/64
832 <<
" is not a quadrilateral face of a pyramid.");
833 const int *pv = pyr_t::FaceVert[i/64];
835 const int *qo = quad_t::Orient[i%64];
838 for (
int j = 0; j < 4; j++)
841 locpm(0, j) = vert.
x;
842 locpm(1, j) = vert.
y;
843 locpm(2, j) = vert.
z;
851 for (
int i = 0; i < geom_factors.Size(); i++)
863 geom_factors.Append(gf);
871 for (
int i = 0; i < face_geom_factors.Size(); i++)
885 face_geom_factors.Append(gf);
889 void Mesh::DeleteGeometricFactors()
891 for (
int i = 0; i < geom_factors.Size(); i++)
893 delete geom_factors[i];
895 geom_factors.SetSize(0);
896 for (
int i = 0; i < face_geom_factors.Size(); i++)
898 delete face_geom_factors[i];
900 face_geom_factors.SetSize(0);
903 void Mesh::GetLocalFaceTransformation(
909 GetLocalPtToSegTransformation(Transf, info);
912 case Element::SEGMENT:
913 if (elem_type == Element::TRIANGLE)
915 GetLocalSegToTriTransformation(Transf, info);
919 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
920 GetLocalSegToQuadTransformation(Transf, info);
924 case Element::TRIANGLE:
925 if (elem_type == Element::TETRAHEDRON)
927 GetLocalTriToTetTransformation(Transf, info);
929 else if (elem_type == Element::WEDGE)
931 GetLocalTriToWdgTransformation(Transf, info);
933 else if (elem_type == Element::PYRAMID)
935 GetLocalTriToPyrTransformation(Transf, info);
939 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for " 940 "face type " << face_type
941 <<
" and element type " << elem_type <<
"\n");
945 case Element::QUADRILATERAL:
946 if (elem_type == Element::HEXAHEDRON)
948 GetLocalQuadToHexTransformation(Transf, info);
950 else if (elem_type == Element::WEDGE)
952 GetLocalQuadToWdgTransformation(Transf, info);
954 else if (elem_type == Element::PYRAMID)
956 GetLocalQuadToPyrTransformation(Transf, info);
960 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for " 961 "face type " << face_type
962 <<
" and element type " << elem_type <<
"\n");
971 FaceInfo &face_info = faces_info[FaceNo];
974 FaceElemTr.SetConfigurationMask(cmask);
975 FaceElemTr.Elem1 = NULL;
976 FaceElemTr.Elem2 = NULL;
980 if (mask & FaceElementTransformations::HAVE_ELEM1)
982 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
983 FaceElemTr.Elem1 = &Transformation;
990 FaceElemTr.Elem2No = face_info.
Elem2No;
991 if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
992 FaceElemTr.Elem2No >= 0)
995 if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
996 { MFEM_ABORT(
"NURBS mesh not supported!"); }
998 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
999 FaceElemTr.Elem2 = &Transformation2;
1004 if (mask & FaceElementTransformations::HAVE_FACE)
1006 GetFaceTransformation(FaceNo, &FaceElemTr);
1011 FaceElemTr.SetGeometryType(GetFaceGeometry(FaceNo));
1015 int face_type = GetFaceElementType(FaceNo);
1016 if (mask & FaceElementTransformations::HAVE_LOC1)
1018 int elem_type = GetElementType(face_info.
Elem1No);
1019 GetLocalFaceTransformation(face_type, elem_type,
1020 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
1023 if ((mask & FaceElementTransformations::HAVE_LOC2) &&
1024 FaceElemTr.Elem2No >= 0)
1026 int elem_type = GetElementType(face_info.
Elem2No);
1027 GetLocalFaceTransformation(face_type, elem_type,
1028 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
1031 if (Nonconforming() && IsSlaveFace(face_info))
1033 ApplyLocalSlaveTransformation(FaceElemTr, face_info,
false);
1038 FaceElemTr.SetConfigurationMask(cmask);
1044 double dist = FaceElemTr.CheckConsistency();
1047 mfem::out <<
"\nInternal error: face id = " << FaceNo
1048 <<
", dist = " << dist <<
'\n';
1049 FaceElemTr.CheckConsistency(1);
1050 MFEM_ABORT(
"internal error");
1060 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
1066 #ifdef MFEM_THREAD_SAFE 1071 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1072 MFEM_ASSERT(nc_faces_info[fi.
NCFace].Slave,
"internal error");
1083 std::swap(composition(0,0), composition(0,1));
1084 std::swap(composition(1,0), composition(1,1));
1105 int fn = GetBdrFace(BdrElemNo);
1108 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
1112 tr = GetFaceElementTransformations(fn, 21);
1113 tr->
Attribute = boundary[BdrElemNo]->GetAttribute();
1115 tr->
ElementType = ElementTransformation::BDR_FACE;
1120 int Mesh::GetBdrFace(
int BdrElemNo)
const 1125 fn = be_to_face[BdrElemNo];
1129 fn = be_to_edge[BdrElemNo];
1133 fn = boundary[BdrElemNo]->GetVertices()[0];
1144 GetFaceElements(
f, &e1, &e2);
1145 GetFaceInfos(
f, &inf1, &inf2, &ncface);
1155 if (
f < GetNumFaces())
1161 face.
tag = FaceInfoTag::LocalConforming;
1162 face.
topology = FaceTopology::Conforming;
1171 face.
tag = FaceInfoTag::LocalSlaveNonconforming;
1172 face.
topology = FaceTopology::Nonconforming;
1177 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1188 face.
tag = FaceInfoTag::Boundary;
1189 face.
topology = FaceTopology::Boundary;
1198 face.
tag = FaceInfoTag::SharedConforming;
1199 face.
topology = FaceTopology::Conforming;
1211 face.
tag = FaceInfoTag::MasterNonconforming;
1212 face.
topology = FaceTopology::Nonconforming;
1221 face.
tag = FaceInfoTag::SharedSlaveNonconforming;
1222 face.
topology = FaceTopology::Nonconforming;
1237 face.
tag = FaceInfoTag::GhostMaster;
1247 face.
tag = FaceInfoTag::GhostSlave;
1248 face.
topology = FaceTopology::Nonconforming;
1265 case FaceInfoTag::LocalConforming:
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::LocalSlaveNonconforming:
1273 res.Elem1No = element[0].index;
1274 res.Elem2No = element[1].index;
1275 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1276 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1277 res.NCFace = ncface;
1279 case FaceInfoTag::Boundary:
1280 res.Elem1No = element[0].index;
1281 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1283 case FaceInfoTag::SharedConforming:
1284 res.Elem1No = element[0].index;
1285 res.Elem2No = -1 - element[1].index;
1286 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1287 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1289 case FaceInfoTag::MasterNonconforming:
1290 res.Elem1No = element[0].index;
1291 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1293 case FaceInfoTag::SharedSlaveNonconforming:
1294 res.Elem1No = element[0].index;
1295 res.Elem2No = -1 - element[1].index;
1296 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1297 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1299 case FaceInfoTag::GhostMaster:
1301 case FaceInfoTag::GhostSlave:
1302 res.Elem1No = element[0].index;
1303 res.Elem2No = -1 - element[1].index;
1304 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1305 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1313 os <<
"face topology=";
1316 case Mesh::FaceTopology::Boundary:
1319 case Mesh::FaceTopology::Conforming:
1322 case Mesh::FaceTopology::Nonconforming:
1323 os <<
"Non-conforming";
1325 case Mesh::FaceTopology::NA:
1330 os <<
"element[0].location=";
1333 case Mesh::ElementLocation::Local:
1336 case Mesh::ElementLocation::FaceNbr:
1339 case Mesh::ElementLocation::NA:
1344 os <<
"element[1].location=";
1347 case Mesh::ElementLocation::Local:
1350 case Mesh::ElementLocation::FaceNbr:
1353 case Mesh::ElementLocation::NA:
1358 os <<
"element[0].conformity=";
1361 case Mesh::ElementConformity::Coincident:
1364 case Mesh::ElementConformity::Superset:
1367 case Mesh::ElementConformity::Subset:
1370 case Mesh::ElementConformity::NA:
1375 os <<
"element[1].conformity=";
1378 case Mesh::ElementConformity::Coincident:
1381 case Mesh::ElementConformity::Superset:
1384 case Mesh::ElementConformity::Subset:
1387 case Mesh::ElementConformity::NA:
1392 os <<
"element[0].index=" << info.
element[0].
index <<
'\n' 1393 <<
"element[1].index=" << info.
element[1].
index <<
'\n' 1398 <<
"ncface=" << info.
ncface << std::endl;
1402 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const 1404 *Elem1 = faces_info[Face].Elem1No;
1405 *Elem2 = faces_info[Face].Elem2No;
1408 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const 1410 *Inf1 = faces_info[Face].Elem1Inf;
1411 *Inf2 = faces_info[Face].Elem2Inf;
1414 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2,
int *NCFace)
const 1416 *Inf1 = faces_info[Face].Elem1Inf;
1417 *Inf2 = faces_info[Face].Elem2Inf;
1418 *NCFace = faces_info[Face].NCFace;
1425 case 1:
return Geometry::POINT;
1426 case 2:
return Geometry::SEGMENT;
1428 if (Face < NumOfFaces)
1430 return faces[Face]->GetGeometryType();
1433 const int nc_face_id = faces_info[Face].NCFace;
1434 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1435 return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
1437 return Geometry::INVALID;
1442 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1447 Array<int> face_to_be(Dim == 2 ? NumOfEdges : NumOfFaces);
1449 for (
int i = 0; i < NumOfBdrElements; i++)
1451 face_to_be[GetBdrElementEdgeIndex(i)] = i;
1461 NumOfElements = NumOfBdrElements = 0;
1462 NumOfEdges = NumOfFaces = 0;
1463 nbInteriorFaces = -1;
1464 nbBoundaryFaces = -1;
1465 meshgen = mesh_geoms = 0;
1471 last_operation = Mesh::NONE;
1474 void Mesh::InitTables()
1477 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1478 face_to_elem = NULL;
1481 void Mesh::SetEmpty()
1487 void Mesh::DestroyTables()
1492 DeleteGeometricFactors();
1502 delete face_to_elem;
1503 face_to_elem = NULL;
1506 void Mesh::DestroyPointers()
1508 if (own_nodes) {
delete Nodes; }
1514 for (
int i = 0; i < NumOfElements; i++)
1516 FreeElement(elements[i]);
1519 for (
int i = 0; i < NumOfBdrElements; i++)
1521 FreeElement(boundary[i]);
1524 for (
int i = 0; i < faces.Size(); i++)
1526 FreeElement(faces[i]);
1532 void Mesh::Destroy()
1536 elements.DeleteAll();
1537 vertices.DeleteAll();
1538 boundary.DeleteAll();
1540 faces_info.DeleteAll();
1541 nc_faces_info.DeleteAll();
1542 be_to_edge.DeleteAll();
1543 be_to_face.DeleteAll();
1551 CoarseFineTr.Clear();
1553 #ifdef MFEM_USE_MEMALLOC 1557 attributes.DeleteAll();
1558 bdr_attributes.DeleteAll();
1561 void Mesh::ResetLazyData()
1563 delete el_to_el; el_to_el = NULL;
1564 delete face_edge; face_edge = NULL;
1565 delete face_to_elem; face_to_elem = NULL;
1566 delete edge_vertex; edge_vertex = NULL;
1567 DeleteGeometricFactors();
1568 nbInteriorFaces = -1;
1569 nbBoundaryFaces = -1;
1572 void Mesh::SetAttributes()
1577 for (
int i = 0; i < attribs.
Size(); i++)
1579 attribs[i] = GetBdrAttribute(i);
1583 attribs.
Copy(bdr_attributes);
1584 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1586 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1590 for (
int i = 0; i < attribs.
Size(); i++)
1592 attribs[i] = GetAttribute(i);
1596 attribs.
Copy(attributes);
1597 if (attributes.Size() > 0 && attributes[0] <= 0)
1599 MFEM_WARNING(
"Non-positive attributes in the domain!");
1603 void Mesh::InitMesh(
int Dim_,
int spaceDim_,
int NVert,
int NElem,
int NBdrElem)
1608 spaceDim = spaceDim_;
1611 vertices.SetSize(NVert);
1614 elements.SetSize(NElem);
1616 NumOfBdrElements = 0;
1617 boundary.SetSize(NBdrElem);
1620 template<
typename T>
1621 static void CheckEnlarge(
Array<T> &array,
int size)
1623 if (size >= array.
Size()) { array.
SetSize(size + 1); }
1626 int Mesh::AddVertex(
double x,
double y,
double z)
1628 CheckEnlarge(vertices, NumOfVertices);
1629 double *v = vertices[NumOfVertices]();
1633 return NumOfVertices++;
1636 int Mesh::AddVertex(
const double *coords)
1638 CheckEnlarge(vertices, NumOfVertices);
1639 vertices[NumOfVertices].SetCoords(spaceDim, coords);
1640 return NumOfVertices++;
1645 MFEM_ASSERT(coords.
Size() >= spaceDim,
1646 "invalid 'coords' size: " << coords.
Size());
1647 return AddVertex(coords.
GetData());
1650 void Mesh::AddVertexParents(
int i,
int p1,
int p2)
1656 if (i < vertices.Size())
1658 double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1659 for (
int j = 0; j < 3; j++)
1661 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1666 int Mesh::AddSegment(
int v1,
int v2,
int attr)
1668 CheckEnlarge(elements, NumOfElements);
1669 elements[NumOfElements] =
new Segment(v1, v2, attr);
1670 return NumOfElements++;
1673 int Mesh::AddSegment(
const int *vi,
int attr)
1675 CheckEnlarge(elements, NumOfElements);
1676 elements[NumOfElements] =
new Segment(vi, attr);
1677 return NumOfElements++;
1680 int Mesh::AddTriangle(
int v1,
int v2,
int v3,
int attr)
1682 CheckEnlarge(elements, NumOfElements);
1683 elements[NumOfElements] =
new Triangle(v1, v2, v3, attr);
1684 return NumOfElements++;
1687 int Mesh::AddTriangle(
const int *vi,
int attr)
1689 CheckEnlarge(elements, NumOfElements);
1690 elements[NumOfElements] =
new Triangle(vi, attr);
1691 return NumOfElements++;
1694 int Mesh::AddQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1696 CheckEnlarge(elements, NumOfElements);
1697 elements[NumOfElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1698 return NumOfElements++;
1701 int Mesh::AddQuad(
const int *vi,
int attr)
1703 CheckEnlarge(elements, NumOfElements);
1705 return NumOfElements++;
1708 int Mesh::AddTet(
int v1,
int v2,
int v3,
int v4,
int attr)
1710 int vi[4] = {v1, v2, v3, v4};
1711 return AddTet(vi, attr);
1714 int Mesh::AddTet(
const int *vi,
int attr)
1716 CheckEnlarge(elements, NumOfElements);
1717 #ifdef MFEM_USE_MEMALLOC 1719 tet = TetMemory.Alloc();
1722 elements[NumOfElements] = tet;
1724 elements[NumOfElements] =
new Tetrahedron(vi, attr);
1726 return NumOfElements++;
1729 int Mesh::AddWedge(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int attr)
1731 CheckEnlarge(elements, NumOfElements);
1732 elements[NumOfElements] =
new Wedge(v1, v2, v3, v4, v5, v6, attr);
1733 return NumOfElements++;
1736 int Mesh::AddWedge(
const int *vi,
int attr)
1738 CheckEnlarge(elements, NumOfElements);
1739 elements[NumOfElements] =
new Wedge(vi, attr);
1740 return NumOfElements++;
1743 int Mesh::AddPyramid(
int v1,
int v2,
int v3,
int v4,
int v5,
int attr)
1745 CheckEnlarge(elements, NumOfElements);
1746 elements[NumOfElements] =
new Pyramid(v1, v2, v3, v4, v5, attr);
1747 return NumOfElements++;
1750 int Mesh::AddPyramid(
const int *vi,
int attr)
1752 CheckEnlarge(elements, NumOfElements);
1753 elements[NumOfElements] =
new Pyramid(vi, attr);
1754 return NumOfElements++;
1757 int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1760 CheckEnlarge(elements, NumOfElements);
1761 elements[NumOfElements] =
1762 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1763 return NumOfElements++;
1766 int Mesh::AddHex(
const int *vi,
int attr)
1768 CheckEnlarge(elements, NumOfElements);
1769 elements[NumOfElements] =
new Hexahedron(vi, attr);
1770 return NumOfElements++;
1773 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1775 static const int hex_to_tet[6][4] =
1777 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1778 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1782 for (
int i = 0; i < 6; i++)
1784 for (
int j = 0; j < 4; j++)
1786 ti[j] = vi[hex_to_tet[i][j]];
1792 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1794 static const int hex_to_wdg[2][6] =
1796 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1800 for (
int i = 0; i < 2; i++)
1802 for (
int j = 0; j < 6; j++)
1804 ti[j] = vi[hex_to_wdg[i][j]];
1810 void Mesh::AddHexAsPyramids(
const int *vi,
int attr)
1812 static const int hex_to_pyr[6][5] =
1814 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
1815 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
1819 for (
int i = 0; i < 6; i++)
1821 for (
int j = 0; j < 5; j++)
1823 ti[j] = vi[hex_to_pyr[i][j]];
1825 AddPyramid(ti, attr);
1831 CheckEnlarge(elements, NumOfElements);
1832 elements[NumOfElements] = elem;
1833 return NumOfElements++;
1838 CheckEnlarge(boundary, NumOfBdrElements);
1839 boundary[NumOfBdrElements] = elem;
1840 return NumOfBdrElements++;
1843 int Mesh::AddBdrSegment(
int v1,
int v2,
int attr)
1845 CheckEnlarge(boundary, NumOfBdrElements);
1846 boundary[NumOfBdrElements] =
new Segment(v1, v2, attr);
1847 return NumOfBdrElements++;
1850 int Mesh::AddBdrSegment(
const int *vi,
int attr)
1852 CheckEnlarge(boundary, NumOfBdrElements);
1853 boundary[NumOfBdrElements] =
new Segment(vi, attr);
1854 return NumOfBdrElements++;
1857 int Mesh::AddBdrTriangle(
int v1,
int v2,
int v3,
int attr)
1859 CheckEnlarge(boundary, NumOfBdrElements);
1860 boundary[NumOfBdrElements] =
new Triangle(v1, v2, v3, attr);
1861 return NumOfBdrElements++;
1864 int Mesh::AddBdrTriangle(
const int *vi,
int attr)
1866 CheckEnlarge(boundary, NumOfBdrElements);
1867 boundary[NumOfBdrElements] =
new Triangle(vi, attr);
1868 return NumOfBdrElements++;
1871 int Mesh::AddBdrQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1873 CheckEnlarge(boundary, NumOfBdrElements);
1874 boundary[NumOfBdrElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1875 return NumOfBdrElements++;
1878 int Mesh::AddBdrQuad(
const int *vi,
int attr)
1880 CheckEnlarge(boundary, NumOfBdrElements);
1882 return NumOfBdrElements++;
1885 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1887 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1890 for (
int i = 0; i < 2; i++)
1892 for (
int j = 0; j < 3; j++)
1894 ti[j] = vi[quad_to_tri[i][j]];
1896 AddBdrTriangle(ti, attr);
1900 int Mesh::AddBdrPoint(
int v,
int attr)
1902 CheckEnlarge(boundary, NumOfBdrElements);
1903 boundary[NumOfBdrElements] =
new Point(&v, attr);
1904 return NumOfBdrElements++;
1907 void Mesh::GenerateBoundaryElements()
1910 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1914 for (i = 0; i < boundary.Size(); i++)
1916 FreeElement(boundary[i]);
1926 NumOfBdrElements = 0;
1927 for (i = 0; i < faces_info.Size(); i++)
1929 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1932 boundary.SetSize(NumOfBdrElements);
1933 be2face.
SetSize(NumOfBdrElements);
1934 for (j = i = 0; i < faces_info.Size(); i++)
1936 if (faces_info[i].Elem2No < 0)
1938 boundary[j] = faces[i]->Duplicate(
this);
1945 void Mesh::FinalizeCheck()
1947 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1948 vertices.Size() == 0,
1949 "incorrect number of vertices: preallocated: " << vertices.Size()
1950 <<
", actually added: " << NumOfVertices);
1951 MFEM_VERIFY(elements.Size() == NumOfElements,
1952 "incorrect number of elements: preallocated: " << elements.Size()
1953 <<
", actually added: " << NumOfElements);
1954 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1955 "incorrect number of boundary elements: preallocated: " 1956 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1959 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1962 CheckElementOrientation(fix_orientation);
1966 MarkTriMeshForRefinement();
1971 el_to_edge =
new Table;
1972 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1974 CheckBdrElementOrientation();
1988 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1989 bool fix_orientation)
1992 if (fix_orientation)
1994 CheckElementOrientation(fix_orientation);
1999 el_to_edge =
new Table;
2000 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2002 CheckBdrElementOrientation();
2022 GeckoProgress(
double limit) : limit(limit) { sw.
Start(); }
2023 virtual bool quit()
const {
return limit > 0 && sw.
UserTime() > limit; }
2026 class GeckoVerboseProgress :
public GeckoProgress
2032 GeckoVerboseProgress(
double limit) : GeckoProgress(limit) {}
2034 virtual void beginorder(
const Graph* graph,
Float cost)
const 2035 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2036 virtual void endorder(
const Graph* graph,
Float cost)
const 2037 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2039 virtual void beginiter(
const Graph* graph,
2042 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window " 2043 << window << std::flush;
2045 virtual void enditer(
const Graph* graph,
Float mincost,
Float cost)
const 2046 {
mfem::out <<
", cost = " << cost << endl; }
2051 int iterations,
int window,
2052 int period,
int seed,
bool verbose,
2058 GeckoProgress progress(time_limit);
2059 GeckoVerboseProgress vprogress(time_limit);
2062 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2069 const Table &my_el_to_el = ElementToElementTable();
2070 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2072 const int *neighid = my_el_to_el.
GetRow(elemid);
2073 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2075 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2080 graph.
order(&functional, iterations, window, period, seed,
2081 verbose ? &vprogress : &progress);
2087 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2090 return graph.
cost();
2101 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
2102 : coord(coord), dir(dir), points(points), mid(mid) {}
2104 bool operator()(
int i)
const 2106 return (points[3*i + coord] < mid) != dir;
2110 static void HilbertSort2D(
int coord1,
2113 const Array<double> &points,
int *beg,
int *end,
2114 double xmin,
double ymin,
double xmax,
double ymax)
2116 if (end - beg <= 1) {
return; }
2118 double xmid = (xmin + xmax)*0.5;
2119 double ymid = (ymin + ymax)*0.5;
2121 int coord2 = (coord1 + 1) % 2;
2124 int *p0 = beg, *p4 = end;
2125 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2126 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2127 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2131 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2132 ymin, xmin, ymid, xmid);
2134 if (p1 != p0 || p2 != p4)
2136 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2137 xmin, ymid, xmid, ymax);
2139 if (p2 != p0 || p3 != p4)
2141 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2142 xmid, ymid, xmax, ymax);
2146 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2147 ymid, xmax, ymin, xmid);
2151 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2152 const Array<double> &points,
int *beg,
int *end,
2153 double xmin,
double ymin,
double zmin,
2154 double xmax,
double ymax,
double zmax)
2156 if (end - beg <= 1) {
return; }
2158 double xmid = (xmin + xmax)*0.5;
2159 double ymid = (ymin + ymax)*0.5;
2160 double zmid = (zmin + zmax)*0.5;
2162 int coord2 = (coord1 + 1) % 3;
2163 int coord3 = (coord1 + 2) % 3;
2166 int *p0 = beg, *p8 = end;
2167 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2168 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2169 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2170 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2171 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2172 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2173 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2177 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2178 zmin, xmin, ymin, zmid, xmid, ymid);
2180 if (p1 != p0 || p2 != p8)
2182 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2183 ymin, zmid, xmin, ymid, zmax, xmid);
2185 if (p2 != p0 || p3 != p8)
2187 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2188 ymid, zmid, xmin, ymax, zmax, xmid);
2190 if (p3 != p0 || p4 != p8)
2192 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2193 xmin, ymax, zmid, xmid, ymid, zmin);
2195 if (p4 != p0 || p5 != p8)
2197 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2198 xmid, ymax, zmid, xmax, ymid, zmin);
2200 if (p5 != p0 || p6 != p8)
2202 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2203 ymax, zmid, xmax, ymid, zmax, xmid);
2205 if (p6 != p0 || p7 != p8)
2207 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2208 ymid, zmid, xmax, ymin, zmax, xmid);
2212 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2213 zmid, xmax, ymin, zmin, xmid, ymid);
2219 MFEM_VERIFY(spaceDim <= 3,
"");
2222 GetBoundingBox(min, max);
2227 if (spaceDim < 3) { points = 0.0; }
2230 for (
int i = 0; i < GetNE(); i++)
2232 GetElementCenter(i, center);
2233 for (
int j = 0; j < spaceDim; j++)
2235 points[3*i + j] = center(j);
2242 indices.
Sort([&](
int a,
int b)
2243 {
return points[3*
a] < points[3*
b]; });
2245 else if (spaceDim == 2)
2248 HilbertSort2D(0,
false,
false,
2249 points, indices.
begin(), indices.
end(),
2250 min(0), min(1), max(0), max(1));
2255 HilbertSort3D(0,
false,
false,
false,
2256 points, indices.
begin(), indices.
end(),
2257 min(0), min(1), min(2), max(0), max(1), max(2));
2262 for (
int i = 0; i < GetNE(); i++)
2264 ordering[indices[i]] = i;
2269 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
2273 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2278 MFEM_WARNING(
"element reordering of non-conforming meshes is not" 2282 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
2312 old_elem_node_vals.SetSize(GetNE());
2313 nodes_fes = Nodes->FESpace();
2316 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2318 nodes_fes->GetElementVDofs(old_elid, old_dofs);
2320 old_elem_node_vals[old_elid] =
new Vector(vals);
2326 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2328 int new_elid = ordering[old_elid];
2329 new_elements[new_elid] = elements[old_elid];
2334 if (reorder_vertices)
2339 vertex_ordering = -1;
2341 int new_vertex_ind = 0;
2342 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2344 int *elem_vert = elements[new_elid]->GetVertices();
2345 int nv = elements[new_elid]->GetNVertices();
2346 for (
int vi = 0; vi < nv; ++vi)
2348 int old_vertex_ind = elem_vert[vi];
2349 if (vertex_ordering[old_vertex_ind] == -1)
2351 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2352 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
2362 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2364 int *elem_vert = elements[new_elid]->GetVertices();
2365 int nv = elements[new_elid]->GetNVertices();
2366 for (
int vi = 0; vi < nv; ++vi)
2368 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2373 for (
int belid = 0; belid < GetNBE(); ++belid)
2375 int *be_vert = boundary[belid]->GetVertices();
2376 int nv = boundary[belid]->GetNVertices();
2377 for (
int vi = 0; vi < nv; ++vi)
2379 be_vert[vi] = vertex_ordering[be_vert[vi]];
2390 el_to_edge =
new Table;
2391 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2396 GetElementToFaceTable();
2406 last_operation = Mesh::NONE;
2407 nodes_fes->Update(
false);
2410 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2412 int new_elid = ordering[old_elid];
2413 nodes_fes->GetElementVDofs(new_elid, new_dofs);
2414 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
2415 delete old_elem_node_vals[old_elid];
2421 void Mesh::MarkForRefinement()
2427 MarkTriMeshForRefinement();
2431 DSTable v_to_v(NumOfVertices);
2432 GetVertexToVertexTable(v_to_v);
2433 MarkTetMeshForRefinement(v_to_v);
2438 void Mesh::MarkTriMeshForRefinement()
2443 for (
int i = 0; i < NumOfElements; i++)
2445 if (elements[i]->GetType() == Element::TRIANGLE)
2447 GetPointMatrix(i, pmat);
2448 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
2459 for (
int i = 0; i < NumOfVertices; i++)
2464 length_idx[j].one = GetLength(i, it.Column());
2465 length_idx[j].two = j;
2472 for (
int i = 0; i < NumOfEdges; i++)
2474 order[length_idx[i].two] = i;
2478 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
2483 GetEdgeOrdering(v_to_v, order);
2485 for (
int i = 0; i < NumOfElements; i++)
2487 if (elements[i]->GetType() == Element::TETRAHEDRON)
2489 elements[i]->MarkEdge(v_to_v, order);
2492 for (
int i = 0; i < NumOfBdrElements; i++)
2494 if (boundary[i]->GetType() == Element::TRIANGLE)
2496 boundary[i]->MarkEdge(v_to_v, order);
2503 if (*old_v_to_v && *old_elem_vert)
2510 if (*old_v_to_v == NULL)
2512 bool need_v_to_v =
false;
2514 for (
int i = 0; i < GetNEdges(); i++)
2520 if (dofs.
Size() > 0)
2528 *old_v_to_v =
new DSTable(NumOfVertices);
2529 GetVertexToVertexTable(*(*old_v_to_v));
2532 if (*old_elem_vert == NULL)
2534 bool need_elem_vert =
false;
2536 for (
int i = 0; i < GetNE(); i++)
2542 if (dofs.
Size() > 1)
2544 need_elem_vert =
true;
2550 *old_elem_vert =
new Table;
2551 (*old_elem_vert)->
MakeI(GetNE());
2552 for (
int i = 0; i < GetNE(); i++)
2554 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2556 (*old_elem_vert)->MakeJ();
2557 for (
int i = 0; i < GetNE(); i++)
2559 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2560 elements[i]->GetNVertices());
2562 (*old_elem_vert)->ShiftUpI();
2575 const int num_edge_dofs = old_dofs.
Size();
2578 const Vector onodes = *Nodes;
2582 int offset = NumOfVertices * old_dofs.
Size();
2586 if (num_edge_dofs > 0)
2588 DSTable new_v_to_v(NumOfVertices);
2589 GetVertexToVertexTable(new_v_to_v);
2591 for (
int i = 0; i < NumOfVertices; i++)
2595 const int old_i = (*old_v_to_v)(i, it.Column());
2596 const int new_i = it.Index();
2597 if (new_i == old_i) {
continue; }
2599 old_dofs.
SetSize(num_edge_dofs);
2600 new_dofs.
SetSize(num_edge_dofs);
2601 for (
int j = 0; j < num_edge_dofs; j++)
2603 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2604 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2608 for (
int j = 0; j < old_dofs.
Size(); j++)
2610 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2614 offset += NumOfEdges * num_edge_dofs;
2622 Table old_face_vertex;
2623 old_face_vertex.
MakeI(NumOfFaces);
2624 for (
int i = 0; i < NumOfFaces; i++)
2628 old_face_vertex.
MakeJ();
2629 for (
int i = 0; i < NumOfFaces; i++)
2631 faces[i]->GetNVertices());
2635 STable3D *faces_tbl = GetElementToFaceTable(1);
2641 for (
int i = 0; i < NumOfFaces; i++)
2643 const int *old_v = old_face_vertex.
GetRow(i);
2645 switch (old_face_vertex.
RowSize(i))
2648 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2652 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2656 new_fdofs[new_i+1] = old_dofs.
Size();
2661 for (
int i = 0; i < NumOfFaces; i++)
2663 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2666 switch (old_face_vertex.
RowSize(i))
2669 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2670 new_v = faces[new_i]->GetVertices();
2671 new_or = GetTriOrientation(old_v, new_v);
2676 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2677 new_v = faces[new_i]->GetVertices();
2678 new_or = GetQuadOrientation(old_v, new_v);
2685 for (
int j = 0; j < old_dofs.
Size(); j++)
2688 const int old_j = dof_ord[j];
2689 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2693 for (
int j = 0; j < old_dofs.
Size(); j++)
2695 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2717 for (
int i = 0; i < GetNE(); i++)
2719 const int *old_v = old_elem_vert->
GetRow(i);
2720 const int *new_v = elements[i]->GetVertices();
2726 case Geometry::SEGMENT:
2727 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2729 case Geometry::TRIANGLE:
2730 new_or = GetTriOrientation(old_v, new_v);
2732 case Geometry::SQUARE:
2733 new_or = GetQuadOrientation(old_v, new_v);
2735 case Geometry::TETRAHEDRON:
2736 new_or = GetTetOrientation(old_v, new_v);
2740 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2741 <<
" FE collection) are not supported yet!");
2745 MFEM_VERIFY(dof_ord != NULL,
2746 "FE collection '" << fec->
Name()
2747 <<
"' does not define reordering for " 2748 << Geometry::Name[geom] <<
" elements!");
2751 for (
int j = 0; j < new_dofs.
Size(); j++)
2754 const int old_j = dof_ord[j];
2755 new_dofs[old_j] = offset + j;
2757 offset += new_dofs.
Size();
2760 for (
int j = 0; j < old_dofs.
Size(); j++)
2762 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2774 GetElementToFaceTable();
2777 CheckBdrElementOrientation();
2782 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2787 CheckBdrElementOrientation();
2792 last_operation = Mesh::NONE;
2797 void Mesh::SetPatchAttribute(
int i,
int attr)
2799 MFEM_ASSERT(NURBSext,
"SetPatchAttribute is only for NURBS meshes");
2800 NURBSext->SetPatchAttribute(i, attr);
2801 const Array<int>& elems = NURBSext->GetPatchElements(i);
2802 for (
auto e : elems)
2804 SetAttribute(e, attr);
2808 int Mesh::GetPatchAttribute(
int i)
const 2810 MFEM_ASSERT(NURBSext,
"GetPatchAttribute is only for NURBS meshes");
2811 return NURBSext->GetPatchAttribute(i);
2814 void Mesh::SetPatchBdrAttribute(
int i,
int attr)
2816 MFEM_ASSERT(NURBSext,
"SetPatchBdrAttribute is only for NURBS meshes");
2817 NURBSext->SetPatchBdrAttribute(i, attr);
2819 const Array<int>& bdryelems = NURBSext->GetPatchBdrElements(i);
2820 for (
auto be : bdryelems)
2822 SetBdrAttribute(be, attr);
2826 int Mesh::GetPatchBdrAttribute(
int i)
const 2828 MFEM_ASSERT(NURBSext,
"GetBdrPatchBdrAttribute is only for NURBS meshes");
2829 return NURBSext->GetPatchBdrAttribute(i);
2832 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2835 CheckElementOrientation(fix_orientation);
2837 if (NumOfBdrElements == 0)
2839 GetElementToFaceTable();
2841 GenerateBoundaryElements();
2846 DSTable v_to_v(NumOfVertices);
2847 GetVertexToVertexTable(v_to_v);
2848 MarkTetMeshForRefinement(v_to_v);
2851 GetElementToFaceTable();
2854 CheckBdrElementOrientation();
2856 if (generate_edges == 1)
2858 el_to_edge =
new Table;
2859 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2873 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2874 bool fix_orientation)
2877 CheckElementOrientation(fix_orientation);
2879 if (NumOfBdrElements == 0)
2881 GetElementToFaceTable();
2883 GenerateBoundaryElements();
2886 GetElementToFaceTable();
2889 CheckBdrElementOrientation();
2891 if (generate_edges == 1)
2893 el_to_edge =
new Table;
2894 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2908 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2911 CheckElementOrientation(fix_orientation);
2913 GetElementToFaceTable();
2916 if (NumOfBdrElements == 0)
2918 GenerateBoundaryElements();
2921 CheckBdrElementOrientation();
2925 el_to_edge =
new Table;
2926 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2938 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2942 Finalize(refine, fix_orientation);
2945 void Mesh::FinalizeTopology(
bool generate_bdr)
2957 bool generate_edges =
true;
2959 if (spaceDim == 0) { spaceDim = Dim; }
2960 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2964 if (tmp_vertex_parents.Size())
2966 MFEM_VERIFY(ncmesh == NULL,
"");
2967 ncmesh =
new NCMesh(
this);
2971 InitFromNCMesh(*ncmesh);
2972 ncmesh->OnMeshUpdated(
this);
2973 GenerateNCFaceInfo();
2977 tmp_vertex_parents.DeleteAll();
2987 GetElementToFaceTable();
2989 if (NumOfBdrElements == 0 && generate_bdr)
2991 GenerateBoundaryElements();
2992 GetElementToFaceTable();
3001 if (Dim > 1 && generate_edges)
3004 if (!el_to_edge) { el_to_edge =
new Table; }
3005 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3009 if (NumOfBdrElements == 0 && generate_bdr)
3011 GenerateBoundaryElements();
3023 if (NumOfBdrElements == 0 && generate_bdr)
3025 GenerateBoundaryElements();
3032 ncmesh->OnMeshUpdated(
this);
3035 GenerateNCFaceInfo();
3042 void Mesh::Finalize(
bool refine,
bool fix_orientation)
3044 if (NURBSext || ncmesh)
3046 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
3047 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
3056 const bool check_orientation =
true;
3057 const bool curved = (Nodes != NULL);
3058 const bool may_change_topology =
3059 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
3060 ( check_orientation && fix_orientation &&
3061 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
3064 Table *old_elem_vert = NULL;
3066 if (curved && may_change_topology)
3068 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3071 if (check_orientation)
3074 CheckElementOrientation(fix_orientation);
3078 MarkForRefinement();
3081 if (may_change_topology)
3085 DoNodeReorder(old_v_to_v, old_elem_vert);
3086 delete old_elem_vert;
3099 CheckBdrElementOrientation();
3104 if (Dim >= 2 && Dim == spaceDim)
3106 const int num_faces = GetNumFaces();
3107 for (
int i = 0; i < num_faces; i++)
3109 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
3110 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology." 3111 " Interior face with incompatible orientations.");
3118 double sx,
double sy,
double sz,
bool sfc_ordering)
3122 int NVert, NElem, NBdrElem;
3124 NVert = (nx+1) * (ny+1) * (nz+1);
3125 NElem = nx * ny * nz;
3126 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3127 if (type == Element::TETRAHEDRON)
3132 else if (type == Element::WEDGE)
3135 NBdrElem += 2*nx*ny;
3137 else if (type == Element::PYRAMID)
3140 NVert += nx * ny * nz;
3143 InitMesh(3, 3, NVert, NElem, NBdrElem);
3149 for (z = 0; z <= nz; z++)
3151 coord[2] = ((double) z / nz) * sz;
3152 for (y = 0; y <= ny; y++)
3154 coord[1] = ((double) y / ny) * sy;
3155 for (x = 0; x <= nx; x++)
3157 coord[0] = ((double) x / nx) * sx;
3162 if (type == Element::PYRAMID)
3164 for (z = 0; z < nz; z++)
3166 coord[2] = (((double) z + 0.5) / nz) * sz;
3167 for (y = 0; y < ny; y++)
3169 coord[1] = (((double) y + 0.5 ) / ny) * sy;
3170 for (x = 0; x < nx; x++)
3172 coord[0] = (((double) x + 0.5 ) / nx) * sx;
3179 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1)) 3180 #define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx) 3183 if (sfc_ordering && type == Element::HEXAHEDRON)
3186 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
3187 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3189 for (
int k = 0; k < nx*ny*nz; k++)
3196 ind[0] = VTX(x , y , z );
3197 ind[1] = VTX(x+1, y , z );
3198 ind[2] = VTX(x+1, y+1, z );
3199 ind[3] = VTX(x , y+1, z );
3200 ind[4] = VTX(x , y , z+1);
3201 ind[5] = VTX(x+1, y , z+1);
3202 ind[6] = VTX(x+1, y+1, z+1);
3203 ind[7] = VTX(x , y+1, z+1);
3211 for (z = 0; z < nz; z++)
3213 for (y = 0; y < ny; y++)
3215 for (x = 0; x < nx; x++)
3218 ind[0] = VTX(x , y , z );
3219 ind[1] = VTX(x+1, y , z );
3220 ind[2] = VTX(x+1, y+1, z );
3221 ind[3] = VTX(x , y+1, z );
3222 ind[4] = VTX(x , y , z+1);
3223 ind[5] = VTX(x+1, y , z+1);
3224 ind[6] = VTX(x+1, y+1, z+1);
3225 ind[7] = VTX( x, y+1, z+1);
3227 if (type == Element::TETRAHEDRON)
3229 AddHexAsTets(ind, 1);
3231 else if (type == Element::WEDGE)
3233 AddHexAsWedges(ind, 1);
3235 else if (type == Element::PYRAMID)
3237 ind[8] = VTXP(x, y, z);
3238 AddHexAsPyramids(ind, 1);
3251 for (y = 0; y < ny; y++)
3253 for (x = 0; x < nx; x++)
3256 ind[0] = VTX(x , y , 0);
3257 ind[1] = VTX(x , y+1, 0);
3258 ind[2] = VTX(x+1, y+1, 0);
3259 ind[3] = VTX(x+1, y , 0);
3261 if (type == Element::TETRAHEDRON)
3263 AddBdrQuadAsTriangles(ind, 1);
3265 else if (type == Element::WEDGE)
3267 AddBdrQuadAsTriangles(ind, 1);
3276 for (y = 0; y < ny; y++)
3278 for (x = 0; x < nx; x++)
3281 ind[0] = VTX(x , y , nz);
3282 ind[1] = VTX(x+1, y , nz);
3283 ind[2] = VTX(x+1, y+1, nz);
3284 ind[3] = VTX(x , y+1, nz);
3286 if (type == Element::TETRAHEDRON)
3288 AddBdrQuadAsTriangles(ind, 6);
3290 else if (type == Element::WEDGE)
3292 AddBdrQuadAsTriangles(ind, 6);
3301 for (z = 0; z < nz; z++)
3303 for (y = 0; y < ny; y++)
3306 ind[0] = VTX(0 , y , z );
3307 ind[1] = VTX(0 , y , z+1);
3308 ind[2] = VTX(0 , y+1, z+1);
3309 ind[3] = VTX(0 , y+1, z );
3311 if (type == Element::TETRAHEDRON)
3313 AddBdrQuadAsTriangles(ind, 5);
3322 for (z = 0; z < nz; z++)
3324 for (y = 0; y < ny; y++)
3327 ind[0] = VTX(nx, y , z );
3328 ind[1] = VTX(nx, y+1, z );
3329 ind[2] = VTX(nx, y+1, z+1);
3330 ind[3] = VTX(nx, y , z+1);
3332 if (type == Element::TETRAHEDRON)
3334 AddBdrQuadAsTriangles(ind, 3);
3343 for (x = 0; x < nx; x++)
3345 for (z = 0; z < nz; z++)
3348 ind[0] = VTX(x , 0, z );
3349 ind[1] = VTX(x+1, 0, z );
3350 ind[2] = VTX(x+1, 0, z+1);
3351 ind[3] = VTX(x , 0, z+1);
3353 if (type == Element::TETRAHEDRON)
3355 AddBdrQuadAsTriangles(ind, 2);
3364 for (x = 0; x < nx; x++)
3366 for (z = 0; z < nz; z++)
3369 ind[0] = VTX(x , ny, z );
3370 ind[1] = VTX(x , ny, z+1);
3371 ind[2] = VTX(x+1, ny, z+1);
3372 ind[3] = VTX(x+1, ny, z );
3374 if (type == Element::TETRAHEDRON)
3376 AddBdrQuadAsTriangles(ind, 4);
3388 ofstream test_stream(
"debug.mesh");
3390 test_stream.close();
3399 double sx,
double sy,
3400 bool generate_edges,
bool sfc_ordering)
3409 if (type == Element::QUADRILATERAL)
3411 NumOfVertices = (nx+1) * (ny+1);
3412 NumOfElements = nx * ny;
3413 NumOfBdrElements = 2 * nx + 2 * ny;
3415 vertices.SetSize(NumOfVertices);
3416 elements.SetSize(NumOfElements);
3417 boundary.SetSize(NumOfBdrElements);
3424 for (j = 0; j < ny+1; j++)
3426 cy = ((double) j / ny) * sy;
3427 for (i = 0; i < nx+1; i++)
3429 cx = ((double) i / nx) * sx;
3430 vertices[k](0) = cx;
3431 vertices[k](1) = cy;
3440 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
3441 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
3443 for (k = 0; k < nx*ny; k++)
3447 ind[0] = i + j*(nx+1);
3448 ind[1] = i + 1 +j*(nx+1);
3449 ind[2] = i + 1 + (j+1)*(nx+1);
3450 ind[3] = i + (j+1)*(nx+1);
3457 for (j = 0; j < ny; j++)
3459 for (i = 0; i < nx; i++)
3461 ind[0] = i + j*(nx+1);
3462 ind[1] = i + 1 +j*(nx+1);
3463 ind[2] = i + 1 + (j+1)*(nx+1);
3464 ind[3] = i + (j+1)*(nx+1);
3473 for (i = 0; i < nx; i++)
3475 boundary[i] =
new Segment(i, i+1, 1);
3476 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3479 for (j = 0; j < ny; j++)
3481 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3482 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3486 else if (type == Element::TRIANGLE)
3488 NumOfVertices = (nx+1) * (ny+1);
3489 NumOfElements = 2 * nx * ny;
3490 NumOfBdrElements = 2 * nx + 2 * ny;
3492 vertices.SetSize(NumOfVertices);
3493 elements.SetSize(NumOfElements);
3494 boundary.SetSize(NumOfBdrElements);
3501 for (j = 0; j < ny+1; j++)
3503 cy = ((double) j / ny) * sy;
3504 for (i = 0; i < nx+1; i++)
3506 cx = ((double) i / nx) * sx;
3507 vertices[k](0) = cx;
3508 vertices[k](1) = cy;
3515 for (j = 0; j < ny; j++)
3517 for (i = 0; i < nx; i++)
3519 ind[0] = i + j*(nx+1);
3520 ind[1] = i + 1 + (j+1)*(nx+1);
3521 ind[2] = i + (j+1)*(nx+1);
3524 ind[1] = i + 1 + j*(nx+1);
3525 ind[2] = i + 1 + (j+1)*(nx+1);
3533 for (i = 0; i < nx; i++)
3535 boundary[i] =
new Segment(i, i+1, 1);
3536 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3539 for (j = 0; j < ny; j++)
3541 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3542 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3549 MFEM_ABORT(
"Unsupported element type.");
3553 CheckElementOrientation();
3555 if (generate_edges == 1)
3557 el_to_edge =
new Table;
3558 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3560 CheckBdrElementOrientation();
3569 attributes.Append(1);
3570 bdr_attributes.Append(1); bdr_attributes.Append(2);
3571 bdr_attributes.Append(3); bdr_attributes.Append(4);
3576 void Mesh::Make1D(
int n,
double sx)
3585 NumOfVertices = n + 1;
3587 NumOfBdrElements = 2;
3588 vertices.SetSize(NumOfVertices);
3589 elements.SetSize(NumOfElements);
3590 boundary.SetSize(NumOfBdrElements);
3593 for (j = 0; j < n+1; j++)
3595 vertices[j](0) = ((
double) j / n) * sx;
3599 for (j = 0; j < n; j++)
3601 elements[j] =
new Segment(j, j+1, 1);
3606 boundary[0] =
new Point(ind, 1);
3608 boundary[1] =
new Point(ind, 2);
3616 attributes.Append(1);
3617 bdr_attributes.Append(1); bdr_attributes.Append(2);
3620 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
3638 last_operation = Mesh::NONE;
3641 elements.SetSize(NumOfElements);
3642 for (
int i = 0; i < NumOfElements; i++)
3644 elements[i] = mesh.
elements[i]->Duplicate(
this);
3651 boundary.SetSize(NumOfBdrElements);
3652 for (
int i = 0; i < NumOfBdrElements; i++)
3654 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
3673 faces.SetSize(mesh.
faces.Size());
3674 for (
int i = 0; i < faces.Size(); i++)
3677 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
3687 face_to_elem = NULL;
3712 if (dynamic_cast<const ParMesh*>(&mesh))
3724 if (mesh.
Nodes && copy_nodes)
3729 FiniteElementCollection::New(fec->
Name());
3733 Nodes->MakeOwner(fec_copy);
3734 *Nodes = *mesh.
Nodes;
3756 int refine,
bool fix_orientation)
3760 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
3761 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
3775 double sx,
double sy,
bool sfc_ordering)
3778 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
3785 double sx,
double sy,
double sz,
bool sfc_ordering)
3788 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
3797 ref_factors = ref_factor;
3810 Mesh::Mesh(
const std::string &filename,
int generate_edges,
int refine,
3811 bool fix_orientation)
3820 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
3824 Load(imesh, generate_edges, refine, fix_orientation);
3829 bool fix_orientation)
3832 Load(input, generate_edges, refine, fix_orientation);
3841 "Not enough vertices in external array : " 3842 "len_vertex_data = "<< len_vertex_data <<
", " 3845 if (vertex_data == (
double *)(
vertices.GetData()))
3847 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
3852 memcpy(vertex_data,
vertices.GetData(),
3861 int *element_attributes,
int num_elements,
3863 int *boundary_attributes,
int num_boundary_elements,
3866 if (space_dimension == -1)
3872 num_boundary_elements);
3875 int boundary_index_stride = num_boundary_elements > 0 ?
3879 vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
3882 for (
int i = 0; i < num_elements; i++)
3885 elements[i]->SetVertices(element_indices + i * element_index_stride);
3886 elements[i]->SetAttribute(element_attributes[i]);
3890 for (
int i = 0; i < num_boundary_elements; i++)
3893 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3894 boundary[i]->SetAttribute(boundary_attributes[i]);
3910 #ifdef MFEM_USE_MEMALLOC 3919 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3932 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3935 for (
int i = 0; i < nv; i++)
3948 for (
int j = 0; j < nv; j++)
4020 MFEM_ABORT(
"invalid element type: " << type);
4027 std::string parse_tag)
4029 int curved = 0, read_gf = 1;
4030 bool finalize_topo =
true;
4034 MFEM_ABORT(
"Input stream is not open");
4041 getline(input, mesh_type);
4045 int mfem_version = 0;
4046 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
4047 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
4051 int mfem_nc_version = 0;
4052 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
4053 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
4061 if (mfem_version == 12 && parse_tag.empty())
4063 parse_tag =
"mfem_mesh_end";
4067 else if (mfem_nc_version)
4069 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4076 MFEM_VERIFY(mfem_nc_version >= 10,
4077 "Legacy nonconforming format (MFEM mesh v1.1) cannot be " 4078 "used to load a parallel nonconforming mesh, sorry.");
4081 input, mfem_nc_version, curved, is_nc);
4086 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
4099 else if (mesh_type ==
"linemesh")
4103 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
4105 if (mesh_type ==
"curved_areamesh2")
4111 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
4115 else if (mesh_type ==
"TrueGrid")
4119 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
4121 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
4123 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4124 "Unsupported VTK format");
4125 ReadVTKMesh(input, curved, read_gf, finalize_topo);
4127 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
4131 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
4135 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
4140 else if (mesh_type ==
"$MeshFormat")
4145 ((mesh_type.size() > 2 &&
4146 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
4147 (mesh_type.size() > 3 &&
4148 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
4153 #ifdef MFEM_USE_NETCDF 4156 MFEM_ABORT(
"NetCDF support requires configuration with" 4157 " MFEM_USE_NETCDF=YES");
4163 MFEM_ABORT(
"Can not determine Cubit mesh filename!" 4164 " Use mfem::named_ifgzstream for input.");
4170 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
4196 bool generate_bdr =
false;
4201 if (curved && read_gf)
4215 if (mfem_version == 12)
4221 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
4222 getline(input, line);
4228 if (line ==
"mfem_mesh_end") {
break; }
4230 while (line != parse_tag);
4232 else if (mfem_nc_version >= 10)
4237 MFEM_VERIFY(ident ==
"mfem_mesh_end",
4238 "invalid mesh: end of file tag not found");
4246 int i, j, ie, ib, iv, *v, nv;
4277 for (i = 0; i < num_pieces; i++)
4284 for (i = 0; i < num_pieces; i++)
4290 for (j = 0; j < m->
GetNE(); j++)
4295 for (j = 0; j < m->
GetNBE(); j++)
4300 for (
int k = 0; k < nv; k++)
4302 v[k] = lvert_vert[v[k]];
4307 for (j = 0; j < m->
GetNV(); j++)
4319 for (i = 0; i < num_pieces; i++)
4330 for (i = 0; i < num_pieces; i++)
4334 for (j = 0; j < m->
GetNE(); j++)
4339 for (
int k = 0; k < nv; k++)
4346 for (j = 0; j < m->
GetNBE(); j++)
4351 for (
int k = 0; k < nv; k++)
4358 for (j = 0; j < m->
GetNV(); j++)
4372 for (i = 0; i < num_pieces; i++)
4374 gf_array[i] = mesh_array[i]->
GetNodes();
4389 ref_factors = ref_factor;
4400 int orig_ne = orig_mesh.
GetNE();
4401 MFEM_VERIFY(ref_factors.
Size() == orig_ne && orig_ne > 0,
4402 "Number of refinement factors must equal number of elements")
4403 MFEM_VERIFY(ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
4406 "Invalid refinement type. Must use closed basis type.");
4408 int min_ref = ref_factors.
Min();
4409 int max_ref = ref_factors.
Max();
4411 bool var_order = (min_ref != max_ref);
4424 for (
int i = 0; i < orig_ne; i++)
4442 for (
int el = 0; el < orig_ne; el++)
4454 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
4455 for (
int i = 0; i < phys_pts.
Width(); i++)
4464 for (
int k = 0; k < nvert; k++)
4467 v[k] = rdofs[c2h_map[cid]];
4480 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
4491 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
4497 for (
int k = 0; k < nvert; k++)
4500 v[k] = rdofs[c2h_map[cid]];
4522 for (
int iel = 0; iel < orig_ne; iel++)
4531 const int *node_map = NULL;
4534 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
4535 const int *vertex_map = vertex_fec.
GetDofMap(geom);
4536 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
4537 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
4540 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
4543 int iv = vertex_map[iv_lex];
4545 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
4547 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
4550 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
4561 using GeomRef = std::pair<Geometry::Type, int>;
4562 std::map<GeomRef, int> point_matrices_offsets;
4564 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4568 GeomRef id(geom, ref_factors[el_coarse]);
4569 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
4575 point_matrices_offsets[id] = n_point_matrices[geom];
4576 n_point_matrices[geom] += nref_el;
4583 int nmatrices = n_point_matrices[geom];
4590 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4593 int ref = ref_factors[el_coarse];
4594 int offset = point_matrices_offsets[GeomRef(geom, ref)];
4600 for (
int k = 0; k < nvert; k++)
4632 "Mesh::MakeSimplicial requires a properly oriented input mesh");
4634 "Mesh::MakeSimplicial does not support non-conforming meshes.")
4641 Mesh copy(orig_mesh);
4646 int nv = orig_mesh.
GetNV();
4647 int ne = orig_mesh.
GetNE();
4648 int nbe = orig_mesh.
GetNBE();
4661 int new_ne = 0, new_nbe = 0;
4662 for (
int i=0; i<ne; ++i)
4666 for (
int i=0; i<nbe; ++i)
4675 for (
int i=0; i<nv; ++i)
4685 if (vglobal == NULL)
4688 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
4689 vglobal = vglobal_id.
GetData();
4692 constexpr
int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
4693 constexpr
int quad_ntris = 2, prism_ntets = 3;
4694 static const int quad_trimap[2][nv_tri*quad_ntris] =
4706 static const int prism_rot[nv_prism*nv_prism] =
4715 static const int prism_f[nv_quad] = {1, 2, 5, 4};
4716 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
4730 static const int hex_rot[nv_hex*nv_hex] =
4732 0, 1, 2, 3, 4, 5, 6, 7,
4733 1, 0, 4, 5, 2, 3, 7, 6,
4734 2, 1, 5, 6, 3, 0, 4, 7,
4735 3, 0, 1, 2, 7, 4, 5, 6,
4736 4, 0, 3, 7, 5, 1, 2, 6,
4737 5, 1, 0, 4, 6, 2, 3, 7,
4738 6, 2, 1, 5, 7, 3, 0, 4,
4739 7, 3, 2, 6, 4, 0, 1, 5
4741 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
4742 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
4743 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
4744 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
4745 static const int hex_tetmap0[nv_tet*5] =
4752 static const int hex_tetmap1[nv_tet*6] =
4759 static const int hex_tetmap2[nv_tet*6] =
4766 static const int hex_tetmap3[nv_tet*6] =
4773 static const int *hex_tetmaps[4] =
4775 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
4778 auto find_min = [](
const int*
a,
int n) {
return std::min_element(
a,
a+n)-
a; };
4780 for (
int i=0; i<ne; ++i)
4782 const int *v = orig_mesh.
elements[i]->GetVertices();
4786 if (num_subdivisions[orig_geom] == 1)
4798 for (
int itri=0; itri<quad_ntris; ++itri)
4803 for (
int iv=0; iv<nv_tri; ++iv)
4805 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
4813 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
4816 int irot = find_min(vg, nv_prism);
4817 for (
int iv=0; iv<nv_prism; ++iv)
4819 int jv = prism_rot[iv + irot*nv_prism];
4824 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
4825 int j = find_min(q, nv_quad);
4826 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
4827 for (
int itet=0; itet<prism_ntets; ++itet)
4832 for (
int iv=0; iv<nv_tet; ++iv)
4834 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
4842 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
4846 int irot = find_min(vg, nv_hex);
4847 for (
int iv=0; iv<nv_hex; ++iv)
4849 int jv = hex_rot[iv + irot*nv_hex];
4859 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
4860 j = find_min(q, nv_quad);
4861 if (j == 0 || j == 2) { bitmask += 4; }
4863 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
4864 j = find_min(q, nv_quad);
4865 if (j == 1 || j == 3) { bitmask += 2; }
4867 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
4868 j = find_min(q, nv_quad);
4869 if (j == 0 || j == 2) { bitmask += 1; }
4872 int nrot = num_rot[bitmask];
4873 for (
int k=0; k<nrot; ++k)
4887 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
4888 int ntets = (ndiags == 0) ? 5 : 6;
4889 const int *tetmap = hex_tetmaps[ndiags];
4890 for (
int itet=0; itet<ntets; ++itet)
4895 for (
int iv=0; iv<nv_tet; ++iv)
4897 v2[iv] = vg[tetmap[itet + iv*ntets]];
4906 for (
int i=0; i<nbe; ++i)
4908 const int *v = orig_mesh.
boundary[i]->GetVertices();
4911 if (num_subdivisions[orig_geom] == 1)
4921 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
4923 int iv_min = find_min(vg, nv_quad);
4924 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
4925 for (
int itri=0; itri<quad_ntris; ++itri)
4930 for (
int iv=0; iv<nv_tri; ++iv)
4932 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
4939 MFEM_ABORT(
"Unreachable");
4953 Mesh periodic_mesh(orig_mesh,
true);
4959 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
4964 for (
int j = 0; j < nv; j++)
4970 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
4975 for (
int j = 0; j < nv; j++)
4982 return periodic_mesh;
4986 const std::vector<Vector> &translations,
double tol)
const 4990 Vector coord(sdim), at(sdim), dx(sdim);
4991 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
4992 xMax = xMin = xDiff = 0.0;
4996 for (
int be = 0; be <
GetNBE(); be++)
5001 for (
int i = 0; i < dofs.
Size(); i++)
5003 bdr_v.insert(dofs[i]);
5006 for (
int j = 0; j < sdim; j++)
5008 xMax[j] = max(xMax[j], coord[j]);
5009 xMin[j] = min(xMin[j], coord[j]);
5013 add(xMax, -1.0, xMin, xDiff);
5014 double dia = xDiff.
Norml2();
5022 map<int, int> replica2primary;
5024 map<int, set<int>> primary2replicas;
5028 for (
int v : bdr_v) { primary2replicas[v]; }
5032 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
5034 if (r ==
p) {
return; }
5035 primary2replicas[
p].insert(r);
5036 replica2primary[r] =
p;
5037 for (
int s : primary2replicas[r])
5039 primary2replicas[
p].insert(
s);
5040 replica2primary[
s] =
p;
5042 primary2replicas.erase(r);
5045 for (
unsigned int i = 0; i < translations.size(); i++)
5047 for (
int vi : bdr_v)
5050 add(coord, translations[i], at);
5052 for (
int vj : bdr_v)
5055 add(at, -1.0, coord, dx);
5056 if (dx.
Norml2() > dia*tol) {
continue; }
5061 bool pi = primary2replicas.find(vi) != primary2replicas.end();
5062 bool pj = primary2replicas.find(vj) != primary2replicas.end();
5068 make_replica(vj, vi);
5073 int owner_of_vj = replica2primary[vj];
5075 make_replica(vi, owner_of_vj);
5081 int owner_of_vi = replica2primary[vi];
5082 make_replica(vj, owner_of_vi);
5089 int owner_of_vi = replica2primary[vi];
5090 int owner_of_vj = replica2primary[vj];
5091 make_replica(owner_of_vj, owner_of_vi);
5098 std::vector<int> v2v(
GetNV());
5099 for (
size_t i = 0; i < v2v.size(); i++)
5103 for (
auto &&r2p : replica2primary)
5105 v2v[r2p.first] = r2p.second;
5112 MFEM_VERIFY(
NURBSext,
"Mesh::RefineNURBSFromFile: Not a NURBS mesh!");
5113 mfem::out<<
"Refining NURBS from refinement file: "<<ref_file<<endl;
5116 ifstream input(ref_file);
5123 mfem::out<<
"Knot vectors in ref_file: "<<nkv<<endl;
5125 MFEM_ABORT(
"Refine file does not have the correct number of knot vectors");
5130 for (
int kv = 0; kv < nkv; kv++)
5132 knotVec[kv] =
new Vector();
5133 knotVec[kv]->
Load(input);
5141 for (
int kv = 0; kv < nkv; kv++)
5151 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5156 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5173 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5178 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5208 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
5232 for (
int i = 0; i <
elements.Size(); i++)
5242 for (
int i = 0; i <
boundary.Size(); i++)
5259 for (
int i = 0; i < vd; i++)
5326 input >> edge_to_knot[j] >> v[0] >> v[1];
5329 edge_to_knot[j] = -1 - edge_to_knot[j];
5349 if (edge_to_knot.
Size() == 0)
5353 constexpr
int notset = -9999999;
5354 edge_to_knot = notset;
5366 edge0[0] = 0; edge1[0] = 2;
5367 edge0[1] = 1; edge1[1] = 3;
5375 edge0[0] = 0; edge1[0] = 2;
5376 edge0[1] = 0; edge1[1] = 4;
5377 edge0[2] = 0; edge1[2] = 6;
5379 edge0[3] = 1; edge1[3] = 3;
5380 edge0[4] = 1; edge1[4] = 5;
5381 edge0[5] = 1; edge1[5] = 7;
5383 edge0[6] = 8; edge1[6] = 9;
5384 edge0[7] = 8; edge1[7] = 10;
5385 edge0[8] = 8; edge1[8] = 11;
5393 int e0, e1, v0, v1, df;
5399 const int *v =
elements[
p]->GetVertices();
5400 for (j = 0; j < edges.
Size(); j++)
5403 const int *e =
elements[
p]->GetEdgeVertices(j);
5416 for (j = 0; j < edge1.
Size(); j++)
5418 e0 = edges[edge0[j]];
5419 e1 = edges[edge1[j]];
5420 v0 = edge_to_knot[e0];
5421 v1 = edge_to_knot[e1];
5422 df = flip*oedge[edge0[j]]*oedge[edge1[j]];
5425 if ((v0 == notset) && (v1 == notset))
5427 edge_to_knot[e0] = knot;
5428 edge_to_knot[e1] = knot;
5434 else if ((v0 != notset) && (v1 == notset))
5436 edge_to_knot[e1] = (df >= 0 ? -v0-1 : v0);
5438 else if ((v0 == notset) && (v1 != notset))
5440 edge_to_knot[e0] = (df >= 0 ? -v1-1 : v1);
5460 for (j = 0; j < edge1.
Size(); j++)
5462 e0 = edges[edge0[j]];
5463 e1 = edges[edge1[j]];
5464 v0 = edge_to_knot[e0];
5465 v1 = edge_to_knot[e1];
5466 v0 = ( v0 >= 0 ? v0 : -v0-1);
5467 v1 = ( v1 >= 0 ? v1 : -v1-1);
5473 edge_to_knot[e1] = (oedge[edge1[j]] >= 0 ? v0 : -v0-1);
5477 edge_to_knot[e0] = (oedge[edge0[j]] >= 0 ? v1 : -v1-1);
5485 while (corrections > 0 && passes <
GetNE() + 1);
5488 if (corrections > 0 )
5490 mfem::err<<
"Edge_to_knot mapping potentially incorrect"<<endl;
5492 mfem::err<<
" corrections = "<<corrections<<endl;
5502 k = edge_to_knot[j];
5503 cnt[(k >= 0 ? k : -k-1)]++;
5507 for (j = 0; j < cnt.
Size(); j++)
5509 cnt[j] = (cnt[j] > 0 ? k++ : -1);
5514 k = edge_to_knot[j];
5515 edge_to_knot[j] = (k >= 0 ? cnt[k]:-cnt[-k-1]-1);
5519 mfem::out<<
"Generated edge to knot mapping:"<<endl;
5523 k = edge_to_knot[j];
5532 mfem::out<<(k >= 0 ? k:-k-1)<<
" "<< v[0] <<
" "<<v[1]<<endl;
5536 if (corrections > 0 ) {
mfem_error(
"Mesh::LoadPatchTopo");}
5542 if (
p.Size() >= v.
Size())
5544 for (
int d = 0; d < v.
Size(); d++)
5552 for (d = 0; d <
p.Size(); d++)
5556 for ( ; d < v.
Size(); d++)
5588 if (dynamic_cast<const H1_FECollection*>(fec)
5589 || dynamic_cast<const L2_FECollection*>(fec))
5598 #ifndef MFEM_USE_MPI 5599 const bool warn =
true;
5602 const bool warn = !pmesh || pmesh->
GetMyRank() == 0;
5606 MFEM_WARNING(
"converting NURBS mesh to order " << order <<
5607 " H1-continuous mesh!\n " 5608 "If this is the desired behavior, you can silence" 5609 " this warning by converting\n " 5610 "the NURBS mesh to high-order mesh in advance by" 5611 " calling the method\n " 5612 "Mesh::SetCurvature().");
5637 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
5656 MFEM_ASSERT(nodes != NULL,
"");
5672 case 1:
return GetNV();
5708 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG)) 5709 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
5714 int i, j, k, wo = 0, fo = 0;
5723 int *vi =
elements[i]->GetVertices();
5726 for (j = 0; j < 3; j++)
5730 for (j = 0; j < 2; j++)
5731 for (k = 0; k < 2; k++)
5733 J(j, k) = v[j+1][k] - v[0][k];
5754 MFEM_ABORT(
"Invalid 2D element type \"" 5771 int *vi =
elements[i]->GetVertices();
5777 for (j = 0; j < 4; j++)
5781 for (j = 0; j < 3; j++)
5782 for (k = 0; k < 3; k++)
5784 J(j, k) = v[j+1][k] - v[0][k];
5843 MFEM_ABORT(
"Invalid 3D element type \"" 5849 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG)) 5852 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / " 5853 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
5857 MFEM_CONTRACT_VAR(fo);
5870 if (test[0] == base[0])
5871 if (test[1] == base[1])
5879 else if (test[0] == base[1])
5880 if (test[1] == base[0])
5889 if (test[1] == base[0])
5900 for (
int j = 0; j < 3; j++)
5901 if (test[aor[j]] != base[j])
5903 mfem::err <<
"Mesh::GetTriOrientation(...)" << endl;
5905 for (
int k = 0; k < 3; k++)
5910 for (
int k = 0; k < 3; k++)
5931 const int oo[6][6] =
5941 int ori_a_c = oo[ori_a_b][ori_b_c];
5947 const int inv_ori[6] = {0, 1, 4, 3, 2, 5};
5948 return inv_ori[ori];
5955 for (i = 0; i < 4; i++)
5956 if (test[i] == base[0])
5963 if (test[(i+1)%4] == base[1])
5972 for (
int j = 0; j < 4; j++)
5973 if (test[aor[j]] != base[j])
5975 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
5977 for (
int k = 0; k < 4; k++)
5982 for (
int k = 0; k < 4; k++)
5991 if (test[(i+1)%4] == base[1])
6008 const int oo[8][8] =
6010 {0, 1, 2, 3, 4, 5, 6, 7},
6011 {1, 0, 3, 2, 5, 4, 7, 6},
6012 {2, 7, 4, 1, 6, 3, 0, 5},
6013 {3, 6, 5, 0, 7, 2, 1, 4},
6014 {4, 5, 6, 7, 0, 1, 2, 3},
6015 {5, 4, 7, 6, 1, 0, 3, 2},
6016 {6, 3, 0, 5, 2, 7, 4, 1},
6017 {7, 2, 1, 4, 3, 6, 5, 0}
6020 int ori_a_c = oo[ori_a_b][ori_b_c];
6026 const int inv_ori[8] = {0, 1, 6, 3, 4, 5, 2, 7};
6027 return inv_ori[ori];
6038 if (test[0] == base[0])
6039 if (test[1] == base[1])
6040 if (test[2] == base[2])
6048 else if (test[2] == base[1])
6049 if (test[3] == base[2])
6058 if (test[1] == base[2])
6066 else if (test[1] == base[0])
6067 if (test[2] == base[1])
6068 if (test[0] == base[2])
6076 else if (test[3] == base[1])
6077 if (test[2] == base[2])
6086 if (test[3] == base[2])
6094 else if (test[2] == base[0])
6095 if (test[3] == base[1])
6096 if (test[0] == base[2])
6104 else if (test[0] == base[1])
6105 if (test[1] == base[2])
6114 if (test[3] == base[2])
6123 if (test[0] == base[1])
6124 if (test[2] == base[2])
6132 else if (test[1] == base[1])
6133 if (test[0] == base[2])
6142 if (test[1] == base[2])
6153 for (
int j = 0; j < 4; j++)
6154 if (test[aor[j]] != base[j])
6179 int *bv =
boundary[i]->GetVertices();
6185 mfem::Swap<int>(bv[0], bv[1]);
6199 if (
faces_info[fi].Elem2No >= 0) {
continue; }
6202 int *bv =
boundary[i]->GetVertices();
6204 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
6205 const int *fv =
faces[fi]->GetVertices();
6222 MFEM_ABORT(
"Invalid 2D boundary element type \"" 6223 << bdr_type <<
"\"");
6228 if (orientation % 2 == 0) {
continue; }
6230 if (!fix_it) {
continue; }
6238 mfem::Swap<int>(bv[0], bv[1]);
6242 mfem::Swap<int>(be[1], be[2]);
6248 mfem::Swap<int>(bv[0], bv[2]);
6252 mfem::Swap<int>(be[0], be[1]);
6253 mfem::Swap<int>(be[2], be[3]);
6266 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / " 6276 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
6287 MFEM_ASSERT(0 <=
dim &&
dim <=
Dim,
"invalid dim: " <<
dim);
6306 mfem_error(
"Mesh::GetElementEdges(...) element to edge table " 6307 "is not generated.");
6310 const int *v =
elements[i]->GetVertices();
6311 const int ne =
elements[i]->GetNEdges();
6313 for (
int j = 0; j < ne; j++)
6315 const int *e =
elements[i]->GetEdgeVertices(j);
6316 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6327 const int *v =
boundary[i]->GetVertices();
6328 cor[0] = (v[0] < v[1]) ? (1) : (-1);
6341 const int *v =
boundary[i]->GetVertices();
6342 const int ne =
boundary[i]->GetNEdges();
6344 for (
int j = 0; j < ne; j++)
6346 const int *e =
boundary[i]->GetEdgeVertices(j);
6347 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6359 const int *v =
faces[i]->GetVertices();
6360 o[0] = (v[0] < v[1]) ? (1) : (-1);
6372 const int *v =
faces[i]->GetVertices();
6373 const int ne =
faces[i]->GetNEdges();
6375 for (
int j = 0; j < ne; j++)
6377 const int *e =
faces[i]->GetEdgeVertices(j);
6378 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6406 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
6457 for (j = 0; j < nv; j++)
6469 for (j = 0; j < nv; j++)
6516 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
6520 int n = el_faces.
Size();
6523 for (
int j = 0; j < n; j++)
6527 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
6531 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
6532 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
6549 for (
auto f : elem_faces)
6584 MFEM_ABORT(
"invalid geometry");
6592 case 1:
return boundary[i]->GetVertices()[0];
6595 default: MFEM_ABORT(
"invalid dimension!");
6605 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6608 const int *bv =
boundary[bdr_el]->GetVertices();
6616 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6627 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6630 const int *bv =
boundary[bdr_el]->GetVertices();
6638 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6665 for (j = 0; j < nv; j++)
6667 pointmat(k, j) =
vertices[v[j]](k);
6682 for (j = 0; j < nv; j++)
6684 pointmat(k, j) =
vertices[v[j]](k);
6696 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
6699 return sqrt(length);
6707 for (
int i = 0; i < elem_array.
Size(); i++)
6712 for (
int i = 0; i < elem_array.
Size(); i++)
6714 const int *v = elem_array[i]->GetVertices();
6715 const int ne = elem_array[i]->GetNEdges();
6716 for (
int j = 0; j < ne; j++)
6718 const int *e = elem_array[i]->GetEdgeVertices(j);
6732 v_to_v.
Push(v[0], v[1]);
6739 const int *v =
elements[i]->GetVertices();
6740 const int ne =
elements[i]->GetNEdges();
6741 for (
int j = 0; j < ne; j++)
6743 const int *e =
elements[i]->GetEdgeVertices(j);
6744 v_to_v.
Push(v[e[0]], v[e[1]]);
6752 int i, NumberOfEdges;
6768 const int *v =
boundary[i]->GetVertices();
6769 be_to_f[i] = v_to_v(v[0], v[1]);
6782 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
6786 return NumberOfEdges;
6845 if (
faces[gf] == NULL)
6877 if (
faces[gf] == NULL)
6887 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. " 6888 "Interior edge found between 2D elements " 6890 <<
" and " << el <<
".");
6891 int *v =
faces[gf]->GetVertices();
6893 if ( v[1] == v0 && v[0] == v1 )
6897 else if ( v[0] == v0 && v[1] == v1 )
6907 MFEM_ABORT(
"internal error");
6913 int v0,
int v1,
int v2)
6915 if (
faces[gf] == NULL)
6925 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. " 6926 "Interior triangular face found connecting elements " 6928 <<
" and " << el <<
".");
6929 int orientation, vv[3] = { v0, v1, v2 };
6936 faces_info[gf].Elem2Inf = 64 * lf + orientation;
6941 int v0,
int v1,
int v2,
int v3)
6953 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. " 6954 "Interior quadrilateral face found connecting elements " 6956 <<
" and " << el <<
".");
6957 int vv[4] = { v0, v1, v2, v3 };
6971 for (i = 0; i <
faces.Size(); i++)
6977 faces.SetSize(nfaces);
6979 for (i = 0; i < nfaces; i++)
6987 const int *v =
elements[i]->GetVertices();
6997 const int ne =
elements[i]->GetNEdges();
6998 for (
int j = 0; j < ne; j++)
7000 const int *e =
elements[i]->GetEdgeVertices(j);
7011 for (
int j = 0; j < 4; j++)
7015 v[fv[0]], v[fv[1]], v[fv[2]]);
7021 for (
int j = 0; j < 2; j++)
7025 v[fv[0]], v[fv[1]], v[fv[2]]);
7027 for (
int j = 2; j < 5; j++)
7031 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7037 for (
int j = 0; j < 1; j++)
7041 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7043 for (
int j = 1; j < 5; j++)
7047 v[fv[0]], v[fv[1]], v[fv[2]]);
7053 for (
int j = 0; j < 6; j++)
7057 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7062 MFEM_ABORT(
"Unexpected type of Element.");
7070 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
7086 for (
int i = 0; i < list.
masters.Size(); i++)
7089 if (master.
index >= nfaces) {
continue; }
7095 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
7096 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
7100 for (
int i = 0; i < list.
slaves.Size(); i++)
7104 if (slave.
index < 0 ||
7105 slave.
index >= nfaces ||
7134 const int *v =
elements[i]->GetVertices();
7139 for (
int j = 0; j < 4; j++)
7142 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
7148 for (
int j = 0; j < 1; j++)
7151 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7153 for (
int j = 1; j < 5; j++)
7156 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
7162 for (
int j = 0; j < 2; j++)
7165 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
7167 for (
int j = 2; j < 5; j++)
7170 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7178 for (
int j = 0; j < 6; j++)
7181 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
7186 MFEM_ABORT(
"Unexpected type of Element.");
7210 for (
int j = 0; j < 4; j++)
7214 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
7220 for (
int j = 0; j < 2; j++)
7224 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
7226 for (
int j = 2; j < 5; j++)
7230 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
7236 for (
int j = 0; j < 1; j++)
7240 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
7242 for (
int j = 1; j < 5; j++)
7246 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
7254 for (
int j = 0; j < 6; j++)
7258 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
7263 MFEM_ABORT(
"Unexpected type of Element.");
7276 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
7281 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
7285 MFEM_ABORT(
"Unexpected type of boundary Element.");
7299 void Rotate3(
int &
a,
int &
b,
int &c)
7331 Table *old_elem_vert = NULL;
7342 int *v =
elements[i]->GetVertices();
7344 Rotate3(v[0], v[1], v[2]);
7347 Rotate3(v[1], v[2], v[3]);
7360 int *v =
boundary[i]->GetVertices();
7362 Rotate3(v[0], v[1], v[2]);
7378 delete old_elem_vert;
7394 if (
p[i] < pmin[i]) { pmin[i] =
p[i]; }
7395 if (
p[i] > pmax[i]) { pmax[i] =
p[i]; }
7409 for (
int i =
spaceDim-1; i >= 0; i--)
7411 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
7412 if (idx < 0) { idx = 0; }
7413 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
7414 part = part * nxyz[i] + idx;
7416 partitioning[el] = part;
7419 return partitioning;
7429 #ifdef MFEM_USE_METIS 7431 int print_messages = 1;
7434 int init_flag, fin_flag;
7435 MPI_Initialized(&init_flag);
7436 MPI_Finalized(&fin_flag);
7437 if (init_flag && !fin_flag)
7441 if (rank != 0) { print_messages = 0; }
7445 int i, *partitioning;
7455 partitioning[i] = 0;
7462 partitioning[i] = i;
7468 #ifndef MFEM_USE_METIS_5 7480 bool freedata =
false;
7482 idx_t *mpartitioning;
7485 if (
sizeof(
idx_t) ==
sizeof(int))
7489 mpartitioning = (
idx_t*) partitioning;
7498 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
7499 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
7500 mpartitioning =
new idx_t[n];
7503 #ifndef MFEM_USE_METIS_5 7506 METIS_SetDefaultOptions(options);
7507 options[METIS_OPTION_CONTIG] = 1;
7515 if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
7520 if (part_method >= 0 && part_method <= 2)
7522 for (i = 0; i < n; i++)
7528 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
7534 if (part_method == 0 || part_method == 3)
7536 #ifndef MFEM_USE_METIS_5 7565 " error in METIS_PartGraphRecursive!");
7572 if (part_method == 1 || part_method == 4)
7574 #ifndef MFEM_USE_METIS_5 7603 " error in METIS_PartGraphKway!");
7610 if (part_method == 2 || part_method == 5)
7612 #ifndef MFEM_USE_METIS_5 7625 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
7642 " error in METIS_PartGraphKway!");
7650 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = " 7654 nparts = (int) mparts;
7655 if (mpartitioning != (
idx_t*)partitioning)
7659 partitioning[k] = mpartitioning[k];
7666 delete[] mpartitioning;
7682 auto count_partition_elements = [&]()
7684 for (i = 0; i < nparts; i++)
7692 psize[partitioning[i]].one++;
7696 for (i = 0; i < nparts; i++)
7698 if (psize[i].one == 0) { empty_parts++; }
7702 count_partition_elements();
7710 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned " 7711 << empty_parts <<
" empty parts!" 7712 <<
" Applying a simple fix ..." << endl;
7715 SortPairs<int,int>(psize, nparts);
7717 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7724 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7726 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
7732 partitioning[j] = psize[nparts-1-i].two;
7739 count_partition_elements();
7743 return partitioning;
7747 mfem_error(
"Mesh::GeneratePartitioning(...): " 7748 "MFEM was compiled without Metis.");
7762 int num_elem, *i_elem_elem, *j_elem_elem;
7764 num_elem = elem_elem.
Size();
7765 i_elem_elem = elem_elem.
GetI();
7766 j_elem_elem = elem_elem.
GetJ();
7771 int stack_p, stack_top_p, elem;
7775 for (i = 0; i < num_elem; i++)
7777 if (partitioning[i] > num_part)
7779 num_part = partitioning[i];
7786 for (i = 0; i < num_part; i++)
7793 for (elem = 0; elem < num_elem; elem++)
7795 if (component[elem] >= 0)
7800 component[elem] = num_comp[partitioning[elem]]++;
7802 elem_stack[stack_top_p++] = elem;
7804 for ( ; stack_p < stack_top_p; stack_p++)
7806 i = elem_stack[stack_p];
7807 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
7810 if (partitioning[k] == partitioning[i])
7812 if (component[k] < 0)
7814 component[k] = component[i];
7815 elem_stack[stack_top_p++] = k;
7817 else if (component[k] != component[i])
7829 int i, n_empty, n_mcomp;
7837 n_empty = n_mcomp = 0;
7838 for (i = 0; i < num_comp.
Size(); i++)
7839 if (num_comp[i] == 0)
7843 else if (num_comp[i] > 1)
7850 mfem::out <<
"Mesh::CheckPartitioning(...) :\n" 7851 <<
"The following subdomains are empty :\n";
7852 for (i = 0; i < num_comp.
Size(); i++)
7853 if (num_comp[i] == 0)
7861 mfem::out <<
"Mesh::CheckPartitioning(...) :\n" 7862 <<
"The following subdomains are NOT connected :\n";
7863 for (i = 0; i < num_comp.
Size(); i++)
7864 if (num_comp[i] > 1)
7870 if (n_empty == 0 && n_mcomp == 0)
7871 mfem::out <<
"Mesh::CheckPartitioning(...) : " 7872 "All subdomains are connected." << endl;
7886 const double *
a = A.
Data();
7887 const double *
b = B.
Data();
7896 c(0) =
a[0]*
a[3]-
a[1]*
a[2];
7897 c(1) =
a[0]*
b[3]-
a[1]*
b[2]+
b[0]*
a[3]-
b[1]*
a[2];
7898 c(2) =
b[0]*
b[3]-
b[1]*
b[2];
7919 c(0) = (
a[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
7920 a[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
7921 a[2] * (
a[3] *
a[7] -
a[4] *
a[6]));
7923 c(1) = (
b[0] * (
a[4] *
a[8] -
a[5] *
a[7]) +
7924 b[1] * (
a[5] *
a[6] -
a[3] *
a[8]) +
7925 b[2] * (
a[3] *
a[7] -
a[4] *
a[6]) +
7927 a[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
7928 a[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
7929 a[2] * (
b[3] *
a[7] -
b[4] *
a[6]) +
7931 a[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
7932 a[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
7933 a[2] * (
a[3] *
b[7] -
a[4] *
b[6]));
7935 c(2) = (
a[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
7936 a[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
7937 a[2] * (
b[3] *
b[7] -
b[4] *
b[6]) +
7939 b[0] * (
a[4] *
b[8] -
a[5] *
b[7]) +
7940 b[1] * (
a[5] *
b[6] -
a[3] *
b[8]) +
7941 b[2] * (
a[3] *
b[7] -
a[4] *
b[6]) +
7943 b[0] * (
b[4] *
a[8] -
b[5] *
a[7]) +
7944 b[1] * (
b[5] *
a[6] -
b[3] *
a[8]) +
7945 b[2] * (
b[3] *
a[7] -
b[4] *
a[6]));
7947 c(3) = (
b[0] * (
b[4] *
b[8] -
b[5] *
b[7]) +
7948 b[1] * (
b[5] *
b[6] -
b[3] *
b[8]) +
7949 b[2] * (
b[3] *
b[7] -
b[4] *
b[6]));
7995 double a = z(2),
b = z(1), c = z(0);
7996 double D =
b*
b-4*
a*c;
8003 x(0) = x(1) = -0.5 *
b /
a;
8008 x(0) = -(x(1) = fabs(0.5 * sqrt(D) /
a));
8016 t = -0.5 * (
b + sqrt(D));
8020 t = -0.5 * (
b - sqrt(D));
8026 Swap<double>(x(0), x(1));
8034 double a = z(2)/z(3),
b = z(1)/z(3), c = z(0)/z(3);
8037 double Q = (
a *
a - 3 *
b) / 9;
8038 double R = (2 *
a *
a *
a - 9 *
a *
b + 27 * c) / 54;
8039 double Q3 = Q * Q * Q;
8046 x(0) = x(1) = x(2) = -
a / 3;
8050 double sqrtQ = sqrt(Q);
8054 x(0) = -2 * sqrtQ -
a / 3;
8055 x(1) = x(2) = sqrtQ -
a / 3;
8059 x(0) = x(1) = - sqrtQ -
a / 3;
8060 x(2) = 2 * sqrtQ -
a / 3;
8067 double theta = acos(R / sqrt(Q3));
8068 double A = -2 * sqrt(Q);
8070 x0 = A * cos(theta / 3) -
a / 3;
8071 x1 = A * cos((theta + 2.0 * M_PI) / 3) -
a / 3;
8072 x2 = A * cos((theta - 2.0 * M_PI) / 3) -
a / 3;
8077 Swap<double>(x0, x1);
8081 Swap<double>(x1, x2);
8084 Swap<double>(x0, x1);
8097 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
8101 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
8103 x(0) = A + Q / A -
a / 3;
8112 const double factor,
const int Dim)
8114 const double c0 = c(0);
8115 c(0) = c0 * (1.0 - pow(factor, -Dim));
8117 for (
int j = 0; j < nr; j++)
8129 c(0) = c0 * (1.0 - pow(factor, Dim));
8131 for (
int j = 0; j < nr; j++)
8150 const double factor = 2.0;
8165 for (
int k = 0; k < nv; k++)
8168 V(j, k) = displacements(v[k]+j*nvs);
8198 for (
int j = 0; j < nv; j++)
8224 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
8227 vertices[i](j) += displacements(j*nv+i);
8235 for (
int i = 0; i < nv; i++)
8238 vert_coord(j*nv+i) =
vertices[i](j);
8244 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
8247 vertices[i](j) = vert_coord(j*nv+i);
8277 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
8294 (*Nodes) += displacements;
8306 node_coord = (*Nodes);
8318 (*Nodes) = node_coord;
8353 mfem::Swap<GridFunction*>(
Nodes, nodes);
8373 for (j = 1; j < n; j++)
8410 int quad_counter = 0;
8425 vertices.SetSize(oelem + quad_counter);
8432 const int attr =
elements[i]->GetAttribute();
8433 int *v =
elements[i]->GetVertices();
8439 for (
int ei = 0; ei < 3; ei++)
8441 for (
int k = 0; k < 2; k++)
8449 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8451 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8453 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8455 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8459 const int qe = quad_counter;
8463 for (
int ei = 0; ei < 4; ei++)
8465 for (
int k = 0; k < 2; k++)
8473 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
8475 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
8477 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
8479 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
8483 MFEM_ABORT(
"unknown element type: " << el_type);
8493 const int attr =
boundary[i]->GetAttribute();
8494 int *v =
boundary[i]->GetVertices();
8503 static const double A = 0.0, B = 0.5, C = 1.0;
8504 static double tri_children[2*3*4] =
8511 static double quad_children[2*4*4] =
8525 for (
int i = 0; i <
elements.Size(); i++)
8546 if (!
Nodes || update_nodes)
8554 static inline double sqr(
const double &x)
8576 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
8579 int NumOfQuadFaces = 0;
8585 for (
int i = 0; i <
faces.Size(); i++)
8589 f2qf[i] = NumOfQuadFaces;
8596 NumOfQuadFaces =
faces.Size();
8600 int hex_counter = 0;
8603 for (
int i = 0; i <
elements.Size(); i++)
8612 int pyr_counter = 0;
8615 for (
int i = 0; i <
elements.Size(); i++)
8633 DSTable *v_to_v_ptr = v_to_v_p;
8649 std::sort(row_start, J_v2v.
end());
8652 for (
int i = 0; i < J_v2v.
Size(); i++)
8654 e2v[J_v2v[i].two] = i;
8667 it.SetIndex(e2v[it.Index()]);
8677 const int oelem = oface + NumOfQuadFaces;
8682 vertices.SetSize(oelem + hex_counter);
8690 const int attr =
elements[i]->GetAttribute();
8691 int *v =
elements[i]->GetVertices();
8698 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8706 for (
int ei = 0; ei < 6; ei++)
8708 for (
int k = 0; k < 2; k++)
8718 const int rt_algo = 1;
8729 double len_sqr, min_len;
8731 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
8732 sqr(J(1,0)-J(1,1)-J(1,2)) +
8733 sqr(J(2,0)-J(2,1)-J(2,2));
8736 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
8737 sqr(J(1,1)-J(1,0)-J(1,2)) +
8738 sqr(J(2,1)-J(2,0)-J(2,2));
8739 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
8741 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
8742 sqr(J(1,2)-J(1,0)-J(1,1)) +
8743 sqr(J(2,2)-J(2,0)-J(2,1));
8744 if (len_sqr < min_len) { rt = 2; }
8749 double Em_data[18], Js_data[9], Jp_data[9];
8752 double ar1, ar2,
kappa, kappa_min;
8754 for (
int s = 0;
s < 3;
s++)
8756 for (
int t = 0;
t < 3;
t++)
8758 Em(
t,
s) = 0.5*J(
t,
s);
8761 for (
int t = 0;
t < 3;
t++)
8763 Em(
t,3) = 0.5*(J(
t,0)+J(
t,1));
8764 Em(
t,4) = 0.5*(J(
t,0)+J(
t,2));
8765 Em(
t,5) = 0.5*(J(
t,1)+J(
t,2));
8769 for (
int t = 0;
t < 3;
t++)
8771 Js(
t,0) = Em(
t,5)-Em(
t,0);
8772 Js(
t,1) = Em(
t,1)-Em(
t,0);
8773 Js(
t,2) = Em(
t,2)-Em(
t,0);
8777 for (
int t = 0;
t < 3;
t++)
8779 Js(
t,0) = Em(
t,5)-Em(
t,0);
8780 Js(
t,1) = Em(
t,2)-Em(
t,0);
8781 Js(
t,2) = Em(
t,4)-Em(
t,0);
8785 kappa_min = std::max(ar1, ar2);
8789 for (
int t = 0;
t < 3;
t++)
8791 Js(
t,0) = Em(
t,0)-Em(
t,1);
8792 Js(
t,1) = Em(
t,4)-Em(
t,1);
8793 Js(
t,2) = Em(
t,2)-Em(
t,1);
8797 for (
int t = 0;
t < 3;
t++)
8799 Js(
t,0) = Em(
t,2)-Em(
t,1);
8800 Js(
t,1) = Em(
t,4)-Em(
t,1);
8801 Js(
t,2) = Em(
t,5)-Em(
t,1);
8805 kappa = std::max(ar1, ar2);
8806 if (
kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
8809 for (
int t = 0;
t < 3;
t++)
8811 Js(
t,0) = Em(
t,0)-Em(
t,2);
8812 Js(
t,1) = Em(
t,1)-Em(
t,2);
8813 Js(
t,2) = Em(
t,3)-Em(
t,2);
8817 for (
int t = 0;
t < 3;
t++)
8819 Js(
t,0) = Em(
t,1)-Em(
t,2);
8820 Js(
t,1) = Em(
t,5)-Em(
t,2);
8821 Js(
t,2) = Em(
t,3)-Em(
t,2);
8825 kappa = std::max(ar1, ar2);
8826 if (
kappa < kappa_min) { rt = 2; }
8829 static const int mv_all[3][4][4] =
8831 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
8832 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
8833 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
8835 const int (&mv)[4][4] = mv_all[rt];
8837 #ifndef MFEM_USE_MEMALLOC 8839 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8841 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8843 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8845 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8847 for (
int k = 0; k < 4; k++)
8849 new_elements[j+4+k] =
8850 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8851 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8855 new_elements[j+0] = tet =
TetMemory.Alloc();
8856 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8858 new_elements[j+1] = tet =
TetMemory.Alloc();
8859 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8861 new_elements[j+2] = tet =
TetMemory.Alloc();
8862 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8864 new_elements[j+3] = tet =
TetMemory.Alloc();
8865 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8867 for (
int k = 0; k < 4; k++)
8869 new_elements[j+4+k] = tet =
TetMemory.Alloc();
8870 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8871 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8874 for (
int k = 0; k < 4; k++)
8879 for (
int k = 0; k < 4; k++)
8893 for (
int fi = 2; fi < 5; fi++)
8895 for (
int k = 0; k < 4; k++)
8902 for (
int ei = 0; ei < 9; ei++)
8904 for (
int k = 0; k < 2; k++)
8911 const int qf2 = f2qf[
f[2]];
8912 const int qf3 = f2qf[
f[3]];
8913 const int qf4 = f2qf[
f[4]];
8916 new Wedge(v[0], oedge+e[0], oedge+e[2],
8917 oedge+e[6], oface+qf2, oface+qf4, attr);
8920 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
8921 oface+qf3, oface+qf4, oface+qf2, attr);
8924 new Wedge(oedge+e[0], v[1], oedge+e[1],
8925 oface+qf2, oedge+e[7], oface+qf3, attr);
8928 new Wedge(oedge+e[2], oedge+e[1], v[2],
8929 oface+qf4, oface+qf3, oedge+e[8], attr);
8932 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
8933 v[3], oedge+e[3], oedge+e[5], attr);
8936 new Wedge(oface+qf3, oface+qf4, oface+qf2,
8937 oedge+e[4], oedge+e[5], oedge+e[3], attr);
8940 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
8941 oedge+e[3], v[4], oedge+e[4], attr);
8944 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
8945 oedge+e[5], oedge+e[4], v[5], attr);
8954 for (
int fi = 0; fi < 1; fi++)
8956 for (
int k = 0; k < 4; k++)
8963 for (
int ei = 0; ei < 8; ei++)
8965 for (
int k = 0; k < 2; k++)
8972 const int qf0 = f2qf[
f[0]];
8975 new Pyramid(v[0], oedge+e[0], oface+qf0,
8976 oedge+e[3], oedge+e[4], attr);
8979 new Pyramid(oedge+e[0], v[1], oedge+e[1],
8980 oface+qf0, oedge+e[5], attr);
8983 new Pyramid(oface+qf0, oedge+e[1], v[2],
8984 oedge+e[2], oedge+e[6], attr);
8987 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
8988 v[3], oedge+e[7], attr);
8991 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
8992 oedge+e[7], v[4], attr);
8995 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
8996 oedge+e[4], oface+qf0, attr);
8998 #ifndef MFEM_USE_MEMALLOC 9000 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
9004 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
9008 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
9012 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
9016 new_elements[j++] = tet =
TetMemory.Alloc();
9017 tet->
Init(oedge+e[0], oedge+e[4], oedge+e[5],
9020 new_elements[j++] = tet =
TetMemory.Alloc();
9021 tet->
Init(oedge+e[1], oedge+e[5], oedge+e[6],
9024 new_elements[j++] = tet =
TetMemory.Alloc();
9025 tet->
Init(oedge+e[2], oedge+e[6], oedge+e[7],
9028 new_elements[j++] = tet =
TetMemory.Alloc();
9029 tet->
Init(oedge+e[3], oedge+e[7], oedge+e[4],
9038 const int he = hex_counter;
9043 if (f2qf.
Size() == 0)
9049 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[
f[k]]; }
9055 for (
int fi = 0; fi < 6; fi++)
9057 for (
int k = 0; k < 4; k++)
9064 for (
int ei = 0; ei < 12; ei++)
9066 for (
int k = 0; k < 2; k++)
9074 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
9075 oedge+e[3], oedge+e[8], oface+qf[1],
9076 oelem+he, oface+qf[4], attr);
9079 oface+qf[0], oface+qf[1], oedge+e[9],
9080 oface+qf[2], oelem+he, attr);
9082 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
9083 oedge+e[2], oelem+he, oface+qf[2],
9084 oedge+e[10], oface+qf[3], attr);
9086 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
9087 v[3], oface+qf[4], oelem+he,
9088 oface+qf[3], oedge+e[11], attr);
9090 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
9091 oface+qf[4], v[4], oedge+e[4],
9092 oface+qf[5], oedge+e[7], attr);
9094 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
9095 oelem+he, oedge+e[4], v[5],
9096 oedge+e[5], oface+qf[5], attr);
9098 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
9099 oface+qf[3], oface+qf[5], oedge+e[5],
9100 v[6], oedge+e[6], attr);
9102 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
9103 oedge+e[11], oedge+e[7], oface+qf[5],
9104 oedge+e[6], v[7], attr);
9109 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
9121 const int attr =
boundary[i]->GetAttribute();
9122 int *v =
boundary[i]->GetVertices();
9129 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
9136 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
9138 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
9140 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
9142 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
9150 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
9152 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
9154 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
9156 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
9160 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
9166 static const double A = 0.0, B = 0.5, C = 1.0, D = -1.0;
9167 static double tet_children[3*4*16] =
9169 A,A,A, B,A,A, A,B,A, A,A,B,
9170 B,A,A, C,A,A, B,B,A, B,A,B,
9171 A,B,A, B,B,A, A,C,A, A,B,B,
9172 A,A,B, B,A,B, A,B,B, A,A,C,
9177 B,A,A, A,B,B, A,B,A, A,A,B,
9178 B,A,A, A,B,B, A,A,B, B,A,B,
9179 B,A,A, A,B,B, B,A,B, B,B,A,
9180 B,A,A, A,B,B, B,B,A, A,B,A,
9182 A,B,A, B,A,A, B,A,B, A,A,B,
9183 A,B,A, A,A,B, B,A,B, A,B,B,
9184 A,B,A, A,B,B, B,A,B, B,B,A,
9185 A,B,A, B,B,A, B,A,B, B,A,A,
9187 A,A,B, B,A,A, A,B,A, B,B,A,
9188 A,A,B, A,B,A, A,B,B, B,B,A,
9189 A,A,B, A,B,B, B,A,B, B,B,A,
9190 A,A,B, B,A,B, B,A,A, B,B,A
9192 static double pyr_children[3*5*10] =
9194 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
9195 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
9196 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
9197 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
9198 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
9199 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
9200 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
9201 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
9202 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
9203 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
9205 static double pri_children[3*6*8] =
9207 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
9208 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
9209 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
9210 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
9211 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
9212 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
9213 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
9214 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
9216 static double hex_children[3*8*8] =
9218 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
9219 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
9220 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
9221 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
9222 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
9223 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
9224 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
9225 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
9237 for (
int i = 0; i <
elements.Size(); i++)
9268 int i, j, ind, nedges;
9275 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
9289 for (j = 0; j < marked_el.
Size(); j++)
9294 int new_v = cnv + j, new_e = cne + j;
9303 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
9305 UseExternalData(seg_children, 1, 2, 3);
9318 int *edge1 =
new int[nedges];
9319 int *edge2 =
new int[nedges];
9320 int *middle =
new int[nedges];
9322 for (i = 0; i < nedges; i++)
9324 edge1[i] = edge2[i] = middle[i] = -1;
9330 for (j = 1; j < v.
Size(); j++)
9332 ind = v_to_v(v[j-1], v[j]);
9333 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
9335 ind = v_to_v(v[0], v[v.
Size()-1]);
9336 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
9340 for (i = 0; i < marked_el.
Size(); i++)
9346 int need_refinement;
9349 need_refinement = 0;
9350 for (i = 0; i < nedges; i++)
9352 if (middle[i] != -1 && edge1[i] != -1)
9354 need_refinement = 1;
9359 while (need_refinement == 1);
9362 int v1[2], v2[2], bisect, temp;
9364 for (i = 0; i < temp; i++)
9367 bisect = v_to_v(v[0], v[1]);
9368 if (middle[bisect] != -1)
9372 v1[0] = v[0]; v1[1] = middle[bisect];
9373 v2[0] = middle[bisect]; v2[1] = v[1];
9379 mfem_error(
"Only bisection of segment is implemented" 9403 MFEM_VERIFY(
GetNE() == 0 ||
9405 "tetrahedral mesh is not marked for refinement:" 9406 " call Finalize(true)");
9413 for (i = 0; i < marked_el.
Size(); i++)
9419 for (i = 0; i < marked_el.
Size(); i++)
9428 for (i = 0; i < marked_el.
Size(); i++)
9445 int need_refinement;
9450 need_refinement = 0;
9458 if (
elements[i]->NeedRefinement(v_to_v))
9460 need_refinement = 1;
9465 while (need_refinement == 1);
9472 need_refinement = 0;
9474 if (
boundary[i]->NeedRefinement(v_to_v))
9476 need_refinement = 1;
9480 while (need_refinement == 1);
9511 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is " 9512 "not supported. Project the NURBS to Nodes first.");
9522 if (!refinements.
Size())
9543 Swap(*mesh2,
false);
9559 const int *fine,
int nfine,
int op)
9561 double error = elem_error[fine[0]];
9563 for (
int i = 1; i < nfine; i++)
9565 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
9567 double err_fine = elem_error[fine[i]];
9570 case 0: error = std::min(error, err_fine);
break;
9571 case 1: error += err_fine;
break;
9572 case 2: error = std::max(error, err_fine);
break;
9579 double threshold,
int nc_limit,
int op)
9581 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
9582 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. " 9583 "Project the NURBS to Nodes first.");
9596 for (
int i = 0; i < dt.
Size(); i++)
9598 if (nc_limit > 0 && !level_ok[i]) {
continue; }
9603 if (error < threshold) { derefs.
Append(i); }
9606 if (!derefs.
Size()) {
return false; }
9613 Swap(*mesh2,
false);
9627 int nc_limit,
int op)
9637 MFEM_ABORT(
"Derefinement is currently supported for non-conforming " 9644 int nc_limit,
int op)
9647 for (
int i = 0; i < tmp.Size(); i++)
9649 tmp[i] = elem_error(i);
9735 #ifdef MFEM_USE_MEMALLOC 9762 for (
int i = 0; i < elem_array.
Size(); i++)
9764 if (elem_array[i]->GetGeometryType() == geom)
9769 elem_vtx.
SetSize(nv*num_elems);
9773 for (
int i = 0; i < elem_array.
Size(); i++)
9779 elem_vtx.
Append(loc_vtx);
9787 for (
int i = 0; i < nelem; i++) { list[i] = i; }
9803 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
9815 default: MFEM_ABORT(
"internal error");
9821 int nonconforming,
int nc_limit)
9831 else if (nonconforming < 0)
9852 for (
int i = 0; i < refinements.
Size(); i++)
9854 el_to_refine[i] = refinements[i].index;
9858 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
9859 if (rt == 1 || rt == 2 || rt == 4)
9863 else if (rt == 3 || rt == 5 || rt == 6)
9881 for (
int i = 0; i < el_to_refine.
Size(); i++)
9883 refinements[i] =
Refinement(el_to_refine[i]);
9890 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. " 9891 "Please project the NURBS to Nodes first, with SetCurvature().");
9894 MFEM_VERIFY(
ncmesh != NULL || dynamic_cast<const ParMesh*>(
this) == NULL,
9895 "Sorry, converting a conforming ParMesh to an NC mesh is " 9903 (simplices_nonconforming && (
meshgen & 0x1)) )
9916 for (
int i = 0; i <
GetNE(); i++)
9918 if ((
double) rand() / RAND_MAX <
prob)
9923 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
9935 for (
int i = 0; i <
GetNE(); i++)
9938 bool refine =
false;
9939 for (
int j = 0; j < v.
Size(); j++)
9944 double d = vert(l) -
vertices[v[j]](l);
9947 if (dist <= eps*eps) { refine =
true;
break; }
9958 int nonconforming,
int nc_limit)
9960 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
9962 for (
int i = 0; i <
GetNE(); i++)
9964 if (elem_error[i] > threshold)
9978 int nonconforming,
int nc_limit)
9982 return RefineByError(tmp, threshold, nonconforming, nc_limit);
9987 int *edge1,
int *edge2,
int *middle)
9990 int v[2][4], v_new, bisect,
t;
10002 bisect = v_to_v(vert[0], vert[1]);
10003 MFEM_ASSERT(bisect >= 0,
"");
10005 if (middle[bisect] == -1)
10008 for (
int d = 0; d <
spaceDim; d++)
10016 if (edge1[bisect] == i)
10018 edge1[bisect] = edge2[bisect];
10021 middle[bisect] = v_new;
10025 v_new = middle[bisect];
10028 edge1[bisect] = -1;
10033 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
10034 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
10055 bisect = v_to_v(v[1][0], v[1][1]);
10056 MFEM_ASSERT(bisect >= 0,
"");
10058 if (edge1[bisect] == i)
10062 else if (edge2[bisect] == i)
10071 MFEM_ABORT(
"Bisection for now works only for triangles.");
10078 int v[2][4], v_new, bisect,
t;
10088 "TETRAHEDRON element is not marked for refinement.");
10093 bisect = v_to_v.
FindId(vert[0], vert[1]);
10097 for (
int j = 0; j < 3; j++)
10110 int type, old_redges[2], flag;
10113 int new_type, new_redges[2][2];
10116 new_redges[0][0] = 2;
10117 new_redges[0][1] = 1;
10118 new_redges[1][0] = 2;
10119 new_redges[1][1] = 1;
10120 int tr1 = -1, tr2 = -1;
10121 switch (old_redges[0])
10124 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
10129 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
10133 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
10136 switch (old_redges[1])
10139 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
10144 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
10148 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
10155 #ifdef MFEM_USE_MEMALLOC 10191 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
10198 int v[2][3], v_new, bisect,
t;
10209 bisect = v_to_v.
FindId(vert[0], vert[1]);
10210 MFEM_ASSERT(bisect >= 0,
"");
10212 MFEM_ASSERT(v_new != -1,
"");
10216 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
10217 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
10227 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for" 10233 int *edge1,
int *edge2,
int *middle)
10236 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
10245 bisect[0] = v_to_v(v[0],v[1]);
10246 bisect[1] = v_to_v(v[1],v[2]);
10247 bisect[2] = v_to_v(v[0],v[2]);
10248 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
10250 for (j = 0; j < 3; j++)
10252 if (middle[bisect[j]] == -1)
10255 for (
int d = 0; d <
spaceDim; d++)
10263 if (edge1[bisect[j]] == i)
10265 edge1[bisect[j]] = edge2[bisect[j]];
10268 middle[bisect[j]] = v_new[j];
10272 v_new[j] = middle[bisect[j]];
10275 edge1[bisect[j]] = -1;
10281 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
10282 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
10283 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
10284 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
10299 tri2->ResetTransform(code);
10300 tri3->ResetTransform(code);
10304 tri2->PushTransform(1);
10305 tri3->PushTransform(2);
10318 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
10354 for (
int i = 0; i < elem_geoms.
Size(); i++)
10362 std::map<unsigned, int> mat_no;
10366 for (
int j = 0; j <
elements.Size(); j++)
10369 unsigned code =
elements[j]->GetTransform();
10372 int &matrix = mat_no[code];
10373 if (!matrix) { matrix = mat_no.size(); }
10383 std::map<unsigned, int>::iterator it;
10384 for (it = mat_no.begin(); it != mat_no.end(); ++it)
10398 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for" 10399 " geom = " << geom);
10409 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
10418 os <<
"areamesh2\n\n";
10422 os <<
"curved_areamesh2\n\n";
10431 os <<
boundary[i]->GetAttribute();
10432 for (j = 0; j < v.
Size(); j++)
10434 os <<
' ' << v[j] + 1;
10446 for (j = 0; j < v.
Size(); j++)
10448 os <<
' ' << v[j] + 1;
10460 for (j = 1; j <
Dim; j++)
10477 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
10485 os <<
"NETGEN_Neutral_Format\n";
10490 for (j = 0; j <
Dim; j++)
10503 os <<
elements[i]->GetAttribute();
10504 for (j = 0; j < nv; j++)
10506 os <<
' ' << ind[j]+1;
10517 os <<
boundary[i]->GetAttribute();
10518 for (j = 0; j < nv; j++)
10520 os <<
' ' << ind[j]+1;
10532 <<
" 0 0 0 0 0 0 0\n" 10533 <<
"0 0 0 1 0 0 0 0 0 0 0\n" 10535 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n" 10536 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
10540 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
10546 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
10547 for (j = 0; j < nv; j++)
10549 os <<
' ' << ind[j]+1;
10558 os <<
boundary[i]->GetAttribute();
10559 for (j = 0; j < nv; j++)
10561 os <<
' ' << ind[j]+1;
10563 os <<
" 1.0 1.0 1.0 1.0\n";
10596 os <<
"\n# mesh curvature GridFunction";
10601 os <<
"\nmfem_mesh_end" << endl;
10606 os << (section_delimiter.empty()
10607 ?
"MFEM mesh v1.0\n" :
"MFEM mesh v1.2\n");
10611 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 10616 "# TETRAHEDRON = 4\n" 10622 os <<
"\ndimension\n" <<
Dim;
10657 if (!section_delimiter.empty())
10659 os << section_delimiter << endl;
10668 os <<
"MFEM NURBS mesh v1.0\n";
10672 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 10678 os <<
"\ndimension\n" <<
Dim 10695 int ki = e_to_k[i];
10700 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
10707 ofstream ofs(fname);
10708 ofs.precision(precision);
10712 #ifdef MFEM_USE_ADIOS2 10722 "# vtk DataFile Version 3.0\n" 10723 "Generated by MFEM\n" 10725 "DATASET UNSTRUCTURED_GRID\n";
10738 for ( ; j < 3; j++)
10754 os << (*Nodes)(vdofs[0]);
10758 os <<
' ' << (*Nodes)(vdofs[j]);
10760 for ( ; j < 3; j++)
10774 size +=
elements[i]->GetNVertices() + 1;
10779 const int *v =
elements[i]->GetVertices();
10780 const int nv =
elements[i]->GetNVertices();
10784 for (
int j = 0; j < nv; j++)
10786 os <<
' ' << v[perm ? perm[j] : j];
10799 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
10800 "Point meshes should have a single dof per element");
10801 size += dofs.
Size() + 1;
10806 if (!strcmp(fec_name,
"Linear") ||
10807 !strcmp(fec_name,
"H1_0D_P1") ||
10808 !strcmp(fec_name,
"H1_1D_P1") ||
10809 !strcmp(fec_name,
"H1_2D_P1") ||
10810 !strcmp(fec_name,
"H1_3D_P1"))
10814 else if (!strcmp(fec_name,
"Quadratic") ||
10815 !strcmp(fec_name,
"H1_1D_P2") ||
10816 !strcmp(fec_name,
"H1_2D_P2") ||
10817 !strcmp(fec_name,
"H1_3D_P2"))
10823 mfem::err <<
"Mesh::PrintVTK : can not save '" 10824 << fec_name <<
"' elements!" << endl;
10833 for (
int j = 0; j < dofs.
Size(); j++)
10835 os <<
' ' << dofs[j];
10838 else if (order == 2)
10840 const int *vtk_mfem;
10841 switch (
elements[i]->GetGeometryType())
10855 for (
int j = 0; j < dofs.
Size(); j++)
10857 os <<
' ' << dofs[vtk_mfem[j]];
10867 int vtk_cell_type = 5;
10871 os << vtk_cell_type <<
'\n';
10876 <<
"SCALARS material int\n" 10877 <<
"LOOKUP_TABLE default\n";
10880 os <<
elements[i]->GetAttribute() <<
'\n';
10887 bool high_order_output,
10888 int compression_level,
10891 int ref = (high_order_output &&
Nodes)
10894 fname = fname +
".vtu";
10896 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
10897 if (compression_level != 0)
10899 os <<
" compressor=\"vtkZLibDataCompressor\"";
10902 os <<
"<UnstructuredGrid>\n";
10903 PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
10904 os <<
"</Piece>\n";
10905 os <<
"</UnstructuredGrid>\n";
10906 os <<
"</VTKFile>" << std::endl;
10913 bool high_order_output,
10914 int compression_level)
10916 PrintVTU(fname, format, high_order_output, compression_level,
true);
10920 bool high_order_output,
int compression_level,
10928 std::vector<char> buf;
10930 auto get_geom = [&](
int i)
10938 int np = 0, nc_ref = 0;
10939 for (
int i = 0; i < ne; i++)
10948 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\"" 10949 << (high_order_output ? ne : nc_ref) <<
"\">\n";
10952 os <<
"<Points>\n";
10953 os <<
"<DataArray type=\"" << type_str
10954 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
10955 for (
int i = 0; i < ne; i++)
10968 for (
int j = 0; j < pmat.
Width(); j++)
10994 os <<
"</DataArray>" << std::endl;
10995 os <<
"</Points>" << std::endl;
10997 os <<
"<Cells>" << std::endl;
10998 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"" 10999 << fmt_str <<
"\">" << std::endl;
11001 std::vector<int> offset;
11004 if (high_order_output)
11007 for (
int iel = 0; iel < ne; iel++)
11011 int nnodes = local_connectivity.
Size();
11012 for (
int i=0; i<nnodes; ++i)
11019 offset.push_back(np);
11025 for (
int i = 0; i < ne; i++)
11031 for (
int j = 0; j < RG.
Size(); )
11034 offset.push_back(coff);
11036 for (
int k = 0; k < nv; k++, j++)
11050 os <<
"</DataArray>" << std::endl;
11052 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\"" 11053 << fmt_str <<
"\">" << std::endl;
11055 for (
size_t ii=0; ii<offset.size(); ii++)
11063 os <<
"</DataArray>" << std::endl;
11064 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\"" 11065 << fmt_str <<
"\">" << std::endl;
11067 const int *vtk_geom_map =
11069 for (
int i = 0; i < ne; i++)
11072 uint8_t vtk_cell_type = 5;
11074 vtk_cell_type = vtk_geom_map[geom];
11076 if (high_order_output)
11085 for (
int j = 0; j < RG.
Size(); j += nv)
11095 os <<
"</DataArray>" << std::endl;
11096 os <<
"</Cells>" << std::endl;
11098 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
11099 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\"" 11100 << fmt_str <<
"\">" << std::endl;
11101 for (
int i = 0; i < ne; i++)
11104 if (high_order_output)
11123 os <<
"</DataArray>" << std::endl;
11124 os <<
"</CellData>" << std::endl;
11135 "# vtk DataFile Version 3.0\n" 11136 "Generated by MFEM\n" 11138 "DATASET UNSTRUCTURED_GRID\n";
11143 os <<
"FIELD FieldData 1\n" 11153 np = nc = size = 0;
11154 for (
int i = 0; i <
GetNE(); i++)
11163 os <<
"POINTS " << np <<
" double\n";
11165 for (
int i = 0; i <
GetNE(); i++)
11172 for (
int j = 0; j < pmat.
Width(); j++)
11174 os << pmat(0, j) <<
' ';
11177 os << pmat(1, j) <<
' ';
11189 os << 0.0 <<
' ' << 0.0;
11196 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
11198 for (
int i = 0; i <
GetNE(); i++)
11205 for (
int j = 0; j < RG.
Size(); )
11208 for (
int k = 0; k < nv; k++, j++)
11210 os <<
' ' << np + RG[j];
11216 os <<
"CELL_TYPES " << nc <<
'\n';
11217 for (
int i = 0; i <
GetNE(); i++)
11225 for (
int j = 0; j < RG.
Size(); j += nv)
11227 os << vtk_cell_type <<
'\n';
11231 os <<
"CELL_DATA " << nc <<
'\n' 11232 <<
"SCALARS material int\n" 11233 <<
"LOOKUP_TABLE default\n";
11234 for (
int i = 0; i <
GetNE(); i++)
11242 os << attr <<
'\n';
11249 srand((
unsigned)time(0));
11250 double a = double(rand()) / (double(RAND_MAX) + 1.);
11251 int el0 = (int)floor(
a *
GetNE());
11253 os <<
"SCALARS element_coloring int\n" 11254 <<
"LOOKUP_TABLE default\n";
11255 for (
int i = 0; i <
GetNE(); i++)
11262 os << coloring[i] + 1 <<
'\n';
11268 os <<
"POINT_DATA " << np <<
'\n' << flush;
11273 int delete_el_to_el = (
el_to_el) ? (0) : (1);
11275 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
11278 const int *i_el_el = el_el.
GetI();
11279 const int *j_el_el = el_el.
GetJ();
11284 stack_p = stack_top_p = 0;
11285 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
11287 if (colors[el] != -2)
11293 el_stack[stack_top_p++] = el;
11295 for ( ; stack_p < stack_top_p; stack_p++)
11297 int i = el_stack[stack_p];
11298 int num_nb = i_el_el[i+1] - i_el_el[i];
11299 if (max_num_col < num_nb + 1)
11301 max_num_col = num_nb + 1;
11303 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
11305 int k = j_el_el[j];
11306 if (colors[k] == -2)
11309 el_stack[stack_top_p++] = k;
11317 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
11319 int i = el_stack[stack_p], col;
11321 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
11323 col = colors[j_el_el[j]];
11326 col_marker[col] = 1;
11330 for (col = 0; col < max_num_col; col++)
11331 if (col_marker[col] == 0)
11339 if (delete_el_to_el)
11347 int elem_attr)
const 11349 if (
Dim != 3 &&
Dim != 2) {
return; }
11351 int i, j, k, l, nv, nbe, *v;
11353 os <<
"MFEM mesh v1.0\n";
11357 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 11362 "# TETRAHEDRON = 4\n" 11367 os <<
"\ndimension\n" <<
Dim 11372 <<
' ' <<
elements[i]->GetGeometryType();
11375 for (j = 0; j < nv; j++)
11387 l = partitioning[l];
11402 os <<
"\nboundary\n" << nbe <<
'\n';
11408 l = partitioning[l];
11411 nv =
faces[i]->GetNVertices();
11412 v =
faces[i]->GetVertices();
11413 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
11414 for (j = 0; j < nv; j++)
11421 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
11422 for (j = nv-1; j >= 0; j--)
11433 nv =
faces[i]->GetNVertices();
11434 v =
faces[i]->GetVertices();
11435 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
11436 for (j = 0; j < nv; j++)
11467 int interior_faces)
11469 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
11470 if (
Dim != 3 &&
Dim != 2) {
return; }
11479 int nv =
elements[i]->GetNVertices();
11480 const int *ind =
elements[i]->GetVertices();
11481 for (
int j = 0; j < nv; j++)
11491 voff[i] = vcount[i-1] + voff[i-1];
11497 vown[i] =
new int[vcount[i]];
11509 int nv =
elements[i]->GetNVertices();
11510 const int *ind =
elements[i]->GetVertices();
11511 for (
int j = 0; j < nv; j++)
11514 vown[ind[j]][vcount[ind[j]]] = i;
11520 vcount[i] = voff[i+1] - voff[i];
11524 for (
int i = 0; i < edge_el.
Size(); i++)
11526 const int *el = edge_el.
GetRow(i);
11529 int k = partitioning[el[0]];
11530 int l = partitioning[el[1]];
11531 if (interior_faces || k != l)
11543 os <<
"areamesh2\n\n" << nbe <<
'\n';
11545 for (
int i = 0; i < edge_el.
Size(); i++)
11547 const int *el = edge_el.
GetRow(i);
11550 int k = partitioning[el[0]];
11551 int l = partitioning[el[1]];
11552 if (interior_faces || k != l)
11557 for (
int j = 0; j < 2; j++)
11558 for (
int s = 0;
s < vcount[ev[j]];
s++)
11559 if (vown[ev[j]][
s] == el[0])
11561 os <<
' ' << voff[ev[j]]+
s+1;
11565 for (
int j = 1; j >= 0; j--)
11566 for (
int s = 0;
s < vcount[ev[j]];
s++)
11567 if (vown[ev[j]][
s] == el[1])
11569 os <<
' ' << voff[ev[j]]+
s+1;
11576 int k = partitioning[el[0]];
11580 for (
int j = 0; j < 2; j++)
11581 for (
int s = 0;
s < vcount[ev[j]];
s++)
11582 if (vown[ev[j]][
s] == el[0])
11584 os <<
' ' << voff[ev[j]]+
s+1;
11594 int nv =
elements[i]->GetNVertices();
11595 const int *ind =
elements[i]->GetVertices();
11596 os << partitioning[i]+1 <<
' ';
11598 for (
int j = 0; j < nv; j++)
11600 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11601 vown[ind[j]][vcount[ind[j]]] = i;
11608 vcount[i] = voff[i+1] - voff[i];
11614 for (
int k = 0; k < vcount[i]; k++)
11616 for (
int j = 0; j <
Dim; j++)
11626 os <<
"NETGEN_Neutral_Format\n";
11630 for (
int k = 0; k < vcount[i]; k++)
11632 for (
int j = 0; j <
Dim; j++)
11643 int nv =
elements[i]->GetNVertices();
11644 const int *ind =
elements[i]->GetVertices();
11645 os << partitioning[i]+1;
11646 for (
int j = 0; j < nv; j++)
11648 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11649 vown[ind[j]][vcount[ind[j]]] = i;
11656 vcount[i] = voff[i+1] - voff[i];
11666 int k = partitioning[
faces_info[i].Elem1No];
11667 l = partitioning[l];
11668 if (interior_faces || k != l)
11685 int k = partitioning[
faces_info[i].Elem1No];
11686 l = partitioning[l];
11687 if (interior_faces || k != l)
11689 int nv =
faces[i]->GetNVertices();
11690 const int *ind =
faces[i]->GetVertices();
11692 for (
int j = 0; j < nv; j++)
11693 for (
int s = 0;
s < vcount[ind[j]];
s++)
11696 os <<
' ' << voff[ind[j]]+
s+1;
11700 for (
int j = nv-1; j >= 0; j--)
11701 for (
int s = 0;
s < vcount[ind[j]];
s++)
11704 os <<
' ' << voff[ind[j]]+
s+1;
11711 int k = partitioning[
faces_info[i].Elem1No];
11712 int nv =
faces[i]->GetNVertices();
11713 const int *ind =
faces[i]->GetVertices();
11715 for (
int j = 0; j < nv; j++)
11716 for (
int s = 0;
s < vcount[ind[j]];
s++)
11719 os <<
' ' << voff[ind[j]]+
s+1;
11735 int k = partitioning[
faces_info[i].Elem1No];
11736 l = partitioning[l];
11737 if (interior_faces || k != l)
11750 <<
" 0 0 0 0 0 0 0\n" 11751 <<
"0 0 0 1 0 0 0 0 0 0 0\n" 11752 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n" 11753 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n" 11754 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11757 for (
int k = 0; k < vcount[i]; k++)
11758 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' ' 11763 int nv =
elements[i]->GetNVertices();
11764 const int *ind =
elements[i]->GetVertices();
11765 os << i+1 <<
' ' << partitioning[i]+1;
11766 for (
int j = 0; j < nv; j++)
11768 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11769 vown[ind[j]][vcount[ind[j]]] = i;
11776 vcount[i] = voff[i+1] - voff[i];
11785 int k = partitioning[
faces_info[i].Elem1No];
11786 l = partitioning[l];
11787 if (interior_faces || k != l)
11789 int nv =
faces[i]->GetNVertices();
11790 const int *ind =
faces[i]->GetVertices();
11792 for (
int j = 0; j < nv; j++)
11793 for (
int s = 0;
s < vcount[ind[j]];
s++)
11796 os <<
' ' << voff[ind[j]]+
s+1;
11798 os <<
" 1.0 1.0 1.0 1.0\n";
11800 for (
int j = nv-1; j >= 0; j--)
11801 for (
int s = 0;
s < vcount[ind[j]];
s++)
11804 os <<
' ' << voff[ind[j]]+
s+1;
11806 os <<
" 1.0 1.0 1.0 1.0\n";
11811 int k = partitioning[
faces_info[i].Elem1No];
11812 int nv =
faces[i]->GetNVertices();
11813 const int *ind =
faces[i]->GetVertices();
11815 for (
int j = 0; j < nv; j++)
11816 for (
int s = 0;
s < vcount[ind[j]];
s++)
11819 os <<
' ' << voff[ind[j]]+
s+1;
11821 os <<
" 1.0 1.0 1.0 1.0\n";
11845 " NURBS mesh is not supported!");
11849 os <<
"MFEM mesh v1.0\n";
11853 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n" 11858 "# TETRAHEDRON = 4\n" 11863 os <<
"\ndimension\n" <<
Dim 11871 const int *
const i_AF_f = Aface_face.
GetI();
11872 const int *
const j_AF_f = Aface_face.
GetJ();
11874 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
11875 for (
const int * iface = j_AF_f + i_AF_f[iAF];
11876 iface < j_AF_f + i_AF_f[iAF+1];
11879 os << iAF+1 <<
' ';
11911 double *cg =
new double[na*
spaceDim];
11912 int *nbea =
new int[na];
11919 for (i = 0; i < na; i++)
11931 for (k = 0; k < vert.
Size(); k++)
11943 for (k = 0; k < vert.
Size(); k++)
11944 if (vn[vert[k]] == 1)
11949 cg[bea*
spaceDim+j] += pointmat(j,k);
11960 for (k = 0; k < vert.
Size(); k++)
11965 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
11981 double *cg =
new double[na*
spaceDim];
11982 int *nbea =
new int[na];
11989 for (i = 0; i < na; i++)
12001 for (k = 0; k < vert.
Size(); k++)
12013 for (k = 0; k < vert.
Size(); k++)
12014 if (vn[vert[k]] == 1)
12019 cg[bea*
spaceDim+j] += pointmat(j,k);
12030 for (k = 0; k < vert.
Size(); k++)
12035 (1-sf)*cg[bea*
spaceDim+j]/nbea[bea];
12051 for (
int i = 0; i <
vertices.Size(); i++)
12053 for (
int j = 0; j <
spaceDim; j++)
12065 xnew.ProjectCoefficient(f_pert);
12074 "incompatible vector dimensions");
12082 for (
int d = 0; d <
spaceDim; d++)
12102 for (
int i = 0; i <
GetNE(); i++)
12107 for (
int j = 0; j < nv; j++)
12112 for (
int i = 0; i <
GetNBE(); i++)
12117 for (
int j = 0; j < nv; j++)
12123 for (
int i = 0; i < v2v.
Size(); i++)
12128 v2v[i] = num_vert++;
12132 if (num_vert == v2v.
Size()) {
return; }
12134 Vector nodes_by_element;
12139 for (
int i = 0; i <
GetNE(); i++)
12146 for (
int i = 0; i <
GetNE(); i++)
12155 for (
int i = 0; i <
GetNE(); i++)
12160 for (
int j = 0; j < nv; j++)
12165 for (
int i = 0; i <
GetNBE(); i++)
12170 for (
int j = 0; j < nv; j++)
12194 for (
int i = 0; i <
GetNE(); i++)
12207 int num_bdr_elem = 0;
12208 int new_bel_to_edge_nnz = 0;
12209 for (
int i = 0; i <
GetNBE(); i++)
12225 if (num_bdr_elem ==
GetNBE()) {
return; }
12229 Table *new_bel_to_edge = NULL;
12233 new_be_to_edge.
Reserve(num_bdr_elem);
12237 new_be_to_face.
Reserve(num_bdr_elem);
12238 new_bel_to_edge =
new Table;
12239 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
12241 for (
int i = 0; i <
GetNBE(); i++)
12252 int row = new_be_to_face.
Size();
12256 int *new_e = new_bel_to_edge->
GetRow(row);
12257 for (
int j = 0; j < ne; j++)
12261 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
12281 for (
int i = 0; i < attribs.
Size(); i++)
12293 #ifdef MFEM_USE_MEMALLOC 12320 const int npts = point_mat.
Width();
12321 if (!npts) {
return 0; }
12322 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
12326 if (!
GetNE()) {
return 0; }
12328 double *data = point_mat.
GetData();
12335 min_dist = std::numeric_limits<double>::max();
12339 for (
int i = 0; i <
GetNE(); i++)
12343 for (
int k = 0; k < npts; k++)
12346 if (dist < min_dist(k))
12348 min_dist(k) = dist;
12357 for (
int k = 0; k < npts; k++)
12361 int res = inv_tr->
Transform(pt, ips[k]);
12364 elem_ids[k] = e_idx[k];
12368 if (pts_found != npts)
12372 for (
int k = 0; k < npts; k++)
12374 if (elem_ids[k] != -1) {
continue; }
12378 for (
int v = 0; v < elvertices.
Size(); v++)
12380 int vv = elvertices[v];
12382 const int* els = vtoel->
GetRow(vv);
12383 for (
int e = 0; e < ne; e++)
12385 if (els[e] == e_idx[k]) {
continue; }
12387 int res = inv_tr->
Transform(pt, ips[k]);
12390 elem_ids[k] = els[e];
12402 for (
int e = 0; e < neigh.
Size(); e++)
12408 int res = inv_tr->
Transform(pt, ips[k]);
12421 if (inv_trans == NULL) {
delete inv_tr; }
12423 if (warn && pts_found != npts)
12425 MFEM_WARNING((npts-pts_found) <<
" points were not found");
12440 MFEM_VERIFY(
Dim == 2 ||
Dim == 3,
"Only 2D/3D meshes supported right now.");
12441 MFEM_VERIFY(
Dim ==
spaceDim,
"Surface meshes not currently supported.");
12458 skew(0) = std::atan2(J.
Det(), col1 * col2);
12461 ori(0) = std::atan2(J(1,0), J(0,0));
12468 Vector col1, col2, col3;
12472 double len1 = col1.
Norml2(),
12479 col1unit *= 1.0/len1;
12480 col2unit *= 1.0/len2;
12481 col3unit *= 1.0/len3;
12487 aspr(0) = len1/std::sqrt(len2*len3),
12488 aspr(1) = len2/std::sqrt(len1*len3);
12491 aspr(2) = std::sqrt(len1/(len2*len3)),
12492 aspr(3) = std::sqrt(len2/(len1*len3));
12495 Vector crosscol12, crosscol13;
12496 col1.
cross3D(col2, crosscol12);
12497 col1.
cross3D(col3, crosscol13);
12498 skew(0) = std::acos(col1unit*col2unit);
12499 skew(1) = std::acos(col1unit*col3unit);
12500 skew(2) = std::atan(len1*volume/(crosscol12*crosscol13));
12506 for (
int d=0; d<
Dim; d++) { rot(d, 0) = col1unit(d); }
12510 rot1 *= col1unit*col2unit;
12512 col1unit.
cross3D(col2unit, rot1);
12514 for (
int d=0; d <
Dim; d++) { rot(d, 1) = rot2(d); }
12517 for (
int d=0; d <
Dim; d++) { rot(d, 2) = rot1(d); }
12518 double delta = sqrt(pow(rot(2,1)-rot(1,2), 2.0) +
12519 pow(rot(0,2)-rot(2,0), 2.0) +
12520 pow(rot(1,0)-rot(0,1), 2.0));
12525 for (
int d = 0; d <
Dim; d++) { Iden(d, d) = 1.0; };
12531 MFEM_ABORT(
"Invalid rotation matrix. Contact TMOP Developers.");
12536 ori(0) = (1./
delta)*(rot(2,1)-rot(1,2));
12537 ori(1) = (1./
delta)*(rot(0,2)-rot(2,0));
12538 ori(2) = (1./
delta)*(rot(1,0)-rot(0,1));
12539 ori(3) = std::acos(0.5*(rot.
Trace()-1.0));
12553 "mixed meshes are not supported!");
12554 MFEM_ASSERT(
mesh->
GetNodes(),
"meshes without nodes are not supported!");
12567 Compute(nodes, d_mt);
12570 void GeometricFactors::Compute(
const GridFunction &nodes,
12577 const int vdim = fespace->
GetVDim();
12578 const int NE = fespace->
GetNE();
12579 const int ND = fe->
GetDof();
12582 unsigned eval_flags = 0;
12607 qi->DisableTensorProducts(!use_tensor_products);
12615 Vector Enodes(vdim*ND*NE, my_d_mt);
12616 elem_restr->Mult(nodes, Enodes);
12617 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
12621 qi->Mult(nodes, eval_flags,
X,
J,
detJ);
12637 const int vdim = fespace->
GetVDim();
12651 face_restr->
Mult(*nodes, Fnodes);
12653 unsigned eval_flags = 0;
12662 J.
SetSize(vdim*vdim*NQ*NF, my_d_mt);
12696 T.Transform(ip, tip);
12700 V(1) = s * ((ip.
y + layer) / n);
12705 V(2) = s * ((ip.
z + layer) / n);
12714 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
12718 int nvy = (closed) ? (ny) : (ny + 1);
12719 int nvt = mesh->
GetNV() * nvy;
12728 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
12733 for (
int i = 0; i < mesh->
GetNV(); i++)
12736 for (
int j = 0; j < nvy; j++)
12738 vc[1] = sy * (double(j) / ny);
12744 for (
int i = 0; i < mesh->
GetNE(); i++)
12749 for (
int j = 0; j < ny; j++)
12752 qv[0] = vert[0] * nvy + j;
12753 qv[1] = vert[1] * nvy + j;
12754 qv[2] = vert[1] * nvy + (j + 1) % nvy;
12755 qv[3] = vert[0] * nvy + (j + 1) % nvy;
12761 for (
int i = 0; i < mesh->
GetNBE(); i++)
12766 for (
int j = 0; j < ny; j++)
12769 sv[0] = vert[0] * nvy + j;
12770 sv[1] = vert[0] * nvy + (j + 1) % nvy;
12774 Swap<int>(sv[0], sv[1]);
12786 for (
int i = 0; i < mesh->
GetNE(); i++)
12792 sv[0] = vert[0] * nvy;
12793 sv[1] = vert[1] * nvy;
12797 sv[0] = vert[1] * nvy + ny;
12798 sv[1] = vert[0] * nvy + ny;
12814 string cname = name;
12815 if (cname ==
"Linear")
12819 else if (cname ==
"Quadratic")
12823 else if (cname ==
"Cubic")
12827 else if (!strncmp(name,
"H1_", 3))
12831 else if (!strncmp(name,
"L2_T", 4))
12835 else if (!strncmp(name,
"L2_", 3))
12842 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : " 12854 for (
int i = 0; i < mesh->
GetNE(); i++)
12857 for (
int j = ny-1; j >= 0; j--)
12874 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
12879 int nvt = mesh->
GetNV() * nvz;
12884 bool wdgMesh =
false;
12885 bool hexMesh =
false;
12889 for (
int i = 0; i < mesh->
GetNV(); i++)
12893 for (
int j = 0; j < nvz; j++)
12895 vc[2] = sz * (double(j) / nz);
12901 for (
int i = 0; i < mesh->
GetNE(); i++)
12911 for (
int j = 0; j < nz; j++)
12914 pv[0] = vert[0] * nvz + j;
12915 pv[1] = vert[1] * nvz + j;
12916 pv[2] = vert[2] * nvz + j;
12917 pv[3] = vert[0] * nvz + (j + 1) % nvz;
12918 pv[4] = vert[1] * nvz + (j + 1) % nvz;
12919 pv[5] = vert[2] * nvz + (j + 1) % nvz;
12926 for (
int j = 0; j < nz; j++)
12929 hv[0] = vert[0] * nvz + j;
12930 hv[1] = vert[1] * nvz + j;
12931 hv[2] = vert[2] * nvz + j;
12932 hv[3] = vert[3] * nvz + j;
12933 hv[4] = vert[0] * nvz + (j + 1) % nvz;
12934 hv[5] = vert[1] * nvz + (j + 1) % nvz;
12935 hv[6] = vert[2] * nvz + (j + 1) % nvz;
12936 hv[7] = vert[3] * nvz + (j + 1) % nvz;
12938 mesh3d->
AddHex(hv, attr);
12942 mfem::err <<
"Extrude2D : Invalid 2D element type \'" 12943 << geom <<
"\'" << endl;
12949 for (
int i = 0; i < mesh->
GetNBE(); i++)
12954 for (
int j = 0; j < nz; j++)
12957 qv[0] = vert[0] * nvz + j;
12958 qv[1] = vert[1] * nvz + j;
12959 qv[2] = vert[1] * nvz + (j + 1) % nvz;
12960 qv[3] = vert[0] * nvz + (j + 1) % nvz;
12969 for (
int i = 0; i < mesh->
GetNE(); i++)
12980 tv[0] = vert[0] * nvz;
12981 tv[1] = vert[2] * nvz;
12982 tv[2] = vert[1] * nvz;
12986 tv[0] = vert[0] * nvz + nz;
12987 tv[1] = vert[1] * nvz + nz;
12988 tv[2] = vert[2] * nvz + nz;
12996 qv[0] = vert[0] * nvz;
12997 qv[1] = vert[3] * nvz;
12998 qv[2] = vert[2] * nvz;
12999 qv[3] = vert[1] * nvz;
13003 qv[0] = vert[0] * nvz + nz;
13004 qv[1] = vert[1] * nvz + nz;
13005 qv[2] = vert[2] * nvz + nz;
13006 qv[3] = vert[3] * nvz + nz;
13012 mfem::err <<
"Extrude2D : Invalid 2D element type \'" 13013 << geom <<
"\'" << endl;
13019 if ( hexMesh && wdgMesh )
13023 else if ( hexMesh )
13027 else if ( wdgMesh )
13040 string cname = name;
13041 if (cname ==
"Linear")
13045 else if (cname ==
"Quadratic")
13049 else if (cname ==
"Cubic")
13053 else if (!strncmp(name,
"H1_", 3))
13057 else if (!strncmp(name,
"L2_T", 4))
13061 else if (!strncmp(name,
"L2_", 3))
13068 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : " 13080 for (
int i = 0; i < mesh->
GetNE(); i++)
13083 for (
int j = nz-1; j >= 0; j--)
13104 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
13105 <<
" 0 0 " << i <<
" -1 0\n";
13112 double mid[3] = {0, 0, 0};
13113 for (
int j = 0; j < 2; j++)
13115 for (
int k = 0; k <
spaceDim; k++)
13121 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" " 13122 << 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
Returns the indices of the degrees of freedom for the interior of the specified edge.
Table * GetEdgeVertexTable() const
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
static int ComposeQuadOrientations(int ori_a_b, int ori_b_c)
const Table & ElementToFaceTable() const
Array< Element * > boundary
int Dimension() const
Dimension of the reference space used within the elements.
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...
static Mesh LoadFromFile(const std::string &filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
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)
Return FiniteElement for reference element of the specified 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...
static int InvertTriOrientation(int ori)
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'. The returned indices are offsets in...
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]
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering f_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
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.
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
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. This is the number of Local Degrees of Freedom.
void GetElementJacobian(int i, DenseMatrix &J, const IntegrationPoint *ip=NULL)
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.
virtual void Save(const std::string &fname, int precision=16) const
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
Return pointer to the i'th element object.
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
Returns the indices of the degrees of freedom for the specified vertices.
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...
double Trace() const
Trace of a square matrix.
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)
std::function< double(const Vector &)> f(double mass_coeff)
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.
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.
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()
double UserTime()
Return the number of user seconds elapsed since the stopwatch was started.
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()
Determine the sets of unique attribute values in domain and boundary elements.
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.
const double * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
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 Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
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)
void ReadLineMesh(std::istream &input)
virtual void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
void GetGeometricParametersFromJacobian(const DenseMatrix &J, double &volume, Vector &aspr, Vector &skew, Vector &ori) const
Computes geometric parameters associated with a Jacobian matrix in 2D/3D. These parameters are (1) Ar...
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.
double FNorm2() const
Compute the square of the Frobenius norm of the matrix.
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
Returns the indices of the degrees of freedom for the interior of the specified face.
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
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 indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
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...
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'. The returned indices are offsets into an ldo...
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
Dimension of the physical space containing the mesh.
static const int FaceVert[NumFaces][MaxFaceVert]
const CoarseFineTransformations & GetRefinementTransforms()
void RefineNURBSFromFile(std::string ref_file)
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]
static int InvertQuadOrientation(int ori)
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
Compute a single vdof corresponding to the index dof and the vector index vd.
Mesh & operator=(Mesh &&mesh)
Move assignment operator.
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)
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
Returns the indices of the degrees of freedom for the interior of the specified element.
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
Return pointer to the i'th boundary element object.
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.
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...
void cross3D(const Vector &vin, Vector &vout) const
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.
static int ComposeTriOrientations(int ori_a_b, int ori_b_c)
MPI_Comm GetGlobalMPI_Comm()
Get MFEM's "global" MPI communicator.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
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.
const T & AsConst(const T &a)
Utility function similar to std::as_const in c++17.
Table * GetVertexToElementTable()
The returned Table should be deleted 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.