15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
17 #include "../general/binaryio.hpp"
18 #include "../general/text.hpp"
19 #include "../general/device.hpp"
20 #include "../general/tic_toc.hpp"
21 #include "../general/gecko.hpp"
22 #include "../fem/quadinterpolator.hpp"
37 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
42 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
47 int*,
int*,
int*,
int*,
int*,
idxtype*);
49 int*,
int*,
int*,
int*,
int*,
idxtype*);
51 int*,
int*,
int*,
int*,
int*,
idxtype*);
68 void Mesh::GetElementCenter(
int i,
Vector ¢er)
71 int geom = GetElementBaseGeometry(i);
86 return pow(fabs(J.
Det()), 1./Dim);
98 double Mesh::GetElementSize(
int i,
int type)
100 return GetElementSize(GetElementTransformation(i), type);
103 double Mesh::GetElementSize(
int i,
const Vector &dir)
107 GetElementJacobian(i, J);
109 return sqrt((d_hat * d_hat) / (dir * dir));
112 double Mesh::GetElementVolume(
int i)
134 for (
int d = 0; d < spaceDim; d++)
143 for (
int i = 0; i < NumOfVertices; i++)
145 coord = GetVertex(i);
146 for (
int d = 0; d < spaceDim; d++)
148 if (coord[d] < min(d)) { min(d) = coord[d]; }
149 if (coord[d] > max(d)) { max(d) = coord[d]; }
155 const bool use_boundary =
false;
156 int ne = use_boundary ? GetNBE() : GetNE();
164 for (
int i = 0; i < ne; i++)
168 GetBdrElementFace(i, &fn, &fo);
170 Tr = GetFaceElementTransformations(fn, 5);
177 T = GetElementTransformation(i);
181 for (
int j = 0; j < pointmat.
Width(); j++)
183 for (
int d = 0; d < pointmat.
Height(); d++)
185 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
186 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
193 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
194 double &kappa_min,
double &kappa_max,
202 sdim = SpaceDimension();
204 if (Vh) { Vh->
SetSize(NumOfElements); }
205 if (Vk) { Vk->
SetSize(NumOfElements); }
208 h_max = kappa_max = -h_min;
209 if (dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
211 for (i = 0; i < NumOfElements; i++)
213 GetElementJacobian(i, J);
214 h = pow(fabs(J.
Weight()), 1.0/
double(dim));
215 kappa = (dim == sdim) ?
217 if (Vh) { (*Vh)(i) = h; }
218 if (Vk) { (*Vk)(i) = kappa; }
220 if (h < h_min) { h_min = h; }
221 if (h > h_max) { h_max = h; }
222 if (kappa < kappa_min) { kappa_min =
kappa; }
223 if (kappa > kappa_max) { kappa_max =
kappa; }
228 void Mesh::PrintElementsByGeometry(
int dim,
232 for (
int g = Geometry::DimStart[dim], first = 1;
233 g < Geometry::DimStart[dim+1]; g++)
235 if (!num_elems_by_geom[g]) {
continue; }
236 if (!first) { out <<
" + "; }
238 out << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
244 double h_min, h_max, kappa_min, kappa_max;
246 out <<
"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, out);
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, out);
312 <<
"Number of elements : " << GetNE() <<
" -- ";
313 PrintElementsByGeometry(Dim, num_elems_by_geom, out);
315 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
316 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, out);
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 out <<
'\n' << std::flush;
331 case Element::POINT :
return &
PointFE;
332 case Element::SEGMENT :
return &
SegmentFE;
337 case Element::WEDGE :
return &
WedgeFE;
339 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
342 MFEM_ABORT(
"Unknown element type");
351 ElTr->
ElementType = ElementTransformation::ELEMENT;
356 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
362 Nodes->FESpace()->GetElementVDofs(i, vdofs);
365 int n = vdofs.
Size()/spaceDim;
367 for (
int k = 0; k < spaceDim; k++)
369 for (
int j = 0; j < n; j++)
371 pm(k,j) = nodes(vdofs[n*k+j]);
374 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
378 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
383 ElTr->
ElementType = ElementTransformation::ELEMENT;
389 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
390 int nv = elements[i]->GetNVertices();
391 const int *v = elements[i]->GetVertices();
392 int n = vertices.Size();
394 for (
int k = 0; k < spaceDim; k++)
396 for (
int j = 0; j < nv; j++)
398 pm(k, j) = nodes(k*n+v[j]);
401 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
405 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
407 Nodes->FESpace()->GetElementVDofs(i, vdofs);
408 int n = vdofs.
Size()/spaceDim;
410 for (
int k = 0; k < spaceDim; k++)
412 for (
int j = 0; j < n; j++)
414 pm(k,j) = nodes(vdofs[n*k+j]);
417 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
423 GetElementTransformation(i, &Transformation);
425 return &Transformation;
430 GetBdrElementTransformation(i, &BdrTransformation);
431 return &BdrTransformation;
438 ElTr->
ElementType = ElementTransformation::BDR_ELEMENT;
443 GetBdrPointMatrix(i, pm);
444 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
454 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
455 int n = vdofs.
Size()/spaceDim;
457 for (
int k = 0; k < spaceDim; k++)
459 for (
int j = 0; j < n; j++)
461 pm(k,j) = nodes(vdofs[n*k+j]);
468 int elem_id, face_info;
469 GetBdrElementAdjacentElement(i, elem_id, face_info);
471 GetLocalFaceTransformation(GetBdrElementType(i),
472 GetElementType(elem_id),
473 FaceElemTr.Loc1.Transf, face_info);
478 Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
479 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
480 "Mesh requires nodal Finite Element.");
482 FaceElemTr.Loc1.Transf.ElementNo = elem_id;
483 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
484 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
485 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
487 ElTr->
SetFE(face_el);
494 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
501 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
502 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
504 for (
int i = 0; i < spaceDim; i++)
506 for (
int j = 0; j < nv; j++)
508 pm(i, j) = vertices[v[j]](i);
511 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
515 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
521 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
522 int n = vdofs.
Size()/spaceDim;
524 for (
int i = 0; i < spaceDim; i++)
526 for (
int j = 0; j < n; j++)
528 pm(i, j) = nodes(vdofs[n*i+j]);
535 FaceInfo &face_info = faces_info[FaceNo];
540 GetLocalFaceTransformation(face_type,
541 GetElementType(face_info.
Elem1No),
542 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
545 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
547 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
548 "Mesh requires nodal Finite Element.");
551 FaceElemTr.Loc1.Transf.ElementNo = face_info.
Elem1No;
552 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
553 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
554 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
563 GetFaceTransformation(FaceNo, &FaceTransformation);
564 return &FaceTransformation;
571 GetFaceTransformation(EdgeNo, EdTr);
576 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
587 GetEdgeVertices(EdgeNo, v);
590 for (
int i = 0; i < spaceDim; i++)
592 for (
int j = 0; j < nv; j++)
594 pm(i, j) = vertices[v[j]](i);
597 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
601 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
607 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
608 int n = vdofs.
Size()/spaceDim;
610 for (
int i = 0; i < spaceDim; i++)
612 for (
int j = 0; j < n; j++)
614 pm(i, j) = nodes(vdofs[n*i+j]);
617 EdTr->
SetFE(edge_el);
621 MFEM_ABORT(
"Not implemented.");
628 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
629 return &EdgeTransformation;
633 void Mesh::GetLocalPtToSegTransformation(
648 void Mesh::GetLocalSegToTriTransformation(
657 tv = tri_t::Edges[i/64];
658 so = seg_t::Orient[i%64];
661 for (
int j = 0; j < 2; j++)
663 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
664 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
668 void Mesh::GetLocalSegToQuadTransformation(
677 qv = quad_t::Edges[i/64];
678 so = seg_t::Orient[i%64];
681 for (
int j = 0; j < 2; j++)
683 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
684 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
688 void Mesh::GetLocalTriToTetTransformation(
696 const int *tv = tet_t::FaceVert[i/64];
699 const int *to = tri_t::Orient[i%64];
703 for (
int j = 0; j < 3; j++)
706 locpm(0, j) = vert.
x;
707 locpm(1, j) = vert.
y;
708 locpm(2, j) = vert.
z;
712 void Mesh::GetLocalTriToWdgTransformation(
720 MFEM_VERIFY(i < 128,
"Local face index " << i/64
721 <<
" is not a triangular face of a wedge.");
722 const int *pv = pri_t::FaceVert[i/64];
725 const int *to = tri_t::Orient[i%64];
729 for (
int j = 0; j < 3; j++)
732 locpm(0, j) = vert.
x;
733 locpm(1, j) = vert.
y;
734 locpm(2, j) = vert.
z;
738 void Mesh::GetLocalQuadToHexTransformation(
746 const int *hv = hex_t::FaceVert[i/64];
748 const int *qo = quad_t::Orient[i%64];
751 for (
int j = 0; j < 4; j++)
754 locpm(0, j) = vert.
x;
755 locpm(1, j) = vert.
y;
756 locpm(2, j) = vert.
z;
760 void Mesh::GetLocalQuadToWdgTransformation(
768 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
769 <<
" is not a quadrilateral face of a wedge.");
770 const int *pv = pri_t::FaceVert[i/64];
772 const int *qo = quad_t::Orient[i%64];
775 for (
int j = 0; j < 4; j++)
778 locpm(0, j) = vert.
x;
779 locpm(1, j) = vert.
y;
780 locpm(2, j) = vert.
z;
788 for (
int i = 0; i < geom_factors.Size(); i++)
800 geom_factors.Append(gf);
808 for (
int i = 0; i < face_geom_factors.Size(); i++)
821 face_geom_factors.Append(gf);
825 void Mesh::DeleteGeometricFactors()
827 for (
int i = 0; i < geom_factors.Size(); i++)
829 delete geom_factors[i];
831 geom_factors.SetSize(0);
832 for (
int i = 0; i < face_geom_factors.Size(); i++)
834 delete face_geom_factors[i];
836 face_geom_factors.SetSize(0);
839 void Mesh::GetLocalFaceTransformation(
845 GetLocalPtToSegTransformation(Transf, info);
848 case Element::SEGMENT:
849 if (elem_type == Element::TRIANGLE)
851 GetLocalSegToTriTransformation(Transf, info);
855 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
856 GetLocalSegToQuadTransformation(Transf, info);
860 case Element::TRIANGLE:
861 if (elem_type == Element::TETRAHEDRON)
863 GetLocalTriToTetTransformation(Transf, info);
867 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
868 GetLocalTriToWdgTransformation(Transf, info);
872 case Element::QUADRILATERAL:
873 if (elem_type == Element::HEXAHEDRON)
875 GetLocalQuadToHexTransformation(Transf, info);
879 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
880 GetLocalQuadToWdgTransformation(Transf, info);
889 FaceInfo &face_info = faces_info[FaceNo];
892 FaceElemTr.SetConfigurationMask(cmask);
893 FaceElemTr.Elem1 = NULL;
894 FaceElemTr.Elem2 = NULL;
898 if (mask & FaceElementTransformations::HAVE_ELEM1)
900 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
901 FaceElemTr.Elem1 = &Transformation;
908 FaceElemTr.Elem2No = face_info.
Elem2No;
909 if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
910 FaceElemTr.Elem2No >= 0)
913 if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
914 { MFEM_ABORT(
"NURBS mesh not supported!"); }
916 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
917 FaceElemTr.Elem2 = &Transformation2;
922 if (mask & FaceElementTransformations::HAVE_FACE)
924 GetFaceTransformation(FaceNo, &FaceElemTr);
929 FaceElemTr.SetGeometryType(GetFaceGeometryType(FaceNo));
933 int face_type = GetFaceElementType(FaceNo);
934 if (mask & FaceElementTransformations::HAVE_LOC1)
936 int elem_type = GetElementType(face_info.
Elem1No);
937 GetLocalFaceTransformation(face_type, elem_type,
938 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
941 if ((mask & FaceElementTransformations::HAVE_LOC2) &&
942 FaceElemTr.Elem2No >= 0)
944 int elem_type = GetElementType(face_info.
Elem2No);
945 GetLocalFaceTransformation(face_type, elem_type,
946 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
949 if (Nonconforming() && IsSlaveFace(face_info))
951 ApplyLocalSlaveTransformation(FaceElemTr, face_info,
false);
956 FaceElemTr.SetConfigurationMask(cmask);
962 double dist = FaceElemTr.CheckConsistency();
965 mfem::out <<
"\nInternal error: face id = " << FaceNo
966 <<
", dist = " << dist <<
'\n';
967 FaceElemTr.CheckConsistency(1);
968 MFEM_ABORT(
"internal error");
978 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
984 #ifdef MFEM_THREAD_SAFE
989 MFEM_ASSERT(fi.
NCFace >= 0,
"");
990 MFEM_ASSERT(nc_faces_info[fi.
NCFace].Slave,
"internal error");
1001 std::swap(composition(0,0), composition(0,1));
1002 std::swap(composition(1,0), composition(1,1));
1023 int fn = GetBdrFace(BdrElemNo);
1026 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
1030 tr = GetFaceElementTransformations(fn, 21);
1031 tr->
Attribute = boundary[BdrElemNo]->GetAttribute();
1033 tr->
ElementType = ElementTransformation::BDR_FACE;
1037 int Mesh::GetBdrFace(
int BdrElemNo)
const
1042 fn = be_to_face[BdrElemNo];
1046 fn = be_to_edge[BdrElemNo];
1050 fn = boundary[BdrElemNo]->GetVertices()[0];
1055 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const
1057 *Elem1 = faces_info[Face].
Elem1No;
1058 *Elem2 = faces_info[Face].Elem2No;
1061 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const
1063 *Inf1 = faces_info[Face].Elem1Inf;
1064 *Inf2 = faces_info[Face].Elem2Inf;
1067 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2,
int *NCFace)
const
1069 *Inf1 = faces_info[Face].Elem1Inf;
1070 *Inf2 = faces_info[Face].Elem2Inf;
1071 *NCFace = faces_info[Face].NCFace;
1078 case 1:
return Geometry::POINT;
1079 case 2:
return Geometry::SEGMENT;
1081 if (Face < NumOfFaces)
1083 return faces[Face]->GetGeometryType();
1086 const int nc_face_id = faces_info[Face].NCFace;
1087 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1088 return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
1090 return Geometry::INVALID;
1095 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1103 NumOfElements = NumOfBdrElements = 0;
1104 NumOfEdges = NumOfFaces = 0;
1105 nbInteriorFaces = -1;
1106 nbBoundaryFaces = -1;
1107 meshgen = mesh_geoms = 0;
1113 last_operation = Mesh::NONE;
1116 void Mesh::InitTables()
1119 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1122 void Mesh::SetEmpty()
1128 void Mesh::DestroyTables()
1133 DeleteGeometricFactors();
1144 void Mesh::DestroyPointers()
1146 if (own_nodes) {
delete Nodes; }
1152 for (
int i = 0; i < NumOfElements; i++)
1154 FreeElement(elements[i]);
1157 for (
int i = 0; i < NumOfBdrElements; i++)
1159 FreeElement(boundary[i]);
1162 for (
int i = 0; i < faces.Size(); i++)
1164 FreeElement(faces[i]);
1170 void Mesh::Destroy()
1174 elements.DeleteAll();
1175 vertices.DeleteAll();
1176 boundary.DeleteAll();
1178 faces_info.DeleteAll();
1179 nc_faces_info.DeleteAll();
1180 be_to_edge.DeleteAll();
1181 be_to_face.DeleteAll();
1189 CoarseFineTr.Clear();
1191 #ifdef MFEM_USE_MEMALLOC
1195 attributes.DeleteAll();
1196 bdr_attributes.DeleteAll();
1199 void Mesh::ResetLazyData()
1201 delete el_to_el; el_to_el = NULL;
1202 delete face_edge; face_edge = NULL;
1203 delete edge_vertex; edge_vertex = NULL;
1204 DeleteGeometricFactors();
1205 nbInteriorFaces = -1;
1206 nbBoundaryFaces = -1;
1209 void Mesh::SetAttributes()
1214 for (
int i = 0; i < attribs.
Size(); i++)
1216 attribs[i] = GetBdrAttribute(i);
1220 attribs.
Copy(bdr_attributes);
1221 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1223 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1227 for (
int i = 0; i < attribs.
Size(); i++)
1229 attribs[i] = GetAttribute(i);
1233 attribs.
Copy(attributes);
1234 if (attributes.Size() > 0 && attributes[0] <= 0)
1236 MFEM_WARNING(
"Non-positive attributes in the domain!");
1240 void Mesh::InitMesh(
int Dim_,
int spaceDim_,
int NVert,
int NElem,
int NBdrElem)
1245 spaceDim = spaceDim_;
1248 vertices.SetSize(NVert);
1251 elements.SetSize(NElem);
1253 NumOfBdrElements = 0;
1254 boundary.SetSize(NBdrElem);
1257 template<
typename T>
1258 static void CheckEnlarge(
Array<T> &array,
int size)
1260 if (size >= array.
Size()) { array.
SetSize(size + 1); }
1263 int Mesh::AddVertex(
double x,
double y,
double z)
1265 CheckEnlarge(vertices, NumOfVertices);
1266 double *v = vertices[NumOfVertices]();
1270 return NumOfVertices++;
1273 int Mesh::AddVertex(
const double *coords)
1275 CheckEnlarge(vertices, NumOfVertices);
1276 vertices[NumOfVertices].SetCoords(spaceDim, coords);
1277 return NumOfVertices++;
1280 void Mesh::AddVertexParents(
int i,
int p1,
int p2)
1286 if (i < vertices.Size())
1288 double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1289 for (
int j = 0; j < 3; j++)
1291 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1296 int Mesh::AddSegment(
int v1,
int v2,
int attr)
1298 CheckEnlarge(elements, NumOfElements);
1299 elements[NumOfElements] =
new Segment(v1, v2, attr);
1300 return NumOfElements++;
1303 int Mesh::AddSegment(
const int *vi,
int attr)
1305 CheckEnlarge(elements, NumOfElements);
1306 elements[NumOfElements] =
new Segment(vi, attr);
1307 return NumOfElements++;
1310 int Mesh::AddTriangle(
int v1,
int v2,
int v3,
int attr)
1312 CheckEnlarge(elements, NumOfElements);
1313 elements[NumOfElements] =
new Triangle(v1, v2, v3, attr);
1314 return NumOfElements++;
1317 int Mesh::AddTriangle(
const int *vi,
int attr)
1319 CheckEnlarge(elements, NumOfElements);
1320 elements[NumOfElements] =
new Triangle(vi, attr);
1321 return NumOfElements++;
1324 int Mesh::AddQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1326 CheckEnlarge(elements, NumOfElements);
1327 elements[NumOfElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1328 return NumOfElements++;
1331 int Mesh::AddQuad(
const int *vi,
int attr)
1333 CheckEnlarge(elements, NumOfElements);
1335 return NumOfElements++;
1338 int Mesh::AddTet(
int v1,
int v2,
int v3,
int v4,
int attr)
1340 int vi[4] = {v1, v2, v3, v4};
1341 return AddTet(vi, attr);
1344 int Mesh::AddTet(
const int *vi,
int attr)
1346 CheckEnlarge(elements, NumOfElements);
1347 #ifdef MFEM_USE_MEMALLOC
1349 tet = TetMemory.Alloc();
1352 elements[NumOfElements] = tet;
1354 elements[NumOfElements] =
new Tetrahedron(vi, attr);
1356 return NumOfElements++;
1359 int Mesh::AddWedge(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int attr)
1361 CheckEnlarge(elements, NumOfElements);
1362 elements[NumOfElements] =
new Wedge(v1, v2, v3, v4, v5, v6, attr);
1363 return NumOfElements++;
1366 int Mesh::AddWedge(
const int *vi,
int attr)
1368 CheckEnlarge(elements, NumOfElements);
1369 elements[NumOfElements] =
new Wedge(vi, attr);
1370 return NumOfElements++;
1373 int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1376 CheckEnlarge(elements, NumOfElements);
1377 elements[NumOfElements] =
1378 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1379 return NumOfElements++;
1382 int Mesh::AddHex(
const int *vi,
int attr)
1384 CheckEnlarge(elements, NumOfElements);
1385 elements[NumOfElements] =
new Hexahedron(vi, attr);
1386 return NumOfElements++;
1389 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1391 static const int hex_to_tet[6][4] =
1393 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1394 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1398 for (
int i = 0; i < 6; i++)
1400 for (
int j = 0; j < 4; j++)
1402 ti[j] = vi[hex_to_tet[i][j]];
1408 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1410 static const int hex_to_wdg[2][6] =
1412 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1416 for (
int i = 0; i < 2; i++)
1418 for (
int j = 0; j < 6; j++)
1420 ti[j] = vi[hex_to_wdg[i][j]];
1428 CheckEnlarge(elements, NumOfElements);
1429 elements[NumOfElements] = elem;
1430 return NumOfElements++;
1435 CheckEnlarge(boundary, NumOfBdrElements);
1436 boundary[NumOfBdrElements] = elem;
1437 return NumOfBdrElements++;
1440 int Mesh::AddBdrSegment(
int v1,
int v2,
int attr)
1442 CheckEnlarge(boundary, NumOfBdrElements);
1443 boundary[NumOfBdrElements] =
new Segment(v1, v2, attr);
1444 return NumOfBdrElements++;
1447 int Mesh::AddBdrSegment(
const int *vi,
int attr)
1449 CheckEnlarge(boundary, NumOfBdrElements);
1450 boundary[NumOfBdrElements] =
new Segment(vi, attr);
1451 return NumOfBdrElements++;
1454 int Mesh::AddBdrTriangle(
int v1,
int v2,
int v3,
int attr)
1456 CheckEnlarge(boundary, NumOfBdrElements);
1457 boundary[NumOfBdrElements] =
new Triangle(v1, v2, v3, attr);
1458 return NumOfBdrElements++;
1461 int Mesh::AddBdrTriangle(
const int *vi,
int attr)
1463 CheckEnlarge(boundary, NumOfBdrElements);
1464 boundary[NumOfBdrElements] =
new Triangle(vi, attr);
1465 return NumOfBdrElements++;
1468 int Mesh::AddBdrQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1470 CheckEnlarge(boundary, NumOfBdrElements);
1471 boundary[NumOfBdrElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1472 return NumOfBdrElements++;
1475 int Mesh::AddBdrQuad(
const int *vi,
int attr)
1477 CheckEnlarge(boundary, NumOfBdrElements);
1479 return NumOfBdrElements++;
1482 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1484 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1487 for (
int i = 0; i < 2; i++)
1489 for (
int j = 0; j < 3; j++)
1491 ti[j] = vi[quad_to_tri[i][j]];
1493 AddBdrTriangle(ti, attr);
1497 int Mesh::AddBdrPoint(
int v,
int attr)
1499 CheckEnlarge(boundary, NumOfBdrElements);
1500 boundary[NumOfBdrElements] =
new Point(&v, attr);
1501 return NumOfBdrElements++;
1504 void Mesh::GenerateBoundaryElements()
1507 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1511 for (i = 0; i < boundary.Size(); i++)
1513 FreeElement(boundary[i]);
1523 NumOfBdrElements = 0;
1524 for (i = 0; i < faces_info.Size(); i++)
1526 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1529 boundary.SetSize(NumOfBdrElements);
1530 be2face.
SetSize(NumOfBdrElements);
1531 for (j = i = 0; i < faces_info.Size(); i++)
1533 if (faces_info[i].Elem2No < 0)
1535 boundary[j] = faces[i]->Duplicate(
this);
1542 void Mesh::FinalizeCheck()
1544 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1545 vertices.Size() == 0,
1546 "incorrect number of vertices: preallocated: " << vertices.Size()
1547 <<
", actually added: " << NumOfVertices);
1548 MFEM_VERIFY(elements.Size() == NumOfElements,
1549 "incorrect number of elements: preallocated: " << elements.Size()
1550 <<
", actually added: " << NumOfElements);
1551 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1552 "incorrect number of boundary elements: preallocated: "
1553 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1556 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1559 CheckElementOrientation(fix_orientation);
1563 MarkTriMeshForRefinement();
1568 el_to_edge =
new Table;
1569 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1571 CheckBdrElementOrientation();
1585 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1586 bool fix_orientation)
1589 if (fix_orientation)
1591 CheckElementOrientation(fix_orientation);
1596 el_to_edge =
new Table;
1597 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1599 CheckBdrElementOrientation();
1619 GeckoProgress(
double limit) : limit(limit) { sw.Start(); }
1620 virtual bool quit()
const {
return limit > 0 && sw.UserTime() > limit; }
1623 class GeckoVerboseProgress :
public GeckoProgress
1629 GeckoVerboseProgress(
double limit) : GeckoProgress(limit) {}
1631 virtual void beginorder(
const Graph* graph,
Float cost)
const
1632 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
1633 virtual void endorder(
const Graph* graph,
Float cost)
const
1634 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
1636 virtual void beginiter(
const Graph* graph,
1639 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
1640 << window << std::flush;
1642 virtual void enditer(
const Graph* graph,
Float mincost,
Float cost)
const
1643 {
mfem::out <<
", cost = " << cost << endl; }
1648 int iterations,
int window,
1649 int period,
int seed,
bool verbose,
1655 GeckoProgress progress(time_limit);
1656 GeckoVerboseProgress vprogress(time_limit);
1659 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1666 const Table &my_el_to_el = ElementToElementTable();
1667 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1669 const int *neighid = my_el_to_el.
GetRow(elemid);
1670 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
1672 graph.
insert_arc(elemid + 1, neighid[i] + 1);
1677 graph.
order(&functional, iterations, window, period, seed,
1678 verbose ? &vprogress : &progress);
1684 ordering[gnodeid - 1] = graph.
rank(gnodeid);
1687 return graph.
cost();
1698 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
1699 : coord(coord), dir(dir), points(points), mid(mid) {}
1701 bool operator()(
int i)
const
1703 return (points[3*i + coord] < mid) != dir;
1707 static void HilbertSort2D(
int coord1,
1710 const Array<double> &points,
int *beg,
int *end,
1711 double xmin,
double ymin,
double xmax,
double ymax)
1713 if (end - beg <= 1) {
return; }
1715 double xmid = (xmin + xmax)*0.5;
1716 double ymid = (ymin + ymax)*0.5;
1718 int coord2 = (coord1 + 1) % 2;
1721 int *p0 = beg, *p4 = end;
1722 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
1723 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
1724 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
1728 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
1729 ymin, xmin, ymid, xmid);
1731 if (p1 != p0 || p2 != p4)
1733 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
1734 xmin, ymid, xmid, ymax);
1736 if (p2 != p0 || p3 != p4)
1738 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
1739 xmid, ymid, xmax, ymax);
1743 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
1744 ymid, xmax, ymin, xmid);
1748 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
1749 const Array<double> &points,
int *beg,
int *end,
1750 double xmin,
double ymin,
double zmin,
1751 double xmax,
double ymax,
double zmax)
1753 if (end - beg <= 1) {
return; }
1755 double xmid = (xmin + xmax)*0.5;
1756 double ymid = (ymin + ymax)*0.5;
1757 double zmid = (zmin + zmax)*0.5;
1759 int coord2 = (coord1 + 1) % 3;
1760 int coord3 = (coord1 + 2) % 3;
1763 int *p0 = beg, *p8 = end;
1764 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
1765 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
1766 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
1767 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
1768 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
1769 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
1770 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
1774 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
1775 zmin, xmin, ymin, zmid, xmid, ymid);
1777 if (p1 != p0 || p2 != p8)
1779 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
1780 ymin, zmid, xmin, ymid, zmax, xmid);
1782 if (p2 != p0 || p3 != p8)
1784 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
1785 ymid, zmid, xmin, ymax, zmax, xmid);
1787 if (p3 != p0 || p4 != p8)
1789 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
1790 xmin, ymax, zmid, xmid, ymid, zmin);
1792 if (p4 != p0 || p5 != p8)
1794 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
1795 xmid, ymax, zmid, xmax, ymid, zmin);
1797 if (p5 != p0 || p6 != p8)
1799 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
1800 ymax, zmid, xmax, ymid, zmax, xmid);
1802 if (p6 != p0 || p7 != p8)
1804 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
1805 ymid, zmid, xmax, ymin, zmax, xmid);
1809 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
1810 zmid, xmax, ymin, zmin, xmid, ymid);
1816 MFEM_VERIFY(spaceDim <= 3,
"");
1819 GetBoundingBox(min, max);
1824 if (spaceDim < 3) { points = 0.0; }
1827 for (
int i = 0; i < GetNE(); i++)
1829 GetElementCenter(i, center);
1830 for (
int j = 0; j < spaceDim; j++)
1832 points[3*i + j] = center(j);
1839 indices.
Sort([&](
int a,
int b)
1840 {
return points[3*
a] < points[3*
b]; });
1842 else if (spaceDim == 2)
1845 HilbertSort2D(0,
false,
false,
1846 points, indices.
begin(), indices.
end(),
1847 min(0), min(1), max(0), max(1));
1852 HilbertSort3D(0,
false,
false,
false,
1853 points, indices.
begin(), indices.
end(),
1854 min(0), min(1), min(2), max(0), max(1), max(2));
1859 for (
int i = 0; i < GetNE(); i++)
1861 ordering[indices[i]] = i;
1866 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
1870 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
1875 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
1879 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
1909 old_elem_node_vals.SetSize(GetNE());
1910 nodes_fes = Nodes->FESpace();
1913 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1915 nodes_fes->GetElementVDofs(old_elid, old_dofs);
1917 old_elem_node_vals[old_elid] =
new Vector(vals);
1923 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
1925 int new_elid = ordering[old_elid];
1926 new_elements[new_elid] = elements[old_elid];
1931 if (reorder_vertices)
1936 vertex_ordering = -1;
1938 int new_vertex_ind = 0;
1939 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1941 int *elem_vert = elements[new_elid]->GetVertices();
1942 int nv = elements[new_elid]->GetNVertices();
1943 for (
int vi = 0; vi < nv; ++vi)
1945 int old_vertex_ind = elem_vert[vi];
1946 if (vertex_ordering[old_vertex_ind] == -1)
1948 vertex_ordering[old_vertex_ind] = new_vertex_ind;
1949 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
1959 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1961 int *elem_vert = elements[new_elid]->GetVertices();
1962 int nv = elements[new_elid]->GetNVertices();
1963 for (
int vi = 0; vi < nv; ++vi)
1965 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
1970 for (
int belid = 0; belid < GetNBE(); ++belid)
1972 int *be_vert = boundary[belid]->GetVertices();
1973 int nv = boundary[belid]->GetNVertices();
1974 for (
int vi = 0; vi < nv; ++vi)
1976 be_vert[vi] = vertex_ordering[be_vert[vi]];
1987 el_to_edge =
new Table;
1988 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1993 GetElementToFaceTable();
2003 last_operation = Mesh::NONE;
2004 nodes_fes->Update(
false);
2007 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
2009 int new_elid = ordering[old_elid];
2010 nodes_fes->GetElementVDofs(new_elid, new_dofs);
2011 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
2012 delete old_elem_node_vals[old_elid];
2018 void Mesh::MarkForRefinement()
2024 MarkTriMeshForRefinement();
2028 DSTable v_to_v(NumOfVertices);
2029 GetVertexToVertexTable(v_to_v);
2030 MarkTetMeshForRefinement(v_to_v);
2035 void Mesh::MarkTriMeshForRefinement()
2040 for (
int i = 0; i < NumOfElements; i++)
2042 if (elements[i]->GetType() == Element::TRIANGLE)
2044 GetPointMatrix(i, pmat);
2045 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
2056 for (
int i = 0; i < NumOfVertices; i++)
2061 length_idx[j].one = GetLength(i, it.Column());
2062 length_idx[j].two = j;
2069 for (
int i = 0; i < NumOfEdges; i++)
2071 order[length_idx[i].two] = i;
2075 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
2080 GetEdgeOrdering(v_to_v, order);
2082 for (
int i = 0; i < NumOfElements; i++)
2084 if (elements[i]->GetType() == Element::TETRAHEDRON)
2086 elements[i]->MarkEdge(v_to_v, order);
2089 for (
int i = 0; i < NumOfBdrElements; i++)
2091 if (boundary[i]->GetType() == Element::TRIANGLE)
2093 boundary[i]->MarkEdge(v_to_v, order);
2100 if (*old_v_to_v && *old_elem_vert)
2107 if (*old_v_to_v == NULL)
2109 bool need_v_to_v =
false;
2111 for (
int i = 0; i < GetNEdges(); i++)
2117 if (dofs.
Size() > 0)
2125 *old_v_to_v =
new DSTable(NumOfVertices);
2126 GetVertexToVertexTable(*(*old_v_to_v));
2129 if (*old_elem_vert == NULL)
2131 bool need_elem_vert =
false;
2133 for (
int i = 0; i < GetNE(); i++)
2139 if (dofs.
Size() > 1)
2141 need_elem_vert =
true;
2147 *old_elem_vert =
new Table;
2148 (*old_elem_vert)->
MakeI(GetNE());
2149 for (
int i = 0; i < GetNE(); i++)
2151 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2153 (*old_elem_vert)->MakeJ();
2154 for (
int i = 0; i < GetNE(); i++)
2156 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2157 elements[i]->GetNVertices());
2159 (*old_elem_vert)->ShiftUpI();
2172 const int num_edge_dofs = old_dofs.
Size();
2175 const Vector onodes = *Nodes;
2179 int offset = NumOfVertices * old_dofs.
Size();
2183 if (num_edge_dofs > 0)
2185 DSTable new_v_to_v(NumOfVertices);
2186 GetVertexToVertexTable(new_v_to_v);
2188 for (
int i = 0; i < NumOfVertices; i++)
2192 const int old_i = (*old_v_to_v)(i, it.Column());
2193 const int new_i = it.Index();
2194 if (new_i == old_i) {
continue; }
2196 old_dofs.
SetSize(num_edge_dofs);
2197 new_dofs.
SetSize(num_edge_dofs);
2198 for (
int j = 0; j < num_edge_dofs; j++)
2200 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2201 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2205 for (
int j = 0; j < old_dofs.
Size(); j++)
2207 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2211 offset += NumOfEdges * num_edge_dofs;
2219 Table old_face_vertex;
2220 old_face_vertex.
MakeI(NumOfFaces);
2221 for (
int i = 0; i < NumOfFaces; i++)
2225 old_face_vertex.
MakeJ();
2226 for (
int i = 0; i < NumOfFaces; i++)
2228 faces[i]->GetNVertices());
2232 STable3D *faces_tbl = GetElementToFaceTable(1);
2238 for (
int i = 0; i < NumOfFaces; i++)
2240 const int *old_v = old_face_vertex.
GetRow(i);
2242 switch (old_face_vertex.
RowSize(i))
2245 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2249 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2253 new_fdofs[new_i+1] = old_dofs.
Size();
2258 for (
int i = 0; i < NumOfFaces; i++)
2260 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2263 switch (old_face_vertex.
RowSize(i))
2266 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2267 new_v = faces[new_i]->GetVertices();
2268 new_or = GetTriOrientation(old_v, new_v);
2273 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2274 new_v = faces[new_i]->GetVertices();
2275 new_or = GetQuadOrientation(old_v, new_v);
2282 for (
int j = 0; j < old_dofs.
Size(); j++)
2285 const int old_j = dof_ord[j];
2286 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2290 for (
int j = 0; j < old_dofs.
Size(); j++)
2292 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2314 for (
int i = 0; i < GetNE(); i++)
2316 const int *old_v = old_elem_vert->
GetRow(i);
2317 const int *new_v = elements[i]->GetVertices();
2323 case Geometry::SEGMENT:
2324 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2326 case Geometry::TRIANGLE:
2327 new_or = GetTriOrientation(old_v, new_v);
2329 case Geometry::SQUARE:
2330 new_or = GetQuadOrientation(old_v, new_v);
2332 case Geometry::TETRAHEDRON:
2333 new_or = GetTetOrientation(old_v, new_v);
2337 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2338 <<
" FE collection) are not supported yet!");
2342 MFEM_VERIFY(dof_ord != NULL,
2343 "FE collection '" << fec->
Name()
2344 <<
"' does not define reordering for "
2345 << Geometry::Name[
geom] <<
" elements!");
2348 for (
int j = 0; j < new_dofs.
Size(); j++)
2351 const int old_j = dof_ord[j];
2352 new_dofs[old_j] = offset + j;
2354 offset += new_dofs.
Size();
2357 for (
int j = 0; j < old_dofs.
Size(); j++)
2359 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2371 GetElementToFaceTable();
2374 CheckBdrElementOrientation();
2379 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2384 CheckBdrElementOrientation();
2389 last_operation = Mesh::NONE;
2394 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2397 CheckElementOrientation(fix_orientation);
2399 if (NumOfBdrElements == 0)
2401 GetElementToFaceTable();
2403 GenerateBoundaryElements();
2408 DSTable v_to_v(NumOfVertices);
2409 GetVertexToVertexTable(v_to_v);
2410 MarkTetMeshForRefinement(v_to_v);
2413 GetElementToFaceTable();
2416 CheckBdrElementOrientation();
2418 if (generate_edges == 1)
2420 el_to_edge =
new Table;
2421 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2435 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2436 bool fix_orientation)
2439 CheckElementOrientation(fix_orientation);
2441 if (NumOfBdrElements == 0)
2443 GetElementToFaceTable();
2445 GenerateBoundaryElements();
2448 GetElementToFaceTable();
2451 CheckBdrElementOrientation();
2453 if (generate_edges == 1)
2455 el_to_edge =
new Table;
2456 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2470 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2473 CheckElementOrientation(fix_orientation);
2475 GetElementToFaceTable();
2478 if (NumOfBdrElements == 0)
2480 GenerateBoundaryElements();
2483 CheckBdrElementOrientation();
2487 el_to_edge =
new Table;
2488 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2500 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2504 Finalize(refine, fix_orientation);
2507 void Mesh::FinalizeTopology(
bool generate_bdr)
2519 bool generate_edges =
true;
2521 if (spaceDim == 0) { spaceDim = Dim; }
2522 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2526 if (tmp_vertex_parents.Size())
2528 MFEM_VERIFY(ncmesh == NULL,
"");
2529 ncmesh =
new NCMesh(
this);
2533 InitFromNCMesh(*ncmesh);
2534 ncmesh->OnMeshUpdated(
this);
2535 GenerateNCFaceInfo();
2539 tmp_vertex_parents.DeleteAll();
2549 GetElementToFaceTable();
2551 if (NumOfBdrElements == 0 && generate_bdr)
2553 GenerateBoundaryElements();
2554 GetElementToFaceTable();
2563 if (Dim > 1 && generate_edges)
2566 if (!el_to_edge) { el_to_edge =
new Table; }
2567 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2571 if (NumOfBdrElements == 0 && generate_bdr)
2573 GenerateBoundaryElements();
2590 ncmesh->OnMeshUpdated(
this);
2593 GenerateNCFaceInfo();
2600 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2602 if (NURBSext || ncmesh)
2604 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
2605 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
2614 const bool check_orientation =
true;
2615 const bool curved = (Nodes != NULL);
2616 const bool may_change_topology =
2617 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
2618 ( check_orientation && fix_orientation &&
2619 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
2622 Table *old_elem_vert = NULL;
2624 if (curved && may_change_topology)
2626 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
2629 if (check_orientation)
2632 CheckElementOrientation(fix_orientation);
2636 MarkForRefinement();
2639 if (may_change_topology)
2643 DoNodeReorder(old_v_to_v, old_elem_vert);
2644 delete old_elem_vert;
2657 CheckBdrElementOrientation();
2662 if (Dim >= 2 && Dim == spaceDim)
2664 const int num_faces = GetNumFaces();
2665 for (
int i = 0; i < num_faces; i++)
2667 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
2668 faces_info[i].Elem2Inf%2 != 0,
"Invalid mesh topology."
2669 " Interior face with incompatible orientations.");
2676 double sx,
double sy,
double sz,
bool sfc_ordering)
2680 int NVert, NElem, NBdrElem;
2682 NVert = (nx+1) * (ny+1) * (nz+1);
2683 NElem = nx * ny * nz;
2684 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
2685 if (type == Element::TETRAHEDRON)
2690 else if (type == Element::WEDGE)
2693 NBdrElem += 2*nx*ny;
2696 InitMesh(3, 3, NVert, NElem, NBdrElem);
2702 for (z = 0; z <= nz; z++)
2704 coord[2] = ((double) z / nz) * sz;
2705 for (y = 0; y <= ny; y++)
2707 coord[1] = ((double) y / ny) * sy;
2708 for (x = 0; x <= nx; x++)
2710 coord[0] = ((double) x / nx) * sx;
2716 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
2719 if (sfc_ordering && type == Element::HEXAHEDRON)
2722 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
2723 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
2725 for (
int k = 0; k < nx*ny*nz; k++)
2731 ind[0] = VTX(x , y , z );
2732 ind[1] = VTX(x+1, y , z );
2733 ind[2] = VTX(x+1, y+1, z );
2734 ind[3] = VTX(x , y+1, z );
2735 ind[4] = VTX(x , y , z+1);
2736 ind[5] = VTX(x+1, y , z+1);
2737 ind[6] = VTX(x+1, y+1, z+1);
2738 ind[7] = VTX(x , y+1, z+1);
2745 for (z = 0; z < nz; z++)
2747 for (y = 0; y < ny; y++)
2749 for (x = 0; x < nx; x++)
2751 ind[0] = VTX(x , y , z );
2752 ind[1] = VTX(x+1, y , z );
2753 ind[2] = VTX(x+1, y+1, z );
2754 ind[3] = VTX(x , y+1, z );
2755 ind[4] = VTX(x , y , z+1);
2756 ind[5] = VTX(x+1, y , z+1);
2757 ind[6] = VTX(x+1, y+1, z+1);
2758 ind[7] = VTX( x, y+1, z+1);
2759 if (type == Element::TETRAHEDRON)
2761 AddHexAsTets(ind, 1);
2763 else if (type == Element::WEDGE)
2765 AddHexAsWedges(ind, 1);
2778 for (y = 0; y < ny; y++)
2780 for (x = 0; x < nx; x++)
2782 ind[0] = VTX(x , y , 0);
2783 ind[1] = VTX(x , y+1, 0);
2784 ind[2] = VTX(x+1, y+1, 0);
2785 ind[3] = VTX(x+1, y , 0);
2786 if (type == Element::TETRAHEDRON)
2788 AddBdrQuadAsTriangles(ind, 1);
2790 else if (type == Element::WEDGE)
2792 AddBdrQuadAsTriangles(ind, 1);
2801 for (y = 0; y < ny; y++)
2803 for (x = 0; x < nx; x++)
2805 ind[0] = VTX(x , y , nz);
2806 ind[1] = VTX(x+1, y , nz);
2807 ind[2] = VTX(x+1, y+1, nz);
2808 ind[3] = VTX(x , y+1, nz);
2809 if (type == Element::TETRAHEDRON)
2811 AddBdrQuadAsTriangles(ind, 6);
2813 else if (type == Element::WEDGE)
2815 AddBdrQuadAsTriangles(ind, 1);
2824 for (z = 0; z < nz; z++)
2826 for (y = 0; y < ny; y++)
2828 ind[0] = VTX(0 , y , z );
2829 ind[1] = VTX(0 , y , z+1);
2830 ind[2] = VTX(0 , y+1, z+1);
2831 ind[3] = VTX(0 , y+1, z );
2832 if (type == Element::TETRAHEDRON)
2834 AddBdrQuadAsTriangles(ind, 5);
2843 for (z = 0; z < nz; z++)
2845 for (y = 0; y < ny; y++)
2847 ind[0] = VTX(nx, y , z );
2848 ind[1] = VTX(nx, y+1, z );
2849 ind[2] = VTX(nx, y+1, z+1);
2850 ind[3] = VTX(nx, y , z+1);
2851 if (type == Element::TETRAHEDRON)
2853 AddBdrQuadAsTriangles(ind, 3);
2862 for (x = 0; x < nx; x++)
2864 for (z = 0; z < nz; z++)
2866 ind[0] = VTX(x , 0, z );
2867 ind[1] = VTX(x+1, 0, z );
2868 ind[2] = VTX(x+1, 0, z+1);
2869 ind[3] = VTX(x , 0, z+1);
2870 if (type == Element::TETRAHEDRON)
2872 AddBdrQuadAsTriangles(ind, 2);
2881 for (x = 0; x < nx; x++)
2883 for (z = 0; z < nz; z++)
2885 ind[0] = VTX(x , ny, z );
2886 ind[1] = VTX(x , ny, z+1);
2887 ind[2] = VTX(x+1, ny, z+1);
2888 ind[3] = VTX(x+1, ny, z );
2889 if (type == Element::TETRAHEDRON)
2891 AddBdrQuadAsTriangles(ind, 4);
2903 ofstream test_stream(
"debug.mesh");
2905 test_stream.close();
2914 double sx,
double sy,
2915 bool generate_edges,
bool sfc_ordering)
2924 if (type == Element::QUADRILATERAL)
2926 NumOfVertices = (nx+1) * (ny+1);
2927 NumOfElements = nx * ny;
2928 NumOfBdrElements = 2 * nx + 2 * ny;
2930 vertices.SetSize(NumOfVertices);
2931 elements.SetSize(NumOfElements);
2932 boundary.SetSize(NumOfBdrElements);
2939 for (j = 0; j < ny+1; j++)
2941 cy = ((double) j / ny) * sy;
2942 for (i = 0; i < nx+1; i++)
2944 cx = ((double) i / nx) * sx;
2945 vertices[k](0) = cx;
2946 vertices[k](1) = cy;
2955 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
2956 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
2958 for (k = 0; k < nx*ny; k++)
2962 ind[0] = i + j*(nx+1);
2963 ind[1] = i + 1 +j*(nx+1);
2964 ind[2] = i + 1 + (j+1)*(nx+1);
2965 ind[3] = i + (j+1)*(nx+1);
2972 for (j = 0; j < ny; j++)
2974 for (i = 0; i < nx; i++)
2976 ind[0] = i + j*(nx+1);
2977 ind[1] = i + 1 +j*(nx+1);
2978 ind[2] = i + 1 + (j+1)*(nx+1);
2979 ind[3] = i + (j+1)*(nx+1);
2988 for (i = 0; i < nx; i++)
2990 boundary[i] =
new Segment(i, i+1, 1);
2991 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
2994 for (j = 0; j < ny; j++)
2996 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
2997 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3001 else if (type == Element::TRIANGLE)
3003 NumOfVertices = (nx+1) * (ny+1);
3004 NumOfElements = 2 * nx * ny;
3005 NumOfBdrElements = 2 * nx + 2 * ny;
3007 vertices.SetSize(NumOfVertices);
3008 elements.SetSize(NumOfElements);
3009 boundary.SetSize(NumOfBdrElements);
3016 for (j = 0; j < ny+1; j++)
3018 cy = ((double) j / ny) * sy;
3019 for (i = 0; i < nx+1; i++)
3021 cx = ((double) i / nx) * sx;
3022 vertices[k](0) = cx;
3023 vertices[k](1) = cy;
3030 for (j = 0; j < ny; j++)
3032 for (i = 0; i < nx; i++)
3034 ind[0] = i + j*(nx+1);
3035 ind[1] = i + 1 + (j+1)*(nx+1);
3036 ind[2] = i + (j+1)*(nx+1);
3039 ind[1] = i + 1 + j*(nx+1);
3040 ind[2] = i + 1 + (j+1)*(nx+1);
3048 for (i = 0; i < nx; i++)
3050 boundary[i] =
new Segment(i, i+1, 1);
3051 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3054 for (j = 0; j < ny; j++)
3056 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3057 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3064 MFEM_ABORT(
"Unsupported element type.");
3068 CheckElementOrientation();
3070 if (generate_edges == 1)
3072 el_to_edge =
new Table;
3073 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3075 CheckBdrElementOrientation();
3084 attributes.Append(1);
3085 bdr_attributes.Append(1); bdr_attributes.Append(2);
3086 bdr_attributes.Append(3); bdr_attributes.Append(4);
3091 void Mesh::Make1D(
int n,
double sx)
3100 NumOfVertices = n + 1;
3102 NumOfBdrElements = 2;
3103 vertices.SetSize(NumOfVertices);
3104 elements.SetSize(NumOfElements);
3105 boundary.SetSize(NumOfBdrElements);
3108 for (j = 0; j < n+1; j++)
3110 vertices[j](0) = ((
double) j / n) * sx;
3114 for (j = 0; j < n; j++)
3116 elements[j] =
new Segment(j, j+1, 1);
3121 boundary[0] =
new Point(ind, 1);
3123 boundary[1] =
new Point(ind, 2);
3131 attributes.Append(1);
3132 bdr_attributes.Append(1); bdr_attributes.Append(2);
3135 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
3153 last_operation = Mesh::NONE;
3156 elements.SetSize(NumOfElements);
3157 for (
int i = 0; i < NumOfElements; i++)
3159 elements[i] = mesh.
elements[i]->Duplicate(
this);
3166 boundary.SetSize(NumOfBdrElements);
3167 for (
int i = 0; i < NumOfBdrElements; i++)
3169 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
3188 faces.SetSize(mesh.
faces.Size());
3189 for (
int i = 0; i < faces.Size(); i++)
3192 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
3226 if (dynamic_cast<const ParMesh*>(&mesh))
3238 if (mesh.
Nodes && copy_nodes)
3243 FiniteElementCollection::New(fec->
Name());
3247 Nodes->MakeOwner(fec_copy);
3248 *Nodes = *mesh.
Nodes;
3270 bool fix_orientation)
3274 if (!imesh) { MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n'); }
3275 else { mesh.
Load(imesh, generate_edges, refine, fix_orientation); }
3289 double sx,
double sy,
bool sfc_ordering)
3292 mesh.
Make2D(nx, ny, type, sx, sy, generate_edges, sfc_ordering);
3299 double sx,
double sy,
double sz,
bool sfc_ordering)
3302 mesh.
Make3D(nx, ny, nz, type, sx, sy, sz, sfc_ordering);
3311 ref_factors = ref_factor;
3325 bool fix_orientation)
3334 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
3338 Load(imesh, generate_edges, refine, fix_orientation);
3343 bool fix_orientation)
3346 Load(input, generate_edges, refine, fix_orientation);
3355 "Not enough vertices in external array : "
3356 "len_vertex_data = "<< len_vertex_data <<
", "
3359 if (vertex_data == (
double *)(
vertices.GetData()))
3361 MFEM_ASSERT(!
vertices.OwnsData(),
"invalid ownership");
3366 memcpy(vertex_data,
vertices.GetData(),
3375 int *element_attributes,
int num_elements,
3377 int *boundary_attributes,
int num_boundary_elements,
3378 int dimension,
int space_dimension)
3380 if (space_dimension == -1)
3382 space_dimension = dimension;
3385 InitMesh(dimension, space_dimension, 0, num_elements,
3386 num_boundary_elements);
3389 int boundary_index_stride = num_boundary_elements > 0 ?
3393 vertices.MakeRef(reinterpret_cast<Vertex*>(vertices_), num_vertices);
3396 for (
int i = 0; i < num_elements; i++)
3399 elements[i]->SetVertices(element_indices + i * element_index_stride);
3400 elements[i]->SetAttribute(element_attributes[i]);
3404 for (
int i = 0; i < num_boundary_elements; i++)
3407 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3408 boundary[i]->SetAttribute(boundary_attributes[i]);
3424 #ifdef MFEM_USE_MEMALLOC
3432 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3445 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3448 for (
int i = 0; i < nv; i++)
3461 for (
int j = 0; j < nv; j++)
3524 MFEM_ABORT(
"invalid element type: " << type);
3531 std::string parse_tag)
3533 int curved = 0, read_gf = 1;
3534 bool finalize_topo =
true;
3538 MFEM_ABORT(
"Input stream is not open");
3545 getline(input, mesh_type);
3549 int mfem_version = 0;
3550 if (mesh_type ==
"MFEM mesh v1.0") { mfem_version = 10; }
3551 else if (mesh_type ==
"MFEM mesh v1.2") { mfem_version = 12; }
3555 int mfem_nc_version = 0;
3556 if (mesh_type ==
"MFEM NC mesh v1.0") { mfem_nc_version = 10; }
3557 else if (mesh_type ==
"MFEM mesh v1.1") { mfem_nc_version = 1 ; }
3565 if (mfem_version == 12 && parse_tag.empty())
3567 parse_tag =
"mfem_mesh_end";
3571 else if (mfem_nc_version)
3573 MFEM_ASSERT(
ncmesh == NULL,
"internal error");
3580 MFEM_VERIFY(mfem_nc_version >= 10,
3581 "Legacy nonconforming format (MFEM mesh v1.1) cannot be "
3582 "used to load a parallel nonconforming mesh, sorry.");
3585 input, mfem_nc_version, curved, is_nc);
3590 ncmesh =
new NCMesh(input, mfem_nc_version, curved, is_nc);
3603 else if (mesh_type ==
"linemesh")
3607 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
3609 if (mesh_type ==
"curved_areamesh2")
3615 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
3619 else if (mesh_type ==
"TrueGrid")
3623 else if (mesh_type.rfind(
"# vtk DataFile Version") == 0)
3625 int major_vtk_version = mesh_type[mesh_type.length()-3] -
'0';
3627 MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
3628 "Unsupported VTK format");
3629 ReadVTKMesh(input, curved, read_gf, finalize_topo);
3631 else if (mesh_type.rfind(
"<VTKFile ") == 0 || mesh_type.rfind(
"<?xml") == 0)
3635 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
3639 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
3644 else if (mesh_type ==
"$MeshFormat")
3649 ((mesh_type.size() > 2 &&
3650 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
3651 (mesh_type.size() > 3 &&
3652 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
3657 #ifdef MFEM_USE_NETCDF
3660 MFEM_ABORT(
"NetCDF support requires configuration with"
3661 " MFEM_USE_NETCDF=YES");
3667 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
3668 " Use mfem::named_ifgzstream for input.");
3674 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
3700 bool generate_bdr =
false;
3705 if (curved && read_gf)
3719 if (mfem_version == 12)
3725 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
3726 getline(input, line);
3732 if (line ==
"mfem_mesh_end") {
break; }
3734 while (line != parse_tag);
3736 else if (mfem_nc_version >= 10)
3741 MFEM_VERIFY(ident ==
"mfem_mesh_end",
3742 "invalid mesh: end of file tag not found");
3750 int i, j, ie, ib, iv, *v, nv;
3781 for (i = 0; i < num_pieces; i++)
3788 for (i = 0; i < num_pieces; i++)
3794 for (j = 0; j < m->
GetNE(); j++)
3799 for (j = 0; j < m->
GetNBE(); j++)
3804 for (
int k = 0; k < nv; k++)
3806 v[k] = lvert_vert[v[k]];
3811 for (j = 0; j < m->
GetNV(); j++)
3823 for (i = 0; i < num_pieces; i++)
3834 for (i = 0; i < num_pieces; i++)
3838 for (j = 0; j < m->
GetNE(); j++)
3843 for (
int k = 0; k < nv; k++)
3850 for (j = 0; j < m->
GetNBE(); j++)
3855 for (
int k = 0; k < nv; k++)
3862 for (j = 0; j < m->
GetNV(); j++)
3876 for (i = 0; i < num_pieces; i++)
3878 gf_array[i] = mesh_array[i]->
GetNodes();
3893 ref_factors = ref_factor;
3904 int orig_ne = orig_mesh.
GetNE();
3905 MFEM_VERIFY(ref_factors.
Size() == orig_ne && orig_ne > 0,
3906 "Number of refinement factors must equal number of elements")
3907 MFEM_VERIFY(ref_factors.
Min() >= 1,
"Refinement factor must be >= 1");
3910 "Invalid refinement type. Must use closed basis type.");
3912 int min_ref = ref_factors.
Min();
3913 int max_ref = ref_factors.
Max();
3915 bool var_order = (min_ref != max_ref);
3928 for (
int i = 0; i < orig_ne; i++)
3946 for (
int el = 0; el < orig_ne; el++)
3958 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[el]);
3959 for (
int i = 0; i < phys_pts.
Width(); i++)
3968 for (
int k = 0; k < nvert; k++)
3971 v[k] = rdofs[c2h_map[cid]];
3978 for (
int el = 0; el < orig_mesh.
GetNBE(); el++)
3989 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[i]);
3995 for (
int k = 0; k < nvert; k++)
3998 v[k] = rdofs[c2h_map[cid]];
4020 for (
int iel = 0; iel < orig_ne; iel++)
4029 const int *node_map = NULL;
4032 if (h1_fec != NULL) { node_map = h1_fec->
GetDofMap(geom); }
4033 const int *vertex_map = vertex_fec.
GetDofMap(geom);
4034 const int *c2h_map = rfec.
GetDofMap(geom, ref_factors[iel]);
4035 for (
int jel = 0; jel < RG.
RefGeoms.
Size()/nvert; jel++)
4038 for (
int iv_lex=0; iv_lex<nvert; ++iv_lex)
4041 int iv = vertex_map[iv_lex];
4043 int pt_idx = c2h_map[RG.
RefGeoms[iv+nvert*jel]];
4045 int node_idx = node_map ? node_map[iv_lex] : iv_lex;
4048 (*Nodes)[dofs[node_idx + d*nvert]] = phys_pts(d,pt_idx);
4059 using GeomRef = std::pair<Geometry::Type, int>;
4060 std::map<GeomRef, int> point_matrices_offsets;
4062 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4066 GeomRef id(geom, ref_factors[el_coarse]);
4067 if (point_matrices_offsets.find(
id) == point_matrices_offsets.end())
4073 point_matrices_offsets[id] = n_point_matrices[
geom];
4074 n_point_matrices[
geom] += nref_el;
4081 int nmatrices = n_point_matrices[
geom];
4088 for (
int el_coarse = 0; el_coarse < orig_ne; ++el_coarse)
4091 int ref = ref_factors[el_coarse];
4092 int offset = point_matrices_offsets[GeomRef(geom, ref)];
4098 for (
int k = 0; k < nvert; k++)
4130 "Mesh::MakeSimplicial requires a properly oriented input mesh");
4132 "Mesh::MakeSimplicial does not support non-conforming meshes.")
4139 Mesh copy(orig_mesh);
4144 int nv = orig_mesh.
GetNV();
4145 int ne = orig_mesh.
GetNE();
4146 int nbe = orig_mesh.
GetNBE();
4159 int new_ne = 0, new_nbe = 0;
4160 for (
int i=0; i<ne; ++i)
4164 for (
int i=0; i<nbe; ++i)
4173 for (
int i=0; i<nv; ++i)
4183 if (vglobal == NULL)
4186 for (
int i=0; i<nv; ++i) { vglobal_id[i] = i; }
4187 vglobal = vglobal_id.
GetData();
4190 constexpr
int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
4191 constexpr
int quad_ntris = 2, prism_ntets = 3;
4192 static const int quad_trimap[2][nv_tri*quad_ntris] =
4204 static const int prism_rot[nv_prism*nv_prism] =
4213 static const int prism_f[nv_quad] = {1, 2, 5, 4};
4214 static const int prism_tetmaps[2][nv_prism*prism_ntets] =
4228 static const int hex_rot[nv_hex*nv_hex] =
4230 0, 1, 2, 3, 4, 5, 6, 7,
4231 1, 0, 4, 5, 2, 3, 7, 6,
4232 2, 1, 5, 6, 3, 0, 4, 7,
4233 3, 0, 1, 2, 7, 4, 5, 6,
4234 4, 0, 3, 7, 5, 1, 2, 6,
4235 5, 1, 0, 4, 6, 2, 3, 7,
4236 6, 2, 1, 5, 7, 3, 0, 4,
4237 7, 3, 2, 6, 4, 0, 1, 5
4239 static const int hex_f0[nv_quad] = {1, 2, 6, 5};
4240 static const int hex_f1[nv_quad] = {2, 3, 7, 6};
4241 static const int hex_f2[nv_quad] = {4, 5, 6, 7};
4242 static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
4243 static const int hex_tetmap0[nv_tet*5] =
4250 static const int hex_tetmap1[nv_tet*6] =
4257 static const int hex_tetmap2[nv_tet*6] =
4264 static const int hex_tetmap3[nv_tet*6] =
4271 static const int *hex_tetmaps[4] =
4273 hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
4276 auto find_min = [](
const int*
a,
int n) {
return std::min_element(a,a+n)-
a; };
4278 for (
int i=0; i<ne; ++i)
4280 const int *v = orig_mesh.
elements[i]->GetVertices();
4284 if (num_subdivisions[orig_geom] == 1)
4296 for (
int itri=0; itri<quad_ntris; ++itri)
4301 for (
int iv=0; iv<nv_tri; ++iv)
4303 v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
4311 for (
int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
4314 int irot = find_min(vg, nv_prism);
4315 for (
int iv=0; iv<nv_prism; ++iv)
4317 int jv = prism_rot[iv + irot*nv_prism];
4322 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
4323 int j = find_min(q, nv_quad);
4324 const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
4325 for (
int itet=0; itet<prism_ntets; ++itet)
4330 for (
int iv=0; iv<nv_tet; ++iv)
4332 v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
4340 for (
int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
4344 int irot = find_min(vg, nv_hex);
4345 for (
int iv=0; iv<nv_hex; ++iv)
4347 int jv = hex_rot[iv + irot*nv_hex];
4357 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
4358 j = find_min(q, nv_quad);
4359 if (j == 0 || j == 2) { bitmask += 4; }
4361 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
4362 j = find_min(q, nv_quad);
4363 if (j == 1 || j == 3) { bitmask += 2; }
4365 for (
int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
4366 j = find_min(q, nv_quad);
4367 if (j == 0 || j == 2) { bitmask += 1; }
4370 int nrot = num_rot[bitmask];
4371 for (
int irot=0; irot<nrot; ++irot)
4385 int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
4386 int ntets = (ndiags == 0) ? 5 : 6;
4387 const int *tetmap = hex_tetmaps[ndiags];
4388 for (
int itet=0; itet<ntets; ++itet)
4393 for (
int iv=0; iv<nv_tet; ++iv)
4395 v2[iv] = vg[tetmap[itet + iv*ntets]];
4404 for (
int i=0; i<nbe; ++i)
4406 const int *v = orig_mesh.
boundary[i]->GetVertices();
4409 if (num_subdivisions[orig_geom] == 1)
4419 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
4421 int iv_min = find_min(vg, nv_quad);
4422 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
4423 for (
int itri=0; itri<quad_ntris; ++itri)
4428 for (
int iv=0; iv<nv_tri; ++iv)
4430 v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
4437 MFEM_ABORT(
"Unreachable");
4451 Mesh periodic_mesh(orig_mesh,
true);
4457 for (
int i = 0; i < periodic_mesh.
GetNE(); i++)
4462 for (
int j = 0; j < nv; j++)
4468 for (
int i = 0; i < periodic_mesh.
GetNBE(); i++)
4473 for (
int j = 0; j < nv; j++)
4480 return periodic_mesh;
4484 const std::vector<Vector> &translations,
double tol)
const
4488 Vector coord(sdim), at(sdim), dx(sdim);
4489 Vector xMax(sdim), xMin(sdim), xDiff(sdim);
4490 xMax = xMin = xDiff = 0.0;
4494 for (
int be = 0; be <
GetNBE(); be++)
4499 for (
int i = 0; i < dofs.
Size(); i++)
4501 bdr_v.insert(dofs[i]);
4504 for (
int j = 0; j < sdim; j++)
4506 xMax[j] = max(xMax[j], coord[j]);
4507 xMin[j] = min(xMin[j], coord[j]);
4511 add(xMax, -1.0, xMin, xDiff);
4512 double dia = xDiff.
Norml2();
4520 map<int, int> replica2primary;
4522 map<int, set<int>> primary2replicas;
4526 for (
int v : bdr_v) { primary2replicas[v]; }
4530 auto make_replica = [&replica2primary, &primary2replicas](
int r,
int p)
4532 if (r ==
p) {
return; }
4533 primary2replicas[
p].insert(r);
4534 replica2primary[r] =
p;
4535 for (
int s : primary2replicas[r])
4537 primary2replicas[
p].insert(
s);
4538 replica2primary[
s] =
p;
4540 primary2replicas.erase(r);
4543 for (
unsigned int i = 0; i < translations.size(); i++)
4545 for (
int vi : bdr_v)
4548 add(coord, translations[i], at);
4550 for (
int vj : bdr_v)
4553 add(at, -1.0, coord, dx);
4554 if (dx.
Norml2() > dia*tol) {
continue; }
4559 bool pi = primary2replicas.find(vi) != primary2replicas.end();
4560 bool pj = primary2replicas.find(vj) != primary2replicas.end();
4566 make_replica(vj, vi);
4571 int owner_of_vj = replica2primary[vj];
4573 make_replica(vi, owner_of_vj);
4579 int owner_of_vi = replica2primary[vi];
4580 make_replica(vj, owner_of_vi);
4587 int owner_of_vi = replica2primary[vi];
4588 int owner_of_vj = replica2primary[vj];
4589 make_replica(owner_of_vj, owner_of_vi);
4596 std::vector<int> v2v(
GetNV());
4597 for (
size_t i = 0; i < v2v.size(); i++)
4601 for (
auto &&r2p : replica2primary)
4603 v2v[r2p.first] = r2p.second;
4612 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
4617 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
4634 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
4639 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
4669 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
4693 for (
int i = 0; i <
elements.Size(); i++)
4703 for (
int i = 0; i <
boundary.Size(); i++)
4720 for (
int i = 0; i < vd; i++)
4774 boundary.SetSize(NumOfBdrElements);
4785 edge_to_knot.
SetSize(NumOfEdges);
4789 input >> edge_to_knot[j] >> v[0] >> v[1];
4792 edge_to_knot[j] = -1 - edge_to_knot[j];
4810 for (
int d = 0; d < v.
Size(); d++)
4818 for (d = 0; d < p.
Size(); d++)
4822 for ( ; d < v.
Size(); d++)
4854 if (dynamic_cast<const H1_FECollection*>(fec)
4855 || dynamic_cast<const L2_FECollection*>(fec))
4884 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
4903 MFEM_ASSERT(nodes != NULL,
"");
4919 case 1:
return GetNV();
4926 static int CountFacesByType(
const Mesh &mesh,
const FaceType type)
4945 if (nf<0) { nf = CountFacesByType(*
this, type); }
4949 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
4950 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
4955 int i, j, k, wo = 0, fo = 0;
4964 int *vi =
elements[i]->GetVertices();
4967 for (j = 0; j < 3; j++)
4971 for (j = 0; j < 2; j++)
4972 for (k = 0; k < 2; k++)
4974 J(j, k) = v[j+1][k] - v[0][k];
4995 MFEM_ABORT(
"Invalid 2D element type \""
5012 int *vi =
elements[i]->GetVertices();
5018 for (j = 0; j < 4; j++)
5022 for (j = 0; j < 3; j++)
5023 for (k = 0; k < 3; k++)
5025 J(j, k) = v[j+1][k] - v[0][k];
5071 MFEM_ABORT(
"Invalid 3D element type \""
5077 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
5080 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
5081 <<
NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
5096 if (test[0] == base[0])
5097 if (test[1] == base[1])
5105 else if (test[0] == base[1])
5106 if (test[1] == base[0])
5115 if (test[1] == base[0])
5126 for (
int j = 0; j < 3; j++)
5127 if (test[aor[j]] != base[j])
5140 for (i = 0; i < 4; i++)
5141 if (test[i] == base[0])
5148 if (test[(i+1)%4] == base[1])
5157 for (
int j = 0; j < 4; j++)
5158 if (test[aor[j]] != base[j])
5160 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
5162 for (
int k = 0; k < 4; k++)
5167 for (
int k = 0; k < 4; k++)
5176 if (test[(i+1)%4] == base[1])
5192 if (test[0] == base[0])
5193 if (test[1] == base[1])
5194 if (test[2] == base[2])
5202 else if (test[2] == base[1])
5203 if (test[3] == base[2])
5212 if (test[1] == base[2])
5220 else if (test[1] == base[0])
5221 if (test[2] == base[1])
5222 if (test[0] == base[2])
5230 else if (test[3] == base[1])
5231 if (test[2] == base[2])
5240 if (test[3] == base[2])
5248 else if (test[2] == base[0])
5249 if (test[3] == base[1])
5250 if (test[0] == base[2])
5258 else if (test[0] == base[1])
5259 if (test[1] == base[2])
5268 if (test[3] == base[2])
5277 if (test[0] == base[1])
5278 if (test[2] == base[2])
5286 else if (test[1] == base[1])
5287 if (test[0] == base[2])
5296 if (test[1] == base[2])
5307 for (
int j = 0; j < 4; j++)
5308 if (test[aor[j]] != base[j])
5333 int *bv =
boundary[i]->GetVertices();
5339 mfem::Swap<int>(bv[0], bv[1]);
5353 if (
faces_info[fi].Elem2No >= 0) {
continue; }
5356 int *bv =
boundary[i]->GetVertices();
5358 MFEM_ASSERT(fi <
faces.Size(),
"internal error");
5359 const int *fv =
faces[fi]->GetVertices();
5376 MFEM_ABORT(
"Invalid 2D boundary element type \""
5377 << bdr_type <<
"\"");
5382 if (orientation % 2 == 0) {
continue; }
5384 if (!fix_it) {
continue; }
5392 mfem::Swap<int>(bv[0], bv[1]);
5396 mfem::Swap<int>(be[1], be[2]);
5402 mfem::Swap<int>(bv[0], bv[2]);
5406 mfem::Swap<int>(be[0], be[1]);
5407 mfem::Swap<int>(be[2], be[3]);
5420 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
5430 MFEM_ASSERT(0 <= dim && dim <=
Dim,
"invalid dim: " << dim);
5441 MFEM_ASSERT(0 <= dim && dim <=
Dim,
"invalid dim: " << dim);
5460 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
5461 "is not generated.");
5464 const int *v =
elements[i]->GetVertices();
5465 const int ne =
elements[i]->GetNEdges();
5467 for (
int j = 0; j < ne; j++)
5469 const int *e =
elements[i]->GetEdgeVertices(j);
5470 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5481 const int *v =
boundary[i]->GetVertices();
5482 cor[0] = (v[0] < v[1]) ? (1) : (-1);
5495 const int *v =
boundary[i]->GetVertices();
5496 const int ne =
boundary[i]->GetNEdges();
5498 for (
int j = 0; j < ne; j++)
5500 const int *e =
boundary[i]->GetEdgeVertices(j);
5501 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5513 const int *v =
faces[i]->GetVertices();
5514 o[0] = (v[0] < v[1]) ? (1) : (-1);
5526 const int *v =
faces[i]->GetVertices();
5527 const int ne =
faces[i]->GetNEdges();
5529 for (
int j = 0; j < ne; j++)
5531 const int *e =
faces[i]->GetEdgeVertices(j);
5532 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
5560 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
5611 for (j = 0; j < nv; j++)
5623 for (j = 0; j < nv; j++)
5670 MFEM_VERIFY(
el_to_face != NULL,
"el_to_face not generated");
5674 int n = faces.
Size();
5677 for (
int j = 0; j < n; j++)
5685 MFEM_ASSERT(
faces_info[faces[j]].Elem2No == i,
"internal error");
5710 MFEM_ABORT(
"invalid geometry");
5718 case 1:
return boundary[i]->GetVertices()[0];
5721 default: MFEM_ABORT(
"invalid dimension!");
5731 MFEM_ASSERT(fi.
Elem1Inf % 64 == 0,
"internal error");
5734 const int *bv =
boundary[bdr_el]->GetVertices();
5742 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
5769 for (j = 0; j < nv; j++)
5771 pointmat(k, j) =
vertices[v[j]](k);
5786 for (j = 0; j < nv; j++)
5788 pointmat(k, j) =
vertices[v[j]](k);
5800 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
5803 return sqrt(length);
5811 for (
int i = 0; i < elem_array.
Size(); i++)
5816 for (
int i = 0; i < elem_array.
Size(); i++)
5818 const int *v = elem_array[i]->GetVertices();
5819 const int ne = elem_array[i]->GetNEdges();
5820 for (
int j = 0; j < ne; j++)
5822 const int *e = elem_array[i]->GetEdgeVertices(j);
5836 v_to_v.
Push(v[0], v[1]);
5843 const int *v =
elements[i]->GetVertices();
5844 const int ne =
elements[i]->GetNEdges();
5845 for (
int j = 0; j < ne; j++)
5847 const int *e =
elements[i]->GetEdgeVertices(j);
5848 v_to_v.
Push(v[e[0]], v[e[1]]);
5856 int i, NumberOfEdges;
5872 const int *v =
boundary[i]->GetVertices();
5873 be_to_f[i] = v_to_v(v[0], v[1]);
5886 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
5890 return NumberOfEdges;
5981 if (
faces[gf] == NULL)
5991 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
5992 "Interior edge found between 2D elements "
5994 <<
" and " << el <<
".");
5995 int *v =
faces[gf]->GetVertices();
5997 if ( v[1] == v0 && v[0] == v1 )
6001 else if ( v[0] == v0 && v[1] == v1 )
6011 MFEM_ABORT(
"internal error");
6017 int v0,
int v1,
int v2)
6019 if (
faces[gf] == NULL)
6029 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6030 "Interior triangular face found connecting elements "
6032 <<
" and " << el <<
".");
6033 int orientation, vv[3] = { v0, v1, v2 };
6040 faces_info[gf].Elem2Inf = 64 * lf + orientation;
6045 int v0,
int v1,
int v2,
int v3)
6057 MFEM_VERIFY(
faces_info[gf].Elem2No < 0,
"Invalid mesh topology. "
6058 "Interior quadrilateral face found connecting elements "
6060 <<
" and " << el <<
".");
6061 int vv[4] = { v0, v1, v2, v3 };
6075 for (i = 0; i <
faces.Size(); i++)
6081 faces.SetSize(nfaces);
6083 for (i = 0; i < nfaces; i++)
6091 const int *v =
elements[i]->GetVertices();
6101 const int ne =
elements[i]->GetNEdges();
6102 for (
int j = 0; j < ne; j++)
6104 const int *e =
elements[i]->GetEdgeVertices(j);
6115 for (
int j = 0; j < 4; j++)
6119 v[fv[0]], v[fv[1]], v[fv[2]]);
6125 for (
int j = 0; j < 2; j++)
6129 v[fv[0]], v[fv[1]], v[fv[2]]);
6131 for (
int j = 2; j < 5; j++)
6135 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6141 for (
int j = 0; j < 6; j++)
6145 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6150 MFEM_ABORT(
"Unexpected type of Element.");
6158 MFEM_VERIFY(
ncmesh,
"missing NCMesh.");
6174 for (
int i = 0; i < list.
masters.Size(); i++)
6177 if (master.
index >= nfaces) {
continue; }
6185 for (
int i = 0; i < list.
slaves.Size(); i++)
6189 if (slave.
index < 0 ||
6190 slave.
index >= nfaces ||
6219 const int *v =
elements[i]->GetVertices();
6224 for (
int j = 0; j < 4; j++)
6227 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6233 for (
int j = 0; j < 2; j++)
6236 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
6238 for (
int j = 2; j < 5; j++)
6241 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6249 for (
int j = 0; j < 6; j++)
6252 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
6257 MFEM_ABORT(
"Unexpected type of Element.");
6281 for (
int j = 0; j < 4; j++)
6285 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6291 for (
int j = 0; j < 2; j++)
6295 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
6297 for (
int j = 2; j < 5; j++)
6301 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6309 for (
int j = 0; j < 6; j++)
6313 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
6318 MFEM_ABORT(
"Unexpected type of Element.");
6331 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
6336 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
6340 MFEM_ABORT(
"Unexpected type of boundary Element.");
6354 void Rotate3(
int &
a,
int &
b,
int &c)
6386 Table *old_elem_vert = NULL;
6397 int *v =
elements[i]->GetVertices();
6399 Rotate3(v[0], v[1], v[2]);
6402 Rotate3(v[1], v[2], v[3]);
6415 int *v =
boundary[i]->GetVertices();
6417 Rotate3(v[0], v[1], v[2]);
6433 delete old_elem_vert;
6449 if (p[i] < pmin[i]) { pmin[i] = p[i]; }
6450 if (p[i] > pmax[i]) { pmax[i] = p[i]; }
6464 for (
int i =
spaceDim-1; i >= 0; i--)
6466 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
6467 if (idx < 0) { idx = 0; }
6468 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
6469 part = part * nxyz[i] + idx;
6471 partitioning[el] = part;
6474 return partitioning;
6479 #ifdef MFEM_USE_METIS
6481 int print_messages = 1;
6484 int init_flag, fin_flag;
6485 MPI_Initialized(&init_flag);
6486 MPI_Finalized(&fin_flag);
6487 if (init_flag && !fin_flag)
6491 if (rank != 0) { print_messages = 0; }
6495 int i, *partitioning;
6505 partitioning[i] = 0;
6512 partitioning[i] = i;
6518 #ifndef MFEM_USE_METIS_5
6530 bool freedata =
false;
6532 idx_t *mpartitioning;
6535 if (
sizeof(
idx_t) ==
sizeof(int))
6539 mpartitioning = (
idx_t*) partitioning;
6548 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
6549 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
6550 mpartitioning =
new idx_t[n];
6553 #ifndef MFEM_USE_METIS_5
6556 METIS_SetDefaultOptions(options);
6557 options[METIS_OPTION_CONTIG] = 1;
6561 if (part_method >= 0 && part_method <= 2)
6563 for (i = 0; i < n; i++)
6569 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
6575 if (part_method == 0 || part_method == 3)
6577 #ifndef MFEM_USE_METIS_5
6606 " error in METIS_PartGraphRecursive!");
6613 if (part_method == 1 || part_method == 4)
6615 #ifndef MFEM_USE_METIS_5
6644 " error in METIS_PartGraphKway!");
6651 if (part_method == 2 || part_method == 5)
6653 #ifndef MFEM_USE_METIS_5
6666 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
6683 " error in METIS_PartGraphKway!");
6691 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
6695 nparts = (int) mparts;
6696 if (mpartitioning != (
idx_t*)partitioning)
6700 partitioning[k] = mpartitioning[k];
6707 delete[] mpartitioning;
6723 auto count_partition_elements = [&]()
6725 for (i = 0; i < nparts; i++)
6733 psize[partitioning[i]].one++;
6737 for (i = 0; i < nparts; i++)
6739 if (psize[i].one == 0) { empty_parts++; }
6743 count_partition_elements();
6751 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
6752 << empty_parts <<
" empty parts!"
6753 <<
" Applying a simple fix ..." << endl;
6756 SortPairs<int,int>(psize, nparts);
6758 for (i = nparts-1; i > nparts-1-empty_parts; i--)
6765 for (i = nparts-1; i > nparts-1-empty_parts; i--)
6767 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
6773 partitioning[j] = psize[nparts-1-i].two;
6780 count_partition_elements();
6784 return partitioning;
6788 mfem_error(
"Mesh::GeneratePartitioning(...): "
6789 "MFEM was compiled without Metis.");
6803 int num_elem, *i_elem_elem, *j_elem_elem;
6805 num_elem = elem_elem.
Size();
6806 i_elem_elem = elem_elem.
GetI();
6807 j_elem_elem = elem_elem.
GetJ();
6812 int stack_p, stack_top_p, elem;
6816 for (i = 0; i < num_elem; i++)
6818 if (partitioning[i] > num_part)
6820 num_part = partitioning[i];
6827 for (i = 0; i < num_part; i++)
6834 for (elem = 0; elem < num_elem; elem++)
6836 if (component[elem] >= 0)
6841 component[elem] = num_comp[partitioning[elem]]++;
6843 elem_stack[stack_top_p++] = elem;
6845 for ( ; stack_p < stack_top_p; stack_p++)
6847 i = elem_stack[stack_p];
6848 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
6851 if (partitioning[k] == partitioning[i])
6853 if (component[k] < 0)
6855 component[k] = component[i];
6856 elem_stack[stack_top_p++] = k;
6858 else if (component[k] != component[i])
6870 int i, n_empty, n_mcomp;
6878 n_empty = n_mcomp = 0;
6879 for (i = 0; i < num_comp.
Size(); i++)
6880 if (num_comp[i] == 0)
6884 else if (num_comp[i] > 1)
6891 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
6892 <<
"The following subdomains are empty :\n";
6893 for (i = 0; i < num_comp.
Size(); i++)
6894 if (num_comp[i] == 0)
6902 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
6903 <<
"The following subdomains are NOT connected :\n";
6904 for (i = 0; i < num_comp.
Size(); i++)
6905 if (num_comp[i] > 1)
6911 if (n_empty == 0 && n_mcomp == 0)
6912 mfem::out <<
"Mesh::CheckPartitioning(...) : "
6913 "All subdomains are connected." << endl;
6927 const double *a = A.
Data();
6928 const double *b = B.
Data();
6937 c(0) = a[0]*a[3]-a[1]*a[2];
6938 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
6939 c(2) = b[0]*b[3]-b[1]*b[2];
6960 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
6961 a[1] * (a[5] * a[6] - a[3] * a[8]) +
6962 a[2] * (a[3] * a[7] - a[4] * a[6]));
6964 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
6965 b[1] * (a[5] * a[6] - a[3] * a[8]) +
6966 b[2] * (a[3] * a[7] - a[4] * a[6]) +
6968 a[0] * (b[4] * a[8] - b[5] * a[7]) +
6969 a[1] * (b[5] * a[6] - b[3] * a[8]) +
6970 a[2] * (b[3] * a[7] - b[4] * a[6]) +
6972 a[0] * (a[4] * b[8] - a[5] * b[7]) +
6973 a[1] * (a[5] * b[6] - a[3] * b[8]) +
6974 a[2] * (a[3] * b[7] - a[4] * b[6]));
6976 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
6977 a[1] * (b[5] * b[6] - b[3] * b[8]) +
6978 a[2] * (b[3] * b[7] - b[4] * b[6]) +
6980 b[0] * (a[4] * b[8] - a[5] * b[7]) +
6981 b[1] * (a[5] * b[6] - a[3] * b[8]) +
6982 b[2] * (a[3] * b[7] - a[4] * b[6]) +
6984 b[0] * (b[4] * a[8] - b[5] * a[7]) +
6985 b[1] * (b[5] * a[6] - b[3] * a[8]) +
6986 b[2] * (b[3] * a[7] - b[4] * a[6]));
6988 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
6989 b[1] * (b[5] * b[6] - b[3] * b[8]) +
6990 b[2] * (b[3] * b[7] - b[4] * b[6]));
7036 double a = z(2), b = z(1), c = z(0);
7037 double D = b*b-4*a*c;
7044 x(0) = x(1) = -0.5 * b /
a;
7049 x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
7057 t = -0.5 * (b + sqrt(D));
7061 t = -0.5 * (b - sqrt(D));
7067 Swap<double>(x(0), x(1));
7075 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
7078 double Q = (a * a - 3 *
b) / 9;
7079 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
7080 double Q3 = Q * Q * Q;
7087 x(0) = x(1) = x(2) = - a / 3;
7091 double sqrtQ = sqrt(Q);
7095 x(0) = -2 * sqrtQ - a / 3;
7096 x(1) = x(2) = sqrtQ - a / 3;
7100 x(0) = x(1) = - sqrtQ - a / 3;
7101 x(2) = 2 * sqrtQ - a / 3;
7108 double theta = acos(R / sqrt(Q3));
7109 double A = -2 * sqrt(Q);
7111 x0 = A * cos(theta / 3) - a / 3;
7112 x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
7113 x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
7118 Swap<double>(x0, x1);
7122 Swap<double>(x1, x2);
7125 Swap<double>(x0, x1);
7138 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
7142 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
7144 x(0) = A + Q / A - a / 3;
7153 const double factor,
const int Dim)
7155 const double c0 = c(0);
7156 c(0) = c0 * (1.0 - pow(factor, -Dim));
7158 for (
int j = 0; j < nr; j++)
7170 c(0) = c0 * (1.0 - pow(factor, Dim));
7172 for (
int j = 0; j < nr; j++)
7191 const double factor = 2.0;
7206 for (
int k = 0; k < nv; k++)
7209 V(j, k) = displacements(v[k]+j*nvs);
7239 for (
int j = 0; j < nv; j++)
7265 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7268 vertices[i](j) += displacements(j*nv+i);
7276 for (
int i = 0; i < nv; i++)
7279 vert_coord(j*nv+i) =
vertices[i](j);
7285 for (
int i = 0, nv =
vertices.Size(); i < nv; i++)
7288 vertices[i](j) = vert_coord(j*nv+i);
7299 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
7318 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
7335 (*Nodes) += displacements;
7347 node_coord = (*Nodes);
7359 (*Nodes) = node_coord;
7388 mfem::Swap<GridFunction*>(
Nodes, nodes);
7405 for (j = 1; j < n; j++)
7439 int quad_counter = 0;
7454 vertices.SetSize(oelem + quad_counter);
7455 new_elements.
SetSize(4 * NumOfElements);
7461 const int attr =
elements[i]->GetAttribute();
7462 int *v =
elements[i]->GetVertices();
7468 for (
int ei = 0; ei < 3; ei++)
7470 for (
int k = 0; k < 2; k++)
7478 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
7480 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
7482 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
7484 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
7488 const int qe = quad_counter;
7492 for (
int ei = 0; ei < 4; ei++)
7494 for (
int k = 0; k < 2; k++)
7502 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
7504 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
7506 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
7508 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
7512 MFEM_ABORT(
"unknown element type: " << el_type);
7522 const int attr =
boundary[i]->GetAttribute();
7523 int *v =
boundary[i]->GetVertices();
7532 static const double A = 0.0, B = 0.5, C = 1.0;
7533 static double tri_children[2*3*4] =
7540 static double quad_children[2*4*4] =
7554 for (
int i = 0; i <
elements.Size(); i++)
7575 if (!
Nodes || update_nodes)
7583 static inline double sqr(
const double &x)
7605 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
7608 int NumOfQuadFaces = 0;
7614 for (
int i = 0; i <
faces.Size(); i++)
7618 f2qf[i] = NumOfQuadFaces;
7625 NumOfQuadFaces =
faces.Size();
7629 int hex_counter = 0;
7632 for (
int i = 0; i <
elements.Size(); i++)
7648 DSTable *v_to_v_ptr = v_to_v_p;
7664 std::sort(row_start, J_v2v.
end());
7667 for (
int i = 0; i < J_v2v.
Size(); i++)
7669 e2v[J_v2v[i].two] = i;
7682 it.SetIndex(e2v[it.Index()]);
7692 const int oelem = oface + NumOfQuadFaces;
7697 vertices.SetSize(oelem + hex_counter);
7705 const int attr =
elements[i]->GetAttribute();
7706 int *v =
elements[i]->GetVertices();
7713 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
7721 for (
int ei = 0; ei < 6; ei++)
7723 for (
int k = 0; k < 2; k++)
7733 const int rt_algo = 1;
7744 double len_sqr, min_len;
7746 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
7747 sqr(J(1,0)-J(1,1)-J(1,2)) +
7748 sqr(J(2,0)-J(2,1)-J(2,2));
7751 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
7752 sqr(J(1,1)-J(1,0)-J(1,2)) +
7753 sqr(J(2,1)-J(2,0)-J(2,2));
7754 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
7756 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
7757 sqr(J(1,2)-J(1,0)-J(1,1)) +
7758 sqr(J(2,2)-J(2,0)-J(2,1));
7759 if (len_sqr < min_len) { rt = 2; }
7764 double Em_data[18], Js_data[9], Jp_data[9];
7767 double ar1, ar2,
kappa, kappa_min;
7769 for (
int s = 0;
s < 3;
s++)
7771 for (
int t = 0;
t < 3;
t++)
7773 Em(
t,
s) = 0.5*J(
t,
s);
7776 for (
int t = 0;
t < 3;
t++)
7778 Em(
t,3) = 0.5*(J(
t,0)+J(
t,1));
7779 Em(
t,4) = 0.5*(J(
t,0)+J(
t,2));
7780 Em(
t,5) = 0.5*(J(
t,1)+J(
t,2));
7784 for (
int t = 0;
t < 3;
t++)
7786 Js(
t,0) = Em(
t,5)-Em(
t,0);
7787 Js(
t,1) = Em(
t,1)-Em(
t,0);
7788 Js(
t,2) = Em(
t,2)-Em(
t,0);
7792 for (
int t = 0;
t < 3;
t++)
7794 Js(
t,0) = Em(
t,5)-Em(
t,0);
7795 Js(
t,1) = Em(
t,2)-Em(
t,0);
7796 Js(
t,2) = Em(
t,4)-Em(
t,0);
7800 kappa_min = std::max(ar1, ar2);
7804 for (
int t = 0;
t < 3;
t++)
7806 Js(
t,0) = Em(
t,0)-Em(
t,1);
7807 Js(
t,1) = Em(
t,4)-Em(
t,1);
7808 Js(
t,2) = Em(
t,2)-Em(
t,1);
7812 for (
int t = 0;
t < 3;
t++)
7814 Js(
t,0) = Em(
t,2)-Em(
t,1);
7815 Js(
t,1) = Em(
t,4)-Em(
t,1);
7816 Js(
t,2) = Em(
t,5)-Em(
t,1);
7820 kappa = std::max(ar1, ar2);
7821 if (kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
7824 for (
int t = 0;
t < 3;
t++)
7826 Js(
t,0) = Em(
t,0)-Em(
t,2);
7827 Js(
t,1) = Em(
t,1)-Em(
t,2);
7828 Js(
t,2) = Em(
t,3)-Em(
t,2);
7832 for (
int t = 0;
t < 3;
t++)
7834 Js(
t,0) = Em(
t,1)-Em(
t,2);
7835 Js(
t,1) = Em(
t,5)-Em(
t,2);
7836 Js(
t,2) = Em(
t,3)-Em(
t,2);
7840 kappa = std::max(ar1, ar2);
7841 if (kappa < kappa_min) { rt = 2; }
7844 static const int mv_all[3][4][4] =
7846 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
7847 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
7848 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
7850 const int (&mv)[4][4] = mv_all[rt];
7852 #ifndef MFEM_USE_MEMALLOC
7854 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
7856 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
7858 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
7860 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
7862 for (
int k = 0; k < 4; k++)
7864 new_elements[j+4+k] =
7865 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
7866 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
7870 new_elements[j+0] = tet =
TetMemory.Alloc();
7871 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
7873 new_elements[j+1] = tet =
TetMemory.Alloc();
7874 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
7876 new_elements[j+2] = tet =
TetMemory.Alloc();
7877 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
7879 new_elements[j+3] = tet =
TetMemory.Alloc();
7880 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
7882 for (
int k = 0; k < 4; k++)
7884 new_elements[j+4+k] = tet =
TetMemory.Alloc();
7885 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
7886 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
7889 for (
int k = 0; k < 4; k++)
7894 for (
int k = 0; k < 4; k++)
7908 for (
int fi = 2; fi < 5; fi++)
7910 for (
int k = 0; k < 4; k++)
7917 for (
int ei = 0; ei < 9; ei++)
7919 for (
int k = 0; k < 2; k++)
7926 const int qf2 = f2qf[f[2]];
7927 const int qf3 = f2qf[f[3]];
7928 const int qf4 = f2qf[f[4]];
7931 new Wedge(v[0], oedge+e[0], oedge+e[2],
7932 oedge+e[6], oface+qf2, oface+qf4, attr);
7935 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
7936 oface+qf3, oface+qf4, oface+qf2, attr);
7939 new Wedge(oedge+e[0], v[1], oedge+e[1],
7940 oface+qf2, oedge+e[7], oface+qf3, attr);
7943 new Wedge(oedge+e[2], oedge+e[1], v[2],
7944 oface+qf4, oface+qf3, oedge+e[8], attr);
7947 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
7948 v[3], oedge+e[3], oedge+e[5], attr);
7951 new Wedge(oface+qf3, oface+qf4, oface+qf2,
7952 oedge+e[4], oedge+e[5], oedge+e[3], attr);
7955 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
7956 oedge+e[3], v[4], oedge+e[4], attr);
7959 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
7960 oedge+e[5], oedge+e[4], v[5], attr);
7967 const int he = hex_counter;
7972 if (f2qf.
Size() == 0)
7978 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[f[k]]; }
7984 for (
int fi = 0; fi < 6; fi++)
7986 for (
int k = 0; k < 4; k++)
7993 for (
int ei = 0; ei < 12; ei++)
7995 for (
int k = 0; k < 2; k++)
8003 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
8004 oedge+e[3], oedge+e[8], oface+qf[1],
8005 oelem+he, oface+qf[4], attr);
8008 oface+qf[0], oface+qf[1], oedge+e[9],
8009 oface+qf[2], oelem+he, attr);
8011 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
8012 oedge+e[2], oelem+he, oface+qf[2],
8013 oedge+e[10], oface+qf[3], attr);
8015 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
8016 v[3], oface+qf[4], oelem+he,
8017 oface+qf[3], oedge+e[11], attr);
8019 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
8020 oface+qf[4], v[4], oedge+e[4],
8021 oface+qf[5], oedge+e[7], attr);
8023 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
8024 oelem+he, oedge+e[4], v[5],
8025 oedge+e[5], oface+qf[5], attr);
8027 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
8028 oface+qf[3], oface+qf[5], oedge+e[5],
8029 v[6], oedge+e[6], attr);
8031 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
8032 oedge+e[11], oedge+e[7], oface+qf[5],
8033 oedge+e[6], v[7], attr);
8038 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
8050 const int attr =
boundary[i]->GetAttribute();
8051 int *v =
boundary[i]->GetVertices();
8058 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
8065 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
8067 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
8069 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
8071 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
8079 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
8081 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
8083 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
8085 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
8089 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
8095 static const double A = 0.0, B = 0.5, C = 1.0;
8096 static double tet_children[3*4*16] =
8098 A,A,A, B,A,A, A,B,A, A,A,B,
8099 B,A,A, C,A,A, B,B,A, B,A,B,
8100 A,B,A, B,B,A, A,C,A, A,B,B,
8101 A,A,B, B,A,B, A,B,B, A,A,C,
8106 B,A,A, A,B,B, A,B,A, A,A,B,
8107 B,A,A, A,B,B, A,A,B, B,A,B,
8108 B,A,A, A,B,B, B,A,B, B,B,A,
8109 B,A,A, A,B,B, B,B,A, A,B,A,
8111 A,B,A, B,A,A, B,A,B, A,A,B,
8112 A,B,A, A,A,B, B,A,B, A,B,B,
8113 A,B,A, A,B,B, B,A,B, B,B,A,
8114 A,B,A, B,B,A, B,A,B, B,A,A,
8116 A,A,B, B,A,A, A,B,A, B,B,A,
8117 A,A,B, A,B,A, A,B,B, B,B,A,
8118 A,A,B, A,B,B, B,A,B, B,B,A,
8119 A,A,B, B,A,B, B,A,A, B,B,A
8121 static double pri_children[3*6*8] =
8123 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
8124 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
8125 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
8126 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
8127 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
8128 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
8129 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
8130 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
8132 static double hex_children[3*8*8] =
8134 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
8135 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
8136 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
8137 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
8138 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
8139 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
8140 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
8141 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
8151 for (
int i = 0; i <
elements.Size(); i++)
8182 int i, j, ind, nedges;
8189 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
8203 for (j = 0; j < marked_el.
Size(); j++)
8208 int new_v = cnv + j, new_e = cne + j;
8217 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
8219 UseExternalData(seg_children, 1, 2, 3);
8232 int *edge1 =
new int[nedges];
8233 int *edge2 =
new int[nedges];
8234 int *middle =
new int[nedges];
8236 for (i = 0; i < nedges; i++)
8238 edge1[i] = edge2[i] = middle[i] = -1;
8244 for (j = 1; j < v.
Size(); j++)
8246 ind = v_to_v(v[j-1], v[j]);
8247 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8249 ind = v_to_v(v[0], v[v.
Size()-1]);
8250 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
8254 for (i = 0; i < marked_el.
Size(); i++)
8260 int need_refinement;
8263 need_refinement = 0;
8264 for (i = 0; i < nedges; i++)
8266 if (middle[i] != -1 && edge1[i] != -1)
8268 need_refinement = 1;
8273 while (need_refinement == 1);
8276 int v1[2], v2[2], bisect, temp;
8278 for (i = 0; i < temp; i++)
8281 bisect = v_to_v(v[0], v[1]);
8282 if (middle[bisect] != -1)
8286 v1[0] = v[0]; v1[1] = middle[bisect];
8287 v2[0] = middle[bisect]; v2[1] = v[1];
8293 mfem_error(
"Only bisection of segment is implemented"
8317 MFEM_VERIFY(
GetNE() == 0 ||
8319 "tetrahedral mesh is not marked for refinement:"
8320 " call Finalize(true)");
8327 for (i = 0; i < marked_el.
Size(); i++)
8333 for (i = 0; i < marked_el.
Size(); i++)
8342 for (i = 0; i < marked_el.
Size(); i++)
8359 int need_refinement;
8364 need_refinement = 0;
8372 if (
elements[i]->NeedRefinement(v_to_v))
8374 need_refinement = 1;
8379 while (need_refinement == 1);
8386 need_refinement = 0;
8388 if (
boundary[i]->NeedRefinement(v_to_v))
8390 need_refinement = 1;
8394 while (need_refinement == 1);
8425 MFEM_VERIFY(!
NURBSext,
"Nonconforming refinement of NURBS meshes is "
8426 "not supported. Project the NURBS to Nodes first.");
8436 if (!refinements.
Size())
8457 Swap(*mesh2,
false);
8473 const int *fine,
int nfine,
int op)
8476 for (
int i = 0; i < nfine; i++)
8478 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
8480 double err_fine = elem_error[fine[i]];
8483 case 0: error = std::min(error, err_fine);
break;
8484 case 1: error += err_fine;
break;
8485 case 2: error = std::max(error, err_fine);
break;
8492 double threshold,
int nc_limit,
int op)
8494 MFEM_VERIFY(
ncmesh,
"Only supported for non-conforming meshes.");
8495 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
8496 "Project the NURBS to Nodes first.");
8509 for (
int i = 0; i < dt.
Size(); i++)
8511 if (nc_limit > 0 && !level_ok[i]) {
continue; }
8516 if (error < threshold) { derefs.
Append(i); }
8519 if (!derefs.
Size()) {
return false; }
8526 Swap(*mesh2,
false);
8540 int nc_limit,
int op)
8550 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
8557 int nc_limit,
int op)
8560 for (
int i = 0; i < tmp.Size(); i++)
8562 tmp[i] = elem_error(i);
8647 #ifdef MFEM_USE_MEMALLOC
8674 for (
int i = 0; i < elem_array.
Size(); i++)
8676 if (elem_array[i]->GetGeometryType() ==
geom)
8681 elem_vtx.
SetSize(nv*num_elems);
8685 for (
int i = 0; i < elem_array.
Size(); i++)
8691 elem_vtx.
Append(loc_vtx);
8699 for (
int i = 0; i < nelem; i++) { list[i] = i; }
8715 else if (ref_algo == 1 &&
meshgen == 1 &&
Dim == 3)
8727 default: MFEM_ABORT(
"internal error");
8733 int nonconforming,
int nc_limit)
8743 else if (nonconforming < 0)
8764 for (
int i = 0; i < refinements.
Size(); i++)
8766 el_to_refine[i] = refinements[i].index;
8770 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
8771 if (rt == 1 || rt == 2 || rt == 4)
8775 else if (rt == 3 || rt == 5 || rt == 6)
8793 for (
int i = 0; i < el_to_refine.
Size(); i++)
8795 refinements[i] =
Refinement(el_to_refine[i]);
8802 MFEM_VERIFY(!
NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
8803 "Please project the NURBS to Nodes first, with SetCurvature().");
8806 MFEM_VERIFY(
ncmesh != NULL || dynamic_cast<const ParMesh*>(
this) == NULL,
8807 "Sorry, converting a conforming ParMesh to an NC mesh is "
8815 (simplices_nonconforming && (
meshgen & 0x1)) )
8828 for (
int i = 0; i <
GetNE(); i++)
8830 if ((
double) rand() / RAND_MAX <
prob)
8835 type = (
Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
8847 for (
int i = 0; i <
GetNE(); i++)
8850 bool refine =
false;
8851 for (
int j = 0; j < v.
Size(); j++)
8856 double d = vert(l) -
vertices[v[j]](l);
8859 if (dist <= eps*eps) { refine =
true;
break; }
8870 int nonconforming,
int nc_limit)
8872 MFEM_VERIFY(elem_error.
Size() ==
GetNE(),
"");
8874 for (
int i = 0; i <
GetNE(); i++)
8876 if (elem_error[i] > threshold)
8890 int nonconforming,
int nc_limit)
8894 return RefineByError(tmp, threshold, nonconforming, nc_limit);
8899 int *edge1,
int *edge2,
int *middle)
8902 int v[2][4], v_new, bisect,
t;
8914 bisect = v_to_v(vert[0], vert[1]);
8915 MFEM_ASSERT(bisect >= 0,
"");
8917 if (middle[bisect] == -1)
8928 if (edge1[bisect] == i)
8930 edge1[bisect] = edge2[bisect];
8933 middle[bisect] = v_new;
8937 v_new = middle[bisect];
8945 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
8946 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
8967 bisect = v_to_v(v[1][0], v[1][1]);
8968 MFEM_ASSERT(bisect >= 0,
"");
8970 if (edge1[bisect] == i)
8974 else if (edge2[bisect] == i)
8983 MFEM_ABORT(
"Bisection for now works only for triangles.");
8990 int v[2][4], v_new, bisect,
t;
8997 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
9001 "TETRAHEDRON element is not marked for refinement.");
9006 bisect = v_to_v.
FindId(vert[0], vert[1]);
9010 for (j = 0; j < 3; j++)
9027 new_redges[0][0] = 2;
9028 new_redges[0][1] = 1;
9029 new_redges[1][0] = 2;
9030 new_redges[1][1] = 1;
9031 int tr1 = -1, tr2 = -1;
9032 switch (old_redges[0])
9035 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
9040 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
9044 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
9047 switch (old_redges[1])
9050 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
9055 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
9059 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
9066 #ifdef MFEM_USE_MEMALLOC
9102 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
9109 int v[2][3], v_new, bisect,
t;
9120 bisect = v_to_v.
FindId(vert[0], vert[1]);
9121 MFEM_ASSERT(bisect >= 0,
"");
9123 MFEM_ASSERT(v_new != -1,
"");
9127 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
9128 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
9138 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
9144 int *edge1,
int *edge2,
int *middle)
9147 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
9156 bisect[0] = v_to_v(v[0],v[1]);
9157 bisect[1] = v_to_v(v[1],v[2]);
9158 bisect[2] = v_to_v(v[0],v[2]);
9159 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
9161 for (j = 0; j < 3; j++)
9163 if (middle[bisect[j]] == -1)
9174 if (edge1[bisect[j]] == i)
9176 edge1[bisect[j]] = edge2[bisect[j]];
9179 middle[bisect[j]] = v_new[j];
9183 v_new[j] = middle[bisect[j]];
9186 edge1[bisect[j]] = -1;
9192 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
9193 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
9194 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
9195 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
9210 tri2->ResetTransform(code);
9211 tri3->ResetTransform(code);
9215 tri2->PushTransform(1);
9216 tri3->PushTransform(2);
9229 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
9265 for (
int i = 0; i < elem_geoms.
Size(); i++)
9273 std::map<unsigned, int> mat_no;
9277 for (
int i = 0; i <
elements.Size(); i++)
9280 unsigned code =
elements[i]->GetTransform();
9283 int &matrix = mat_no[code];
9284 if (!matrix) { matrix = mat_no.size(); }
9294 std::map<unsigned, int>::iterator it;
9295 for (it = mat_no.begin(); it != mat_no.end(); ++it)
9309 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
9310 " geom = " << geom);
9320 MFEM_ASSERT(
Dim==
spaceDim,
"2D Manifold meshes not supported");
9329 out <<
"areamesh2\n\n";
9333 out <<
"curved_areamesh2\n\n";
9342 out <<
boundary[i]->GetAttribute();
9343 for (j = 0; j < v.
Size(); j++)
9345 out <<
' ' << v[j] + 1;
9357 for (j = 0; j < v.
Size(); j++)
9359 out <<
' ' << v[j] + 1;
9371 for (j = 1; j <
Dim; j++)
9388 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
9396 out <<
"NETGEN_Neutral_Format\n";
9401 for (j = 0; j <
Dim; j++)
9414 out <<
elements[i]->GetAttribute();
9415 for (j = 0; j < nv; j++)
9417 out <<
' ' << ind[j]+1;
9428 out <<
boundary[i]->GetAttribute();
9429 for (j = 0; j < nv; j++)
9431 out <<
' ' << ind[j]+1;
9443 <<
" 0 0 0 0 0 0 0\n"
9444 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
9446 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
9447 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
9451 <<
' ' <<
vertices[i](2) <<
" 0.0\n";
9457 out << i+1 <<
' ' <<
elements[i]->GetAttribute();
9458 for (j = 0; j < nv; j++)
9460 out <<
' ' << ind[j]+1;
9469 out <<
boundary[i]->GetAttribute();
9470 for (j = 0; j < nv; j++)
9472 out <<
' ' << ind[j]+1;
9474 out <<
" 1.0 1.0 1.0 1.0\n";
9507 out <<
"\n# mesh curvature GridFunction";
9512 out <<
"\nmfem_mesh_end" << endl;
9517 out << (section_delimiter.empty()
9518 ?
"MFEM mesh v1.0\n" :
"MFEM mesh v1.2\n");
9522 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
9527 "# TETRAHEDRON = 4\n"
9532 out <<
"\ndimension\n" <<
Dim;
9567 if (!section_delimiter.empty())
9569 out << section_delimiter << endl;
9578 out <<
"MFEM NURBS mesh v1.0\n";
9582 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
9588 out <<
"\ndimension\n" <<
Dim
9610 out << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
9617 ofstream ofs(fname);
9618 ofs.precision(precision);
9622 #ifdef MFEM_USE_ADIOS2
9632 "# vtk DataFile Version 3.0\n"
9633 "Generated by MFEM\n"
9635 "DATASET UNSTRUCTURED_GRID\n";
9664 out << (*Nodes)(vdofs[0]);
9668 out <<
' ' << (*Nodes)(vdofs[j]);
9684 size +=
elements[i]->GetNVertices() + 1;
9686 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
9689 const int *v =
elements[i]->GetVertices();
9690 const int nv =
elements[i]->GetNVertices();
9694 for (
int j = 0; j < nv; j++)
9696 out <<
' ' << v[perm ? perm[j] : j];
9709 MFEM_ASSERT(
Dim != 0 || dofs.
Size() == 1,
9710 "Point meshes should have a single dof per element");
9711 size += dofs.
Size() + 1;
9713 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
9716 if (!strcmp(fec_name,
"Linear") ||
9717 !strcmp(fec_name,
"H1_0D_P1") ||
9718 !strcmp(fec_name,
"H1_1D_P1") ||
9719 !strcmp(fec_name,
"H1_2D_P1") ||
9720 !strcmp(fec_name,
"H1_3D_P1"))
9724 else if (!strcmp(fec_name,
"Quadratic") ||
9725 !strcmp(fec_name,
"H1_1D_P2") ||
9726 !strcmp(fec_name,
"H1_2D_P2") ||
9727 !strcmp(fec_name,
"H1_3D_P2"))
9733 mfem::err <<
"Mesh::PrintVTK : can not save '"
9734 << fec_name <<
"' elements!" << endl;
9743 for (
int j = 0; j < dofs.
Size(); j++)
9745 out <<
' ' << dofs[j];
9748 else if (order == 2)
9750 const int *vtk_mfem;
9751 switch (
elements[i]->GetGeometryType())
9765 for (
int j = 0; j < dofs.
Size(); j++)
9767 out <<
' ' << dofs[vtk_mfem[j]];
9777 int vtk_cell_type = 5;
9781 out << vtk_cell_type <<
'\n';
9785 out <<
"CELL_DATA " << NumOfElements <<
'\n'
9786 <<
"SCALARS material int\n"
9787 <<
"LOOKUP_TABLE default\n";
9790 out <<
elements[i]->GetAttribute() <<
'\n';
9797 bool high_order_output,
9798 int compression_level,
9801 int ref = (high_order_output &&
Nodes)
9804 fname = fname +
".vtu";
9806 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
9807 if (compression_level != 0)
9809 out <<
" compressor=\"vtkZLibDataCompressor\"";
9812 out <<
"<UnstructuredGrid>\n";
9813 PrintVTU(out, ref, format, high_order_output, compression_level, bdr);
9814 out <<
"</Piece>\n";
9815 out <<
"</UnstructuredGrid>\n";
9816 out <<
"</VTKFile>" << std::endl;
9823 bool high_order_output,
9824 int compression_level)
9826 PrintVTU(fname, format, high_order_output, compression_level,
true);
9829 template <
typename T>
9840 const uint8_t &val,
const char *suffix,
9843 if (format ==
VTKFormat::ASCII) { out << static_cast<int>(val) << suffix; }
9849 const double &val,
const char *suffix,
9854 bin_io::AppendBytes<float>(buf, float(val));
9862 out << val << suffix;
9868 const float &val,
const char *suffix,
9873 else {
out << val << suffix; }
9877 int compression_level)
9885 bool high_order_output,
int compression_level,
9893 std::vector<char> buf;
9895 auto get_geom = [&](
int i)
9903 int np = 0, nc_ref = 0, size = 0;
9904 for (
int i = 0; i < ne; i++)
9914 out <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
9915 << (high_order_output ? ne : nc_ref) <<
"\">\n";
9918 out <<
"<Points>\n";
9919 out <<
"<DataArray type=\"" << type_str
9920 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
9921 for (
int i = 0; i < ne; i++)
9934 for (
int j = 0; j < pmat.
Width(); j++)
9960 out <<
"</DataArray>" << std::endl;
9961 out <<
"</Points>" << std::endl;
9963 out <<
"<Cells>" << std::endl;
9964 out <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
9965 << fmt_str <<
"\">" << std::endl;
9967 std::vector<int> offset;
9970 if (high_order_output)
9973 for (
int iel = 0; iel < ne; iel++)
9977 int nnodes = local_connectivity.
Size();
9978 for (
int i=0; i<nnodes; ++i)
9984 offset.push_back(np);
9990 for (
int i = 0; i < ne; i++)
9996 for (
int j = 0; j < RG.
Size(); )
9999 offset.push_back(coff);
10001 for (
int k = 0; k < nv; k++, j++)
10014 out <<
"</DataArray>" << std::endl;
10016 out <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
10017 << fmt_str <<
"\">" << std::endl;
10019 for (
size_t ii=0; ii<offset.size(); ii++)
10027 out <<
"</DataArray>" << std::endl;
10028 out <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
10029 << fmt_str <<
"\">" << std::endl;
10031 const int *vtk_geom_map =
10033 for (
int i = 0; i < ne; i++)
10036 uint8_t vtk_cell_type = 5;
10038 vtk_cell_type = vtk_geom_map[
geom];
10040 if (high_order_output)
10049 for (
int j = 0; j < RG.
Size(); j += nv)
10059 out <<
"</DataArray>" << std::endl;
10060 out <<
"</Cells>" << std::endl;
10062 out <<
"<CellData Scalars=\"attribute\">" << std::endl;
10063 out <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
10064 << fmt_str <<
"\">" << std::endl;
10065 for (
int i = 0; i < ne; i++)
10068 if (high_order_output)
10087 out <<
"</DataArray>" << std::endl;
10088 out <<
"</CellData>" << std::endl;
10099 "# vtk DataFile Version 3.0\n"
10100 "Generated by MFEM\n"
10102 "DATASET UNSTRUCTURED_GRID\n";
10107 out <<
"FIELD FieldData 1\n"
10117 np = nc = size = 0;
10118 for (
int i = 0; i <
GetNE(); i++)
10127 out <<
"POINTS " << np <<
" double\n";
10129 for (
int i = 0; i <
GetNE(); i++)
10136 for (
int j = 0; j < pmat.
Width(); j++)
10138 out << pmat(0, j) <<
' ';
10141 out << pmat(1, j) <<
' ';
10153 out << 0.0 <<
' ' << 0.0;
10160 out <<
"CELLS " << nc <<
' ' << size <<
'\n';
10162 for (
int i = 0; i <
GetNE(); i++)
10169 for (
int j = 0; j < RG.
Size(); )
10172 for (
int k = 0; k < nv; k++, j++)
10174 out <<
' ' << np + RG[j];
10180 out <<
"CELL_TYPES " << nc <<
'\n';
10181 for (
int i = 0; i <
GetNE(); i++)
10189 for (
int j = 0; j < RG.
Size(); j += nv)
10191 out << vtk_cell_type <<
'\n';
10195 out <<
"CELL_DATA " << nc <<
'\n'
10196 <<
"SCALARS material int\n"
10197 <<
"LOOKUP_TABLE default\n";
10198 for (
int i = 0; i <
GetNE(); i++)
10206 out << attr <<
'\n';
10213 srand((
unsigned)time(0));
10214 double a = double(rand()) / (double(RAND_MAX) + 1.);
10215 int el0 = (int)floor(a *
GetNE());
10217 out <<
"SCALARS element_coloring int\n"
10218 <<
"LOOKUP_TABLE default\n";
10219 for (
int i = 0; i <
GetNE(); i++)
10226 out << coloring[i] + 1 <<
'\n';
10232 out <<
"POINT_DATA " << np <<
'\n' << flush;
10237 int delete_el_to_el = (
el_to_el) ? (0) : (1);
10239 int num_el =
GetNE(), stack_p, stack_top_p, max_num_col;
10242 const int *i_el_el = el_el.
GetI();
10243 const int *j_el_el = el_el.
GetJ();
10248 stack_p = stack_top_p = 0;
10249 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
10251 if (colors[el] != -2)
10257 el_stack[stack_top_p++] = el;
10259 for ( ; stack_p < stack_top_p; stack_p++)
10261 int i = el_stack[stack_p];
10262 int num_nb = i_el_el[i+1] - i_el_el[i];
10263 if (max_num_col < num_nb + 1)
10265 max_num_col = num_nb + 1;
10267 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10269 int k = j_el_el[j];
10270 if (colors[k] == -2)
10273 el_stack[stack_top_p++] = k;
10281 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
10283 int i = el_stack[stack_p], col;
10285 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
10287 col = colors[j_el_el[j]];
10290 col_marker[col] = 1;
10294 for (col = 0; col < max_num_col; col++)
10295 if (col_marker[col] == 0)
10303 if (delete_el_to_el)
10311 int elem_attr)
const
10313 if (
Dim != 3 &&
Dim != 2) {
return; }
10315 int i, j, k, l, nv, nbe, *v;
10317 out <<
"MFEM mesh v1.0\n";
10321 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10326 "# TETRAHEDRON = 4\n"
10331 out <<
"\ndimension\n" <<
Dim
10336 <<
' ' <<
elements[i]->GetGeometryType();
10339 for (j = 0; j < nv; j++)
10341 out <<
' ' << v[j];
10351 l = partitioning[l];
10366 out <<
"\nboundary\n" << nbe <<
'\n';
10372 l = partitioning[l];
10375 nv =
faces[i]->GetNVertices();
10376 v =
faces[i]->GetVertices();
10377 out << k+1 <<
' ' <<
faces[i]->GetGeometryType();
10378 for (j = 0; j < nv; j++)
10380 out <<
' ' << v[j];
10385 out << l+1 <<
' ' <<
faces[i]->GetGeometryType();
10386 for (j = nv-1; j >= 0; j--)
10388 out <<
' ' << v[j];
10397 nv =
faces[i]->GetNVertices();
10398 v =
faces[i]->GetVertices();
10399 out << k+1 <<
' ' <<
faces[i]->GetGeometryType();
10400 for (j = 0; j < nv; j++)
10402 out <<
' ' << v[j];
10424 out <<
"\nnodes\n";
10431 int interior_faces)
10433 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported\n");
10434 if (
Dim != 3 &&
Dim != 2) {
return; }
10450 for (j = 0; j < nv; j++)
10456 int *voff =
new int[NumOfVertices+1];
10460 voff[i] = vcount[i-1] + voff[i-1];
10466 vown[i] =
new int[vcount[i]];
10483 for (j = 0; j < nv; j++)
10486 vown[ind[j]][vcount[ind[j]]] = i;
10492 vcount[i] = voff[i+1] - voff[i];
10496 for (i = 0; i < edge_el.
Size(); i++)
10498 const int *el = edge_el.
GetRow(i);
10501 k = partitioning[el[0]];
10502 l = partitioning[el[1]];
10503 if (interior_faces || k != l)
10515 out <<
"areamesh2\n\n" << nbe <<
'\n';
10517 for (i = 0; i < edge_el.
Size(); i++)
10519 const int *el = edge_el.
GetRow(i);
10522 k = partitioning[el[0]];
10523 l = partitioning[el[1]];
10524 if (interior_faces || k != l)
10529 for (j = 0; j < 2; j++)
10530 for (s = 0; s < vcount[ev[j]]; s++)
10531 if (vown[ev[j]][s] == el[0])
10533 out <<
' ' << voff[ev[j]]+s+1;
10537 for (j = 1; j >= 0; j--)
10538 for (s = 0; s < vcount[ev[j]]; s++)
10539 if (vown[ev[j]][s] == el[1])
10541 out <<
' ' << voff[ev[j]]+s+1;
10548 k = partitioning[el[0]];
10552 for (j = 0; j < 2; j++)
10553 for (s = 0; s < vcount[ev[j]]; s++)
10554 if (vown[ev[j]][s] == el[0])
10556 out <<
' ' << voff[ev[j]]+s+1;
10563 out << NumOfElements <<
'\n';
10568 out << partitioning[i]+1 <<
' ';
10570 for (j = 0; j < nv; j++)
10572 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
10573 vown[ind[j]][vcount[ind[j]]] = i;
10580 vcount[i] = voff[i+1] - voff[i];
10586 for (k = 0; k < vcount[i]; k++)
10588 for (j = 0; j <
Dim; j++)
10598 out <<
"NETGEN_Neutral_Format\n";
10602 for (k = 0; k < vcount[i]; k++)
10604 for (j = 0; j <
Dim; j++)
10612 out << NumOfElements <<
'\n';
10617 out << partitioning[i]+1;
10618 for (j = 0; j < nv; j++)
10620 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
10621 vown[ind[j]][vcount[ind[j]]] = i;
10628 vcount[i] = voff[i+1] - voff[i];
10638 l = partitioning[l];
10639 if (interior_faces || k != l)
10649 out << nbe <<
'\n';
10654 l = partitioning[l];
10655 if (interior_faces || k != l)
10657 nv =
faces[i]->GetNVertices();
10658 ind =
faces[i]->GetVertices();
10660 for (j = 0; j < nv; j++)
10661 for (s = 0; s < vcount[ind[j]]; s++)
10662 if (vown[ind[j]][s] == faces_info[i].Elem1No)
10664 out <<
' ' << voff[ind[j]]+s+1;
10668 for (j = nv-1; j >= 0; j--)
10669 for (s = 0; s < vcount[ind[j]]; s++)
10670 if (vown[ind[j]][s] == faces_info[i].Elem2No)
10672 out <<
' ' << voff[ind[j]]+s+1;
10680 nv =
faces[i]->GetNVertices();
10681 ind =
faces[i]->GetVertices();
10683 for (j = 0; j < nv; j++)
10684 for (s = 0; s < vcount[ind[j]]; s++)
10685 if (vown[ind[j]][s] == faces_info[i].Elem1No)
10687 out <<
' ' << voff[ind[j]]+s+1;
10702 l = partitioning[l];
10703 if (interior_faces || k != l)
10714 out <<
"TrueGrid\n"
10716 <<
" 0 0 0 0 0 0 0\n"
10717 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
10718 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
10719 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
10720 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
10723 for (k = 0; k < vcount[i]; k++)
10724 out << voff[i]+k <<
" 0.0 " <<
vertices[i](0) <<
' '
10731 out << i+1 <<
' ' << partitioning[i]+1;
10732 for (j = 0; j < nv; j++)
10734 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
10735 vown[ind[j]][vcount[ind[j]]] = i;
10742 vcount[i] = voff[i+1] - voff[i];
10750 l = partitioning[l];
10751 if (interior_faces || k != l)
10753 nv =
faces[i]->GetNVertices();
10754 ind =
faces[i]->GetVertices();
10756 for (j = 0; j < nv; j++)
10757 for (s = 0; s < vcount[ind[j]]; s++)
10758 if (vown[ind[j]][s] == faces_info[i].Elem1No)
10760 out <<
' ' << voff[ind[j]]+s+1;
10762 out <<
" 1.0 1.0 1.0 1.0\n";
10764 for (j = nv-1; j >= 0; j--)
10765 for (s = 0; s < vcount[ind[j]]; s++)
10766 if (vown[ind[j]][s] == faces_info[i].Elem2No)
10768 out <<
' ' << voff[ind[j]]+s+1;
10770 out <<
" 1.0 1.0 1.0 1.0\n";
10776 nv =
faces[i]->GetNVertices();
10777 ind =
faces[i]->GetVertices();
10779 for (j = 0; j < nv; j++)
10780 for (s = 0; s < vcount[ind[j]]; s++)
10781 if (vown[ind[j]][s] == faces_info[i].Elem1No)
10783 out <<
' ' << voff[ind[j]]+s+1;
10785 out <<
" 1.0 1.0 1.0 1.0\n";
10808 " NURBS mesh is not supported!");
10812 out <<
"MFEM mesh v1.0\n";
10816 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10821 "# TETRAHEDRON = 4\n"
10826 out <<
"\ndimension\n" <<
Dim
10834 const int *
const i_AF_f = Aface_face.
GetI();
10835 const int *
const j_AF_f = Aface_face.
GetJ();
10837 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
10838 for (
const int * iface = j_AF_f + i_AF_f[iAF];
10839 iface < j_AF_f + i_AF_f[iAF+1];
10842 out << iAF+1 <<
' ';
10863 out <<
"\nnodes\n";
10874 double *cg =
new double[na*
spaceDim];
10875 int *nbea =
new int[na];
10882 for (i = 0; i < na; i++)
10886 cg[i*spaceDim+j] = 0.0;
10894 for (k = 0; k < vert.
Size(); k++)
10906 for (k = 0; k < vert.
Size(); k++)
10907 if (vn[vert[k]] == 1)
10912 cg[bea*spaceDim+j] += pointmat(j,k);
10923 for (k = 0; k < vert.
Size(); k++)
10928 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
10944 double *cg =
new double[na*
spaceDim];
10945 int *nbea =
new int[na];
10952 for (i = 0; i < na; i++)
10956 cg[i*spaceDim+j] = 0.0;
10964 for (k = 0; k < vert.
Size(); k++)
10976 for (k = 0; k < vert.
Size(); k++)
10977 if (vn[vert[k]] == 1)
10982 cg[bea*spaceDim+j] += pointmat(j,k);
10993 for (k = 0; k < vert.
Size(); k++)
10998 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
11014 for (
int i = 0; i <
vertices.Size(); i++)
11016 for (
int j = 0; j <
spaceDim; j++)
11028 xnew.ProjectCoefficient(f_pert);
11036 "incompatible vector dimensions");
11044 for (
int d = 0; d <
spaceDim; d++)
11046 vertices[i](d) = xnew(d + spaceDim*i);
11063 for (
int i = 0; i <
GetNE(); i++)
11068 for (
int j = 0; j < nv; j++)
11073 for (
int i = 0; i <
GetNBE(); i++)
11078 for (
int j = 0; j < nv; j++)
11084 for (
int i = 0; i < v2v.
Size(); i++)
11089 v2v[i] = num_vert++;
11093 if (num_vert == v2v.
Size()) {
return; }
11095 Vector nodes_by_element;
11100 for (
int i = 0; i <
GetNE(); i++)
11107 for (
int i = 0; i <
GetNE(); i++)
11116 for (
int i = 0; i <
GetNE(); i++)
11121 for (
int j = 0; j < nv; j++)
11126 for (
int i = 0; i <
GetNBE(); i++)
11131 for (
int j = 0; j < nv; j++)
11155 for (
int i = 0; i <
GetNE(); i++)
11168 int num_bdr_elem = 0;
11169 int new_bel_to_edge_nnz = 0;
11170 for (
int i = 0; i <
GetNBE(); i++)
11186 if (num_bdr_elem ==
GetNBE()) {
return; }
11190 Table *new_bel_to_edge = NULL;
11194 new_be_to_edge.
Reserve(num_bdr_elem);
11198 new_be_to_face.
Reserve(num_bdr_elem);
11199 new_bel_to_edge =
new Table;
11200 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
11202 for (
int i = 0; i <
GetNBE(); i++)
11213 int row = new_be_to_face.
Size();
11217 int *new_e = new_bel_to_edge->
GetRow(row);
11218 for (
int j = 0; j < ne; j++)
11222 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
11242 for (
int i = 0; i < attribs.
Size(); i++)
11254 #ifdef MFEM_USE_MEMALLOC
11281 const int npts = point_mat.
Width();
11282 if (!npts) {
return 0; }
11283 MFEM_VERIFY(point_mat.
Height() ==
spaceDim,
"Invalid points matrix");
11287 if (!
GetNE()) {
return 0; }
11289 double *data = point_mat.
GetData();
11296 min_dist = std::numeric_limits<double>::max();
11300 for (
int i = 0; i <
GetNE(); i++)
11304 for (
int k = 0; k < npts; k++)
11307 if (dist < min_dist(k))
11309 min_dist(k) = dist;
11318 for (
int k = 0; k < npts; k++)
11322 int res = inv_tr->
Transform(pt, ips[k]);
11325 elem_ids[k] = e_idx[k];
11329 if (pts_found != npts)
11333 for (
int k = 0; k < npts; k++)
11335 if (elem_ids[k] != -1) {
continue; }
11339 for (
int v = 0; v < vertices.
Size(); v++)
11341 int vv = vertices[v];
11343 const int* els = vtoel->
GetRow(vv);
11344 for (
int e = 0; e < ne; e++)
11346 if (els[e] == e_idx[k]) {
continue; }
11348 int res = inv_tr->
Transform(pt, ips[k]);
11351 elem_ids[k] = els[e];
11363 for (
int e = 0; e < neigh.
Size(); e++)
11369 int res = inv_tr->
Transform(pt, ips[k]);
11382 if (inv_trans == NULL) {
delete inv_tr; }
11384 if (warn && pts_found != npts)
11386 MFEM_WARNING((npts-pts_found) <<
" points were not found");
11400 "mixed meshes are not supported!");
11401 MFEM_ASSERT(mesh->
GetNodes(),
"meshes without nodes are not supported!");
11414 Compute(nodes, d_mt);
11417 void GeometricFactors::Compute(
const GridFunction &nodes,
11424 const int vdim = fespace->
GetVDim();
11425 const int NE = fespace->
GetNE();
11426 const int ND = fe->
GetDof();
11429 unsigned eval_flags = 0;
11439 J.
SetSize(dim*vdim*NQ*NE, my_d_mt);
11454 qi->DisableTensorProducts(!use_tensor_products);
11462 Vector Enodes(vdim*ND*NE, my_d_mt);
11463 elem_restr->Mult(nodes, Enodes);
11464 qi->Mult(Enodes, eval_flags,
X,
J,
detJ);
11468 qi->Mult(nodes, eval_flags,
X,
J,
detJ);
11483 const int vdim = fespace->
GetVDim();
11492 face_restr->
Mult(*nodes, Fnodes);
11494 unsigned eval_flags = 0;
11535 V(1) = s * ((ip.
y + layer) / n);
11540 V(2) = s * ((ip.
z + layer) / n);
11549 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
11553 int nvy = (closed) ? (ny) : (ny + 1);
11554 int nvt = mesh->
GetNV() * nvy;
11563 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
11568 for (
int i = 0; i < mesh->
GetNV(); i++)
11571 for (
int j = 0; j < nvy; j++)
11573 vc[1] = sy * (double(j) / ny);
11579 for (
int i = 0; i < mesh->
GetNE(); i++)
11584 for (
int j = 0; j < ny; j++)
11587 qv[0] = vert[0] * nvy + j;
11588 qv[1] = vert[1] * nvy + j;
11589 qv[2] = vert[1] * nvy + (j + 1) % nvy;
11590 qv[3] = vert[0] * nvy + (j + 1) % nvy;
11596 for (
int i = 0; i < mesh->
GetNBE(); i++)
11601 for (
int j = 0; j < ny; j++)
11604 sv[0] = vert[0] * nvy + j;
11605 sv[1] = vert[0] * nvy + (j + 1) % nvy;
11609 Swap<int>(sv[0], sv[1]);
11621 for (
int i = 0; i < mesh->
GetNE(); i++)
11627 sv[0] = vert[0] * nvy;
11628 sv[1] = vert[1] * nvy;
11632 sv[0] = vert[1] * nvy + ny;
11633 sv[1] = vert[0] * nvy + ny;
11649 string cname = name;
11650 if (cname ==
"Linear")
11654 else if (cname ==
"Quadratic")
11658 else if (cname ==
"Cubic")
11662 else if (!strncmp(name,
"H1_", 3))
11666 else if (!strncmp(name,
"L2_T", 4))
11670 else if (!strncmp(name,
"L2_", 3))
11677 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
11689 for (
int i = 0; i < mesh->
GetNE(); i++)
11692 for (
int j = ny-1; j >= 0; j--)
11709 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
11714 int nvt = mesh->
GetNV() * nvz;
11719 bool wdgMesh =
false;
11720 bool hexMesh =
false;
11724 for (
int i = 0; i < mesh->
GetNV(); i++)
11728 for (
int j = 0; j < nvz; j++)
11730 vc[2] = sz * (double(j) / nz);
11736 for (
int i = 0; i < mesh->
GetNE(); i++)
11746 for (
int j = 0; j < nz; j++)
11749 pv[0] = vert[0] * nvz + j;
11750 pv[1] = vert[1] * nvz + j;
11751 pv[2] = vert[2] * nvz + j;
11752 pv[3] = vert[0] * nvz + (j + 1) % nvz;
11753 pv[4] = vert[1] * nvz + (j + 1) % nvz;
11754 pv[5] = vert[2] * nvz + (j + 1) % nvz;
11761 for (
int j = 0; j < nz; j++)
11764 hv[0] = vert[0] * nvz + j;
11765 hv[1] = vert[1] * nvz + j;
11766 hv[2] = vert[2] * nvz + j;
11767 hv[3] = vert[3] * nvz + j;
11768 hv[4] = vert[0] * nvz + (j + 1) % nvz;
11769 hv[5] = vert[1] * nvz + (j + 1) % nvz;
11770 hv[6] = vert[2] * nvz + (j + 1) % nvz;
11771 hv[7] = vert[3] * nvz + (j + 1) % nvz;
11773 mesh3d->
AddHex(hv, attr);
11777 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
11778 << geom <<
"\'" << endl;
11784 for (
int i = 0; i < mesh->
GetNBE(); i++)
11789 for (
int j = 0; j < nz; j++)
11792 qv[0] = vert[0] * nvz + j;
11793 qv[1] = vert[1] * nvz + j;
11794 qv[2] = vert[1] * nvz + (j + 1) % nvz;
11795 qv[3] = vert[0] * nvz + (j + 1) % nvz;
11804 for (
int i = 0; i < mesh->
GetNE(); i++)
11815 tv[0] = vert[0] * nvz;
11816 tv[1] = vert[2] * nvz;
11817 tv[2] = vert[1] * nvz;
11821 tv[0] = vert[0] * nvz + nz;
11822 tv[1] = vert[1] * nvz + nz;
11823 tv[2] = vert[2] * nvz + nz;
11831 qv[0] = vert[0] * nvz;
11832 qv[1] = vert[3] * nvz;
11833 qv[2] = vert[2] * nvz;
11834 qv[3] = vert[1] * nvz;
11838 qv[0] = vert[0] * nvz + nz;
11839 qv[1] = vert[1] * nvz + nz;
11840 qv[2] = vert[2] * nvz + nz;
11841 qv[3] = vert[3] * nvz + nz;
11847 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
11848 << geom <<
"\'" << endl;
11854 if ( hexMesh && wdgMesh )
11858 else if ( hexMesh )
11862 else if ( wdgMesh )
11875 string cname = name;
11876 if (cname ==
"Linear")
11880 else if (cname ==
"Quadratic")
11884 else if (cname ==
"Cubic")
11888 else if (!strncmp(name,
"H1_", 3))
11892 else if (!strncmp(name,
"L2_T", 4))
11896 else if (!strncmp(name,
"L2_", 3))
11903 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
11915 for (
int i = 0; i < mesh->
GetNE(); i++)
11918 for (
int j = nz-1; j >= 0; j--)
11939 out << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
11940 <<
" 0 0 " << i <<
" -1 0\n";
11947 double mid[3] = {0, 0, 0};
11948 for (
int j = 0; j < 2; j++)
11950 for (
int k = 0; k <
spaceDim; k++)
11955 out << NumOfVertices+i <<
" "
11956 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
11957 << 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]
virtual void Print(std::ostream &out=mfem::out) const
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)
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]
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.
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.
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void PrintWithPartitioning(int *partitioning, std::ostream &out, int elem_attr=0) const
Prints the mesh with boundary elements given by the boundary of the subdomains, so that the boundary ...
static FiniteElement * GetTransformationFEforElementType(Element::Type)
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.
void WriteBinaryOrASCII< float >(std::ostream &out, std::vector< char > &buf, const float &val, const char *suffix, VTKFormat format)
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)
void PrintVTU(std::ostream &out, int ref=1, VTKFormat format=VTKFormat::ASCII, bool high_order_output=false, int compression_level=0, bool bdr_elements=false)
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 GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
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)
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 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.
void WriteBase64WithSizeAndClear(std::ostream &out, std::vector< char > &buf, int compression_level)
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.
static const int Edges[NumEdges][2]
void MoveNodes(const Vector &displacements)
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()
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]
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)
void WriteBinaryOrASCII< uint8_t >(std::ostream &out, std::vector< char > &buf, const uint8_t &val, const char *suffix, VTKFormat format)
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)
NURBSExtension * StealNURBSext()
Array< DenseMatrix * > point_matrices[Geometry::NumGeom]
List of unique point matrices for each slave geometry.
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.
void WriteBinaryOrASCII(std::ostream &out, std::vector< char > &buf, const T &val, const char *suffix, VTKFormat format)
const Element * GetElement(int i) const
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
void GetElementTopo(Array< Element * > &elements) const
signed char local
local number within 'element'
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
int GetVDim()
Returns dimension of the vector.
virtual void GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'.
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 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]
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.
class H1_WedgeElement WedgeFE
void FindPartitioningComponents(Table &elem_elem, const Array< int > &partitioning, Array< int > &component, Array< int > &num_comp)
Data type tetrahedron element.
void GetElementInteriorDofs(int i, Array< int > &dofs) const
List of mesh geometries stored as Array<Geometry::Type>.
A general vector function coefficient.
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 GetFaceElements(int Face, int *Elem1, int *Elem2) const
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.
void PrintVTK(std::ostream &out)
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
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
int SpaceDimension() const
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.
virtual void PrintXG(std::ostream &out=mfem::out) const
Print the mesh to the given stream using Netgen/Truegrid format.
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
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)
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]
void WriteBinaryOrASCII< double >(std::ostream &out, std::vector< char > &buf, const double &val, const char *suffix, VTKFormat format)
const FiniteElementSpace * GetNodalFESpace() const
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Array< int > leaf_elements
finest elements, in Mesh ordering (+ ghosts)
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
static const int Orient[NumOrient][NumVert]
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
void GetFaceEdges(int i, Array< int > &edges, Array< int > &o) const
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
const char * VTKByteOrder()
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]
void WriteVTKEncodedCompressed(std::ostream &out, const void *bytes, uint32_t nbytes, int compression_level)
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
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.
void AppendBytes(std::vector< char > &vec, const T &val)
Append the binary representation of val to the byte buffer vec.
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
NCMesh * ncmesh
Optional non-conforming 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.
Evaluate the derivatives at quadrature points.
void ReadTrueGridMesh(std::istream &input)
const FiniteElementCollection * FEColl() const
void GetNodes(Vector &node_coord) const
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
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)
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
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]
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)
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)
int matrix
Index into the DenseTensor corresponding to the parent Geometry::Type stored in CoarseFineTransformat...
Arbitrary order "L2-conforming" discontinuous finite elements.
Evaluate the values at quadrature points.
void Printer(std::ostream &out=mfem::out, std::string section_delimiter="") const
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...