15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
17 #include "../general/binaryio.hpp"
18 #include "../general/text.hpp"
19 #include "../general/device.hpp"
20 #include "../general/tic_toc.hpp"
21 #include "../general/gecko.hpp"
22 #include "../fem/quadinterpolator.hpp"
37 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
42 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
47 int*,
int*,
int*,
int*,
int*,
idxtype*);
49 int*,
int*,
int*,
int*,
int*,
idxtype*);
51 int*,
int*,
int*,
int*,
int*,
idxtype*);
68 void Mesh::GetElementCenter(
int i,
Vector ¢er)
71 int geom = GetElementBaseGeometry(i);
98 double Mesh::GetElementSize(
int i,
int type)
100 return GetElementSize(GetElementTransformation(i), type);
103 double Mesh::GetElementSize(
int i,
const Vector &dir)
107 GetElementJacobian(i, J);
109 return sqrt((d_hat * d_hat) / (dir * dir));
112 double Mesh::GetElementVolume(
int i)
134 for (
int d = 0; d < spaceDim; d++)
143 for (
int i = 0; i < NumOfVertices; i++)
145 coord = GetVertex(i);
146 for (
int d = 0; d < spaceDim; d++)
148 if (coord[d] < min(d)) { min(d) = coord[d]; }
149 if (coord[d] > max(d)) { max(d) = coord[d]; }
155 const bool use_boundary =
false;
156 int ne = use_boundary ? GetNBE() : GetNE();
164 for (
int i = 0; i < ne; i++)
168 GetBdrElementFace(i, &fn, &fo);
170 Tr = GetFaceElementTransformations(fn, 5);
177 T = GetElementTransformation(i);
181 for (
int j = 0; j < pointmat.
Width(); j++)
183 for (
int d = 0; d < pointmat.
Height(); d++)
185 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
186 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
193 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
194 double &kappa_min,
double &kappa_max,
202 sdim = SpaceDimension();
204 if (Vh) { Vh->
SetSize(NumOfElements); }
205 if (Vk) { Vk->
SetSize(NumOfElements); }
208 h_max = kappa_max = -h_min;
209 if (dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
211 for (i = 0; i < NumOfElements; i++)
213 GetElementJacobian(i, J);
214 h =
pow(fabs(J.
Weight()), 1.0/
double(dim));
215 kappa = (dim == sdim) ?
217 if (Vh) { (*Vh)(i) = h; }
218 if (Vk) { (*Vk)(i) = kappa; }
220 if (h < h_min) { h_min = h; }
221 if (h > h_max) { h_max = h; }
222 if (kappa < kappa_min) { kappa_min =
kappa; }
223 if (kappa > kappa_max) { kappa_max =
kappa; }
228 void Mesh::PrintElementsByGeometry(
int dim,
232 for (
int g = Geometry::DimStart[dim], first = 1;
233 g < Geometry::DimStart[dim+1]; g++)
235 if (!num_elems_by_geom[g]) {
continue; }
236 if (!first) { os <<
" + "; }
238 os << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
242 void Mesh::PrintCharacteristics(
Vector *Vh,
Vector *Vk, std::ostream &os)
244 double h_min, h_max, kappa_min, kappa_max;
246 os <<
"Mesh Characteristics:";
248 this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
250 Array<int> num_elems_by_geom(Geometry::NumGeom);
251 num_elems_by_geom = 0;
252 for (
int i = 0; i < GetNE(); i++)
254 num_elems_by_geom[GetElementBaseGeometry(i)]++;
258 <<
"Dimension : " << Dimension() <<
'\n'
259 <<
"Space dimension : " << SpaceDimension();
263 <<
"Number of vertices : " << GetNV() <<
'\n'
264 <<
"Number of elements : " << GetNE() <<
'\n'
265 <<
"Number of bdr elem : " << GetNBE() <<
'\n';
270 <<
"Number of vertices : " << GetNV() <<
'\n'
271 <<
"Number of elements : " << GetNE() <<
'\n'
272 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
273 <<
"h_min : " << h_min <<
'\n'
274 <<
"h_max : " << h_max <<
'\n';
279 <<
"Number of vertices : " << GetNV() <<
'\n'
280 <<
"Number of edges : " << GetNEdges() <<
'\n'
281 <<
"Number of elements : " << GetNE() <<
" -- ";
282 PrintElementsByGeometry(2, num_elems_by_geom, os);
284 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
285 <<
"Euler Number : " << EulerNumber2D() <<
'\n'
286 <<
"h_min : " << h_min <<
'\n'
287 <<
"h_max : " << h_max <<
'\n'
288 <<
"kappa_min : " << kappa_min <<
'\n'
289 <<
"kappa_max : " << kappa_max <<
'\n';
293 Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
294 num_bdr_elems_by_geom = 0;
295 for (
int i = 0; i < GetNBE(); i++)
297 num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
299 Array<int> num_faces_by_geom(Geometry::NumGeom);
300 num_faces_by_geom = 0;
301 for (
int i = 0; i < GetNFaces(); i++)
303 num_faces_by_geom[GetFaceBaseGeometry(i)]++;
307 <<
"Number of vertices : " << GetNV() <<
'\n'
308 <<
"Number of edges : " << GetNEdges() <<
'\n'
309 <<
"Number of faces : " << GetNFaces() <<
" -- ";
310 PrintElementsByGeometry(Dim-1, num_faces_by_geom, os);
312 <<
"Number of elements : " << GetNE() <<
" -- ";
313 PrintElementsByGeometry(Dim, num_elems_by_geom, os);
315 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
316 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, os);
318 <<
"Euler Number : " << EulerNumber() <<
'\n'
319 <<
"h_min : " << h_min <<
'\n'
320 <<
"h_max : " << h_max <<
'\n'
321 <<
"kappa_min : " << kappa_min <<
'\n'
322 <<
"kappa_max : " << kappa_max <<
'\n';
324 os <<
'\n' << std::flush;
331 case Element::POINT :
return &
PointFE;
332 case Element::SEGMENT :
return &
SegmentFE;
337 case Element::WEDGE :
return &
WedgeFE;
338 case Element::PYRAMID :
return &
PyramidFE;
340 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
343 MFEM_ABORT(
"Unknown element type");
352 ElTr->
ElementType = ElementTransformation::ELEMENT;
358 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
364 Nodes->FESpace()->GetElementVDofs(i, vdofs);
367 int n = vdofs.
Size()/spaceDim;
369 for (
int k = 0; k < spaceDim; k++)
371 for (
int j = 0; j < n; j++)
373 pm(k,j) = nodes(vdofs[n*k+j]);
376 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
380 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
385 ElTr->
ElementType = ElementTransformation::ELEMENT;
392 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
393 int nv = elements[i]->GetNVertices();
394 const int *v = elements[i]->GetVertices();
395 int n = vertices.Size();
397 for (
int k = 0; k < spaceDim; k++)
399 for (
int j = 0; j < nv; j++)
401 pm(k, j) = nodes(k*n+v[j]);
404 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
408 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
410 Nodes->FESpace()->GetElementVDofs(i, vdofs);
411 int n = vdofs.
Size()/spaceDim;
413 for (
int k = 0; k < spaceDim; k++)
415 for (
int j = 0; j < n; j++)
417 pm(k,j) = nodes(vdofs[n*k+j]);
420 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
426 GetElementTransformation(i, &Transformation);
428 return &Transformation;
433 GetBdrElementTransformation(i, &BdrTransformation);
434 return &BdrTransformation;
441 ElTr->
ElementType = ElementTransformation::BDR_ELEMENT;
447 GetBdrPointMatrix(i, pm);
448 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
458 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
459 int n = vdofs.
Size()/spaceDim;
461 for (
int k = 0; k < spaceDim; k++)
463 for (
int j = 0; j < n; j++)
465 pm(k,j) = nodes(vdofs[n*k+j]);
472 int elem_id, face_info;
473 GetBdrElementAdjacentElement(i, elem_id, face_info);
475 GetLocalFaceTransformation(GetBdrElementType(i),
476 GetElementType(elem_id),
477 FaceElemTr.Loc1.Transf, face_info);
482 Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
483 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
484 "Mesh requires nodal Finite Element.");
486 FaceElemTr.Loc1.Transf.ElementNo = elem_id;
487 FaceElemTr.Loc1.Transf.mesh =
this;
488 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
489 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
490 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
492 ElTr->
SetFE(face_el);
499 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
507 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
508 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
510 for (
int i = 0; i < spaceDim; i++)
512 for (
int j = 0; j < nv; j++)
514 pm(i, j) = vertices[v[j]](i);
517 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
521 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
527 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
528 int n = vdofs.
Size()/spaceDim;
530 for (
int i = 0; i < spaceDim; i++)
532 for (
int j = 0; j < n; j++)
534 pm(i, j) = nodes(vdofs[n*i+j]);
541 FaceInfo &face_info = faces_info[FaceNo];
546 GetLocalFaceTransformation(face_type,
547 GetElementType(face_info.
Elem1No),
548 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
551 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
553 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
554 "Mesh requires nodal Finite Element.");
557 FaceElemTr.Loc1.Transf.ElementNo = face_info.
Elem1No;
558 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
559 FaceElemTr.Loc1.Transf.mesh =
this;
560 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
561 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
570 GetFaceTransformation(FaceNo, &FaceTransformation);
571 return &FaceTransformation;
578 GetFaceTransformation(EdgeNo, EdTr);
583 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
595 GetEdgeVertices(EdgeNo, v);
598 for (
int i = 0; i < spaceDim; i++)
600 for (
int j = 0; j < nv; j++)
602 pm(i, j) = vertices[v[j]](i);
605 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
609 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
615 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
616 int n = vdofs.
Size()/spaceDim;
618 for (
int i = 0; i < spaceDim; i++)
620 for (
int j = 0; j < n; j++)
622 pm(i, j) = nodes(vdofs[n*i+j]);
625 EdTr->
SetFE(edge_el);
629 MFEM_ABORT(
"Not implemented.");
636 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
637 return &EdgeTransformation;
641 void Mesh::GetLocalPtToSegTransformation(
656 void Mesh::GetLocalSegToTriTransformation(
665 tv = tri_t::Edges[i/64];
666 so = seg_t::Orient[i%64];
669 for (
int j = 0; j < 2; j++)
671 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
672 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
676 void Mesh::GetLocalSegToQuadTransformation(
685 qv = quad_t::Edges[i/64];
686 so = seg_t::Orient[i%64];
689 for (
int j = 0; j < 2; j++)
691 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
692 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
696 void Mesh::GetLocalTriToTetTransformation(
704 const int *tv = tet_t::FaceVert[i/64];
707 const int *to = tri_t::Orient[i%64];
711 for (
int j = 0; j < 3; j++)
714 locpm(0, j) = vert.
x;
715 locpm(1, j) = vert.
y;
716 locpm(2, j) = vert.
z;
720 void Mesh::GetLocalTriToWdgTransformation(
728 MFEM_VERIFY(i < 128,
"Local face index " << i/64
729 <<
" is not a triangular face of a wedge.");
730 const int *pv = pri_t::FaceVert[i/64];
733 const int *to = tri_t::Orient[i%64];
737 for (
int j = 0; j < 3; j++)
740 locpm(0, j) = vert.
x;
741 locpm(1, j) = vert.
y;
742 locpm(2, j) = vert.
z;
746 void Mesh::GetLocalTriToPyrTransformation(
753 MFEM_VERIFY(i >= 64,
"Local face index " << i/64
754 <<
" is not a triangular face of a pyramid.");
755 const int *pv = pyr_t::FaceVert[i/64];
758 const int *to = tri_t::Orient[i%64];
762 for (
int j = 0; j < 3; j++)
765 locpm(0, j) = vert.
x;
766 locpm(1, j) = vert.
y;
767 locpm(2, j) = vert.
z;
771 void Mesh::GetLocalQuadToHexTransformation(
779 const int *hv = hex_t::FaceVert[i/64];
781 const int *qo = quad_t::Orient[i%64];
784 for (
int j = 0; j < 4; j++)
787 locpm(0, j) = vert.
x;
788 locpm(1, j) = vert.
y;
789 locpm(2, j) = vert.
z;
793 void Mesh::GetLocalQuadToWdgTransformation(
801 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
802 <<
" is not a quadrilateral face of a wedge.");
803 const int *pv = pri_t::FaceVert[i/64];
805 const int *qo = quad_t::Orient[i%64];
808 for (
int j = 0; j < 4; j++)
811 locpm(0, j) = vert.
x;
812 locpm(1, j) = vert.
y;
813 locpm(2, j) = vert.
z;
817 void Mesh::GetLocalQuadToPyrTransformation(
824 MFEM_VERIFY(i < 64,
"Local face index " << i/64
825 <<
" is not a quadrilateral face of a pyramid.");
826 const int *pv = pyr_t::FaceVert[i/64];
828 const int *qo = quad_t::Orient[i%64];
831 for (
int j = 0; j < 4; j++)
834 locpm(0, j) = vert.
x;
835 locpm(1, j) = vert.
y;
836 locpm(2, j) = vert.
z;
844 for (
int i = 0; i < geom_factors.Size(); i++)
856 geom_factors.Append(gf);
864 for (
int i = 0; i < face_geom_factors.Size(); i++)
877 face_geom_factors.Append(gf);
881 void Mesh::DeleteGeometricFactors()
883 for (
int i = 0; i < geom_factors.Size(); i++)
885 delete geom_factors[i];
887 geom_factors.SetSize(0);
888 for (
int i = 0; i < face_geom_factors.Size(); i++)
890 delete face_geom_factors[i];
892 face_geom_factors.SetSize(0);
895 void Mesh::GetLocalFaceTransformation(
901 GetLocalPtToSegTransformation(Transf, info);
904 case Element::SEGMENT:
905 if (elem_type == Element::TRIANGLE)
907 GetLocalSegToTriTransformation(Transf, info);
911 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
912 GetLocalSegToQuadTransformation(Transf, info);
916 case Element::TRIANGLE:
917 if (elem_type == Element::TETRAHEDRON)
919 GetLocalTriToTetTransformation(Transf, info);
921 else if (elem_type == Element::WEDGE)
923 GetLocalTriToWdgTransformation(Transf, info);
925 else if (elem_type == Element::PYRAMID)
927 GetLocalTriToPyrTransformation(Transf, info);
931 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
932 "face type " << face_type
933 <<
" and element type " << elem_type <<
"\n");
937 case Element::QUADRILATERAL:
938 if (elem_type == Element::HEXAHEDRON)
940 GetLocalQuadToHexTransformation(Transf, info);
942 else if (elem_type == Element::WEDGE)
944 GetLocalQuadToWdgTransformation(Transf, info);
946 else if (elem_type == Element::PYRAMID)
948 GetLocalQuadToPyrTransformation(Transf, info);
952 MFEM_ABORT(
"Mesh::GetLocalFaceTransformation not defined for "
953 "face type " << face_type
954 <<
" and element type " << elem_type <<
"\n");
963 FaceInfo &face_info = faces_info[FaceNo];
966 FaceElemTr.SetConfigurationMask(cmask);
967 FaceElemTr.Elem1 = NULL;
968 FaceElemTr.Elem2 = NULL;
972 if (mask & FaceElementTransformations::HAVE_ELEM1)
974 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
975 FaceElemTr.Elem1 = &Transformation;
982 FaceElemTr.Elem2No = face_info.
Elem2No;
983 if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
984 FaceElemTr.Elem2No >= 0)
987 if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
988 { MFEM_ABORT(
"NURBS mesh not supported!"); }
990 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
991 FaceElemTr.Elem2 = &Transformation2;
996 if (mask & FaceElementTransformations::HAVE_FACE)
998 GetFaceTransformation(FaceNo, &FaceElemTr);
1003 FaceElemTr.SetGeometryType(GetFaceGeometryType(FaceNo));
1007 int face_type = GetFaceElementType(FaceNo);
1008 if (mask & FaceElementTransformations::HAVE_LOC1)
1010 int elem_type = GetElementType(face_info.
Elem1No);
1011 GetLocalFaceTransformation(face_type, elem_type,
1012 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
1015 if ((mask & FaceElementTransformations::HAVE_LOC2) &&
1016 FaceElemTr.Elem2No >= 0)
1018 int elem_type = GetElementType(face_info.
Elem2No);
1019 GetLocalFaceTransformation(face_type, elem_type,
1020 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
1023 if (Nonconforming() && IsSlaveFace(face_info))
1025 ApplyLocalSlaveTransformation(FaceElemTr, face_info,
false);
1030 FaceElemTr.SetConfigurationMask(cmask);
1036 double dist = FaceElemTr.CheckConsistency();
1039 mfem::out <<
"\nInternal error: face id = " << FaceNo
1040 <<
", dist = " << dist <<
'\n';
1041 FaceElemTr.CheckConsistency(1);
1042 MFEM_ABORT(
"internal error");
1052 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
1058 #ifdef MFEM_THREAD_SAFE
1063 MFEM_ASSERT(fi.
NCFace >= 0,
"");
1064 MFEM_ASSERT(nc_faces_info[fi.
NCFace].Slave,
"internal error");
1075 std::swap(composition(0,0), composition(0,1));
1076 std::swap(composition(1,0), composition(1,1));
1097 int fn = GetBdrFace(BdrElemNo);
1100 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
1104 tr = GetFaceElementTransformations(fn, 21);
1105 tr->
Attribute = boundary[BdrElemNo]->GetAttribute();
1107 tr->
ElementType = ElementTransformation::BDR_FACE;
1112 int Mesh::GetBdrFace(
int BdrElemNo)
const
1117 fn = be_to_face[BdrElemNo];
1121 fn = be_to_edge[BdrElemNo];
1125 fn = boundary[BdrElemNo]->GetVertices()[0];
1136 GetFaceElements(f, &e1, &e2);
1137 GetFaceInfos(f, &inf1, &inf2, &ncface);
1147 if (f < GetNumFaces())
1153 face.
tag = FaceInfoTag::LocalConforming;
1154 face.
topology = FaceTopology::Conforming;
1163 face.
tag = FaceInfoTag::LocalSlaveNonconforming;
1164 face.
topology = FaceTopology::Nonconforming;
1169 MFEM_ASSERT(inf2%64==0,
"unexpected slave face orientation.");
1180 face.
tag = FaceInfoTag::Boundary;
1181 face.
topology = FaceTopology::Boundary;
1190 face.
tag = FaceInfoTag::SharedConforming;
1191 face.
topology = FaceTopology::Conforming;
1203 face.
tag = FaceInfoTag::MasterNonconforming;
1204 face.
topology = FaceTopology::Nonconforming;
1213 face.
tag = FaceInfoTag::SharedSlaveNonconforming;
1214 face.
topology = FaceTopology::Nonconforming;
1229 face.
tag = FaceInfoTag::GhostMaster;
1239 face.
tag = FaceInfoTag::GhostSlave;
1240 face.
topology = FaceTopology::Nonconforming;
1257 case FaceInfoTag::LocalConforming:
1258 res.
Elem1No = element[0].index;
1259 res.Elem2No = element[1].index;
1260 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1261 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1262 res.NCFace = ncface;
1264 case FaceInfoTag::LocalSlaveNonconforming:
1265 res.Elem1No = element[0].index;
1266 res.Elem2No = element[1].index;
1267 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1268 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1269 res.NCFace = ncface;
1271 case FaceInfoTag::Boundary:
1272 res.Elem1No = element[0].index;
1273 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1275 case FaceInfoTag::SharedConforming:
1276 res.Elem1No = element[0].index;
1277 res.Elem2No = -1 - element[1].index;
1278 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1279 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1281 case FaceInfoTag::MasterNonconforming:
1282 res.Elem1No = element[0].index;
1283 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1285 case FaceInfoTag::SharedSlaveNonconforming:
1286 res.Elem1No = element[0].index;
1287 res.Elem2No = -1 - element[1].index;
1288 res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
1289 res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
1291 case FaceInfoTag::GhostMaster:
1293 case FaceInfoTag::GhostSlave:
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;
1305 os <<
"face topology=";
1308 case Mesh::FaceTopology::Boundary:
1311 case Mesh::FaceTopology::Conforming:
1314 case Mesh::FaceTopology::Nonconforming:
1315 os <<
"Non-conforming";
1317 case Mesh::FaceTopology::NA:
1321 os <<
"element[0].location=";
1324 case Mesh::ElementLocation::Local:
1327 case Mesh::ElementLocation::FaceNbr:
1330 case Mesh::ElementLocation::NA:
1335 os <<
"element[1].location=";
1338 case Mesh::ElementLocation::Local:
1341 case Mesh::ElementLocation::FaceNbr:
1344 case Mesh::ElementLocation::NA:
1349 os <<
"element[0].conformity=";
1352 case Mesh::ElementConformity::Coincident:
1355 case Mesh::ElementConformity::Superset:
1358 case Mesh::ElementConformity::Subset:
1361 case Mesh::ElementConformity::NA:
1366 os <<
"element[1].conformity=";
1369 case Mesh::ElementConformity::Coincident:
1372 case Mesh::ElementConformity::Superset:
1375 case Mesh::ElementConformity::Subset:
1378 case Mesh::ElementConformity::NA:
1383 os <<
"element[0].index=" << info.
element[0].
index << std::endl
1384 <<
"element[1].index=" << info.
element[1].
index << std::endl
1389 <<
"ncface=" << info.
ncface << std::endl;
1393 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const
1395 *Elem1 = faces_info[Face].Elem1No;
1396 *Elem2 = faces_info[Face].Elem2No;
1399 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const
1401 *Inf1 = faces_info[Face].Elem1Inf;
1402 *Inf2 = faces_info[Face].Elem2Inf;
1405 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2,
int *NCFace)
const
1407 *Inf1 = faces_info[Face].Elem1Inf;
1408 *Inf2 = faces_info[Face].Elem2Inf;
1409 *NCFace = faces_info[Face].NCFace;
1416 case 1:
return Geometry::POINT;
1417 case 2:
return Geometry::SEGMENT;
1419 if (Face < NumOfFaces)
1421 return faces[Face]->GetGeometryType();
1424 const int nc_face_id = faces_info[Face].NCFace;
1425 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1426 return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
1428 return Geometry::INVALID;
1433 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1441 NumOfElements = NumOfBdrElements = 0;
1442 NumOfEdges = NumOfFaces = 0;
1443 nbInteriorFaces = -1;
1444 nbBoundaryFaces = -1;
1445 meshgen = mesh_geoms = 0;
1451 last_operation = Mesh::NONE;
1454 void Mesh::InitTables()
1457 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1460 void Mesh::SetEmpty()
1466 void Mesh::DestroyTables()
1471 DeleteGeometricFactors();
1482 void Mesh::DestroyPointers()
1484 if (own_nodes) {
delete Nodes; }
1490 for (
int i = 0; i < NumOfElements; i++)
1492 FreeElement(elements[i]);
1495 for (
int i = 0; i < NumOfBdrElements; i++)
1497 FreeElement(boundary[i]);
1500 for (
int i = 0; i < faces.Size(); i++)
1502 FreeElement(faces[i]);
1508 void Mesh::Destroy()
1512 elements.DeleteAll();
1513 vertices.DeleteAll();
1514 boundary.DeleteAll();
1516 faces_info.DeleteAll();
1517 nc_faces_info.DeleteAll();
1518 be_to_edge.DeleteAll();
1519 be_to_face.DeleteAll();
1527 CoarseFineTr.Clear();
1529 #ifdef MFEM_USE_MEMALLOC
1533 attributes.DeleteAll();
1534 bdr_attributes.DeleteAll();
1537 void Mesh::ResetLazyData()
1539 delete el_to_el; el_to_el = NULL;
1540 delete face_edge; face_edge = NULL;
1541 delete edge_vertex; edge_vertex = NULL;
1542 DeleteGeometricFactors();
1543 nbInteriorFaces = -1;
1544 nbBoundaryFaces = -1;
1547 void Mesh::SetAttributes()
1552 for (
int i = 0; i < attribs.
Size(); i++)
1554 attribs[i] = GetBdrAttribute(i);
1558 attribs.
Copy(bdr_attributes);
1559 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1561 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1565 for (
int i = 0; i < attribs.
Size(); i++)
1567 attribs[i] = GetAttribute(i);
1571 attribs.
Copy(attributes);
1572 if (attributes.Size() > 0 && attributes[0] <= 0)
1574 MFEM_WARNING(
"Non-positive attributes in the domain!");
1578 void Mesh::InitMesh(
int Dim_,
int spaceDim_,
int NVert,
int NElem,
int NBdrElem)
1583 spaceDim = spaceDim_;
1586 vertices.SetSize(NVert);
1589 elements.SetSize(NElem);
1591 NumOfBdrElements = 0;
1592 boundary.SetSize(NBdrElem);
1595 template<
typename T>
1596 static void CheckEnlarge(
Array<T> &array,
int size)
1598 if (size >= array.
Size()) { array.
SetSize(size + 1); }
1601 int Mesh::AddVertex(
double x,
double y,
double z)
1603 CheckEnlarge(vertices, NumOfVertices);
1604 double *v = vertices[NumOfVertices]();
1608 return NumOfVertices++;
1611 int Mesh::AddVertex(
const double *coords)
1613 CheckEnlarge(vertices, NumOfVertices);
1614 vertices[NumOfVertices].SetCoords(spaceDim, coords);
1615 return NumOfVertices++;
1618 void Mesh::AddVertexParents(
int i,
int p1,
int p2)
1624 if (i < vertices.Size())
1626 double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1627 for (
int j = 0; j < 3; j++)
1629 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1634 int Mesh::AddSegment(
int v1,
int v2,
int attr)
1636 CheckEnlarge(elements, NumOfElements);
1637 elements[NumOfElements] =
new Segment(v1, v2, attr);
1638 return NumOfElements++;
1641 int Mesh::AddSegment(
const int *vi,
int attr)
1643 CheckEnlarge(elements, NumOfElements);
1644 elements[NumOfElements] =
new Segment(vi, attr);
1645 return NumOfElements++;
1648 int Mesh::AddTriangle(
int v1,
int v2,
int v3,
int attr)
1650 CheckEnlarge(elements, NumOfElements);
1651 elements[NumOfElements] =
new Triangle(v1, v2, v3, attr);
1652 return NumOfElements++;
1655 int Mesh::AddTriangle(
const int *vi,
int attr)
1657 CheckEnlarge(elements, NumOfElements);
1658 elements[NumOfElements] =
new Triangle(vi, attr);
1659 return NumOfElements++;
1662 int Mesh::AddQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1664 CheckEnlarge(elements, NumOfElements);
1665 elements[NumOfElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1666 return NumOfElements++;
1669 int Mesh::AddQuad(
const int *vi,
int attr)
1671 CheckEnlarge(elements, NumOfElements);
1673 return NumOfElements++;
1676 int Mesh::AddTet(
int v1,
int v2,
int v3,
int v4,
int attr)
1678 int vi[4] = {v1, v2, v3, v4};
1679 return AddTet(vi, attr);
1682 int Mesh::AddTet(
const int *vi,
int attr)
1684 CheckEnlarge(elements, NumOfElements);
1685 #ifdef MFEM_USE_MEMALLOC
1687 tet = TetMemory.Alloc();
1690 elements[NumOfElements] = tet;
1692 elements[NumOfElements] =
new Tetrahedron(vi, attr);
1694 return NumOfElements++;
1697 int Mesh::AddWedge(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int attr)
1699 CheckEnlarge(elements, NumOfElements);
1700 elements[NumOfElements] =
new Wedge(v1, v2, v3, v4, v5, v6, attr);
1701 return NumOfElements++;
1704 int Mesh::AddWedge(
const int *vi,
int attr)
1706 CheckEnlarge(elements, NumOfElements);
1707 elements[NumOfElements] =
new Wedge(vi, attr);
1708 return NumOfElements++;
1711 int Mesh::AddPyramid(
int v1,
int v2,
int v3,
int v4,
int v5,
int attr)
1713 CheckEnlarge(elements, NumOfElements);
1714 elements[NumOfElements] =
new Pyramid(v1, v2, v3, v4, v5, attr);
1715 return NumOfElements++;
1718 int Mesh::AddPyramid(
const int *vi,
int attr)
1720 CheckEnlarge(elements, NumOfElements);
1721 elements[NumOfElements] =
new Pyramid(vi, attr);
1722 return NumOfElements++;
1725 int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1728 CheckEnlarge(elements, NumOfElements);
1729 elements[NumOfElements] =
1730 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1731 return NumOfElements++;
1734 int Mesh::AddHex(
const int *vi,
int attr)
1736 CheckEnlarge(elements, NumOfElements);
1737 elements[NumOfElements] =
new Hexahedron(vi, attr);
1738 return NumOfElements++;
1741 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1743 static const int hex_to_tet[6][4] =
1745 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1746 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1750 for (
int i = 0; i < 6; i++)
1752 for (
int j = 0; j < 4; j++)
1754 ti[j] = vi[hex_to_tet[i][j]];
1760 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1762 static const int hex_to_wdg[2][6] =
1764 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1768 for (
int i = 0; i < 2; i++)
1770 for (
int j = 0; j < 6; j++)
1772 ti[j] = vi[hex_to_wdg[i][j]];
1778 void Mesh::AddHexAsPyramids(
const int *vi,
int attr)
1780 static const int hex_to_pyr[6][5] =
1782 { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
1783 { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
1787 for (
int i = 0; i < 6; i++)
1789 for (
int j = 0; j < 5; j++)
1791 ti[j] = vi[hex_to_pyr[i][j]];
1793 AddPyramid(ti, attr);
1799 CheckEnlarge(elements, NumOfElements);
1800 elements[NumOfElements] = elem;
1801 return NumOfElements++;
1806 CheckEnlarge(boundary, NumOfBdrElements);
1807 boundary[NumOfBdrElements] = elem;
1808 return NumOfBdrElements++;
1811 int Mesh::AddBdrSegment(
int v1,
int v2,
int attr)
1813 CheckEnlarge(boundary, NumOfBdrElements);
1814 boundary[NumOfBdrElements] =
new Segment(v1, v2, attr);
1815 return NumOfBdrElements++;
1818 int Mesh::AddBdrSegment(
const int *vi,
int attr)
1820 CheckEnlarge(boundary, NumOfBdrElements);
1821 boundary[NumOfBdrElements] =
new Segment(vi, attr);
1822 return NumOfBdrElements++;
1825 int Mesh::AddBdrTriangle(
int v1,
int v2,
int v3,
int attr)
1827 CheckEnlarge(boundary, NumOfBdrElements);
1828 boundary[NumOfBdrElements] =
new Triangle(v1, v2, v3, attr);
1829 return NumOfBdrElements++;
1832 int Mesh::AddBdrTriangle(
const int *vi,
int attr)
1834 CheckEnlarge(boundary, NumOfBdrElements);
1835 boundary[NumOfBdrElements] =
new Triangle(vi, attr);
1836 return NumOfBdrElements++;
1839 int Mesh::AddBdrQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1841 CheckEnlarge(boundary, NumOfBdrElements);
1842 boundary[NumOfBdrElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1843 return NumOfBdrElements++;
1846 int Mesh::AddBdrQuad(
const int *vi,
int attr)
1848 CheckEnlarge(boundary, NumOfBdrElements);
1850 return NumOfBdrElements++;
1853 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1855 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1858 for (
int i = 0; i < 2; i++)
1860 for (
int j = 0; j < 3; j++)
1862 ti[j] = vi[quad_to_tri[i][j]];
1864 AddBdrTriangle(ti, attr);
1868 int Mesh::AddBdrPoint(
int v,
int attr)
1870 CheckEnlarge(boundary, NumOfBdrElements);
1871 boundary[NumOfBdrElements] =
new Point(&v, attr);
1872 return NumOfBdrElements++;
1875 void Mesh::GenerateBoundaryElements()
1878 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1882 for (i = 0; i < boundary.Size(); i++)
1884 FreeElement(boundary[i]);
1894 NumOfBdrElements = 0;
1895 for (i = 0; i < faces_info.Size(); i++)
1897 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1900 boundary.SetSize(NumOfBdrElements);
1901 be2face.
SetSize(NumOfBdrElements);
1902 for (j = i = 0; i < faces_info.Size(); i++)
1904 if (faces_info[i].Elem2No < 0)
1906 boundary[j] = faces[i]->Duplicate(
this);
1913 void Mesh::FinalizeCheck()
1915 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1916 vertices.Size() == 0,
1917 "incorrect number of vertices: preallocated: " << vertices.Size()
1918 <<
", actually added: " << NumOfVertices);
1919 MFEM_VERIFY(elements.Size() == NumOfElements,
1920 "incorrect number of elements: preallocated: " << elements.Size()
1921 <<
", actually added: " << NumOfElements);
1922 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1923 "incorrect number of boundary elements: preallocated: "
1924 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1927 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1930 CheckElementOrientation(fix_orientation);
1934 MarkTriMeshForRefinement();
1939 el_to_edge =
new Table;
1940 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1942 CheckBdrElementOrientation();
1956 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1957 bool fix_orientation)
1960 if (fix_orientation)
1962 CheckElementOrientation(fix_orientation);
1967 el_to_edge =
new Table;
1968 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1970 CheckBdrElementOrientation();
1990 GeckoProgress(
double limit) : limit(limit) { sw.Start(); }
1991 virtual bool quit()
const {
return limit > 0 && sw.UserTime() > limit; }
1994 class GeckoVerboseProgress :
public GeckoProgress
2000 GeckoVerboseProgress(
double limit) : GeckoProgress(limit) {}
2002 virtual void beginorder(
const Graph* graph,
Float cost)
const
2003 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
2004 virtual void endorder(
const Graph* graph,
Float cost)
const
2005 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
2007 virtual void beginiter(
const Graph* graph,
2010 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
2011 << window << std::flush;
2013 virtual void enditer(
const Graph* graph,
Float mincost,
Float cost)
const
2014 {
mfem::out <<
", cost = " << cost << endl; }
2019 int iterations,
int window,
2020 int period,
int seed,
bool verbose,
2026 GeckoProgress progress(time_limit);
2027 GeckoVerboseProgress vprogress(time_limit);
2030 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2037 const Table &my_el_to_el = ElementToElementTable();
2038 for (
int elemid = 0; elemid < GetNE(); ++elemid)
2040 const int *neighid = my_el_to_el.
GetRow(elemid);
2041 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
2043 graph.
insert_arc(elemid + 1, neighid[i] + 1);
2048 graph.
order(&functional, iterations, window, period, seed,
2049 verbose ? &vprogress : &progress);
2055 ordering[gnodeid - 1] = graph.
rank(gnodeid);
2058 return graph.
cost();
2069 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
2070 : coord(coord), dir(dir), points(points), mid(mid) {}
2072 bool operator()(
int i)
const
2074 return (points[3*i + coord] < mid) != dir;
2078 static void HilbertSort2D(
int coord1,
2081 const Array<double> &points,
int *beg,
int *end,
2082 double xmin,
double ymin,
double xmax,
double ymax)
2084 if (end - beg <= 1) {
return; }
2086 double xmid = (xmin + xmax)*0.5;
2087 double ymid = (ymin + ymax)*0.5;
2089 int coord2 = (coord1 + 1) % 2;
2092 int *p0 = beg, *p4 = end;
2093 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
2094 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
2095 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
2099 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
2100 ymin, xmin, ymid, xmid);
2102 if (p1 != p0 || p2 != p4)
2104 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
2105 xmin, ymid, xmid, ymax);
2107 if (p2 != p0 || p3 != p4)
2109 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
2110 xmid, ymid, xmax, ymax);
2114 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
2115 ymid, xmax, ymin, xmid);
2119 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
2120 const Array<double> &points,
int *beg,
int *end,
2121 double xmin,
double ymin,
double zmin,
2122 double xmax,
double ymax,
double zmax)
2124 if (end - beg <= 1) {
return; }
2126 double xmid = (xmin + xmax)*0.5;
2127 double ymid = (ymin + ymax)*0.5;
2128 double zmid = (zmin + zmax)*0.5;
2130 int coord2 = (coord1 + 1) % 3;
2131 int coord3 = (coord1 + 2) % 3;
2134 int *p0 = beg, *p8 = end;
2135 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
2136 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
2137 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
2138 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
2139 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
2140 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
2141 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
2145 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
2146 zmin, xmin, ymin, zmid, xmid, ymid);
2148 if (p1 != p0 || p2 != p8)
2150 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
2151 ymin, zmid, xmin, ymid, zmax, xmid);
2153 if (p2 != p0 || p3 != p8)
2155 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
2156 ymid, zmid, xmin, ymax, zmax, xmid);
2158 if (p3 != p0 || p4 != p8)
2160 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
2161 xmin, ymax, zmid, xmid, ymid, zmin);
2163 if (p4 != p0 || p5 != p8)
2165 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
2166 xmid, ymax, zmid, xmax, ymid, zmin);
2168 if (p5 != p0 || p6 != p8)
2170 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
2171 ymax, zmid, xmax, ymid, zmax, xmid);
2173 if (p6 != p0 || p7 != p8)
2175 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
2176 ymid, zmid, xmax, ymin, zmax, xmid);
2180 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
2181 zmid, xmax, ymin, zmin, xmid, ymid);
2187 MFEM_VERIFY(spaceDim <= 3,
"");
2190 GetBoundingBox(min, max);
2195 if (spaceDim < 3) { points = 0.0; }
2198 for (
int i = 0; i < GetNE(); i++)
2200 GetElementCenter(i, center);
2201 for (
int j = 0; j < spaceDim; j++)
2203 points[3*i + j] = center(j);
2210 indices.
Sort([&](
int a,
int b)
2211 {
return points[3*
a] < points[3*
b]; });
2213 else if (spaceDim == 2)
2216 HilbertSort2D(0,
false,
false,
2217 points, indices.
begin(), indices.
end(),
2218 min(0), min(1), max(0), max(1));
2223 HilbertSort3D(0,
false,
false,
false,
2224 points, indices.
begin(), indices.
end(),
2225 min(0), min(1), min(2), max(0), max(1), max(2));
2230 for (
int i = 0; i < GetNE(); i++)
2232 ordering[indices[i]] = i;
2237 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
2241 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
2246 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
2250 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
2280 old_elem_node_vals.SetSize(GetNE());
2281 nodes_fes = Nodes->FESpace();
2284 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2286 nodes_fes->GetElementVDofs(old_elid, old_dofs);
2288 old_elem_node_vals[old_elid] =
new Vector(vals);
2294 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
2296 int new_elid = ordering[old_elid];
2297 new_elements[new_elid] = elements[old_elid];
2302 if (reorder_vertices)
2307 vertex_ordering = -1;
2309 int new_vertex_ind = 0;
2310 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2312 int *elem_vert = elements[new_elid]->GetVertices();
2313 int nv = elements[new_elid]->GetNVertices();
2314 for (
int vi = 0; vi < nv; ++vi)
2316 int old_vertex_ind = elem_vert[vi];
2317 if (vertex_ordering[old_vertex_ind] == -1)
2319 vertex_ordering[old_vertex_ind] = new_vertex_ind;
2320 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
2330 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
2332 int *elem_vert = elements[new_elid]->GetVertices();
2333 int nv = elements[new_elid]->GetNVertices();
2334 for (
int vi = 0; vi < nv; ++vi)
2336 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
2341 for (
int belid = 0; belid < GetNBE(); ++belid)
2343 int *be_vert = boundary[belid]->GetVertices();
2344 int nv = boundary[belid]->GetNVertices();
2345 for (
int vi = 0; vi < nv; ++vi)
2347 be_vert[vi] = vertex_ordering[be_vert[vi]];
2358 el_to_edge =
new Table;
2359 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2364 GetElementToFaceTable();
2374 last_operation = Mesh::NONE;
2375 nodes_fes->Update(
false);
2378 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2380 int new_elid = ordering[old_elid];
2381 nodes_fes->GetElementVDofs(new_elid, new_dofs);
2382 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
2383 delete old_elem_node_vals[old_elid];
2389 void Mesh::MarkForRefinement()
2395 MarkTriMeshForRefinement();
2399 DSTable v_to_v(NumOfVertices);
2400 GetVertexToVertexTable(v_to_v);
2401 MarkTetMeshForRefinement(v_to_v);
2406 void Mesh::MarkTriMeshForRefinement()
2411 for (
int i = 0; i < NumOfElements; i++)
2413 if (elements[i]->GetType() == Element::TRIANGLE)
2415 GetPointMatrix(i, pmat);
2416 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
2427 for (
int i = 0; i < NumOfVertices; i++)
2432 length_idx[j].one = GetLength(i, it.Column());
2433 length_idx[j].two = j;
2440 for (
int i = 0; i < NumOfEdges; i++)
2442 order[length_idx[i].two] = i;
2446 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
2451 GetEdgeOrdering(v_to_v, order);
2453 for (
int i = 0; i < NumOfElements; i++)
2455 if (elements[i]->GetType() == Element::TETRAHEDRON)
2457 elements[i]->MarkEdge(v_to_v, order);
2460 for (
int i = 0; i < NumOfBdrElements; i++)
2462 if (boundary[i]->GetType() == Element::TRIANGLE)
2464 boundary[i]->MarkEdge(v_to_v, order);
2471 if (*old_v_to_v && *old_elem_vert)
2478 if (*old_v_to_v == NULL)
2480 bool need_v_to_v =
false;
2482 for (
int i = 0; i < GetNEdges(); i++)
2488 if (dofs.
Size() > 0)
2496 *old_v_to_v =
new DSTable(NumOfVertices);
2497 GetVertexToVertexTable(*(*old_v_to_v));
2500 if (*old_elem_vert == NULL)
2502 bool need_elem_vert =
false;
2504 for (
int i = 0; i < GetNE(); i++)
2510 if (dofs.
Size() > 1)
2512 need_elem_vert =
true;
2518 *old_elem_vert =
new Table;
2519 (*old_elem_vert)->
MakeI(GetNE());
2520 for (
int i = 0; i < GetNE(); i++)
2522 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2524 (*old_elem_vert)->MakeJ();
2525 for (
int i = 0; i < GetNE(); i++)
2527 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2528 elements[i]->GetNVertices());
2530 (*old_elem_vert)->ShiftUpI();
2543 const int num_edge_dofs = old_dofs.
Size();
2546 const Vector onodes = *Nodes;
2550 int offset = NumOfVertices * old_dofs.
Size();
2554 if (num_edge_dofs > 0)
2556 DSTable new_v_to_v(NumOfVertices);
2557 GetVertexToVertexTable(new_v_to_v);
2559 for (
int i = 0; i < NumOfVertices; i++)
2563 const int old_i = (*old_v_to_v)(i, it.Column());
2564 const int new_i = it.Index();
2565 if (new_i == old_i) {
continue; }
2567 old_dofs.
SetSize(num_edge_dofs);
2568 new_dofs.
SetSize(num_edge_dofs);
2569 for (
int j = 0; j < num_edge_dofs; j++)
2571 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2572 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2576 for (
int j = 0; j < old_dofs.
Size(); j++)
2578 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2582 offset += NumOfEdges * num_edge_dofs;
2590 Table old_face_vertex;
2591 old_face_vertex.
MakeI(NumOfFaces);
2592 for (
int i = 0; i < NumOfFaces; i++)
2596 old_face_vertex.
MakeJ();
2597 for (
int i = 0; i < NumOfFaces; i++)
2599 faces[i]->GetNVertices());
2603 STable3D *faces_tbl = GetElementToFaceTable(1);
2609 for (
int i = 0; i < NumOfFaces; i++)
2611 const int *old_v = old_face_vertex.
GetRow(i);
2613 switch (old_face_vertex.
RowSize(i))
2616 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2620 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2624 new_fdofs[new_i+1] = old_dofs.
Size();
2629 for (
int i = 0; i < NumOfFaces; i++)
2631 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2634 switch (old_face_vertex.
RowSize(i))
2637 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2638 new_v = faces[new_i]->GetVertices();
2639 new_or = GetTriOrientation(old_v, new_v);
2644 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2645 new_v = faces[new_i]->GetVertices();
2646 new_or = GetQuadOrientation(old_v, new_v);
2653 for (
int j = 0; j < old_dofs.
Size(); j++)
2656 const int old_j = dof_ord[j];
2657 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2661 for (
int j = 0; j < old_dofs.
Size(); j++)
2663 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2685 for (
int i = 0; i < GetNE(); i++)
2687 const int *old_v = old_elem_vert->
GetRow(i);
2688 const int *new_v = elements[i]->GetVertices();
2694 case Geometry::SEGMENT:
2695 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2697 case Geometry::TRIANGLE:
2698 new_or = GetTriOrientation(old_v, new_v);
2700 case Geometry::SQUARE:
2701 new_or = GetQuadOrientation(old_v, new_v);
2703 case Geometry::TETRAHEDRON:
2704 new_or = GetTetOrientation(old_v, new_v);
2708 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2709 <<
" FE collection) are not supported yet!");
2713 MFEM_VERIFY(dof_ord != NULL,
2714 "FE collection '" << fec->
Name()
2715 <<
"' does not define reordering for "
2716 << Geometry::Name[geom] <<
" elements!");
2719 for (
int j = 0; j < new_dofs.
Size(); j++)
2722 const int old_j = dof_ord[j];
2723 new_dofs[old_j] = offset + j;
2725 offset += new_dofs.
Size();
2728 for (
int j = 0; j < old_dofs.
Size(); j++)
2730 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2742 GetElementToFaceTable();
2745 CheckBdrElementOrientation();
2750 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2755 CheckBdrElementOrientation();
2760 last_operation = Mesh::NONE;
2765 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2768 CheckElementOrientation(fix_orientation);
2770 if (NumOfBdrElements == 0)
2772 GetElementToFaceTable();
2774 GenerateBoundaryElements();
2779 DSTable v_to_v(NumOfVertices);
2780 GetVertexToVertexTable(v_to_v);
2781 MarkTetMeshForRefinement(v_to_v);
2784 GetElementToFaceTable();
2787 CheckBdrElementOrientation();
2789 if (generate_edges == 1)
2791 el_to_edge =
new Table;
2792 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2806 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2807 bool fix_orientation)
2810 CheckElementOrientation(fix_orientation);
2812 if (NumOfBdrElements == 0)
2814 GetElementToFaceTable();
2816 GenerateBoundaryElements();
2819 GetElementToFaceTable();
2822 CheckBdrElementOrientation();
2824 if (generate_edges == 1)
2826 el_to_edge =
new Table;
2827 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2841 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2844 CheckElementOrientation(fix_orientation);
2846 GetElementToFaceTable();
2849 if (NumOfBdrElements == 0)
2851 GenerateBoundaryElements();
2854 CheckBdrElementOrientation();
2858 el_to_edge =
new Table;
2859 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2871 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2875 Finalize(refine, fix_orientation);
2878 void Mesh::FinalizeTopology(
bool generate_bdr)
2890 bool generate_edges =
true;
2892 if (spaceDim == 0) { spaceDim = Dim; }
2893 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2897 if (tmp_vertex_parents.Size())
2899 MFEM_VERIFY(ncmesh == NULL,
"");
2900 ncmesh =
new NCMesh(
this);
2904 InitFromNCMesh(*ncmesh);
2905 ncmesh->OnMeshUpdated(
this);
2906 GenerateNCFaceInfo();
2910 tmp_vertex_parents.DeleteAll();
2920 GetElementToFaceTable();
2922 if (NumOfBdrElements == 0 && generate_bdr)
2924 GenerateBoundaryElements();
2925 GetElementToFaceTable();
2934 if (Dim > 1 && generate_edges)
2937 if (!el_to_edge) { el_to_edge =
new Table; }
2938 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2942 if (NumOfBdrElements == 0 && generate_bdr)
2944 GenerateBoundaryElements();
2961 ncmesh->OnMeshUpdated(
this);
2964 GenerateNCFaceInfo();
2971 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2973 if (NURBSext || ncmesh)
2975 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
2976 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
2985 const bool check_orientation =
true;
2986 const bool curved = (Nodes != NULL);
2987 const bool may_change_topology =
2988 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
2989 ( check_orientation && fix_orientation &&
2990 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
2993 Table *old_elem_vert = NULL;
2995 if (curved && may_change_topology)
2997 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
3000 if (check_orientation)
3003 CheckElementOrientation(fix_orientation);
3007 MarkForRefinement();
3010 if (may_change_topology)
3014 DoNodeReorder(old_v_to_v, old_elem_vert);
3015 delete old_elem_vert;
3028 CheckBdrElementOrientation();
3033 if (Dim >= 2 && Dim == spaceDim)
3035 const int num_faces = GetNumFaces();
3036 for (
int i = 0; i < num_faces; i++)
3038 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
3039 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
3040 " Interior face with incompatible orientations.");
3047 double sx,
double sy,
double sz,
bool sfc_ordering)
3051 int NVert, NElem, NBdrElem;
3053 NVert = (nx+1) * (ny+1) * (nz+1);
3054 NElem = nx * ny * nz;
3055 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
3056 if (type == Element::TETRAHEDRON)
3061 else if (type == Element::WEDGE)
3064 NBdrElem += 2*nx*ny;
3066 else if (type == Element::PYRAMID)
3069 NVert += nx * ny * nz;
3072 InitMesh(3, 3, NVert, NElem, NBdrElem);
3078 for (z = 0; z <= nz; z++)
3080 coord[2] = ((double) z / nz) * sz;
3081 for (y = 0; y <= ny; y++)
3083 coord[1] = ((double) y / ny) * sy;
3084 for (x = 0; x <= nx; x++)
3086 coord[0] = ((double) x / nx) * sx;
3091 if (type == Element::PYRAMID)
3093 for (z = 0; z < nz; z++)
3095 coord[2] = (((double) z + 0.5) / nz) * sz;
3096 for (y = 0; y < ny; y++)
3098 coord[1] = (((double) y + 0.5 ) / ny) * sy;
3099 for (x = 0; x < nx; x++)
3101 coord[0] = (((double) x + 0.5 ) / nx) * sx;
3108 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
3109 #define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
3112 if (sfc_ordering && type == Element::HEXAHEDRON)
3115 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
3116 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
3118 for (
int k = 0; k < nx*ny*nz; k++)
3125 ind[0] = VTX(x , y , z );
3126 ind[1] = VTX(x+1, y , z );
3127 ind[2] = VTX(x+1, y+1, z );
3128 ind[3] = VTX(x , y+1, z );
3129 ind[4] = VTX(x , y , z+1);
3130 ind[5] = VTX(x+1, y , z+1);
3131 ind[6] = VTX(x+1, y+1, z+1);
3132 ind[7] = VTX(x , y+1, z+1);
3140 for (z = 0; z < nz; z++)
3142 for (y = 0; y < ny; y++)
3144 for (x = 0; x < nx; x++)
3147 ind[0] = VTX(x , y , z );
3148 ind[1] = VTX(x+1, y , z );
3149 ind[2] = VTX(x+1, y+1, z );
3150 ind[3] = VTX(x , y+1, z );
3151 ind[4] = VTX(x , y , z+1);
3152 ind[5] = VTX(x+1, y , z+1);
3153 ind[6] = VTX(x+1, y+1, z+1);
3154 ind[7] = VTX( x, y+1, z+1);
3156 if (type == Element::TETRAHEDRON)
3158 AddHexAsTets(ind, 1);
3160 else if (type == Element::WEDGE)
3162 AddHexAsWedges(ind, 1);
3164 else if (type == Element::PYRAMID)
3166 ind[8] = VTXP(x, y, z);
3167 AddHexAsPyramids(ind, 1);
3180 for (y = 0; y < ny; y++)
3182 for (x = 0; x < nx; x++)
3185 ind[0] = VTX(x , y , 0);
3186 ind[1] = VTX(x , y+1, 0);
3187 ind[2] = VTX(x+1, y+1, 0);
3188 ind[3] = VTX(x+1, y , 0);
3190 if (type == Element::TETRAHEDRON)
3192 AddBdrQuadAsTriangles(ind, 1);
3194 else if (type == Element::WEDGE)
3196 AddBdrQuadAsTriangles(ind, 1);
3205 for (y = 0; y < ny; y++)
3207 for (x = 0; x < nx; x++)
3210 ind[0] = VTX(x , y , nz);
3211 ind[1] = VTX(x+1, y , nz);
3212 ind[2] = VTX(x+1, y+1, nz);
3213 ind[3] = VTX(x , y+1, nz);
3215 if (type == Element::TETRAHEDRON)
3217 AddBdrQuadAsTriangles(ind, 6);
3219 else if (type == Element::WEDGE)
3221 AddBdrQuadAsTriangles(ind, 6);
3230 for (z = 0; z < nz; z++)
3232 for (y = 0; y < ny; y++)
3235 ind[0] = VTX(0 , y , z );
3236 ind[1] = VTX(0 , y , z+1);
3237 ind[2] = VTX(0 , y+1, z+1);
3238 ind[3] = VTX(0 , y+1, z );
3240 if (type == Element::TETRAHEDRON)
3242 AddBdrQuadAsTriangles(ind, 5);
3251 for (z = 0; z < nz; z++)
3253 for (y = 0; y < ny; y++)
3256 ind[0] = VTX(nx, y , z );
3257 ind[1] = VTX(nx, y+1, z );
3258 ind[2] = VTX(nx, y+1, z+1);
3259 ind[3] = VTX(nx, y , z+1);
3261 if (type == Element::TETRAHEDRON)
3263 AddBdrQuadAsTriangles(ind, 3);
3272 for (x = 0; x < nx; x++)
3274 for (z = 0; z < nz; z++)
3277 ind[0] = VTX(x , 0, z );
3278 ind[1] = VTX(x+1, 0, z );
3279 ind[2] = VTX(x+1, 0, z+1);
3280 ind[3] = VTX(x , 0, z+1);
3282 if (type == Element::TETRAHEDRON)
3284 AddBdrQuadAsTriangles(ind, 2);
3293 for (x = 0; x < nx; x++)
3295 for (z = 0; z < nz; z++)
3298 ind[0] = VTX(x , ny, z );
3299 ind[1] = VTX(x , ny, z+1);
3300 ind[2] = VTX(x+1, ny, z+1);
3301 ind[3] = VTX(x+1, ny, z );
3303 if (type == Element::TETRAHEDRON)
3305 AddBdrQuadAsTriangles(ind, 4);
3317 ofstream test_stream(
"debug.mesh");
3319 test_stream.close();
3328 double sx,
double sy,
3329 bool generate_edges,
bool sfc_ordering)
3338 if (type == Element::QUADRILATERAL)
3340 NumOfVertices = (nx+1) * (ny+1);
3341 NumOfElements = nx * ny;
3342 NumOfBdrElements = 2 * nx + 2 * ny;
3344 vertices.SetSize(NumOfVertices);
3345 elements.SetSize(NumOfElements);
3346 boundary.SetSize(NumOfBdrElements);
3353 for (j = 0; j < ny+1; j++)
3355 cy = ((double) j / ny) * sy;
3356 for (i = 0; i < nx+1; i++)
3358 cx = ((double) i / nx) * sx;
3359 vertices[k](0) = cx;
3360 vertices[k](1) = cy;
3369 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
3370 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
3372 for (k = 0; k < nx*ny; k++)
3376 ind[0] = i + j*(nx+1);
3377 ind[1] = i + 1 +j*(nx+1);
3378 ind[2] = i + 1 + (j+1)*(nx+1);
3379 ind[3] = i + (j+1)*(nx+1);
3386 for (j = 0; j < ny; j++)
3388 for (i = 0; i < nx; i++)
3390 ind[0] = i + j*(nx+1);
3391 ind[1] = i + 1 +j*(nx+1);
3392 ind[2] = i + 1 + (j+1)*(nx+1);
3393 ind[3] = i + (j+1)*(nx+1);
3402 for (i = 0; i < nx; i++)
3404 boundary[i] =
new Segment(i, i+1, 1);
3405 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3408 for (j = 0; j < ny; j++)
3410 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3411 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3415 else if (type == Element::TRIANGLE)
3417 NumOfVertices = (nx+1) * (ny+1);
3418 NumOfElements = 2 * nx * ny;
3419 NumOfBdrElements = 2 * nx + 2 * ny;
3421 vertices.SetSize(NumOfVertices);
3422 elements.SetSize(NumOfElements);
3423 boundary.SetSize(NumOfBdrElements);
3430 for (j = 0; j < ny+1; j++)
3432 cy = ((double) j / ny) * sy;
3433 for (i = 0; i < nx+1; i++)
3435 cx = ((double) i / nx) * sx;
3436 vertices[k](0) = cx;
3437 vertices[k](1) = cy;
3444 for (j = 0; j < ny; j++)
3446 for (i = 0; i < nx; i++)
3448 ind[0] = i + j*(nx+1);
3449 ind[1] = i + 1 + (j+1)*(nx+1);
3450 ind[2] = i + (j+1)*(nx+1);
3453 ind[1] = i + 1 + j*(nx+1);
3454 ind[2] = i + 1 + (j+1)*(nx+1);
3462 for (i = 0; i < nx; i++)
3464 boundary[i] =
new Segment(i, i+1, 1);
3465 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3468 for (j = 0; j < ny; j++)
3470 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3471 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3478 MFEM_ABORT(
"Unsupported element type.");
3482 CheckElementOrientation();
3484 if (generate_edges == 1)
3486 el_to_edge =
new Table;
3487 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3489 CheckBdrElementOrientation();
3498 attributes.Append(1);
3499 bdr_attributes.Append(1); bdr_attributes.Append(2);
3500 bdr_attributes.Append(3); bdr_attributes.Append(4);
3505 void Mesh::Make1D(
int n,
double sx)
3514 NumOfVertices = n + 1;
3516 NumOfBdrElements = 2;
3517 vertices.SetSize(NumOfVertices);
3518 elements.SetSize(NumOfElements);
3519 boundary.SetSize(NumOfBdrElements);
3522 for (j = 0; j < n+1; j++)
3524 vertices[j](0) = ((
double) j / n) * sx;
3528 for (j = 0; j < n; j++)
3530 elements[j] =
new Segment(j, j+1, 1);
3535 boundary[0] =
new Point(ind, 1);
3537 boundary[1] =
new Point(ind, 2);
3545 attributes.Append(1);
3546 bdr_attributes.Append(1); bdr_attributes.Append(2);
3549 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
3567 last_operation = Mesh::NONE;
3570 elements.SetSize(NumOfElements);
3571 for (
int i = 0; i < NumOfElements; i++)
3573 elements[i] = mesh.
elements[i]->Duplicate(
this);
3580 boundary.SetSize(NumOfBdrElements);
3581 for (
int i = 0; i < NumOfBdrElements; i++)
3583 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
3602 faces.SetSize(mesh.
faces.Size());
3603 for (
int i = 0; i < faces.Size(); i++)
3606 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
3640 if (dynamic_cast<const ParMesh*>(&mesh))
3652 if (mesh.
Nodes && copy_nodes)
3657 FiniteElementCollection::New(fec->
Name());
3661 Nodes->MakeOwner(fec_copy);
3662 *Nodes = *mesh.
Nodes;
3684 bool fix_orientation)
3688 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
3689 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
3703 double sx,
double sy,
bool sfc_ordering)
3706 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
3713 double sx,
double sy,
double sz,
bool sfc_ordering)
3716 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
3725 ref_factors = ref_factor;
3739 bool fix_orientation)
3748 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
3752 Load(imesh, generate_edges, refine, fix_orientation);
3757 bool fix_orientation)
3760 Load(input, generate_edges, refine, fix_orientation);
3769 "Not enough vertices in external array : "
3770 "len_vertex_data = "<< len_vertex_data <<
", "
3773 if (vertex_data == (
double *)(
vertices.GetData()))
3775 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
3780 memcpy(vertex_data,
vertices.GetData(),
3789 int *element_attributes,
int num_elements,
3791 int *boundary_attributes,
int num_boundary_elements,
3792 int dimension,
int space_dimension)
3794 if (space_dimension == -1)
3796 space_dimension = dimension;
3799 InitMesh(dimension, space_dimension, 0, num_elements,
3800 num_boundary_elements);
3803 int boundary_index_stride = num_boundary_elements > 0 ?
3807 vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
3810 for (
int i = 0; i < num_elements; i++)
3813 elements[i]->SetVertices(element_indices + i * element_index_stride);
3814 elements[i]->SetAttribute(element_attributes[i]);
3818 for (
int i = 0; i < num_boundary_elements; i++)
3821 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3822 boundary[i]->SetAttribute(boundary_attributes[i]);
3838 #ifdef MFEM_USE_MEMALLOC
3847 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3860 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3863 for (
int i = 0; i < nv; i++)
3876 for (
int j = 0; j < nv; j++)
3948 MFEM_ABORT(
"invalid element type: " << type);
3955 std::string parse_tag)
3957 int curved = 0, read_gf = 1;
3958 bool finalize_topo =
true;
3962 MFEM_ABORT(
"Input stream is not open");
3969 getline(input, mesh_type);
3973 int mfem_version = 0;
3974 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
3975 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
3979 int mfem_nc_version = 0;
3980 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
3981 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
3989 if (mfem_version == 12 && parse_tag.empty())
3991 parse_tag =
"mfem_mesh_end";
3995 else if (mfem_nc_version)
3997 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
4004 MFEM_VERIFY(mfem_nc_version >= 10,
4005 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
4006 "used to load a parallel nonconforming mesh, sorry.");
4009 input, mfem_nc_version, curved, is_nc);
4014 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
4027 else if (mesh_type ==
"linemesh")
4031 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
4033 if (mesh_type ==
"curved_areamesh2")
4039 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
4043 else if (mesh_type ==
"TrueGrid")
4047 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
4049 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
4051 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
4052 "Unsupported VTK format");
4053 ReadVTKMesh(input, curved, read_gf, finalize_topo);
4055 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
4059 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
4063 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
4068 else if (mesh_type ==
"$MeshFormat")
4073 ((mesh_type.size() > 2 &&
4074 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
4075 (mesh_type.size() > 3 &&
4076 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
4081 #ifdef MFEM_USE_NETCDF
4084 MFEM_ABORT(
"NetCDF support requires configuration with"
4085 " MFEM_USE_NETCDF=YES");
4091 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
4092 " Use mfem::named_ifgzstream for input.");
4098 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
4124 bool generate_bdr =
false;
4129 if (curved && read_gf)
4143 if (mfem_version == 12)
4149 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
4150 getline(input, line);
4156 if (line ==
"mfem_mesh_end") {
break; }
4158 while (line != parse_tag);
4160 else if (mfem_nc_version >= 10)
4165 MFEM_VERIFY(ident ==
"mfem_mesh_end",
4166 "invalid mesh: end of file tag not found");
4174 int i, j, ie, ib, iv, *v, nv;
4205 for (i = 0; i < num_pieces; i++)
4212 for (i = 0; i < num_pieces; i++)
4218 for (j = 0; j < m->
GetNE(); j++)
4223 for (j = 0; j < m->
GetNBE(); j++)
4228 for (
int k = 0; k < nv; k++)
4230 v[k] = lvert_vert[v[k]];
4235 for (j = 0; j < m->
GetNV(); j++)
4247 for (i = 0; i < num_pieces; i++)
4258 for (i = 0; i < num_pieces; i++)
4262 for (j = 0; j < m->
GetNE(); j++)
4267 for (
int k = 0; k < nv; k++)
4274 for (j = 0; j < m->
GetNBE(); j++)
4279 for (
int k = 0; k < nv; k++)
4286 for (j = 0; j < m->
GetNV(); j++)
4300 for (i = 0; i < num_pieces; i++)
4302 gf_array[i] = mesh_array[i]->
GetNodes();
4317 ref_factors = ref_factor;
4328 int orig_ne = orig_mesh.
GetNE();
4329 MFEM_VERIFY(ref_factors.
Size() == orig_ne && orig_ne > 0,
4330 "Number of refinement factors must equal number of elements")
4331 MFEM_VERIFY(ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
4334 "Invalid refinement type. Must use closed basis type.");
4336 int min_ref = ref_factors.
Min();
4337 int max_ref = ref_factors.
Max();
4339 bool var_order = (min_ref != max_ref);
4352 for (
int i = 0; i < orig_ne; i++)
4370 for (
int el = 0; el < orig_ne; el++)
4382 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
4383 for (
int i = 0; i < phys_pts.
Width(); i++)
4392 for (
int k = 0; k < nvert; k++)
4395 v[k] = rdofs[c2h_map[cid]];
4408 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
4419 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
4425 for (
int k = 0; k < nvert; k++)
4428 v[k] = rdofs[c2h_map[cid]];
4450 for (
int iel = 0; iel < orig_ne; iel++)
4459 const int *node_map = NULL;
4462 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
4463 const int *vertex_map = vertex_fec.
GetDofMap(geom);
4464 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
4465 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
4468 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
4471 int iv = vertex_map[iv_lex];
4473 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
4475 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
4478 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
4489 using GeomRef = std::pair<Geometry::Type, int>;
4490 std::map<GeomRef, int> point_matrices_offsets;
4492 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4496 GeomRef id(geom, ref_factors[el_coarse]);
4497 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
4503 point_matrices_offsets[id] = n_point_matrices[geom];
4504 n_point_matrices[geom] += nref_el;
4511 int nmatrices = n_point_matrices[geom];
4518 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4521 int ref = ref_factors[el_coarse];
4522 int offset = point_matrices_offsets[GeomRef(geom, ref)];
4528 for (
int k = 0; k < nvert; k++)
4560 "Mesh::MakeSimplicial requires a properly oriented input mesh");
4562 "Mesh::MakeSimplicial does not support non-conforming meshes.")
4569 Mesh copy(orig_mesh);
4574 int nv = orig_mesh.
GetNV();
4575 int ne = orig_mesh.
GetNE();
4576 int nbe = orig_mesh.
GetNBE();
4589 int new_ne = 0, new_nbe = 0;
4590 for (
int i=0; i<ne; ++i)
4594 for (
int i=0; i<nbe; ++i)
4603 for (
int i=0; i<nv; ++i)
4613 if (vglobal == NULL)
4616 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
4617 vglobal = vglobal_id.
GetData();
4620 constexpr
int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
4621 constexpr
int quad_ntris = 2, prism_ntets = 3;
4622 static const int quad_trimap[2][nv_tri*quad_ntris] =
4634 static const int prism_rot[nv_prism*nv_prism] =
4643 static const int prism_f[nv_quad] = {1, 2, 5, 4};
4644 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
4658 static const int hex_rot[nv_hex*nv_hex] =
4660 0, 1, 2, 3, 4, 5, 6, 7,
4661 1, 0, 4, 5, 2, 3, 7, 6,
4662 2, 1, 5, 6, 3, 0, 4, 7,
4663 3, 0, 1, 2, 7, 4, 5, 6,
4664 4, 0, 3, 7, 5, 1, 2, 6,
4665 5, 1, 0, 4, 6, 2, 3, 7,
4666 6, 2, 1, 5, 7, 3, 0, 4,
4667 7, 3, 2, 6, 4, 0, 1, 5
4669 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
4670 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
4671 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
4672 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
4673 static const int hex_tetmap0[nv_tet*5] =
4680 static const int hex_tetmap1[nv_tet*6] =
4687 static const int hex_tetmap2[nv_tet*6] =
4694 static const int hex_tetmap3[nv_tet*6] =
4701 static const int *hex_tetmaps[4] =
4703 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
4706 auto find_min = [](
const int*
a,
int n) {
return std::min_element(a,a+n)-
a; };
4708 for (
int i=0; i<ne; ++i)
4710 const int *v = orig_mesh.
elements[i]->GetVertices();
4714 if (num_subdivisions[orig_geom] == 1)
4726 for (
int itri=0; itri<quad_ntris; ++itri)
4731 for (
int iv=0; iv<nv_tri; ++iv)
4733 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
4741 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
4744 int irot = find_min(vg, nv_prism);
4745 for (
int iv=0; iv<nv_prism; ++iv)
4747 int jv = prism_rot[iv + irot*nv_prism];
4752 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
4753 int j = find_min(q, nv_quad);
4754 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
4755 for (
int itet=0; itet<prism_ntets; ++itet)
4760 for (
int iv=0; iv<nv_tet; ++iv)
4762 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
4770 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
4774 int irot = find_min(vg, nv_hex);
4775 for (
int iv=0; iv<nv_hex; ++iv)
4777 int jv = hex_rot[iv + irot*nv_hex];
4787 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
4788 j = find_min(q, nv_quad);
4789 if (j == 0 || j == 2) { bitmask += 4; }
4791 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
4792 j = find_min(q, nv_quad);
4793 if (j == 1 || j == 3) { bitmask += 2; }
4795 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
4796 j = find_min(q, nv_quad);
4797 if (j == 0 || j == 2) { bitmask += 1; }
4800 int nrot = num_rot[bitmask];
4801 for (
int k=0; k<nrot; ++k)
4815 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
4816 int ntets = (ndiags == 0) ? 5 : 6;
4817 const int *tetmap = hex_tetmaps[ndiags];
4818 for (
int itet=0; itet<ntets; ++itet)
4823 for (
int iv=0; iv<nv_tet; ++iv)
4825 v2[iv] = vg[tetmap[itet + iv*ntets]];
4834 for (
int i=0; i<nbe; ++i)
4836 const int *v = orig_mesh.
boundary[i]->GetVertices();
4839 if (num_subdivisions[orig_geom] == 1)
4849 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
4851 int iv_min = find_min(vg, nv_quad);
4852 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
4853 for (
int itri=0; itri<quad_ntris; ++itri)
4858 for (
int iv=0; iv<nv_tri; ++iv)
4860 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
4867 MFEM_ABORT(
"Unreachable");
4881 Mesh periodic_mesh(orig_mesh,
true);
4887 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
4892 for (
int j = 0; j < nv; j++)
4898 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
4903 for (
int j = 0; j < nv; j++)
4910 return periodic_mesh;
4914 const std::vector<Vector> &translations,
double tol)
const
4918 Vector coord(sdim), at(sdim), dx(sdim);
4919 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
4920 xMax = xMin = xDiff = 0.0;
4924 for (
int be = 0; be <
GetNBE(); be++)
4929 for (
int i = 0; i < dofs.
Size(); i++)
4931 bdr_v.insert(dofs[i]);
4934 for (
int j = 0; j < sdim; j++)
4936 xMax[j] = max(xMax[j], coord[j]);
4937 xMin[j] = min(xMin[j], coord[j]);
4941 add(xMax, -1.0, xMin, xDiff);
4942 double dia = xDiff.
Norml2();
4950 map<int, int> replica2primary;
4952 map<int, set<int>> primary2replicas;
4956 for (
int v : bdr_v) { primary2replicas[v]; }
4960 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
4962 if (r ==
p) {
return; }
4963 primary2replicas[
p].insert(r);
4964 replica2primary[r] =
p;
4965 for (
int s : primary2replicas[r])
4967 primary2replicas[
p].insert(
s);
4968 replica2primary[
s] =
p;
4970 primary2replicas.erase(r);
4973 for (
unsigned int i = 0; i < translations.size(); i++)
4975 for (
int vi : bdr_v)
4978 add(coord, translations[i], at);
4980 for (
int vj : bdr_v)
4983 add(at, -1.0, coord, dx);
4984 if (dx.
Norml2() > dia*tol) {
continue; }
4989 bool pi = primary2replicas.find(vi) != primary2replicas.end();
4990 bool pj = primary2replicas.find(vj) != primary2replicas.end();
4996 make_replica(vj, vi);
5001 int owner_of_vj = replica2primary[vj];
5003 make_replica(vi, owner_of_vj);
5009 int owner_of_vi = replica2primary[vi];
5010 make_replica(vj, owner_of_vi);
5017 int owner_of_vi = replica2primary[vi];
5018 int owner_of_vj = replica2primary[vj];
5019 make_replica(owner_of_vj, owner_of_vi);
5026 std::vector<int> v2v(
GetNV());
5027 for (
size_t i = 0; i < v2v.size(); i++)
5031 for (
auto &&r2p : replica2primary)
5033 v2v[r2p.first] = r2p.second;
5042 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5047 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5064 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
5069 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
5099 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
5123 for (
int i = 0; i <
elements.Size(); i++)
5133 for (
int i = 0; i <
boundary.Size(); i++)
5150 for (
int i = 0; i < vd; i++)
5204 boundary.SetSize(NumOfBdrElements);
5215 edge_to_knot.
SetSize(NumOfEdges);
5219 input >> edge_to_knot[j] >> v[0] >> v[1];
5222 edge_to_knot[j] = -1 - edge_to_knot[j];
5240 for (
int d = 0; d < v.
Size(); d++)
5248 for (d = 0; d < p.
Size(); d++)
5252 for ( ; d < v.
Size(); d++)
5284 if (dynamic_cast<const H1_FECollection*>(fec)
5285 || dynamic_cast<const L2_FECollection*>(fec))
5314 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
5333 MFEM_ASSERT(nodes != NULL,
"");
5349 case 1:
return GetNV();
5385 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5386 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
5391 int i, j, k, wo = 0, fo = 0;
5400 int *vi =
elements[i]->GetVertices();
5403 for (j = 0; j < 3; j++)
5407 for (j = 0; j < 2; j++)
5408 for (k = 0; k < 2; k++)
5410 J(j, k) = v[j+1][k] - v[0][k];
5431 MFEM_ABORT(
"Invalid 2D element type \""
5448 int *vi =
elements[i]->GetVertices();
5454 for (j = 0; j < 4; j++)
5458 for (j = 0; j < 3; j++)
5459 for (k = 0; k < 3; k++)
5461 J(j, k) = v[j+1][k] - v[0][k];
5520 MFEM_ABORT(
"Invalid 3D element type \""
5526 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5529 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
5530 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
5545 if (test[0] == base[0])
5546 if (test[1] == base[1])
5554 else if (test[0] == base[1])
5555 if (test[1] == base[0])
5564 if (test[1] == base[0])
5575 for (
int j = 0; j < 3; j++)
5576 if (test[aor[j]] != base[j])
5589 for (i = 0; i < 4; i++)
5590 if (test[i] == base[0])
5597 if (test[(i+1)%4] == base[1])
5606 for (
int j = 0; j < 4; j++)
5607 if (test[aor[j]] != base[j])
5609 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
5611 for (
int k = 0; k < 4; k++)
5616 for (
int k = 0; k < 4; k++)
5625 if (test[(i+1)%4] == base[1])
5641 if (test[0] == base[0])
5642 if (test[1] == base[1])
5643 if (test[2] == base[2])
5651 else if (test[2] == base[1])
5652 if (test[3] == base[2])
5661 if (test[1] == base[2])
5669 else if (test[1] == base[0])
5670 if (test[2] == base[1])
5671 if (test[0] == base[2])
5679 else if (test[3] == base[1])
5680 if (test[2] == base[2])
5689 if (test[3] == base[2])
5697 else if (test[2] == base[0])
5698 if (test[3] == base[1])
5699 if (test[0] == base[2])
5707 else if (test[0] == base[1])
5708 if (test[1] == base[2])
5717 if (test[3] == base[2])
5726 if (test[0] == base[1])
5727 if (test[2] == base[2])
5735 else if (test[1] == base[1])
5736 if (test[0] == base[2])
5745 if (test[1] == base[2])
5756 for (
int j = 0; j < 4; j++)
5757 if (test[aor[j]] != base[j])
5782 int *bv =
boundary[i]->GetVertices();
5788 mfem::Swap<int>(bv[0], bv[1]);
5802 if (
faces_info[fi].Elem2No >= 0) {
continue; }
5805 int *bv =
boundary[i]->GetVertices();
5807 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
5808 const int *fv =
faces[fi]->GetVertices();
5825 MFEM_ABORT(
"Invalid 2D boundary element type \""
5826 << bdr_type <<
"\"");
5831 if (orientation % 2 == 0) {
continue; }
5833 if (!fix_it) {
continue; }
5841 mfem::Swap<int>(bv[0], bv[1]);
5845 mfem::Swap<int>(be[1], be[2]);
5851 mfem::Swap<int>(bv[0], bv[2]);
5855 mfem::Swap<int>(be[0], be[1]);
5856 mfem::Swap<int>(be[2], be[3]);
5869 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
5879 MFEM_ASSERT(0 <= dim && dim <=
Dim,
"invalid dim: " << dim);
5890 MFEM_ASSERT(0 <= dim && dim <=
Dim,
"invalid dim: " << dim);
5909 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
5910 "is not generated.");
5913 const int *v =
elements[i]->GetVertices();
5914 const int ne =
elements[i]->GetNEdges();
5916 for (
int j = 0; j < ne; j++)
5918 const int *e =
elements[i]->GetEdgeVertices(j);
5919 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5930 const int *v =
boundary[i]->GetVertices();
5931 cor[0] = (v[0] < v[1]) ? (1) : (-1);
5944 const int *v =
boundary[i]->GetVertices();
5945 const int ne =
boundary[i]->GetNEdges();
5947 for (
int j = 0; j < ne; j++)
5949 const int *e =
boundary[i]->GetEdgeVertices(j);
5950 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5962 const int *v =
faces[i]->GetVertices();
5963 o[0] = (v[0] < v[1]) ? (1) : (-1);
5975 const int *v =
faces[i]->GetVertices();
5976 const int ne =
faces[i]->GetNEdges();
5978 for (
int j = 0; j < ne; j++)
5980 const int *e =
faces[i]->GetEdgeVertices(j);
5981 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
6009 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
6060 for (j = 0; j < nv; j++)
6072 for (j = 0; j < nv; j++)
6119 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
6123 int n = el_faces.
Size();
6126 for (
int j = 0; j < n; j++)
6130 ori[j] =
faces_info[el_faces[j]].Elem1Inf % 64;
6134 MFEM_ASSERT(
faces_info[el_faces[j]].Elem2No == i,
"internal error");
6135 ori[j] =
faces_info[el_faces[j]].Elem2Inf % 64;
6159 MFEM_ABORT(
"invalid geometry");
6167 case 1:
return boundary[i]->GetVertices()[0];
6170 default: MFEM_ABORT(
"invalid dimension!");
6180 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
6183 const int *bv =
boundary[bdr_el]->GetVertices();
6191 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
6218 for (j = 0; j < nv; j++)
6220 pointmat(k, j) =
vertices[v[j]](k);
6235 for (j = 0; j < nv; j++)
6237 pointmat(k, j) =
vertices[v[j]](k);
6249 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
6252 return sqrt(length);
6260 for (
int i = 0; i < elem_array.
Size(); i++)
6265 for (
int i = 0; i < elem_array.
Size(); i++)
6267 const int *v = elem_array[i]->GetVertices();
6268 const int ne = elem_array[i]->GetNEdges();
6269 for (
int j = 0; j < ne; j++)
6271 const int *e = elem_array[i]->GetEdgeVertices(j);
6285 v_to_v.
Push(v[0], v[1]);
6292 const int *v =
elements[i]->GetVertices();
6293 const int ne =
elements[i]->GetNEdges();
6294 for (
int j = 0; j < ne; j++)
6296 const int *e =
elements[i]->GetEdgeVertices(j);
6297 v_to_v.
Push(v[e[0]], v[e[1]]);
6305 int i, NumberOfEdges;
6321 const int *v =
boundary[i]->GetVertices();
6322 be_to_f[i] = v_to_v(v[0], v[1]);
6335 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
6339 return NumberOfEdges;
6430 if (
faces[gf] == NULL)
6440 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6441 "Interior edge found between 2D elements "
6443 <<
" and " << el <<
".");
6444 int *v =
faces[gf]->GetVertices();
6446 if ( v[1] == v0 && v[0] == v1 )
6450 else if ( v[0] == v0 && v[1] == v1 )
6460 MFEM_ABORT(
"internal error");
6466 int v0,
int v1,
int v2)
6468 if (
faces[gf] == NULL)
6478 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6479 "Interior triangular face found connecting elements "
6481 <<
" and " << el <<
".");
6482 int orientation, vv[3] = { v0, v1, v2 };
6489 faces_info[gf].Elem2Inf = 64 * lf + orientation;
6494 int v0,
int v1,
int v2,
int v3)
6506 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6507 "Interior quadrilateral face found connecting elements "
6509 <<
" and " << el <<
".");
6510 int vv[4] = { v0, v1, v2, v3 };
6524 for (i = 0; i <
faces.Size(); i++)
6530 faces.SetSize(nfaces);
6532 for (i = 0; i < nfaces; i++)
6540 const int *v =
elements[i]->GetVertices();
6550 const int ne =
elements[i]->GetNEdges();
6551 for (
int j = 0; j < ne; j++)
6553 const int *e =
elements[i]->GetEdgeVertices(j);
6564 for (
int j = 0; j < 4; j++)
6568 v[fv[0]], v[fv[1]], v[fv[2]]);
6574 for (
int j = 0; j < 2; j++)
6578 v[fv[0]], v[fv[1]], v[fv[2]]);
6580 for (
int j = 2; j < 5; j++)
6584 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6590 for (
int j = 0; j < 1; j++)
6594 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6596 for (
int j = 1; j < 5; j++)
6600 v[fv[0]], v[fv[1]], v[fv[2]]);
6606 for (
int j = 0; j < 6; j++)
6610 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6615 MFEM_ABORT(
"Unexpected type of Element.");
6623 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
6639 for (
int i = 0; i < list.
masters.Size(); i++)
6642 if (master.
index >= nfaces) {
continue; }
6648 MFEM_ASSERT(master_fi.
Elem2No == -1,
"internal error");
6649 MFEM_ASSERT(master_fi.
Elem2Inf == -1,
"internal error");
6653 for (
int i = 0; i < list.
slaves.Size(); i++)
6657 if (slave.
index < 0 ||
6658 slave.
index >= nfaces ||
6687 const int *v =
elements[i]->GetVertices();
6692 for (
int j = 0; j < 4; j++)
6695 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6701 for (
int j = 0; j < 1; j++)
6704 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6706 for (
int j = 1; j < 5; j++)
6709 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6715 for (
int j = 0; j < 2; j++)
6718 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6720 for (
int j = 2; j < 5; j++)
6723 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6731 for (
int j = 0; j < 6; j++)
6734 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6739 MFEM_ABORT(
"Unexpected type of Element.");
6763 for (
int j = 0; j < 4; j++)
6767 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6773 for (
int j = 0; j < 2; j++)
6777 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6779 for (
int j = 2; j < 5; j++)
6783 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6789 for (
int j = 0; j < 1; j++)
6793 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6795 for (
int j = 1; j < 5; j++)
6799 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6807 for (
int j = 0; j < 6; j++)
6811 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6816 MFEM_ABORT(
"Unexpected type of Element.");
6829 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
6834 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
6838 MFEM_ABORT(
"Unexpected type of boundary Element.");
6852 void Rotate3(
int &
a,
int &
b,
int &c)
6884 Table *old_elem_vert = NULL;
6895 int *v =
elements[i]->GetVertices();
6897 Rotate3(v[0], v[1], v[2]);
6900 Rotate3(v[1], v[2], v[3]);
6913 int *v =
boundary[i]->GetVertices();
6915 Rotate3(v[0], v[1], v[2]);
6931 delete old_elem_vert;
6947 if (p[i] < pmin[i]) { pmin[i] = p[i]; }
6948 if (p[i] > pmax[i]) { pmax[i] = p[i]; }
6962 for (
int i =
spaceDim-1; i >= 0; i--)
6964 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
6965 if (idx < 0) { idx = 0; }
6966 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
6967 part = part * nxyz[i] + idx;
6969 partitioning[el] = part;
6972 return partitioning;
6977 #ifdef MFEM_USE_METIS
6979 int print_messages = 1;
6982 int init_flag, fin_flag;
6983 MPI_Initialized(&init_flag);
6984 MPI_Finalized(&fin_flag);
6985 if (init_flag && !fin_flag)
6989 if (rank != 0) { print_messages = 0; }
6993 int i, *partitioning;
7003 partitioning[i] = 0;
7010 partitioning[i] = i;
7016 #ifndef MFEM_USE_METIS_5
7028 bool freedata =
false;
7030 idx_t *mpartitioning;
7033 if (
sizeof(
idx_t) ==
sizeof(int))
7037 mpartitioning = (
idx_t*) partitioning;
7046 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
7047 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
7048 mpartitioning =
new idx_t[n];
7051 #ifndef MFEM_USE_METIS_5
7054 METIS_SetDefaultOptions(options);
7055 options[METIS_OPTION_CONTIG] = 1;
7059 if (part_method >= 0 && part_method <= 2)
7061 for (i = 0; i < n; i++)
7067 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
7073 if (part_method == 0 || part_method == 3)
7075 #ifndef MFEM_USE_METIS_5
7104 " error in METIS_PartGraphRecursive!");
7111 if (part_method == 1 || part_method == 4)
7113 #ifndef MFEM_USE_METIS_5
7142 " error in METIS_PartGraphKway!");
7149 if (part_method == 2 || part_method == 5)
7151 #ifndef MFEM_USE_METIS_5
7164 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
7181 " error in METIS_PartGraphKway!");
7189 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
7193 nparts = (int) mparts;
7194 if (mpartitioning != (
idx_t*)partitioning)
7198 partitioning[k] = mpartitioning[k];
7205 delete[] mpartitioning;
7221 auto count_partition_elements = [&]()
7223 for (i = 0; i < nparts; i++)
7231 psize[partitioning[i]].one++;
7235 for (i = 0; i < nparts; i++)
7237 if (psize[i].one == 0) { empty_parts++; }
7241 count_partition_elements();
7249 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
7250 << empty_parts <<
" empty parts!"
7251 <<
" Applying a simple fix ..." << endl;
7254 SortPairs<int,int>(psize, nparts);
7256 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7263 for (i = nparts-1; i > nparts-1-empty_parts; i--)
7265 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
7271 partitioning[j] = psize[nparts-1-i].two;
7278 count_partition_elements();
7282 return partitioning;
7286 mfem_error(
"Mesh::GeneratePartitioning(...): "
7287 "MFEM was compiled without Metis.");
7301 int num_elem, *i_elem_elem, *j_elem_elem;
7303 num_elem = elem_elem.
Size();
7304 i_elem_elem = elem_elem.
GetI();
7305 j_elem_elem = elem_elem.
GetJ();
7310 int stack_p, stack_top_p, elem;
7314 for (i = 0; i < num_elem; i++)
7316 if (partitioning[i] > num_part)
7318 num_part = partitioning[i];
7325 for (i = 0; i < num_part; i++)
7332 for (elem = 0; elem < num_elem; elem++)
7334 if (component[elem] >= 0)
7339 component[elem] = num_comp[partitioning[elem]]++;
7341 elem_stack[stack_top_p++] = elem;
7343 for ( ; stack_p < stack_top_p; stack_p++)
7345 i = elem_stack[stack_p];
7346 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
7349 if (partitioning[k] == partitioning[i])
7351 if (component[k] < 0)
7353 component[k] = component[i];
7354 elem_stack[stack_top_p++] = k;
7356 else if (component[k] != component[i])
7368 int i, n_empty, n_mcomp;
7376 n_empty = n_mcomp = 0;
7377 for (i = 0; i < num_comp.
Size(); i++)
7378 if (num_comp[i] == 0)
7382 else if (num_comp[i] > 1)
7389 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
7390 <<
"The following subdomains are empty :\n";
7391 for (i = 0; i < num_comp.
Size(); i++)
7392 if (num_comp[i] == 0)
7400 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
7401 <<
"The following subdomains are NOT connected :\n";
7402 for (i = 0; i < num_comp.
Size(); i++)
7403 if (num_comp[i] > 1)
7409 if (n_empty == 0 && n_mcomp == 0)
7410 mfem::out <<
"Mesh::CheckPartitioning(...) : "
7411 "All subdomains are connected." << endl;
7425 const double *a = A.
Data();
7426 const double *b = B.
Data();
7435 c(0) = a[0]*a[3]-a[1]*a[2];
7436 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
7437 c(2) = b[0]*b[3]-b[1]*b[2];
7458 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
7459 a[1] * (a[5] * a[6] - a[3] * a[8]) +
7460 a[2] * (a[3] * a[7] - a[4] * a[6]));
7462 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
7463 b[1] * (a[5] * a[6] - a[3] * a[8]) +
7464 b[2] * (a[3] * a[7] - a[4] * a[6]) +
7466 a[0] * (b[4] * a[8] - b[5] * a[7]) +
7467 a[1] * (b[5] * a[6] - b[3] * a[8]) +
7468 a[2] * (b[3] * a[7] - b[4] * a[6]) +
7470 a[0] * (a[4] * b[8] - a[5] * b[7]) +
7471 a[1] * (a[5] * b[6] - a[3] * b[8]) +
7472 a[2] * (a[3] * b[7] - a[4] * b[6]));
7474 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
7475 a[1] * (b[5] * b[6] - b[3] * b[8]) +
7476 a[2] * (b[3] * b[7] - b[4] * b[6]) +
7478 b[0] * (a[4] * b[8] - a[5] * b[7]) +
7479 b[1] * (a[5] * b[6] - a[3] * b[8]) +
7480 b[2] * (a[3] * b[7] - a[4] * b[6]) +
7482 b[0] * (b[4] * a[8] - b[5] * a[7]) +
7483 b[1] * (b[5] * a[6] - b[3] * a[8]) +
7484 b[2] * (b[3] * a[7] - b[4] * a[6]));
7486 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
7487 b[1] * (b[5] * b[6] - b[3] * b[8]) +
7488 b[2] * (b[3] * b[7] - b[4] * b[6]));
7534 double a = z(2), b = z(1), c = z(0);
7535 double D = b*b-4*a*c;
7542 x(0) = x(1) = -0.5 * b /
a;
7547 x(0) = -(x(1) = fabs(0.5 *
sqrt(D) / a));
7555 t = -0.5 * (b +
sqrt(D));
7559 t = -0.5 * (b -
sqrt(D));
7565 Swap<double>(x(0), x(1));
7573 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
7576 double Q = (a * a - 3 *
b) / 9;
7577 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
7578 double Q3 = Q * Q * Q;
7585 x(0) = x(1) = x(2) = - a / 3;
7589 double sqrtQ =
sqrt(Q);
7593 x(0) = -2 * sqrtQ - a / 3;
7594 x(1) = x(2) = sqrtQ - a / 3;
7598 x(0) = x(1) = - sqrtQ - a / 3;
7599 x(2) = 2 * sqrtQ - a / 3;
7607 double A = -2 *
sqrt(Q);
7609 x0 = A *
cos(theta / 3) - a / 3;
7610 x1 = A *
cos((theta + 2.0 * M_PI) / 3) - a / 3;
7611 x2 = A *
cos((theta - 2.0 * M_PI) / 3) - a / 3;
7616 Swap<double>(x0, x1);
7620 Swap<double>(x1, x2);
7623 Swap<double>(x0, x1);
7636 A = -
pow(
sqrt(R2 - Q3) + R, 1.0/3.0);
7640 A =
pow(
sqrt(R2 - Q3) - R, 1.0/3.0);
7642 x(0) = A + Q / A - a / 3;
7651 const double factor,
const int Dim)
7653 const double c0 = c(0);
7654 c(0) = c0 * (1.0 -
pow(factor, -Dim));
7656 for (
int j = 0; j < nr; j++)
7668 c(0) = c0 * (1.0 -
pow(factor, Dim));
7670 for (
int j = 0; j < nr; j++)
7689 const double factor = 2.0;
7704 for (
int k = 0; k < nv; k++)
7707 V(j, k) = displacements(v[k]+j*nvs);
7737 for (
int j = 0; j < nv; j++)
7763 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7766 vertices[i](j) += displacements(j*nv+i);
7774 for (
int i = 0; i < nv; i++)
7777 vert_coord(j*nv+i) =
vertices[i](j);
7783 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7786 vertices[i](j) = vert_coord(j*nv+i);
7797 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
7816 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
7833 (*Nodes) += displacements;
7845 node_coord = (*Nodes);
7857 (*Nodes) = node_coord;
7886 mfem::Swap<GridFunction*>(
Nodes, nodes);
7903 for (j = 1; j < n; j++)
7937 int quad_counter = 0;
7952 vertices.SetSize(oelem + quad_counter);
7953 new_elements.
SetSize(4 * NumOfElements);
7959 const int attr =
elements[i]->GetAttribute();
7960 int *v =
elements[i]->GetVertices();
7966 for (
int ei = 0; ei < 3; ei++)
7968 for (
int k = 0; k < 2; k++)
7976 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
7978 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
7980 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
7982 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
7986 const int qe = quad_counter;
7990 for (
int ei = 0; ei < 4; ei++)
7992 for (
int k = 0; k < 2; k++)
8000 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
8002 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
8004 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
8006 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
8010 MFEM_ABORT(
"unknown element type: " << el_type);
8020 const int attr =
boundary[i]->GetAttribute();
8021 int *v =
boundary[i]->GetVertices();
8030 static const double A = 0.0, B = 0.5, C = 1.0;
8031 static double tri_children[2*3*4] =
8038 static double quad_children[2*4*4] =
8052 for (
int i = 0; i <
elements.Size(); i++)
8073 if (!
Nodes || update_nodes)
8081 static inline double sqr(
const double &x)
8103 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
8106 int NumOfQuadFaces = 0;
8112 for (
int i = 0; i <
faces.Size(); i++)
8116 f2qf[i] = NumOfQuadFaces;
8123 NumOfQuadFaces =
faces.Size();
8127 int hex_counter = 0;
8130 for (
int i = 0; i <
elements.Size(); i++)
8139 int pyr_counter = 0;
8142 for (
int i = 0; i <
elements.Size(); i++)
8160 DSTable *v_to_v_ptr = v_to_v_p;
8176 std::sort(row_start, J_v2v.
end());
8179 for (
int i = 0; i < J_v2v.
Size(); i++)
8181 e2v[J_v2v[i].two] = i;
8194 it.SetIndex(e2v[it.Index()]);
8204 const int oelem = oface + NumOfQuadFaces;
8209 vertices.SetSize(oelem + hex_counter);
8217 const int attr =
elements[i]->GetAttribute();
8218 int *v =
elements[i]->GetVertices();
8225 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8233 for (
int ei = 0; ei < 6; ei++)
8235 for (
int k = 0; k < 2; k++)
8245 const int rt_algo = 1;
8256 double len_sqr, min_len;
8258 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
8259 sqr(J(1,0)-J(1,1)-J(1,2)) +
8260 sqr(J(2,0)-J(2,1)-J(2,2));
8263 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
8264 sqr(J(1,1)-J(1,0)-J(1,2)) +
8265 sqr(J(2,1)-J(2,0)-J(2,2));
8266 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
8268 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
8269 sqr(J(1,2)-J(1,0)-J(1,1)) +
8270 sqr(J(2,2)-J(2,0)-J(2,1));
8271 if (len_sqr < min_len) { rt = 2; }
8276 double Em_data[18], Js_data[9], Jp_data[9];
8279 double ar1, ar2,
kappa, kappa_min;
8281 for (
int s = 0;
s < 3;
s++)
8283 for (
int t = 0;
t < 3;
t++)
8285 Em(
t,
s) = 0.5*J(
t,
s);
8288 for (
int t = 0;
t < 3;
t++)
8290 Em(
t,3) = 0.5*(J(
t,0)+J(
t,1));
8291 Em(
t,4) = 0.5*(J(
t,0)+J(
t,2));
8292 Em(
t,5) = 0.5*(J(
t,1)+J(
t,2));
8296 for (
int t = 0;
t < 3;
t++)
8298 Js(
t,0) = Em(
t,5)-Em(
t,0);
8299 Js(
t,1) = Em(
t,1)-Em(
t,0);
8300 Js(
t,2) = Em(
t,2)-Em(
t,0);
8304 for (
int t = 0;
t < 3;
t++)
8306 Js(
t,0) = Em(
t,5)-Em(
t,0);
8307 Js(
t,1) = Em(
t,2)-Em(
t,0);
8308 Js(
t,2) = Em(
t,4)-Em(
t,0);
8312 kappa_min = std::max(ar1, ar2);
8316 for (
int t = 0;
t < 3;
t++)
8318 Js(
t,0) = Em(
t,0)-Em(
t,1);
8319 Js(
t,1) = Em(
t,4)-Em(
t,1);
8320 Js(
t,2) = Em(
t,2)-Em(
t,1);
8324 for (
int t = 0;
t < 3;
t++)
8326 Js(
t,0) = Em(
t,2)-Em(
t,1);
8327 Js(
t,1) = Em(
t,4)-Em(
t,1);
8328 Js(
t,2) = Em(
t,5)-Em(
t,1);
8332 kappa = std::max(ar1, ar2);
8333 if (kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
8336 for (
int t = 0;
t < 3;
t++)
8338 Js(
t,0) = Em(
t,0)-Em(
t,2);
8339 Js(
t,1) = Em(
t,1)-Em(
t,2);
8340 Js(
t,2) = Em(
t,3)-Em(
t,2);
8344 for (
int t = 0;
t < 3;
t++)
8346 Js(
t,0) = Em(
t,1)-Em(
t,2);
8347 Js(
t,1) = Em(
t,5)-Em(
t,2);
8348 Js(
t,2) = Em(
t,3)-Em(
t,2);
8352 kappa = std::max(ar1, ar2);
8353 if (kappa < kappa_min) { rt = 2; }
8356 static const int mv_all[3][4][4] =
8358 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
8359 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
8360 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
8362 const int (&mv)[4][4] = mv_all[rt];
8364 #ifndef MFEM_USE_MEMALLOC
8366 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8368 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8370 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8372 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8374 for (
int k = 0; k < 4; k++)
8376 new_elements[j+4+k] =
8377 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8378 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8382 new_elements[j+0] = tet =
TetMemory.Alloc();
8383 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
8385 new_elements[j+1] = tet =
TetMemory.Alloc();
8386 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
8388 new_elements[j+2] = tet =
TetMemory.Alloc();
8389 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
8391 new_elements[j+3] = tet =
TetMemory.Alloc();
8392 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
8394 for (
int k = 0; k < 4; k++)
8396 new_elements[j+4+k] = tet =
TetMemory.Alloc();
8397 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
8398 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
8401 for (
int k = 0; k < 4; k++)
8406 for (
int k = 0; k < 4; k++)
8420 for (
int fi = 2; fi < 5; fi++)
8422 for (
int k = 0; k < 4; k++)
8429 for (
int ei = 0; ei < 9; ei++)
8431 for (
int k = 0; k < 2; k++)
8438 const int qf2 = f2qf[f[2]];
8439 const int qf3 = f2qf[f[3]];
8440 const int qf4 = f2qf[f[4]];
8443 new Wedge(v[0], oedge+e[0], oedge+e[2],
8444 oedge+e[6], oface+qf2, oface+qf4, attr);
8447 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
8448 oface+qf3, oface+qf4, oface+qf2, attr);
8451 new Wedge(oedge+e[0], v[1], oedge+e[1],
8452 oface+qf2, oedge+e[7], oface+qf3, attr);
8455 new Wedge(oedge+e[2], oedge+e[1], v[2],
8456 oface+qf4, oface+qf3, oedge+e[8], attr);
8459 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
8460 v[3], oedge+e[3], oedge+e[5], attr);
8463 new Wedge(oface+qf3, oface+qf4, oface+qf2,
8464 oedge+e[4], oedge+e[5], oedge+e[3], attr);
8467 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
8468 oedge+e[3], v[4], oedge+e[4], attr);
8471 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
8472 oedge+e[5], oedge+e[4], v[5], attr);
8481 for (
int fi = 0; fi < 1; fi++)
8483 for (
int k = 0; k < 4; k++)
8490 for (
int ei = 0; ei < 8; ei++)
8492 for (
int k = 0; k < 2; k++)
8499 const int qf0 = f2qf[f[0]];
8502 new Pyramid(v[0], oedge+e[0], oface+qf0,
8503 oedge+e[3], oedge+e[4], attr);
8506 new Pyramid(oedge+e[0], v[1], oedge+e[1],
8507 oface+qf0, oedge+e[5], attr);
8510 new Pyramid(oface+qf0, oedge+e[1], v[2],
8511 oedge+e[2], oedge+e[6], attr);
8514 new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
8515 v[3], oedge+e[7], attr);
8518 new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
8519 oedge+e[7], v[4], attr);
8522 new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
8523 oedge+e[4], oface+qf0, attr);
8526 new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
8530 new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
8534 new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
8538 new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
8546 const int he = hex_counter;
8551 if (f2qf.
Size() == 0)
8557 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[f[k]]; }
8563 for (
int fi = 0; fi < 6; fi++)
8565 for (
int k = 0; k < 4; k++)
8572 for (
int ei = 0; ei < 12; ei++)
8574 for (
int k = 0; k < 2; k++)
8582 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
8583 oedge+e[3], oedge+e[8], oface+qf[1],
8584 oelem+he, oface+qf[4], attr);
8587 oface+qf[0], oface+qf[1], oedge+e[9],
8588 oface+qf[2], oelem+he, attr);
8590 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
8591 oedge+e[2], oelem+he, oface+qf[2],
8592 oedge+e[10], oface+qf[3], attr);
8594 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
8595 v[3], oface+qf[4], oelem+he,
8596 oface+qf[3], oedge+e[11], attr);
8598 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
8599 oface+qf[4], v[4], oedge+e[4],
8600 oface+qf[5], oedge+e[7], attr);
8602 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
8603 oelem+he, oedge+e[4], v[5],
8604 oedge+e[5], oface+qf[5], attr);
8606 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
8607 oface+qf[3], oface+qf[5], oedge+e[5],
8608 v[6], oedge+e[6], attr);
8610 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
8611 oedge+e[11], oedge+e[7], oface+qf[5],
8612 oedge+e[6], v[7], attr);
8617 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
8629 const int attr =
boundary[i]->GetAttribute();
8630 int *v =
boundary[i]->GetVertices();
8637 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8644 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8646 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8648 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8650 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8658 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
8660 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
8662 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
8664 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
8668 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
8674 static const double A = 0.0, B = 0.5, C = 1.0, D = -1.0;
8675 static double tet_children[3*4*16] =
8677 A,A,A, B,A,A, A,B,A, A,A,B,
8678 B,A,A, C,A,A, B,B,A, B,A,B,
8679 A,B,A, B,B,A, A,C,A, A,B,B,
8680 A,A,B, B,A,B, A,B,B, A,A,C,
8685 B,A,A, A,B,B, A,B,A, A,A,B,
8686 B,A,A, A,B,B, A,A,B, B,A,B,
8687 B,A,A, A,B,B, B,A,B, B,B,A,
8688 B,A,A, A,B,B, B,B,A, A,B,A,
8690 A,B,A, B,A,A, B,A,B, A,A,B,
8691 A,B,A, A,A,B, B,A,B, A,B,B,
8692 A,B,A, A,B,B, B,A,B, B,B,A,
8693 A,B,A, B,B,A, B,A,B, B,A,A,
8695 A,A,B, B,A,A, A,B,A, B,B,A,
8696 A,A,B, A,B,A, A,B,B, B,B,A,
8697 A,A,B, A,B,B, B,A,B, B,B,A,
8698 A,A,B, B,A,B, B,A,A, B,B,A
8700 static double pyr_children[3*5*10] =
8702 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
8703 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
8704 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
8705 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
8706 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
8707 A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
8708 B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
8709 C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
8710 B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
8711 A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
8713 static double pri_children[3*6*8] =
8715 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
8716 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
8717 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
8718 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
8719 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
8720 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
8721 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
8722 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
8724 static double hex_children[3*8*8] =
8726 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
8727 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
8728 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
8729 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
8730 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
8731 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
8732 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
8733 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
8745 for (
int i = 0; i <
elements.Size(); i++)
8756 NumOfElements = 8 * NumOfElements + 2 * pyr_counter;
8776 int i, j, ind, nedges;
8783 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
8797 for (j = 0; j < marked_el.
Size(); j++)
8802 int new_v = cnv + j, new_e = cne + j;
8811 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
8813 UseExternalData(seg_children, 1, 2, 3);
8826 int *edge1 =
new int[nedges];
8827 int *edge2 =
new int[nedges];
8828 int *middle =
new int[nedges];
8830 for (i = 0; i < nedges; i++)
8832 edge1[i] = edge2[i] = middle[i] = -1;
8838 for (j = 1; j < v.
Size(); j++)
8840 ind = v_to_v(v[j-1], v[j]);
8841 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8843 ind = v_to_v(v[0], v[v.
Size()-1]);
8844 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8848 for (i = 0; i < marked_el.
Size(); i++)
8854 int need_refinement;
8857 need_refinement = 0;
8858 for (i = 0; i < nedges; i++)
8860 if (middle[i] != -1 && edge1[i] != -1)
8862 need_refinement = 1;
8867 while (need_refinement == 1);
8870 int v1[2], v2[2], bisect, temp;
8872 for (i = 0; i < temp; i++)
8875 bisect = v_to_v(v[0], v[1]);
8876 if (middle[bisect] != -1)
8880 v1[0] = v[0]; v1[1] = middle[bisect];
8881 v2[0] = middle[bisect]; v2[1] = v[1];
8887 mfem_error(
"Only bisection of segment is implemented"
8911 MFEM_VERIFY(
GetNE() == 0 ||
8913 "tetrahedral mesh is not marked for refinement:"
8914 " call Finalize(true)");
8921 for (i = 0; i < marked_el.
Size(); i++)
8927 for (i = 0; i < marked_el.
Size(); i++)
8936 for (i = 0; i < marked_el.
Size(); i++)
8953 int need_refinement;
8958 need_refinement = 0;
8966 if (
elements[i]->NeedRefinement(v_to_v))
8968 need_refinement = 1;
8973 while (need_refinement == 1);
8980 need_refinement = 0;
8982 if (
boundary[i]->NeedRefinement(v_to_v))
8984 need_refinement = 1;
8988 while (need_refinement == 1);
9019 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
9020 "not supported. Project the NURBS to Nodes first.");
9030 if (!refinements.
Size())
9051 Swap(*mesh2,
false);
9067 const int *fine,
int nfine,
int op)
9070 for (
int i = 0; i < nfine; i++)
9072 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
9074 double err_fine = elem_error[fine[i]];
9077 case 0: error = std::min(error, err_fine);
break;
9078 case 1: error += err_fine;
break;
9079 case 2: error = std::max(error, err_fine);
break;
9086 double threshold,
int nc_limit,
int op)
9088 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
9089 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
9090 "Project the NURBS to Nodes first.");
9103 for (
int i = 0; i < dt.
Size(); i++)
9105 if (nc_limit > 0 && !level_ok[i]) {
continue; }
9110 if (error < threshold) { derefs.
Append(i); }
9113 if (!derefs.
Size()) {
return false; }
9120 Swap(*mesh2,
false);
9134 int nc_limit,
int op)
9144 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
9151 int nc_limit,
int op)
9154 for (
int i = 0; i < tmp.Size(); i++)
9156 tmp[i] = elem_error(i);
9241 #ifdef MFEM_USE_MEMALLOC
9268 for (
int i = 0; i < elem_array.
Size(); i++)
9270 if (elem_array[i]->GetGeometryType() == geom)
9275 elem_vtx.
SetSize(nv*num_elems);
9279 for (
int i = 0; i < elem_array.
Size(); i++)
9285 elem_vtx.
Append(loc_vtx);
9293 for (
int i = 0; i < nelem; i++) { list[i] = i; }
9309 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
9321 default: MFEM_ABORT(
"internal error");
9327 int nonconforming,
int nc_limit)
9337 else if (nonconforming < 0)
9358 for (
int i = 0; i < refinements.
Size(); i++)
9360 el_to_refine[i] = refinements[i].index;
9364 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
9365 if (rt == 1 || rt == 2 || rt == 4)
9369 else if (rt == 3 || rt == 5 || rt == 6)
9387 for (
int i = 0; i < el_to_refine.
Size(); i++)
9389 refinements[i] =
Refinement(el_to_refine[i]);
9396 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
9397 "Please project the NURBS to Nodes first, with SetCurvature().");
9400 MFEM_VERIFY(
ncmesh != NULL || dynamic_cast<const ParMesh*>(
this) == NULL,
9401 "Sorry, converting a conforming ParMesh to an NC mesh is "
9409 (simplices_nonconforming && (
meshgen & 0x1)) )
9422 for (
int i = 0; i <
GetNE(); i++)
9424 if ((
double) rand() / RAND_MAX <
prob)
9429 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
9441 for (
int i = 0; i <
GetNE(); i++)
9444 bool refine =
false;
9445 for (
int j = 0; j < v.
Size(); j++)
9450 double d = vert(l) -
vertices[v[j]](l);
9453 if (dist <= eps*eps) { refine =
true;
break; }
9464 int nonconforming,
int nc_limit)
9466 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
9468 for (
int i = 0; i <
GetNE(); i++)
9470 if (elem_error[i] > threshold)
9484 int nonconforming,
int nc_limit)
9488 return RefineByError(tmp, threshold, nonconforming, nc_limit);
9493 int *edge1,
int *edge2,
int *middle)
9496 int v[2][4], v_new, bisect,
t;
9508 bisect = v_to_v(vert[0], vert[1]);
9509 MFEM_ASSERT(bisect >= 0,
"");
9511 if (middle[bisect] == -1)
9522 if (edge1[bisect] == i)
9524 edge1[bisect] = edge2[bisect];
9527 middle[bisect] = v_new;
9531 v_new = middle[bisect];
9539 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9540 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9561 bisect = v_to_v(v[1][0], v[1][1]);
9562 MFEM_ASSERT(bisect >= 0,
"");
9564 if (edge1[bisect] == i)
9568 else if (edge2[bisect] == i)
9577 MFEM_ABORT(
"Bisection for now works only for triangles.");
9584 int v[2][4], v_new, bisect,
t;
9591 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
9595 "TETRAHEDRON element is not marked for refinement.");
9600 bisect = v_to_v.
FindId(vert[0], vert[1]);
9604 for (j = 0; j < 3; j++)
9621 new_redges[0][0] = 2;
9622 new_redges[0][1] = 1;
9623 new_redges[1][0] = 2;
9624 new_redges[1][1] = 1;
9625 int tr1 = -1, tr2 = -1;
9626 switch (old_redges[0])
9629 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
9634 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
9638 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
9641 switch (old_redges[1])
9644 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
9649 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
9653 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
9660 #ifdef MFEM_USE_MEMALLOC
9696 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
9703 int v[2][3], v_new, bisect,
t;
9714 bisect = v_to_v.
FindId(vert[0], vert[1]);
9715 MFEM_ASSERT(bisect >= 0,
"");
9717 MFEM_ASSERT(v_new != -1,
"");
9721 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9722 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9732 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
9738 int *edge1,
int *edge2,
int *middle)
9741 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
9750 bisect[0] = v_to_v(v[0],v[1]);
9751 bisect[1] = v_to_v(v[1],v[2]);
9752 bisect[2] = v_to_v(v[0],v[2]);
9753 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
9755 for (j = 0; j < 3; j++)
9757 if (middle[bisect[j]] == -1)
9768 if (edge1[bisect[j]] == i)
9770 edge1[bisect[j]] = edge2[bisect[j]];
9773 middle[bisect[j]] = v_new[j];
9777 v_new[j] = middle[bisect[j]];
9780 edge1[bisect[j]] = -1;
9786 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
9787 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
9788 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
9789 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
9804 tri2->ResetTransform(code);
9805 tri3->ResetTransform(code);
9809 tri2->PushTransform(1);
9810 tri3->PushTransform(2);
9823 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
9859 for (
int i = 0; i < elem_geoms.
Size(); i++)
9867 std::map<unsigned, int> mat_no;
9871 for (
int j = 0; j <
elements.Size(); j++)
9874 unsigned code =
elements[j]->GetTransform();
9877 int &matrix = mat_no[code];
9878 if (!matrix) { matrix = mat_no.size(); }
9888 std::map<unsigned, int>::iterator it;
9889 for (it = mat_no.begin(); it != mat_no.end(); ++it)
9903 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
9904 " geom = " << geom);
9914 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
9923 os <<
"areamesh2\n\n";
9927 os <<
"curved_areamesh2\n\n";
9937 for (j = 0; j < v.
Size(); j++)
9939 os <<
' ' << v[j] + 1;
9951 for (j = 0; j < v.
Size(); j++)
9953 os <<
' ' << v[j] + 1;
9965 for (j = 1; j <
Dim; j++)
9982 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
9990 os <<
"NETGEN_Neutral_Format\n";
9995 for (j = 0; j <
Dim; j++)
10008 os <<
elements[i]->GetAttribute();
10009 for (j = 0; j < nv; j++)
10011 os <<
' ' << ind[j]+1;
10022 os <<
boundary[i]->GetAttribute();
10023 for (j = 0; j < nv; j++)
10025 os <<
' ' << ind[j]+1;
10037 <<
" 0 0 0 0 0 0 0\n"
10038 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
10040 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
10041 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
10045 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
10051 os << i+1 <<
' ' <<
elements[i]->GetAttribute();
10052 for (j = 0; j < nv; j++)
10054 os <<
' ' << ind[j]+1;
10063 os <<
boundary[i]->GetAttribute();
10064 for (j = 0; j < nv; j++)
10066 os <<
' ' << ind[j]+1;
10068 os <<
" 1.0 1.0 1.0 1.0\n";
10101 os <<
"\n# mesh curvature GridFunction";
10106 os <<
"\nmfem_mesh_end" << endl;
10111 os << (section_delimiter.empty()
10112 ?
"MFEM mesh v1.0\n" :
"MFEM mesh v1.2\n");
10116 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10121 "# TETRAHEDRON = 4\n"
10127 os <<
"\ndimension\n" <<
Dim;
10162 if (!section_delimiter.empty())
10164 os << section_delimiter << endl;
10173 os <<
"MFEM NURBS mesh v1.0\n";
10177 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10183 os <<
"\ndimension\n" <<
Dim
10200 int ki = e_to_k[i];
10205 os << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
10212 ofstream ofs(fname);
10213 ofs.precision(precision);
10217 #ifdef MFEM_USE_ADIOS2
10227 "# vtk DataFile Version 3.0\n"
10228 "Generated by MFEM\n"
10230 "DATASET UNSTRUCTURED_GRID\n";
10243 for ( ; j < 3; j++)
10259 os << (*Nodes)(vdofs[0]);
10263 os <<
' ' << (*Nodes)(vdofs[j]);
10265 for ( ; j < 3; j++)
10279 size +=
elements[i]->GetNVertices() + 1;
10281 os <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
10284 const int *v =
elements[i]->GetVertices();
10285 const int nv =
elements[i]->GetNVertices();
10289 for (
int j = 0; j < nv; j++)
10291 os <<
' ' << v[perm ? perm[j] : j];
10304 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
10305 "Point meshes should have a single dof per element");
10306 size += dofs.
Size() + 1;
10308 os <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
10311 if (!strcmp(fec_name,
"Linear") ||
10312 !strcmp(fec_name,
"H1_0D_P1") ||
10313 !strcmp(fec_name,
"H1_1D_P1") ||
10314 !strcmp(fec_name,
"H1_2D_P1") ||
10315 !strcmp(fec_name,
"H1_3D_P1"))
10319 else if (!strcmp(fec_name,
"Quadratic") ||
10320 !strcmp(fec_name,
"H1_1D_P2") ||
10321 !strcmp(fec_name,
"H1_2D_P2") ||
10322 !strcmp(fec_name,
"H1_3D_P2"))
10328 mfem::err <<
"Mesh::PrintVTK : can not save '"
10329 << fec_name <<
"' elements!" << endl;
10338 for (
int j = 0; j < dofs.
Size(); j++)
10340 os <<
' ' << dofs[j];
10343 else if (order == 2)
10345 const int *vtk_mfem;
10346 switch (
elements[i]->GetGeometryType())
10360 for (
int j = 0; j < dofs.
Size(); j++)
10362 os <<
' ' << dofs[vtk_mfem[j]];
10372 int vtk_cell_type = 5;
10376 os << vtk_cell_type <<
'\n';
10380 os <<
"CELL_DATA " << NumOfElements <<
'\n'
10381 <<
"SCALARS material int\n"
10382 <<
"LOOKUP_TABLE default\n";
10385 os <<
elements[i]->GetAttribute() <<
'\n';
10392 bool high_order_output,
10393 int compression_level,
10396 int ref = (high_order_output &&
Nodes)
10399 fname = fname +
".vtu";
10401 os <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
10402 if (compression_level != 0)
10404 os <<
" compressor=\"vtkZLibDataCompressor\"";
10407 os <<
"<UnstructuredGrid>\n";
10408 PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
10409 os <<
"</Piece>\n";
10410 os <<
"</UnstructuredGrid>\n";
10411 os <<
"</VTKFile>" << std::endl;
10418 bool high_order_output,
10419 int compression_level)
10421 PrintVTU(fname, format, high_order_output, compression_level,
true);
10425 bool high_order_output,
int compression_level,
10433 std::vector<char> buf;
10435 auto get_geom = [&](
int i)
10443 int np = 0, nc_ref = 0, size = 0;
10444 for (
int i = 0; i < ne; i++)
10454 os <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
10455 << (high_order_output ? ne : nc_ref) <<
"\">\n";
10458 os <<
"<Points>\n";
10459 os <<
"<DataArray type=\"" << type_str
10460 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
10461 for (
int i = 0; i < ne; i++)
10474 for (
int j = 0; j < pmat.
Width(); j++)
10500 os <<
"</DataArray>" << std::endl;
10501 os <<
"</Points>" << std::endl;
10503 os <<
"<Cells>" << std::endl;
10504 os <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
10505 << fmt_str <<
"\">" << std::endl;
10507 std::vector<int> offset;
10510 if (high_order_output)
10513 for (
int iel = 0; iel < ne; iel++)
10517 int nnodes = local_connectivity.
Size();
10518 for (
int i=0; i<nnodes; ++i)
10525 offset.push_back(np);
10531 for (
int i = 0; i < ne; i++)
10537 for (
int j = 0; j < RG.
Size(); )
10540 offset.push_back(coff);
10542 for (
int k = 0; k < nv; k++, j++)
10556 os <<
"</DataArray>" << std::endl;
10558 os <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
10559 << fmt_str <<
"\">" << std::endl;
10561 for (
size_t ii=0; ii<offset.size(); ii++)
10569 os <<
"</DataArray>" << std::endl;
10570 os <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
10571 << fmt_str <<
"\">" << std::endl;
10573 const int *vtk_geom_map =
10575 for (
int i = 0; i < ne; i++)
10578 uint8_t vtk_cell_type = 5;
10580 vtk_cell_type = vtk_geom_map[geom];
10582 if (high_order_output)
10591 for (
int j = 0; j < RG.
Size(); j += nv)
10601 os <<
"</DataArray>" << std::endl;
10602 os <<
"</Cells>" << std::endl;
10604 os <<
"<CellData Scalars=\"attribute\">" << std::endl;
10605 os <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
10606 << fmt_str <<
"\">" << std::endl;
10607 for (
int i = 0; i < ne; i++)
10610 if (high_order_output)
10629 os <<
"</DataArray>" << std::endl;
10630 os <<
"</CellData>" << std::endl;
10641 "# vtk DataFile Version 3.0\n"
10642 "Generated by MFEM\n"
10644 "DATASET UNSTRUCTURED_GRID\n";
10649 os <<
"FIELD FieldData 1\n"
10659 np = nc = size = 0;
10660 for (
int i = 0; i <
GetNE(); i++)
10669 os <<
"POINTS " << np <<
" double\n";
10671 for (
int i = 0; i <
GetNE(); i++)
10678 for (
int j = 0; j < pmat.
Width(); j++)
10680 os << pmat(0, j) <<
' ';
10683 os << pmat(1, j) <<
' ';
10695 os << 0.0 <<
' ' << 0.0;
10702 os <<
"CELLS " << nc <<
' ' << size <<
'\n';
10704 for (
int i = 0; i <
GetNE(); i++)
10711 for (
int j = 0; j < RG.
Size(); )
10714 for (
int k = 0; k < nv; k++, j++)
10716 os <<
' ' << np + RG[j];
10722 os <<
"CELL_TYPES " << nc <<
'\n';
10723 for (
int i = 0; i <
GetNE(); i++)
10731 for (
int j = 0; j < RG.
Size(); j += nv)
10733 os << vtk_cell_type <<
'\n';
10737 os <<
"CELL_DATA " << nc <<
'\n'
10738 <<
"SCALARS material int\n"
10739 <<
"LOOKUP_TABLE default\n";
10740 for (
int i = 0; i <
GetNE(); i++)
10748 os << attr <<
'\n';
10755 srand((
unsigned)time(0));
10756 double a = double(rand()) / (double(RAND_MAX) + 1.);
10757 int el0 = (int)floor(a *
GetNE());
10759 os <<
"SCALARS element_coloring int\n"
10760 <<
"LOOKUP_TABLE default\n";
10761 for (
int i = 0; i <
GetNE(); i++)
10768 os << coloring[i] + 1 <<
'\n';
10774 os <<
"POINT_DATA " << np <<
'\n' << flush;
10779 int delete_el_to_el = (
el_to_el) ? (0) : (1);
10781 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
10784 const int *i_el_el = el_el.
GetI();
10785 const int *j_el_el = el_el.
GetJ();
10790 stack_p = stack_top_p = 0;
10791 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
10793 if (colors[el] != -2)
10799 el_stack[stack_top_p++] = el;
10801 for ( ; stack_p < stack_top_p; stack_p++)
10803 int i = el_stack[stack_p];
10804 int num_nb = i_el_el[i+1] - i_el_el[i];
10805 if (max_num_col < num_nb + 1)
10807 max_num_col = num_nb + 1;
10809 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10811 int k = j_el_el[j];
10812 if (colors[k] == -2)
10815 el_stack[stack_top_p++] = k;
10823 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
10825 int i = el_stack[stack_p], col;
10827 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10829 col = colors[j_el_el[j]];
10832 col_marker[col] = 1;
10836 for (col = 0; col < max_num_col; col++)
10837 if (col_marker[col] == 0)
10845 if (delete_el_to_el)
10853 int elem_attr)
const
10855 if (
Dim != 3 &&
Dim != 2) {
return; }
10857 int i, j, k, l, nv, nbe, *v;
10859 os <<
"MFEM mesh v1.0\n";
10863 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10868 "# TETRAHEDRON = 4\n"
10873 os <<
"\ndimension\n" <<
Dim
10878 <<
' ' <<
elements[i]->GetGeometryType();
10881 for (j = 0; j < nv; j++)
10893 l = partitioning[l];
10908 os <<
"\nboundary\n" << nbe <<
'\n';
10914 l = partitioning[l];
10917 nv =
faces[i]->GetNVertices();
10918 v =
faces[i]->GetVertices();
10919 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
10920 for (j = 0; j < nv; j++)
10927 os << l+1 <<
' ' <<
faces[i]->GetGeometryType();
10928 for (j = nv-1; j >= 0; j--)
10939 nv =
faces[i]->GetNVertices();
10940 v =
faces[i]->GetVertices();
10941 os << k+1 <<
' ' <<
faces[i]->GetGeometryType();
10942 for (j = 0; j < nv; j++)
10973 int interior_faces)
10975 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
10976 if (
Dim != 3 &&
Dim != 2) {
return; }
10985 int nv =
elements[i]->GetNVertices();
10986 const int *ind =
elements[i]->GetVertices();
10987 for (
int j = 0; j < nv; j++)
10993 int *voff =
new int[NumOfVertices+1];
10997 voff[i] = vcount[i-1] + voff[i-1];
11003 vown[i] =
new int[vcount[i]];
11015 int nv =
elements[i]->GetNVertices();
11016 const int *ind =
elements[i]->GetVertices();
11017 for (
int j = 0; j < nv; j++)
11020 vown[ind[j]][vcount[ind[j]]] = i;
11026 vcount[i] = voff[i+1] - voff[i];
11030 for (
int i = 0; i < edge_el.
Size(); i++)
11032 const int *el = edge_el.
GetRow(i);
11035 int k = partitioning[el[0]];
11036 int l = partitioning[el[1]];
11037 if (interior_faces || k != l)
11049 os <<
"areamesh2\n\n" << nbe <<
'\n';
11051 for (
int i = 0; i < edge_el.
Size(); i++)
11053 const int *el = edge_el.
GetRow(i);
11056 int k = partitioning[el[0]];
11057 int l = partitioning[el[1]];
11058 if (interior_faces || k != l)
11063 for (
int j = 0; j < 2; j++)
11064 for (
int s = 0;
s < vcount[ev[j]];
s++)
11065 if (vown[ev[j]][
s] == el[0])
11067 os <<
' ' << voff[ev[j]]+
s+1;
11071 for (
int j = 1; j >= 0; j--)
11072 for (
int s = 0;
s < vcount[ev[j]];
s++)
11073 if (vown[ev[j]][
s] == el[1])
11075 os <<
' ' << voff[ev[j]]+
s+1;
11082 int k = partitioning[el[0]];
11086 for (
int j = 0; j < 2; j++)
11087 for (
int s = 0;
s < vcount[ev[j]];
s++)
11088 if (vown[ev[j]][
s] == el[0])
11090 os <<
' ' << voff[ev[j]]+
s+1;
11097 os << NumOfElements <<
'\n';
11100 int nv =
elements[i]->GetNVertices();
11101 const int *ind =
elements[i]->GetVertices();
11102 os << partitioning[i]+1 <<
' ';
11104 for (
int j = 0; j < nv; j++)
11106 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11107 vown[ind[j]][vcount[ind[j]]] = i;
11114 vcount[i] = voff[i+1] - voff[i];
11120 for (
int k = 0; k < vcount[i]; k++)
11122 for (
int j = 0; j <
Dim; j++)
11132 os <<
"NETGEN_Neutral_Format\n";
11136 for (
int k = 0; k < vcount[i]; k++)
11138 for (
int j = 0; j <
Dim; j++)
11146 os << NumOfElements <<
'\n';
11149 int nv =
elements[i]->GetNVertices();
11150 const int *ind =
elements[i]->GetVertices();
11151 os << partitioning[i]+1;
11152 for (
int j = 0; j < nv; j++)
11154 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11155 vown[ind[j]][vcount[ind[j]]] = i;
11162 vcount[i] = voff[i+1] - voff[i];
11172 int k = partitioning[
faces_info[i].Elem1No];
11173 l = partitioning[l];
11174 if (interior_faces || k != l)
11191 int k = partitioning[
faces_info[i].Elem1No];
11192 l = partitioning[l];
11193 if (interior_faces || k != l)
11195 int nv =
faces[i]->GetNVertices();
11196 const int *ind =
faces[i]->GetVertices();
11198 for (
int j = 0; j < nv; j++)
11199 for (
int s = 0;
s < vcount[ind[j]];
s++)
11200 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11202 os <<
' ' << voff[ind[j]]+
s+1;
11206 for (
int j = nv-1; j >= 0; j--)
11207 for (
int s = 0;
s < vcount[ind[j]];
s++)
11208 if (vown[ind[j]][
s] == faces_info[i].Elem2No)
11210 os <<
' ' << voff[ind[j]]+
s+1;
11217 int k = partitioning[
faces_info[i].Elem1No];
11218 int nv =
faces[i]->GetNVertices();
11219 const int *ind =
faces[i]->GetVertices();
11221 for (
int j = 0; j < nv; j++)
11222 for (
int s = 0;
s < vcount[ind[j]];
s++)
11223 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11225 os <<
' ' << voff[ind[j]]+
s+1;
11241 int k = partitioning[
faces_info[i].Elem1No];
11242 l = partitioning[l];
11243 if (interior_faces || k != l)
11256 <<
" 0 0 0 0 0 0 0\n"
11257 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
11258 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
11259 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
11260 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
11263 for (
int k = 0; k < vcount[i]; k++)
11264 os << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
11269 int nv =
elements[i]->GetNVertices();
11270 const int *ind =
elements[i]->GetVertices();
11271 os << i+1 <<
' ' << partitioning[i]+1;
11272 for (
int j = 0; j < nv; j++)
11274 os <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
11275 vown[ind[j]][vcount[ind[j]]] = i;
11282 vcount[i] = voff[i+1] - voff[i];
11291 int k = partitioning[
faces_info[i].Elem1No];
11292 l = partitioning[l];
11293 if (interior_faces || k != l)
11295 int nv =
faces[i]->GetNVertices();
11296 const int *ind =
faces[i]->GetVertices();
11298 for (
int j = 0; j < nv; j++)
11299 for (
int s = 0;
s < vcount[ind[j]];
s++)
11300 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11302 os <<
' ' << voff[ind[j]]+
s+1;
11304 os <<
" 1.0 1.0 1.0 1.0\n";
11306 for (
int j = nv-1; j >= 0; j--)
11307 for (
int s = 0;
s < vcount[ind[j]];
s++)
11308 if (vown[ind[j]][
s] == faces_info[i].Elem2No)
11310 os <<
' ' << voff[ind[j]]+
s+1;
11312 os <<
" 1.0 1.0 1.0 1.0\n";
11317 int k = partitioning[
faces_info[i].Elem1No];
11318 int nv =
faces[i]->GetNVertices();
11319 const int *ind =
faces[i]->GetVertices();
11321 for (
int j = 0; j < nv; j++)
11322 for (
int s = 0;
s < vcount[ind[j]];
s++)
11323 if (vown[ind[j]][
s] == faces_info[i].Elem1No)
11325 os <<
' ' << voff[ind[j]]+
s+1;
11327 os <<
" 1.0 1.0 1.0 1.0\n";
11351 " NURBS mesh is not supported!");
11355 os <<
"MFEM mesh v1.0\n";
11359 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
11364 "# TETRAHEDRON = 4\n"
11369 os <<
"\ndimension\n" <<
Dim
11377 const int *
const i_AF_f = Aface_face.
GetI();
11378 const int *
const j_AF_f = Aface_face.
GetJ();
11380 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
11381 for (
const int * iface = j_AF_f + i_AF_f[iAF];
11382 iface < j_AF_f + i_AF_f[iAF+1];
11385 os << iAF+1 <<
' ';
11417 double *cg =
new double[na*
spaceDim];
11418 int *nbea =
new int[na];
11425 for (i = 0; i < na; i++)
11429 cg[i*spaceDim+j] = 0.0;
11437 for (k = 0; k < vert.
Size(); k++)
11449 for (k = 0; k < vert.
Size(); k++)
11450 if (vn[vert[k]] == 1)
11455 cg[bea*spaceDim+j] += pointmat(j,k);
11466 for (k = 0; k < vert.
Size(); k++)
11471 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
11487 double *cg =
new double[na*
spaceDim];
11488 int *nbea =
new int[na];
11495 for (i = 0; i < na; i++)
11499 cg[i*spaceDim+j] = 0.0;
11507 for (k = 0; k < vert.
Size(); k++)
11519 for (k = 0; k < vert.
Size(); k++)
11520 if (vn[vert[k]] == 1)
11525 cg[bea*spaceDim+j] += pointmat(j,k);
11536 for (k = 0; k < vert.
Size(); k++)
11541 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
11557 for (
int i = 0; i <
vertices.Size(); i++)
11559 for (
int j = 0; j <
spaceDim; j++)
11571 xnew.ProjectCoefficient(f_pert);
11579 "incompatible vector dimensions");
11587 for (
int d = 0; d <
spaceDim; d++)
11589 vertices[i](d) = xnew(d + spaceDim*i);
11606 for (
int i = 0; i <
GetNE(); i++)
11611 for (
int j = 0; j < nv; j++)
11616 for (
int i = 0; i <
GetNBE(); i++)
11621 for (
int j = 0; j < nv; j++)
11627 for (
int i = 0; i < v2v.
Size(); i++)
11632 v2v[i] = num_vert++;
11636 if (num_vert == v2v.
Size()) {
return; }
11638 Vector nodes_by_element;
11643 for (
int i = 0; i <
GetNE(); i++)
11650 for (
int i = 0; i <
GetNE(); i++)
11659 for (
int i = 0; i <
GetNE(); i++)
11664 for (
int j = 0; j < nv; j++)
11669 for (
int i = 0; i <
GetNBE(); i++)
11674 for (
int j = 0; j < nv; j++)
11698 for (
int i = 0; i <
GetNE(); i++)
11711 int num_bdr_elem = 0;
11712 int new_bel_to_edge_nnz = 0;
11713 for (
int i = 0; i <
GetNBE(); i++)
11729 if (num_bdr_elem ==
GetNBE()) {
return; }
11733 Table *new_bel_to_edge = NULL;
11737 new_be_to_edge.
Reserve(num_bdr_elem);
11741 new_be_to_face.
Reserve(num_bdr_elem);
11742 new_bel_to_edge =
new Table;
11743 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
11745 for (
int i = 0; i <
GetNBE(); i++)
11756 int row = new_be_to_face.
Size();
11760 int *new_e = new_bel_to_edge->
GetRow(row);
11761 for (
int j = 0; j < ne; j++)
11765 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
11785 for (
int i = 0; i < attribs.
Size(); i++)
11797 #ifdef MFEM_USE_MEMALLOC
11824 const int npts = point_mat.
Width();
11825 if (!npts) {
return 0; }
11826 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
11830 if (!
GetNE()) {
return 0; }
11832 double *data = point_mat.
GetData();
11839 min_dist = std::numeric_limits<double>::max();
11843 for (
int i = 0; i <
GetNE(); i++)
11847 for (
int k = 0; k < npts; k++)
11850 if (dist < min_dist(k))
11852 min_dist(k) = dist;
11861 for (
int k = 0; k < npts; k++)
11865 int res = inv_tr->
Transform(pt, ips[k]);
11868 elem_ids[k] = e_idx[k];
11872 if (pts_found != npts)
11876 for (
int k = 0; k < npts; k++)
11878 if (elem_ids[k] != -1) {
continue; }
11882 for (
int v = 0; v < elvertices.
Size(); v++)
11884 int vv = elvertices[v];
11886 const int* els = vtoel->
GetRow(vv);
11887 for (
int e = 0; e < ne; e++)
11889 if (els[e] == e_idx[k]) {
continue; }
11891 int res = inv_tr->
Transform(pt, ips[k]);
11894 elem_ids[k] = els[e];
11906 for (
int e = 0; e < neigh.
Size(); e++)
11912 int res = inv_tr->
Transform(pt, ips[k]);
11925 if (inv_trans == NULL) {
delete inv_tr; }
11927 if (warn && pts_found != npts)
11929 MFEM_WARNING((npts-pts_found) <<
" points were not found");
11943 "mixed meshes are not supported!");
11944 MFEM_ASSERT(mesh->
GetNodes(),
"meshes without nodes are not supported!");
11957 Compute(nodes, d_mt);
11960 void GeometricFactors::Compute(
const GridFunction &nodes,
11967 const int vdim = fespace->
GetVDim();
11968 const int NE = fespace->
GetNE();
11969 const int ND = fe->
GetDof();
11972 unsigned eval_flags = 0;
11982 J.
SetSize(dim*vdim*NQ*NE, my_d_mt);
11997 qi->DisableTensorProducts(!use_tensor_products);
12005 Vector Enodes(vdim*ND*NE, my_d_mt);
12006 elem_restr->Mult(nodes, Enodes);
12007 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
12011 qi->Mult(nodes, eval_flags,
X,
J,
detJ);
12026 const int vdim = fespace->
GetVDim();
12035 face_restr->
Mult(*nodes, Fnodes);
12037 unsigned eval_flags = 0;
12078 V(1) = s * ((ip.
y + layer) / n);
12083 V(2) = s * ((ip.
z + layer) / n);
12092 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
12096 int nvy = (closed) ? (ny) : (ny + 1);
12097 int nvt = mesh->
GetNV() * nvy;
12106 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
12111 for (
int i = 0; i < mesh->
GetNV(); i++)
12114 for (
int j = 0; j < nvy; j++)
12116 vc[1] = sy * (double(j) / ny);
12122 for (
int i = 0; i < mesh->
GetNE(); i++)
12127 for (
int j = 0; j < ny; j++)
12130 qv[0] = vert[0] * nvy + j;
12131 qv[1] = vert[1] * nvy + j;
12132 qv[2] = vert[1] * nvy + (j + 1) % nvy;
12133 qv[3] = vert[0] * nvy + (j + 1) % nvy;
12139 for (
int i = 0; i < mesh->
GetNBE(); i++)
12144 for (
int j = 0; j < ny; j++)
12147 sv[0] = vert[0] * nvy + j;
12148 sv[1] = vert[0] * nvy + (j + 1) % nvy;
12152 Swap<int>(sv[0], sv[1]);
12164 for (
int i = 0; i < mesh->
GetNE(); i++)
12170 sv[0] = vert[0] * nvy;
12171 sv[1] = vert[1] * nvy;
12175 sv[0] = vert[1] * nvy + ny;
12176 sv[1] = vert[0] * nvy + ny;
12192 string cname = name;
12193 if (cname ==
"Linear")
12197 else if (cname ==
"Quadratic")
12201 else if (cname ==
"Cubic")
12205 else if (!strncmp(name,
"H1_", 3))
12209 else if (!strncmp(name,
"L2_T", 4))
12213 else if (!strncmp(name,
"L2_", 3))
12220 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
12232 for (
int i = 0; i < mesh->
GetNE(); i++)
12235 for (
int j = ny-1; j >= 0; j--)
12252 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
12257 int nvt = mesh->
GetNV() * nvz;
12262 bool wdgMesh =
false;
12263 bool hexMesh =
false;
12267 for (
int i = 0; i < mesh->
GetNV(); i++)
12271 for (
int j = 0; j < nvz; j++)
12273 vc[2] = sz * (double(j) / nz);
12279 for (
int i = 0; i < mesh->
GetNE(); i++)
12289 for (
int j = 0; j < nz; j++)
12292 pv[0] = vert[0] * nvz + j;
12293 pv[1] = vert[1] * nvz + j;
12294 pv[2] = vert[2] * nvz + j;
12295 pv[3] = vert[0] * nvz + (j + 1) % nvz;
12296 pv[4] = vert[1] * nvz + (j + 1) % nvz;
12297 pv[5] = vert[2] * nvz + (j + 1) % nvz;
12304 for (
int j = 0; j < nz; j++)
12307 hv[0] = vert[0] * nvz + j;
12308 hv[1] = vert[1] * nvz + j;
12309 hv[2] = vert[2] * nvz + j;
12310 hv[3] = vert[3] * nvz + j;
12311 hv[4] = vert[0] * nvz + (j + 1) % nvz;
12312 hv[5] = vert[1] * nvz + (j + 1) % nvz;
12313 hv[6] = vert[2] * nvz + (j + 1) % nvz;
12314 hv[7] = vert[3] * nvz + (j + 1) % nvz;
12316 mesh3d->
AddHex(hv, attr);
12320 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
12321 << geom <<
"\'" << endl;
12327 for (
int i = 0; i < mesh->
GetNBE(); i++)
12332 for (
int j = 0; j < nz; j++)
12335 qv[0] = vert[0] * nvz + j;
12336 qv[1] = vert[1] * nvz + j;
12337 qv[2] = vert[1] * nvz + (j + 1) % nvz;
12338 qv[3] = vert[0] * nvz + (j + 1) % nvz;
12347 for (
int i = 0; i < mesh->
GetNE(); i++)
12358 tv[0] = vert[0] * nvz;
12359 tv[1] = vert[2] * nvz;
12360 tv[2] = vert[1] * nvz;
12364 tv[0] = vert[0] * nvz + nz;
12365 tv[1] = vert[1] * nvz + nz;
12366 tv[2] = vert[2] * nvz + nz;
12374 qv[0] = vert[0] * nvz;
12375 qv[1] = vert[3] * nvz;
12376 qv[2] = vert[2] * nvz;
12377 qv[3] = vert[1] * nvz;
12381 qv[0] = vert[0] * nvz + nz;
12382 qv[1] = vert[1] * nvz + nz;
12383 qv[2] = vert[2] * nvz + nz;
12384 qv[3] = vert[3] * nvz + nz;
12390 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
12391 << geom <<
"\'" << endl;
12397 if ( hexMesh && wdgMesh )
12401 else if ( hexMesh )
12405 else if ( wdgMesh )
12418 string cname = name;
12419 if (cname ==
"Linear")
12423 else if (cname ==
"Quadratic")
12427 else if (cname ==
"Cubic")
12431 else if (!strncmp(name,
"H1_", 3))
12435 else if (!strncmp(name,
"L2_T", 4))
12439 else if (!strncmp(name,
"L2_", 3))
12446 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
12458 for (
int i = 0; i < mesh->
GetNE(); i++)
12461 for (
int j = nz-1; j >= 0; j--)
12482 os << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
12483 <<
" 0 0 " << i <<
" -1 0\n";
12490 double mid[3] = {0, 0, 0};
12491 for (
int j = 0; j < 2; j++)
12493 for (
int k = 0; k <
spaceDim; k++)
12498 os << NumOfVertices+i <<
" "
12499 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
12500 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
int GetNPoints() const
Returns the number of the points in the integration rule.
static Mesh MakePeriodic(const Mesh &orig_mesh, const std::vector< int > &v2v)
Create a periodic mesh by identifying vertices of orig_mesh.
Abstract class for all finite elements.
const IntegrationRule * IntRule
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Geometry::Type GetGeometryType() const
void Loader(std::istream &input, int generate_edges=0, std::string parse_tag="")
Mesh * Make3D(int nsteps, double rstep, double aspect, int order, bool sfc)
Ordering::Type GetOrdering() const
Return the ordering method.
int Size() const
Return the logical size of the array.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
void SetCoordsFromPatches(Vector &Nodes)
void PrintSurfaces(const Table &Aface_face, std::ostream &out) const
Print set of disjoint surfaces:
void Get(double *p, const int dim) const
virtual void Save(const char *fname, int precision=16) const
void GetPointMatrix(int i, DenseMatrix &pointmat) const
static const int vtk_quadratic_hex[27]
int * CartesianPartitioning(int nxyz[])
void METIS_PartGraphVKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
int GetDim() const
Returns the reference space dimension for the finite element.
int GetNDofs() const
Returns number of degrees of freedom.
Arc::Index insert_arc(Node::Index i, Node::Index j, Float w=1, Float b=1)
FaceInformation GetFaceInformation(int f) const
Class for an integration rule - an Array of IntegrationPoint.
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
void Init(int ind1, int ind2, int ind3, int ind4, int attr=1, int ref_flag=0)
Initialize the vertex indices and the attribute of a Tetrahedron.
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
void ScaleElements(double sf)
static const int HighOrderMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to arbitrary-order Lagrange VTK geometries.
Class for grid function - Vector with associated FE space.
int DofToVDof(int dof, int vd, int ndofs=-1) const
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T...
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
T * end()
STL-like end. Returns pointer after the last element of the array.
void ReadVTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo)
void UseExternalData(double *ext_data, int i, int j, int k)
void FreeElement(Element *E)
void Mult(const Vector &x, Vector &y) const override=0
Extract the face degrees of freedom from x into y.
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
int CheckElementOrientation(bool fix_it=true)
Check (and optionally attempt to fix) the orientation of the elements.
void SetVertices(const Vector &vert_coord)
virtual Element * Duplicate(Mesh *m) const =0
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
static const int vtk_quadratic_tet[10]
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
void Make2D(int nx, int ny, Element::Type type, double sx, double sy, bool generate_edges, bool sfc_ordering)
Vector J
Jacobians of the element transformations at all quadrature points.
void AddColumnsInRow(int r, int ncol)
Vector X
Mapped (physical) coordinates of all quadrature points.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
static Mesh MakeSimplicial(const Mesh &orig_mesh)
Base class for vector Coefficients that optionally depend on time and space.
void UniformRefinement3D_base(Array< int > *f2qf=NULL, DSTable *v_to_v_p=NULL, bool update_nodes=true)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
class LinearPyramidFiniteElement PyramidFE
Array< Element * > boundary
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
int * GeneratePartitioning(int nparts, int part_method=1)
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
CoarseFineTransformations CoarseFineTr
virtual void LimitNCLevel(int max_nc_level)
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
void MoveVertices(const Vector &displacements)
void SetSize(int s)
Resize the vector to size s.
int GetNBE() const
Returns number of boundary elements.
static FiniteElement * GetTransformationFEforElementType(Element::Type)
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'.
void SetElementOrder(int i, int p)
Sets the order of the i'th finite element.
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, double sx=1.0, double sy=1.0, double sz=1.0, bool sfc_ordering=true)
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Lists all edges/faces in the nonconforming mesh.
virtual void UniformRefinement2D()
Refine a mixed 2D mesh uniformly.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Element::Type GetElementType(int i) const
Returns the type of element i.
double Norml2() const
Returns the l2 norm of the vector.
void ReadNetgen2DMesh(std::istream &input, int &curved)
void ShiftRight(int &a, int &b, int &c)
void SetDims(int rows, int nnz)
static const int FaceVert[NumFaces][MaxFaceVert]
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
static Mesh MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
Create a refined (by any factor) version of orig_mesh.
unsigned matrix
index into NCList::point_matrices[geom]
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
T * GetData()
Returns the data.
Vector detJ
Determinants of the Jacobians at all quadrature points.
static const int FaceVert[NumFaces][MaxFaceVert]
static const int Edges[NumEdges][2]
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
int Push(int r, int c, int f)
Check to see if this entry is in the table and add it to the table if it is not there. Returns the number assigned to the table entry.
Piecewise-(bi/tri)linear continuous finite elements.
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
void Print(std::ostream &out) const
I/O: Print the mesh in "MFEM NC mesh v1.0" format.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Mesh * Extrude1D(Mesh *mesh, const int ny, const double sy, const bool closed)
Extrude a 1D mesh.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
void Transform(void(*f)(const Vector &, Vector &))
int GetBdrElementEdgeIndex(int i) const
void AverageVertices(const int *indexes, int n, int result)
Averages the vertices with given indexes and saves the result in vertices[result].
int GetNE() const
Returns number of elements.
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
void FinalizeWedgeMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a wedge Mesh.
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
FDualNumber< tbase > acos(const FDualNumber< tbase > &f)
acos([dual number])
Evaluate the derivatives at quadrature points.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
NodeExtrudeCoefficient(const int dim, const int n_, const double s_)
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
void GetElementLocalToGlobal(Array< int > &lelem_elem)
void SetMeshGen()
Determine the mesh generator bitmask meshgen, see MeshGenerator().
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
Data type quadrilateral element.
void GetVertices(Vector &vert_coord) const
void ReadNetgen3DMesh(std::istream &input)
bool UsesTensorBasis(const FiniteElementSpace &fes)
Data arrays will be written in ASCII format.
double * GetData() const
Returns the matrix data array.
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
void ReadInlineMesh(std::istream &input, bool generate_edges=false)
double * GetData() const
Return a pointer to the beginning of the Vector data.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
void RemoveInternalBoundaries()
void DebugDump(std::ostream &out) const
Output an NCMesh-compatible debug dump.
Mesh * Extrude2D(Mesh *mesh, const int nz, const double sz)
Extrude a 2D mesh.
Geometry::Type GetElementBaseGeometry(int i) const
void Print(const Mesh &mesh, const adios2stream::mode print_mode=mode::sync)
int Size_of_connections() const
int spaceDim
dimensions of the elements and the vertex coordinates
virtual long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
friend class NURBSExtension
Geometry::Type GetBdrElementGeometry(int i) const
void GetVertexToVertexTable(DSTable &) const
void RedRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
void WriteBinaryOrASCII(std::ostream &os, std::vector< char > &buf, const T &val, const char *suffix, VTKFormat format)
Write either ASCII data to the stream or binary data to the buffer depending on the given format...
void ReadCubit(const char *filename, int &curved, int &read_gf)
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
void add(const Vector &v1, const Vector &v2, Vector &v)
void DeleteAll()
Delete the whole array.
void AddConnections(int r, const int *c, int nc)
int master
master number (in Mesh numbering)
void InitRefinementTransforms()
const NCList & GetFaceList()
Return the current list of conforming and nonconforming faces.
void UniformRefinement2D_base(bool update_nodes=true)
Array< NCFaceInfo > nc_faces_info
Element * ReadElement(std::istream &)
void KnotInsert(Array< KnotVector * > &kv)
void OnMeshUpdated(Mesh *mesh)
std::vector< int > CreatePeriodicVertexMapping(const std::vector< Vector > &translations, double tol=1e-8) const
Creates a mapping v2v from the vertex indices of the mesh such that coincident vertices under the giv...
Element * NewElement(int geom)
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
void SetVerticesFromNodes(const GridFunction *nodes)
Helper to set vertex coordinates given a high-order curvature function.
Structure for storing face geometric factors: coordinates, Jacobians, determinants of the Jacobians...
ElementTransformation * GetBdrElementTransformation(int i)
Returns the transformation defining the i-th boundary element.
virtual void UpdateMeshPointer(Mesh *new_mesh)
int AddBdrTriangle(int v1, int v2, int v3, int attr=1)
Mesh * Make2D(int nsteps, double rstep, double phi, double aspect, int order, bool sfc)
void PrintTopo(std::ostream &out, const Array< int > &e_to_k) const
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void MakeRefined_(Mesh &orig_mesh, const Array< int > ref_factors, int ref_type)
Internal function used in Mesh::MakeRefined.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Piecewise-(bi)cubic continuous finite elements.
int GetNE() const
Returns number of elements in the mesh.
uint rank(Node::Index i) const
PointFiniteElement PointFE
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void GetElementJacobian(int i, DenseMatrix &J)
Data type hexahedron element.
static int GetTetOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void GetVertexDofs(int i, Array< int > &dofs) const
int AddVertex(double x, double y=0.0, double z=0.0)
bool Nonconforming() const
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Vector J
Jacobians of the element transformations at all quadrature points.
Data type Pyramid element.
static const int Edges[NumEdges][2]
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void MoveNodes(const Vector &displacements)
void PrintWithPartitioning(int *partitioning, std::ostream &os, int elem_attr=0) const
Prints the mesh with boundary elements given by the boundary of the subdomains, so that the boundary ...
double f(const Vector &xvec)
Symmetric 3D Table stored as an array of rows each of which has a stack of column, floor, number nodes. The number of the node is assigned by counting the nodes from zero as they are pushed into the table. Diagonals of any kind are not allowed so the row, column and floor must all be different for each node. Only one node is stored for all 6 symmetric entries that are indexable by unique triplets of row, column, and floor.
const Table & GetDerefinementTable()
void PrintVTK(std::ostream &os)
Native ordering as defined by the FiniteElement.
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void Make3D(int nx, int ny, int nz, Element::Type type, double sx, double sy, double sz, bool sfc_ordering)
void AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
virtual void SetAttributes()
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
void EnsureNCMesh(bool simplices_nonconforming=false)
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
Mesh * GetMesh() const
Returns the mesh.
Vector detJ
Determinants of the Jacobians at all quadrature points.
int FindCoarseElement(int i)
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
virtual void Derefine(const Array< int > &derefs)
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
void FindNeighbors(int elem, Array< int > &neighbors, const Array< int > *search_set=NULL)
void Mult(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
Interpolate the E-vector e_vec to quadrature points.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
static const int Map[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to linear VTK geometries.
void AddQuadFaceElement(int lf, int gf, int el, int v0, int v1, int v2, int v3)
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
static const int NumVerts[NumGeom]
void CheckPartitioning(int *partitioning_)
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void AddConnection(int r, int c)
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
void ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
STable3D * GetElementToFaceTable(int ret_ftbl=0)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Type
Constants for the classes derived from Element.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const
void CheckDisplacements(const Vector &displacements, double &tmax)
This structure stores the low level information necessary to interpret the configuration of elements ...
NURBSExtension * StealNURBSext()
Array< DenseMatrix * > point_matrices[Geometry::NumGeom]
List of unique point matrices for each slave geometry.
VTKFormat
Data array format for VTK and VTU files.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
int AddBdrSegment(int v1, int v2, int attr=1)
int AddElement(Element *elem)
The parameter elem should be allocated using the NewElement() method.
int GetMaxElementOrder() const
Return the maximum polynomial order.
void ReadLineMesh(std::istream &input)
void GetBdrElementTopo(Array< Element * > &boundary) const
virtual void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
GeometryRefiner GlobGeometryRefiner
const CoarseFineTransformations & GetRefinementTransforms()
void SetLayer(const int l)
int GetAttribute() const
Return element's attribute.
int GetElementToEdgeTable(Table &, Array< int > &)
T * begin()
STL-like begin. Returns pointer to the first element of the array.
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i...
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
void SetKnotsFromPatches()
void RandomRefinement(double prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Data type triangle element.
const Element * GetElement(int i) const
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
void GetElementTopo(Array< Element * > &elements) const
signed char local
local number within 'element'
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
int GetVDim()
Returns dimension of the vector.
void SetNodalFESpace(FiniteElementSpace *nfes)
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
int Size() const
Returns the number of TYPE I elements.
int GetVDim() const
Returns vector dimension.
Vector normal
Normals at all quadrature points.
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
static const int Edges[NumEdges][2]
FiniteElementSpace * FESpace()
virtual MFEM_DEPRECATED void ReorientTetMesh()
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
static const int QuadraticMap[Geometry::NUM_GEOMETRIES]
Map from MFEM's Geometry::Type to legacy quadratic VTK geometries/.
void GetColumn(int c, Vector &col) const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
void ReadMFEMMesh(std::istream &input, int version, int &curved)
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
This function is not public anymore. Use GeneralRefinement instead.
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
void FindPartitioningComponents(Table &elem_elem, const Array< int > &partitioning, Array< int > &component, Array< int > &num_comp)
Data type tetrahedron element.
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
void GetElementInteriorDofs(int i, Array< int > &dofs) const
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
List of mesh geometries stored as Array<Geometry::Type>.
A general vector function coefficient.
void PrintVTU(std::ostream &os, int ref=1, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0, bool bdr_elements=false)
int CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
int SpaceDimension() const
GeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, MemoryType d_mt=MemoryType::DEFAULT)
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
void PrintElementsWithPartitioning(int *partitioning, std::ostream &out, int interior_faces=0)
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
static const int Edges[NumEdges][2]
void AddPointFaceElement(int lf, int gf, int el)
Used in GenerateFaces()
void Make1D(int n, double sx=1.0)
Creates a 1D mesh for the interval [0,sx] divided into n equal intervals.
double * Data() const
Returns the matrix data array.
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
void Swap(Array< T > &, Array< T > &)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
virtual void NonconformingRefinement(const Array< Refinement > &refinements, int nc_limit=0)
This function is not public anymore. Use GeneralRefinement instead.
virtual void PushTransform(int tr)
Add 'tr' to the current chain of coarse-fine transformations.
void FinalizeMesh(int refine=0, bool fix_orientation=true)
Finalize the construction of any type of Mesh.
double p(const Vector &x, double t)
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
const NURBSExtension * GetNURBSext() const
int SpaceDimension() const
Return the space dimension of the NCMesh.
Nonconforming edge/face within a bigger edge/face.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
static Mesh LoadFromFile(const char *filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
static void PrintElement(const Element *, std::ostream &)
virtual void NURBSUniformRefinement()
Refine NURBS mesh.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
const IntegrationRule * IntRule
MemoryType
Memory types supported by MFEM.
Vector X
Mapped (physical) coordinates of all quadrature points.
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
static void GetElementArrayEdgeTable(const Array< Element * > &elem_array, const DSTable &v_to_v, Table &el_to_edge)
void AddAColumnInRow(int r)
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
int AddBdrQuad(int v1, int v2, int v3, int v4, int attr=1)
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
int FindRoots(const Vector &z, Vector &x)
bool IsSlaveFace(const FaceInfo &fi) const
int Push4(int r, int c, int f, int t)
Check to see if this entry is in the table and add it to the table if it is not there. The entry is addressed by the three smallest values of (r,c,f,t). Returns the number assigned to the table entry.
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, double sx=1.0, double sy=1.0, bool sfc_ordering=true)
void InitMesh(int Dim_, int spaceDim_, int NVert, int NElem, int NBdrElem)
Begin construction of a mesh.
void ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf, bool &finalize_topo, const std::string &xml_prefix="")
Helper struct for defining a connectivity table, see Table::MakeFromList.
Geometry::Type GetElementGeometry(int i) const
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
class Linear3DFiniteElement TetrahedronFE
A class to initialize the size of a Tensor.
static Mesh MakeCartesian1D(int n, double sx=1.0)
Array< Element * > elements
const Table & ElementToElementTable()
double GetLength(int i, int j) const
Return the length of the segment from node i to node j.
void SetRelaxedHpConformity(bool relaxed=true)
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Linear2DFiniteElement TriangleFE
static const int FaceVert[NumFaces][MaxFaceVert]
const CoarseFineTransformations & GetRefinementTransforms()
bool RefineByError(const Array< double > &elem_error, double threshold, int nonconforming=-1, int nc_limit=0)
Evaluate the values at quadrature points.
Operation GetLastOperation() const
Return type of last modification of the mesh.
Table * GetFaceToElementTable() const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
NURBSExtension * NURBSext
Optional NURBS mesh extension.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void PrintBdrVTU(std::string fname, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0)
void GetNode(int i, double *coord) const
virtual const char * Name() const
const Table & ElementToFaceTable() const
static const int Edges[NumEdges][2]
Class for integration point with weight.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
Element * ReadElementWithoutAttr(std::istream &)
void Swap(Mesh &other, bool non_geometry)
FaceGeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags, FaceType type)
int GetNumFacesWithGhost() const
Return the number of faces (3D), edges (2D) or vertices (1D) including ghost faces.
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
static const int vtk_quadratic_wedge[18]
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
static const int DimStart[MaxDim+2]
const FiniteElementSpace * GetNodalFESpace() const
class LinearWedgeFiniteElement WedgeFE
virtual void Print(std::ostream &os=mfem::out) const
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Array< int > leaf_elements
finest elements, in Mesh ordering (+ ghosts)
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
static const int Orient[NumOrient][NumVert]
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
void DegreeElevate(int rel_degree, int degree=16)
signed char geom
Geometry::Type (faces only) (char to save RAM)
void FindTMax(Vector &c, Vector &x, double &tmax, const double factor, const int Dim)
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int index(int i, int j, int nx, int ny)
static const int Edges[NumEdges][2]
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
Lexicographic ordering for tensor-product FiniteElements.
Array< FaceInfo > faces_info
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
virtual int GetNVertices() const =0
int parent
Coarse Element index in the coarse mesh.
void SetAttribute(const int attr)
Set element's attribute.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
static const int Orient[NumOrient][NumVert]
STable3D * GetFacesTable()
Piecewise-(bi)quadratic continuous finite elements.
int Dimension() const
Return the dimension of the NCMesh.
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
NCMesh * ncmesh
Optional nonconforming mesh extension.
Geometry::Type GetBdrElementBaseGeometry(int i) const
virtual void CheckDerefinementNCLevel(const Table &deref_table, Array< int > &level_ok, int max_nc_level)
virtual void UniformRefinement3D()
Refine a mixed 3D mesh uniformly.
int AddBdrElement(Element *elem)
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
static void PrintElementWithoutAttr(const Element *, std::ostream &)
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
BlockArray< Element > elements
virtual bool NonconformingDerefinement(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
NC version of GeneralDerefinement.
void Clear()
Clear the contents of the Mesh.
int GetNEdges() const
Return the number of edges.
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
double AggregateError(const Array< double > &elem_error, const int *fine, int nfine, int op)
Derefinement helper.
int GetNFaces() const
Return the number of faces in a 3D mesh.
BiLinear2DFiniteElement QuadrilateralFE
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
Mesh & operator=(Mesh &&mesh)
Move assignment operstor.
virtual void PrintXG(std::ostream &os=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
Evaluate the derivatives at quadrature points.
void ReadTrueGridMesh(std::istream &input)
const FiniteElementCollection * FEColl() const
void GetNodes(Vector &node_coord) const
virtual int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type, does not count master nonconforming face...
int NumberOfEntries() const
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
void SetNode(int i, const double *coord)
virtual void PushTransform(int tr)
Add 'tr' to the current chain of coarse-fine transformations.
void CreateVTKElementConnectivity(Array< int > &con, Geometry::Type geom, int ref)
Create the VTK element connectivity array for a given element geometry and refinement level...
int GetId(int p1, int p2)
Get id of item whose parents are p1, p2... Create it if it doesn't exist.
Arbitrary order H1-conforming (continuous) finite elements.
void XYZ_VectorFunction(const Vector &p, Vector &v)
void GetElementColoring(Array< int > &colors, int el0=0)
void ChangeVertexDataOwnership(double *vertices, int len_vertices, bool zerocopy=false)
Set the internal Vertex array to point to the given vertices array without assuming ownership of the ...
const std::string filename
const NCList & GetEdgeList()
Return the current list of conforming and nonconforming edges.
MemAlloc< Tetrahedron, 1024 > TetMemory
Table * GetFaceEdgeTable() const
Returns the face-to-edge Table (3D)
bool IsGhost(const Element &el) const
Return true if the Element el is a ghost element.
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
Node::Index insert_node(Float length=1)
TriLinear3DFiniteElement HexahedronFE
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
void SetNodes(const Vector &node_coord)
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
static const int * VertexPermutation[Geometry::NUM_GEOMETRIES]
Permutation from MFEM's vertex ordering to VTK's vertex ordering.
void ScaleSubdomains(double sf)
void DegreeElevate(int rel_degree, int degree=16)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void InitFromNCMesh(const NCMesh &ncmesh)
Initialize vertices/elements/boundary/tables from a nonconforming mesh.
Base class for operators that extracts Face degrees of freedom.
void Bisection(int i, const DSTable &, int *, int *, int *)
Bisect a triangle: element with index i is bisected.
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn't exist.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Rank 3 tensor (array of matrices)
Class for parallel meshes.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Abstract data type element.
int GetAttribute(int i) const
Return the attribute of element i.
void AddTriangleFaceElement(int lf, int gf, int el, int v0, int v1, int v2)
void WriteBase64WithSizeAndClear(std::ostream &os, std::vector< char > &buf, int compression_level)
Encode in base 64 (and potentially compress) the given data, write it to the output stream (with a he...
const Table & ElementToEdgeTable() const
Data type line segment element.
Linear1DFiniteElement SegmentFE
Array< GeometricFactors * > geom_factors
Optional geometric factors.
int GetNFDofs() const
Number of all scalar face-interior dofs.
void BdrBisection(int i, const HashTable< Hashed2 > &)
Bisect a boundary triangle: boundary element with index i is bisected.
MPI_Comm GetGlobalMPI_Comm()
Get MFEM's "global" MPI communicator.
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
int NumberOfElements()
Return the number of elements added to the table.
void GenerateNCFaceInfo()
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
virtual Type GetType() const =0
Returns element's type.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Table * GetVertexToElementTable()
The returned Table must be destroyed by the caller.
const Element * GetBdrElement(int i) const
void ConvertToPatches(const Vector &Nodes)
static const int Orient[NumOrient][NumVert]
void GetElementData(const Array< Element * > &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void UpdateNodes()
Update the nodes of a curved mesh after refinement.
Defines the position of a fine element within a coarse element.
void KnotInsert(Array< KnotVector * > &kv)
void Print(std::ostream &out) const
virtual void Refine(const Array< Refinement > &refinements)
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
void GreenRefinement(int i, const DSTable &v_to_v, int *edge1, int *edge2, int *middle)
Arbitrary order "L2-conforming" discontinuous finite elements.
Evaluate the values at quadrature points.
double f(const Vector &p)
void Printer(std::ostream &out=mfem::out, std::string section_delimiter="") const
static const int FaceVert[NumFaces][MaxFaceVert]
Class used to extrude the nodes of a mesh.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...