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"
35 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
40 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
45 int*,
int*,
int*,
int*,
int*,
idxtype*);
47 int*,
int*,
int*,
int*,
int*,
idxtype*);
49 int*,
int*,
int*,
int*,
int*,
idxtype*);
66 void Mesh::GetElementCenter(
int i,
Vector ¢er)
69 int geom = GetElementBaseGeometry(i);
74 double Mesh::GetElementSize(
int i,
int type)
77 GetElementJacobian(i, J);
80 return pow(fabs(J.
Det()), 1./Dim);
92 double Mesh::GetElementSize(
int i,
const Vector &dir)
96 GetElementJacobian(i, J);
98 return sqrt((d_hat * d_hat) / (dir * dir));
101 double Mesh::GetElementVolume(
int i)
123 for (
int d = 0; d < spaceDim; d++)
132 for (
int i = 0; i < NumOfVertices; i++)
134 coord = GetVertex(i);
135 for (
int d = 0; d < spaceDim; d++)
137 if (coord[d] < min(d)) { min(d) = coord[d]; }
138 if (coord[d] > max(d)) { max(d) = coord[d]; }
144 const bool use_boundary =
false;
145 int ne = use_boundary ? GetNBE() : GetNE();
153 for (
int i = 0; i < ne; i++)
157 GetBdrElementFace(i, &fn, &fo);
159 Tr = GetFaceElementTransformations(fn, 5);
166 T = GetElementTransformation(i);
170 for (
int j = 0; j < pointmat.
Width(); j++)
172 for (
int d = 0; d < pointmat.
Height(); d++)
174 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
175 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
182 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
183 double &kappa_min,
double &kappa_max,
191 sdim = SpaceDimension();
193 if (Vh) { Vh->
SetSize(NumOfElements); }
194 if (Vk) { Vk->
SetSize(NumOfElements); }
197 h_max = kappa_max = -h_min;
198 if (dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
200 for (i = 0; i < NumOfElements; i++)
202 GetElementJacobian(i, J);
203 h = pow(fabs(J.
Weight()), 1.0/
double(dim));
204 kappa = (dim == sdim) ?
206 if (Vh) { (*Vh)(i) = h; }
207 if (Vk) { (*Vk)(i) = kappa; }
209 if (h < h_min) { h_min = h; }
210 if (h > h_max) { h_max = h; }
211 if (kappa < kappa_min) { kappa_min =
kappa; }
212 if (kappa > kappa_max) { kappa_max =
kappa; }
217 void Mesh::PrintElementsByGeometry(
int dim,
221 for (
int g = Geometry::DimStart[dim], first = 1;
222 g < Geometry::DimStart[dim+1]; g++)
224 if (!num_elems_by_geom[g]) {
continue; }
225 if (!first) { out <<
" + "; }
227 out << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
233 double h_min, h_max, kappa_min, kappa_max;
235 out <<
"Mesh Characteristics:";
237 this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
239 Array<int> num_elems_by_geom(Geometry::NumGeom);
240 num_elems_by_geom = 0;
241 for (
int i = 0; i < GetNE(); i++)
243 num_elems_by_geom[GetElementBaseGeometry(i)]++;
247 <<
"Dimension : " << Dimension() <<
'\n'
248 <<
"Space dimension : " << SpaceDimension();
252 <<
"Number of vertices : " << GetNV() <<
'\n'
253 <<
"Number of elements : " << GetNE() <<
'\n'
254 <<
"Number of bdr elem : " << GetNBE() <<
'\n';
259 <<
"Number of vertices : " << GetNV() <<
'\n'
260 <<
"Number of elements : " << GetNE() <<
'\n'
261 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
262 <<
"h_min : " << h_min <<
'\n'
263 <<
"h_max : " << h_max <<
'\n';
268 <<
"Number of vertices : " << GetNV() <<
'\n'
269 <<
"Number of edges : " << GetNEdges() <<
'\n'
270 <<
"Number of elements : " << GetNE() <<
" -- ";
271 PrintElementsByGeometry(2, num_elems_by_geom, out);
273 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
274 <<
"Euler Number : " << EulerNumber2D() <<
'\n'
275 <<
"h_min : " << h_min <<
'\n'
276 <<
"h_max : " << h_max <<
'\n'
277 <<
"kappa_min : " << kappa_min <<
'\n'
278 <<
"kappa_max : " << kappa_max <<
'\n';
282 Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
283 num_bdr_elems_by_geom = 0;
284 for (
int i = 0; i < GetNBE(); i++)
286 num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
288 Array<int> num_faces_by_geom(Geometry::NumGeom);
289 num_faces_by_geom = 0;
290 for (
int i = 0; i < GetNFaces(); i++)
292 num_faces_by_geom[GetFaceBaseGeometry(i)]++;
296 <<
"Number of vertices : " << GetNV() <<
'\n'
297 <<
"Number of edges : " << GetNEdges() <<
'\n'
298 <<
"Number of faces : " << GetNFaces() <<
" -- ";
299 PrintElementsByGeometry(Dim-1, num_faces_by_geom, out);
301 <<
"Number of elements : " << GetNE() <<
" -- ";
302 PrintElementsByGeometry(Dim, num_elems_by_geom, out);
304 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
305 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, out);
307 <<
"Euler Number : " << EulerNumber() <<
'\n'
308 <<
"h_min : " << h_min <<
'\n'
309 <<
"h_max : " << h_max <<
'\n'
310 <<
"kappa_min : " << kappa_min <<
'\n'
311 <<
"kappa_max : " << kappa_max <<
'\n';
313 out <<
'\n' << std::flush;
320 case Element::POINT :
return &
PointFE;
321 case Element::SEGMENT :
return &
SegmentFE;
326 case Element::WEDGE :
return &
WedgeFE;
328 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
331 MFEM_ABORT(
"Unknown element type");
340 ElTr->
ElementType = ElementTransformation::ELEMENT;
345 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
351 Nodes->FESpace()->GetElementVDofs(i, vdofs);
354 int n = vdofs.
Size()/spaceDim;
356 for (
int k = 0; k < spaceDim; k++)
358 for (
int j = 0; j < n; j++)
360 pm(k,j) = nodes(vdofs[n*k+j]);
363 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
367 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
372 ElTr->
ElementType = ElementTransformation::ELEMENT;
378 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
379 int nv = elements[i]->GetNVertices();
380 const int *v = elements[i]->GetVertices();
381 int n = vertices.Size();
383 for (
int k = 0; k < spaceDim; k++)
385 for (
int j = 0; j < nv; j++)
387 pm(k, j) = nodes(k*n+v[j]);
390 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
394 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
396 Nodes->FESpace()->GetElementVDofs(i, vdofs);
397 int n = vdofs.
Size()/spaceDim;
399 for (
int k = 0; k < spaceDim; k++)
401 for (
int j = 0; j < n; j++)
403 pm(k,j) = nodes(vdofs[n*k+j]);
406 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
412 GetElementTransformation(i, &Transformation);
414 return &Transformation;
419 GetBdrElementTransformation(i, &BdrTransformation);
420 return &BdrTransformation;
427 ElTr->
ElementType = ElementTransformation::BDR_ELEMENT;
432 GetBdrPointMatrix(i, pm);
433 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
443 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
444 int n = vdofs.
Size()/spaceDim;
446 for (
int k = 0; k < spaceDim; k++)
448 for (
int j = 0; j < n; j++)
450 pm(k,j) = nodes(vdofs[n*k+j]);
457 int elem_id, face_info;
458 GetBdrElementAdjacentElement(i, elem_id, face_info);
460 GetLocalFaceTransformation(GetBdrElementType(i),
461 GetElementType(elem_id),
462 FaceElemTr.Loc1.Transf, face_info);
466 Nodes->FESpace()->GetTraceElement(elem_id,
467 GetBdrElementBaseGeometry(i));
470 FaceElemTr.Loc1.Transf.ElementNo = elem_id;
471 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
472 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
473 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
475 ElTr->
SetFE(face_el);
482 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
489 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
490 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
492 for (
int i = 0; i < spaceDim; i++)
494 for (
int j = 0; j < nv; j++)
496 pm(i, j) = vertices[v[j]](i);
499 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
503 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
509 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
510 int n = vdofs.
Size()/spaceDim;
512 for (
int i = 0; i < spaceDim; i++)
514 for (
int j = 0; j < n; j++)
516 pm(i, j) = nodes(vdofs[n*i+j]);
523 FaceInfo &face_info = faces_info[FaceNo];
528 GetLocalFaceTransformation(face_type,
529 GetElementType(face_info.
Elem1No),
530 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
533 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
537 FaceElemTr.Loc1.Transf.ElementNo = face_info.
Elem1No;
538 FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
539 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
540 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
549 GetFaceTransformation(FaceNo, &FaceTransformation);
550 return &FaceTransformation;
557 GetFaceTransformation(EdgeNo, EdTr);
562 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
573 GetEdgeVertices(EdgeNo, v);
576 for (
int i = 0; i < spaceDim; i++)
578 for (
int j = 0; j < nv; j++)
580 pm(i, j) = vertices[v[j]](i);
583 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
587 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
591 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
592 int n = vdofs.
Size()/spaceDim;
594 for (
int i = 0; i < spaceDim; i++)
596 for (
int j = 0; j < n; j++)
598 pm(i, j) = (*Nodes)(vdofs[n*i+j]);
601 EdTr->
SetFE(edge_el);
605 MFEM_ABORT(
"Not implemented.");
612 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
613 return &EdgeTransformation;
617 void Mesh::GetLocalPtToSegTransformation(
632 void Mesh::GetLocalSegToTriTransformation(
641 tv = tri_t::Edges[i/64];
642 so = seg_t::Orient[i%64];
645 for (
int j = 0; j < 2; j++)
647 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
648 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
652 void Mesh::GetLocalSegToQuadTransformation(
661 qv = quad_t::Edges[i/64];
662 so = seg_t::Orient[i%64];
665 for (
int j = 0; j < 2; j++)
667 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
668 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
672 void Mesh::GetLocalTriToTetTransformation(
680 const int *tv = tet_t::FaceVert[i/64];
683 const int *to = tri_t::Orient[i%64];
687 for (
int j = 0; j < 3; j++)
690 locpm(0, j) = vert.
x;
691 locpm(1, j) = vert.
y;
692 locpm(2, j) = vert.
z;
696 void Mesh::GetLocalTriToWdgTransformation(
704 MFEM_VERIFY(i < 128,
"Local face index " << i/64
705 <<
" is not a triangular face of a wedge.");
706 const int *pv = pri_t::FaceVert[i/64];
709 const int *to = tri_t::Orient[i%64];
713 for (
int j = 0; j < 3; j++)
716 locpm(0, j) = vert.
x;
717 locpm(1, j) = vert.
y;
718 locpm(2, j) = vert.
z;
722 void Mesh::GetLocalQuadToHexTransformation(
730 const int *hv = hex_t::FaceVert[i/64];
732 const int *qo = quad_t::Orient[i%64];
735 for (
int j = 0; j < 4; j++)
738 locpm(0, j) = vert.
x;
739 locpm(1, j) = vert.
y;
740 locpm(2, j) = vert.
z;
744 void Mesh::GetLocalQuadToWdgTransformation(
752 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
753 <<
" is not a quadrilateral face of a wedge.");
754 const int *pv = pri_t::FaceVert[i/64];
756 const int *qo = quad_t::Orient[i%64];
759 for (
int j = 0; j < 4; j++)
762 locpm(0, j) = vert.
x;
763 locpm(1, j) = vert.
y;
764 locpm(2, j) = vert.
z;
771 for (
int i = 0; i < geom_factors.Size(); i++)
783 geom_factors.Append(gf);
791 for (
int i = 0; i < face_geom_factors.Size(); i++)
804 face_geom_factors.Append(gf);
808 void Mesh::DeleteGeometricFactors()
810 for (
int i = 0; i < geom_factors.Size(); i++)
812 delete geom_factors[i];
814 geom_factors.SetSize(0);
815 for (
int i = 0; i < face_geom_factors.Size(); i++)
817 delete face_geom_factors[i];
819 face_geom_factors.SetSize(0);
822 void Mesh::GetLocalFaceTransformation(
828 GetLocalPtToSegTransformation(Transf, info);
831 case Element::SEGMENT:
832 if (elem_type == Element::TRIANGLE)
834 GetLocalSegToTriTransformation(Transf, info);
838 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
839 GetLocalSegToQuadTransformation(Transf, info);
843 case Element::TRIANGLE:
844 if (elem_type == Element::TETRAHEDRON)
846 GetLocalTriToTetTransformation(Transf, info);
850 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
851 GetLocalTriToWdgTransformation(Transf, info);
855 case Element::QUADRILATERAL:
856 if (elem_type == Element::HEXAHEDRON)
858 GetLocalQuadToHexTransformation(Transf, info);
862 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
863 GetLocalQuadToWdgTransformation(Transf, info);
872 FaceInfo &face_info = faces_info[FaceNo];
875 FaceElemTr.SetConfigurationMask(cmask);
876 FaceElemTr.Elem1 = NULL;
877 FaceElemTr.Elem2 = NULL;
881 if (mask & FaceElementTransformations::HAVE_ELEM1)
883 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
884 FaceElemTr.Elem1 = &Transformation;
891 FaceElemTr.Elem2No = face_info.
Elem2No;
892 if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
893 FaceElemTr.Elem2No >= 0)
896 if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
897 { MFEM_ABORT(
"NURBS mesh not supported!"); }
899 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
900 FaceElemTr.Elem2 = &Transformation2;
905 if (mask & FaceElementTransformations::HAVE_FACE)
907 GetFaceTransformation(FaceNo, &FaceElemTr);
912 FaceElemTr.SetGeometryType(GetFaceGeometryType(FaceNo));
916 int face_type = GetFaceElementType(FaceNo);
917 if (mask & FaceElementTransformations::HAVE_LOC1)
919 int elem_type = GetElementType(face_info.
Elem1No);
920 GetLocalFaceTransformation(face_type, elem_type,
921 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
924 if ((mask & FaceElementTransformations::HAVE_LOC2) &&
925 FaceElemTr.Elem2No >= 0)
927 int elem_type = GetElementType(face_info.
Elem2No);
928 GetLocalFaceTransformation(face_type, elem_type,
929 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
932 if (Nonconforming() && IsSlaveFace(face_info))
934 ApplyLocalSlaveTransformation(FaceElemTr, face_info,
false);
939 FaceElemTr.SetConfigurationMask(cmask);
945 double dist = FaceElemTr.CheckConsistency();
948 mfem::out <<
"\nInternal error: face id = " << FaceNo
949 <<
", dist = " << dist <<
'\n';
950 FaceElemTr.CheckConsistency(1);
951 MFEM_ABORT(
"internal error");
961 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
967 #ifdef MFEM_THREAD_SAFE
972 MFEM_ASSERT(fi.
NCFace >= 0,
"");
973 MFEM_ASSERT(nc_faces_info[fi.
NCFace].Slave,
"internal error");
984 std::swap(composition(0,0), composition(0,1));
985 std::swap(composition(1,0), composition(1,1));
1009 fn = be_to_face[BdrElemNo];
1013 fn = be_to_edge[BdrElemNo];
1017 fn = boundary[BdrElemNo]->GetVertices()[0];
1020 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
1024 tr = GetFaceElementTransformations(fn, 21);
1025 tr->
Attribute = boundary[BdrElemNo]->GetAttribute();
1027 tr->
ElementType = ElementTransformation::BDR_FACE;
1031 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const
1033 *Elem1 = faces_info[Face].
Elem1No;
1034 *Elem2 = faces_info[Face].Elem2No;
1037 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const
1039 *Inf1 = faces_info[Face].Elem1Inf;
1040 *Inf2 = faces_info[Face].Elem2Inf;
1047 case 1:
return Geometry::POINT;
1048 case 2:
return Geometry::SEGMENT;
1050 if (Face < NumOfFaces)
1052 return faces[Face]->GetGeometryType();
1055 const int nc_face_id = faces_info[Face].NCFace;
1056 MFEM_ASSERT(nc_face_id >= 0,
"parent ghost faces are not supported");
1057 return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
1059 return Geometry::INVALID;
1064 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1072 NumOfElements = NumOfBdrElements = 0;
1073 NumOfEdges = NumOfFaces = 0;
1074 nbInteriorFaces = -1;
1075 nbBoundaryFaces = -1;
1076 meshgen = mesh_geoms = 0;
1082 last_operation = Mesh::NONE;
1085 void Mesh::InitTables()
1088 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1091 void Mesh::SetEmpty()
1097 void Mesh::DestroyTables()
1102 DeleteGeometricFactors();
1113 void Mesh::DestroyPointers()
1115 if (own_nodes) {
delete Nodes; }
1121 for (
int i = 0; i < NumOfElements; i++)
1123 FreeElement(elements[i]);
1126 for (
int i = 0; i < NumOfBdrElements; i++)
1128 FreeElement(boundary[i]);
1131 for (
int i = 0; i < faces.Size(); i++)
1133 FreeElement(faces[i]);
1139 void Mesh::Destroy()
1143 elements.DeleteAll();
1144 vertices.DeleteAll();
1145 boundary.DeleteAll();
1147 faces_info.DeleteAll();
1148 nc_faces_info.DeleteAll();
1149 be_to_edge.DeleteAll();
1150 be_to_face.DeleteAll();
1158 CoarseFineTr.Clear();
1160 #ifdef MFEM_USE_MEMALLOC
1164 attributes.DeleteAll();
1165 bdr_attributes.DeleteAll();
1168 void Mesh::ResetLazyData()
1170 delete el_to_el; el_to_el = NULL;
1171 delete face_edge; face_edge = NULL;
1172 delete edge_vertex; edge_vertex = NULL;
1173 DeleteGeometricFactors();
1174 nbInteriorFaces = -1;
1175 nbBoundaryFaces = -1;
1178 void Mesh::SetAttributes()
1183 for (
int i = 0; i < attribs.
Size(); i++)
1185 attribs[i] = GetBdrAttribute(i);
1189 attribs.
Copy(bdr_attributes);
1190 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1192 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1196 for (
int i = 0; i < attribs.
Size(); i++)
1198 attribs[i] = GetAttribute(i);
1202 attribs.
Copy(attributes);
1203 if (attributes.Size() > 0 && attributes[0] <= 0)
1205 MFEM_WARNING(
"Non-positive attributes in the domain!");
1209 void Mesh::InitMesh(
int _Dim,
int _spaceDim,
int NVert,
int NElem,
int NBdrElem)
1214 spaceDim = _spaceDim;
1217 vertices.SetSize(NVert);
1220 elements.SetSize(NElem);
1222 NumOfBdrElements = 0;
1223 boundary.SetSize(NBdrElem);
1226 template<
typename T>
1227 static void CheckEnlarge(
Array<T> &array,
int size)
1229 if (size >= array.
Size()) { array.
SetSize(size + 1); }
1232 int Mesh::AddVertex(
double x,
double y,
double z)
1234 CheckEnlarge(vertices, NumOfVertices);
1235 double *v = vertices[NumOfVertices]();
1239 return NumOfVertices++;
1242 int Mesh::AddVertex(
const double *coords)
1244 CheckEnlarge(vertices, NumOfVertices);
1245 vertices[NumOfVertices].SetCoords(spaceDim, coords);
1246 return NumOfVertices++;
1249 void Mesh::AddVertexParents(
int i,
int p1,
int p2)
1255 if (i < vertices.Size())
1257 double *vi = vertices[i](), *vp1 = vertices[p1](), *vp2 = vertices[p2]();
1258 for (
int j = 0; j < 3; j++)
1260 vi[j] = (vp1[j] + vp2[j]) * 0.5;
1265 int Mesh::AddSegment(
int v1,
int v2,
int attr)
1267 CheckEnlarge(elements, NumOfElements);
1268 elements[NumOfElements] =
new Segment(v1, v2, attr);
1269 return NumOfElements++;
1272 int Mesh::AddSegment(
const int *vi,
int attr)
1274 CheckEnlarge(elements, NumOfElements);
1275 elements[NumOfElements] =
new Segment(vi, attr);
1276 return NumOfElements++;
1279 int Mesh::AddTriangle(
int v1,
int v2,
int v3,
int attr)
1281 CheckEnlarge(elements, NumOfElements);
1282 elements[NumOfElements] =
new Triangle(v1, v2, v3, attr);
1283 return NumOfElements++;
1286 int Mesh::AddTriangle(
const int *vi,
int attr)
1288 CheckEnlarge(elements, NumOfElements);
1289 elements[NumOfElements] =
new Triangle(vi, attr);
1290 return NumOfElements++;
1293 int Mesh::AddQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1295 CheckEnlarge(elements, NumOfElements);
1296 elements[NumOfElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1297 return NumOfElements++;
1300 int Mesh::AddQuad(
const int *vi,
int attr)
1302 CheckEnlarge(elements, NumOfElements);
1304 return NumOfElements++;
1307 int Mesh::AddTet(
int v1,
int v2,
int v3,
int v4,
int attr)
1309 int vi[4] = {v1, v2, v3, v4};
1310 return AddTet(vi, attr);
1313 int Mesh::AddTet(
const int *vi,
int attr)
1315 CheckEnlarge(elements, NumOfElements);
1316 #ifdef MFEM_USE_MEMALLOC
1318 tet = TetMemory.Alloc();
1321 elements[NumOfElements] = tet;
1323 elements[NumOfElements] =
new Tetrahedron(vi, attr);
1325 return NumOfElements++;
1328 int Mesh::AddWedge(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int attr)
1330 CheckEnlarge(elements, NumOfElements);
1331 elements[NumOfElements] =
new Wedge(v1, v2, v3, v4, v5, v6, attr);
1332 return NumOfElements++;
1335 int Mesh::AddWedge(
const int *vi,
int attr)
1337 CheckEnlarge(elements, NumOfElements);
1338 elements[NumOfElements] =
new Wedge(vi, attr);
1339 return NumOfElements++;
1342 int Mesh::AddHex(
int v1,
int v2,
int v3,
int v4,
int v5,
int v6,
int v7,
int v8,
1345 CheckEnlarge(elements, NumOfElements);
1346 elements[NumOfElements] =
1347 new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8, attr);
1348 return NumOfElements++;
1351 int Mesh::AddHex(
const int *vi,
int attr)
1353 CheckEnlarge(elements, NumOfElements);
1354 elements[NumOfElements] =
new Hexahedron(vi, attr);
1355 return NumOfElements++;
1358 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1360 static const int hex_to_tet[6][4] =
1362 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1363 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1367 for (
int i = 0; i < 6; i++)
1369 for (
int j = 0; j < 4; j++)
1371 ti[j] = vi[hex_to_tet[i][j]];
1377 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1379 static const int hex_to_wdg[2][6] =
1381 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1385 for (
int i = 0; i < 2; i++)
1387 for (
int j = 0; j < 6; j++)
1389 ti[j] = vi[hex_to_wdg[i][j]];
1397 CheckEnlarge(elements, NumOfElements);
1398 elements[NumOfElements] = elem;
1399 return NumOfElements++;
1404 CheckEnlarge(boundary, NumOfBdrElements);
1405 boundary[NumOfBdrElements] = elem;
1406 return NumOfBdrElements++;
1409 int Mesh::AddBdrSegment(
int v1,
int v2,
int attr)
1411 CheckEnlarge(boundary, NumOfBdrElements);
1412 boundary[NumOfBdrElements] =
new Segment(v1, v2, attr);
1413 return NumOfBdrElements++;
1416 int Mesh::AddBdrSegment(
const int *vi,
int attr)
1418 CheckEnlarge(boundary, NumOfBdrElements);
1419 boundary[NumOfBdrElements] =
new Segment(vi, attr);
1420 return NumOfBdrElements++;
1423 int Mesh::AddBdrTriangle(
int v1,
int v2,
int v3,
int attr)
1425 CheckEnlarge(boundary, NumOfBdrElements);
1426 boundary[NumOfBdrElements] =
new Triangle(v1, v2, v3, attr);
1427 return NumOfBdrElements++;
1430 int Mesh::AddBdrTriangle(
const int *vi,
int attr)
1432 CheckEnlarge(boundary, NumOfBdrElements);
1433 boundary[NumOfBdrElements] =
new Triangle(vi, attr);
1434 return NumOfBdrElements++;
1437 int Mesh::AddBdrQuad(
int v1,
int v2,
int v3,
int v4,
int attr)
1439 CheckEnlarge(boundary, NumOfBdrElements);
1440 boundary[NumOfBdrElements] =
new Quadrilateral(v1, v2, v3, v4, attr);
1441 return NumOfBdrElements++;
1444 int Mesh::AddBdrQuad(
const int *vi,
int attr)
1446 CheckEnlarge(boundary, NumOfBdrElements);
1448 return NumOfBdrElements++;
1451 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1453 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1456 for (
int i = 0; i < 2; i++)
1458 for (
int j = 0; j < 3; j++)
1460 ti[j] = vi[quad_to_tri[i][j]];
1462 AddBdrTriangle(ti, attr);
1466 void Mesh::GenerateBoundaryElements()
1469 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1473 for (i = 0; i < boundary.Size(); i++)
1475 FreeElement(boundary[i]);
1485 NumOfBdrElements = 0;
1486 for (i = 0; i < faces_info.Size(); i++)
1488 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1491 boundary.SetSize(NumOfBdrElements);
1492 be2face.
SetSize(NumOfBdrElements);
1493 for (j = i = 0; i < faces_info.Size(); i++)
1495 if (faces_info[i].Elem2No < 0)
1497 boundary[j] = faces[i]->Duplicate(
this);
1504 void Mesh::FinalizeCheck()
1506 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1507 vertices.Size() == 0,
1508 "incorrect number of vertices: preallocated: " << vertices.Size()
1509 <<
", actually added: " << NumOfVertices);
1510 MFEM_VERIFY(elements.Size() == NumOfElements,
1511 "incorrect number of elements: preallocated: " << elements.Size()
1512 <<
", actually added: " << NumOfElements);
1513 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1514 "incorrect number of boundary elements: preallocated: "
1515 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1518 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1521 CheckElementOrientation(fix_orientation);
1525 MarkTriMeshForRefinement();
1530 el_to_edge =
new Table;
1531 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1533 CheckBdrElementOrientation();
1547 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1548 bool fix_orientation)
1551 if (fix_orientation)
1553 CheckElementOrientation(fix_orientation);
1558 el_to_edge =
new Table;
1559 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1561 CheckBdrElementOrientation();
1581 GeckoProgress(
double limit) : limit(limit) { sw.Start(); }
1582 virtual bool quit()
const {
return limit > 0 && sw.UserTime() > limit; }
1585 class GeckoVerboseProgress :
public GeckoProgress
1591 GeckoVerboseProgress(
double limit) : GeckoProgress(limit) {}
1593 virtual void beginorder(
const Graph* graph,
Float cost)
const
1594 {
mfem::out <<
"Begin Gecko ordering, cost = " << cost << std::endl; }
1595 virtual void endorder(
const Graph* graph,
Float cost)
const
1596 {
mfem::out <<
"End ordering, cost = " << cost << std::endl; }
1598 virtual void beginiter(
const Graph* graph,
1601 mfem::out <<
"Iteration " << iter <<
"/" << maxiter <<
", window "
1602 << window << std::flush;
1604 virtual void enditer(
const Graph* graph,
Float mincost,
Float cost)
const
1605 {
mfem::out <<
", cost = " << cost << endl; }
1610 int iterations,
int window,
1611 int period,
int seed,
bool verbose,
1617 GeckoProgress progress(time_limit);
1618 GeckoVerboseProgress vprogress(time_limit);
1621 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1628 const Table &my_el_to_el = ElementToElementTable();
1629 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1631 const int *neighid = my_el_to_el.
GetRow(elemid);
1632 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
1634 graph.
insert_arc(elemid + 1, neighid[i] + 1);
1639 graph.
order(&functional, iterations, window, period, seed,
1640 verbose ? &vprogress : &progress);
1646 ordering[gnodeid - 1] = graph.
rank(gnodeid);
1649 return graph.
cost();
1660 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
1661 : coord(coord), dir(dir), points(points), mid(mid) {}
1663 bool operator()(
int i)
const
1665 return (points[3*i + coord] < mid) != dir;
1669 static void HilbertSort2D(
int coord1,
1672 const Array<double> &points,
int *beg,
int *end,
1673 double xmin,
double ymin,
double xmax,
double ymax)
1675 if (end - beg <= 1) {
return; }
1677 double xmid = (xmin + xmax)*0.5;
1678 double ymid = (ymin + ymax)*0.5;
1680 int coord2 = (coord1 + 1) % 2;
1683 int *p0 = beg, *p4 = end;
1684 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
1685 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
1686 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
1690 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
1691 ymin, xmin, ymid, xmid);
1693 if (p1 != p0 || p2 != p4)
1695 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
1696 xmin, ymid, xmid, ymax);
1698 if (p2 != p0 || p3 != p4)
1700 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
1701 xmid, ymid, xmax, ymax);
1705 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
1706 ymid, xmax, ymin, xmid);
1710 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
1711 const Array<double> &points,
int *beg,
int *end,
1712 double xmin,
double ymin,
double zmin,
1713 double xmax,
double ymax,
double zmax)
1715 if (end - beg <= 1) {
return; }
1717 double xmid = (xmin + xmax)*0.5;
1718 double ymid = (ymin + ymax)*0.5;
1719 double zmid = (zmin + zmax)*0.5;
1721 int coord2 = (coord1 + 1) % 3;
1722 int coord3 = (coord1 + 2) % 3;
1725 int *p0 = beg, *p8 = end;
1726 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
1727 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
1728 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
1729 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
1730 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
1731 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
1732 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
1736 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
1737 zmin, xmin, ymin, zmid, xmid, ymid);
1739 if (p1 != p0 || p2 != p8)
1741 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
1742 ymin, zmid, xmin, ymid, zmax, xmid);
1744 if (p2 != p0 || p3 != p8)
1746 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
1747 ymid, zmid, xmin, ymax, zmax, xmid);
1749 if (p3 != p0 || p4 != p8)
1751 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
1752 xmin, ymax, zmid, xmid, ymid, zmin);
1754 if (p4 != p0 || p5 != p8)
1756 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
1757 xmid, ymax, zmid, xmax, ymid, zmin);
1759 if (p5 != p0 || p6 != p8)
1761 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
1762 ymax, zmid, xmax, ymid, zmax, xmid);
1764 if (p6 != p0 || p7 != p8)
1766 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
1767 ymid, zmid, xmax, ymin, zmax, xmid);
1771 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
1772 zmid, xmax, ymin, zmin, xmid, ymid);
1778 MFEM_VERIFY(spaceDim <= 3,
"");
1781 GetBoundingBox(min, max);
1786 if (spaceDim < 3) { points = 0.0; }
1789 for (
int i = 0; i < GetNE(); i++)
1791 GetElementCenter(i, center);
1792 for (
int j = 0; j < spaceDim; j++)
1794 points[3*i + j] = center(j);
1801 indices.
Sort([&](
int a,
int b)
1802 {
return points[3*
a] < points[3*
b]; });
1804 else if (spaceDim == 2)
1807 HilbertSort2D(0,
false,
false,
1808 points, indices.
begin(), indices.
end(),
1809 min(0), min(1), max(0), max(1));
1814 HilbertSort3D(0,
false,
false,
false,
1815 points, indices.
begin(), indices.
end(),
1816 min(0), min(1), min(2), max(0), max(1), max(2));
1821 for (
int i = 0; i < GetNE(); i++)
1823 ordering[indices[i]] = i;
1828 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
1832 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
1837 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
1841 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
1871 old_elem_node_vals.SetSize(GetNE());
1872 nodes_fes = Nodes->FESpace();
1875 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1877 nodes_fes->GetElementVDofs(old_elid, old_dofs);
1879 old_elem_node_vals[old_elid] =
new Vector(vals);
1885 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
1887 int new_elid = ordering[old_elid];
1888 new_elements[new_elid] = elements[old_elid];
1893 if (reorder_vertices)
1898 vertex_ordering = -1;
1900 int new_vertex_ind = 0;
1901 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1903 int *elem_vert = elements[new_elid]->GetVertices();
1904 int nv = elements[new_elid]->GetNVertices();
1905 for (
int vi = 0; vi < nv; ++vi)
1907 int old_vertex_ind = elem_vert[vi];
1908 if (vertex_ordering[old_vertex_ind] == -1)
1910 vertex_ordering[old_vertex_ind] = new_vertex_ind;
1911 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
1921 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1923 int *elem_vert = elements[new_elid]->GetVertices();
1924 int nv = elements[new_elid]->GetNVertices();
1925 for (
int vi = 0; vi < nv; ++vi)
1927 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
1932 for (
int belid = 0; belid < GetNBE(); ++belid)
1934 int *be_vert = boundary[belid]->GetVertices();
1935 int nv = boundary[belid]->GetNVertices();
1936 for (
int vi = 0; vi < nv; ++vi)
1938 be_vert[vi] = vertex_ordering[be_vert[vi]];
1949 el_to_edge =
new Table;
1950 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1955 GetElementToFaceTable();
1965 last_operation = Mesh::NONE;
1966 nodes_fes->Update(
false);
1969 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1971 int new_elid = ordering[old_elid];
1972 nodes_fes->GetElementVDofs(new_elid, new_dofs);
1973 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
1974 delete old_elem_node_vals[old_elid];
1980 void Mesh::MarkForRefinement()
1986 MarkTriMeshForRefinement();
1990 DSTable v_to_v(NumOfVertices);
1991 GetVertexToVertexTable(v_to_v);
1992 MarkTetMeshForRefinement(v_to_v);
1997 void Mesh::MarkTriMeshForRefinement()
2002 for (
int i = 0; i < NumOfElements; i++)
2004 if (elements[i]->GetType() == Element::TRIANGLE)
2006 GetPointMatrix(i, pmat);
2007 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
2018 for (
int i = 0; i < NumOfVertices; i++)
2023 length_idx[j].one = GetLength(i, it.Column());
2024 length_idx[j].two = j;
2031 for (
int i = 0; i < NumOfEdges; i++)
2033 order[length_idx[i].two] = i;
2037 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
2042 GetEdgeOrdering(v_to_v, order);
2044 for (
int i = 0; i < NumOfElements; i++)
2046 if (elements[i]->GetType() == Element::TETRAHEDRON)
2048 elements[i]->MarkEdge(v_to_v, order);
2051 for (
int i = 0; i < NumOfBdrElements; i++)
2053 if (boundary[i]->GetType() == Element::TRIANGLE)
2055 boundary[i]->MarkEdge(v_to_v, order);
2062 if (*old_v_to_v && *old_elem_vert)
2069 if (*old_v_to_v == NULL)
2071 bool need_v_to_v =
false;
2073 for (
int i = 0; i < GetNEdges(); i++)
2079 if (dofs.
Size() > 0)
2087 *old_v_to_v =
new DSTable(NumOfVertices);
2088 GetVertexToVertexTable(*(*old_v_to_v));
2091 if (*old_elem_vert == NULL)
2093 bool need_elem_vert =
false;
2095 for (
int i = 0; i < GetNE(); i++)
2101 if (dofs.
Size() > 1)
2103 need_elem_vert =
true;
2109 *old_elem_vert =
new Table;
2110 (*old_elem_vert)->
MakeI(GetNE());
2111 for (
int i = 0; i < GetNE(); i++)
2113 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
2115 (*old_elem_vert)->MakeJ();
2116 for (
int i = 0; i < GetNE(); i++)
2118 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
2119 elements[i]->GetNVertices());
2121 (*old_elem_vert)->ShiftUpI();
2134 const int num_edge_dofs = old_dofs.
Size();
2137 const Vector onodes = *Nodes;
2141 int offset = NumOfVertices * old_dofs.
Size();
2145 if (num_edge_dofs > 0)
2147 DSTable new_v_to_v(NumOfVertices);
2148 GetVertexToVertexTable(new_v_to_v);
2150 for (
int i = 0; i < NumOfVertices; i++)
2154 const int old_i = (*old_v_to_v)(i, it.Column());
2155 const int new_i = it.Index();
2156 if (new_i == old_i) {
continue; }
2158 old_dofs.
SetSize(num_edge_dofs);
2159 new_dofs.
SetSize(num_edge_dofs);
2160 for (
int j = 0; j < num_edge_dofs; j++)
2162 old_dofs[j] = offset + old_i * num_edge_dofs + j;
2163 new_dofs[j] = offset + new_i * num_edge_dofs + j;
2167 for (
int j = 0; j < old_dofs.
Size(); j++)
2169 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2173 offset += NumOfEdges * num_edge_dofs;
2181 Table old_face_vertex;
2182 old_face_vertex.
MakeI(NumOfFaces);
2183 for (
int i = 0; i < NumOfFaces; i++)
2187 old_face_vertex.
MakeJ();
2188 for (
int i = 0; i < NumOfFaces; i++)
2190 faces[i]->GetNVertices());
2194 STable3D *faces_tbl = GetElementToFaceTable(1);
2200 for (
int i = 0; i < NumOfFaces; i++)
2202 const int *old_v = old_face_vertex.
GetRow(i);
2204 switch (old_face_vertex.
RowSize(i))
2207 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2211 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2215 new_fdofs[new_i+1] = old_dofs.
Size();
2220 for (
int i = 0; i < NumOfFaces; i++)
2222 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2225 switch (old_face_vertex.
RowSize(i))
2228 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2229 new_v = faces[new_i]->GetVertices();
2230 new_or = GetTriOrientation(old_v, new_v);
2235 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2236 new_v = faces[new_i]->GetVertices();
2237 new_or = GetQuadOrientation(old_v, new_v);
2244 for (
int j = 0; j < old_dofs.
Size(); j++)
2247 const int old_j = dof_ord[j];
2248 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2252 for (
int j = 0; j < old_dofs.
Size(); j++)
2254 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2276 for (
int i = 0; i < GetNE(); i++)
2278 const int *old_v = old_elem_vert->
GetRow(i);
2279 const int *new_v = elements[i]->GetVertices();
2285 case Geometry::SEGMENT:
2286 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2288 case Geometry::TRIANGLE:
2289 new_or = GetTriOrientation(old_v, new_v);
2291 case Geometry::SQUARE:
2292 new_or = GetQuadOrientation(old_v, new_v);
2294 case Geometry::TETRAHEDRON:
2295 new_or = GetTetOrientation(old_v, new_v);
2299 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2300 <<
" FE collection) are not supported yet!");
2304 MFEM_VERIFY(dof_ord != NULL,
2305 "FE collection '" << fec->
Name()
2306 <<
"' does not define reordering for "
2307 << Geometry::Name[
geom] <<
" elements!");
2310 for (
int j = 0; j < new_dofs.
Size(); j++)
2313 const int old_j = dof_ord[j];
2314 new_dofs[old_j] = offset + j;
2316 offset += new_dofs.
Size();
2319 for (
int j = 0; j < old_dofs.
Size(); j++)
2321 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2333 GetElementToFaceTable();
2336 CheckBdrElementOrientation();
2341 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2346 CheckBdrElementOrientation();
2351 last_operation = Mesh::NONE;
2356 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2359 CheckElementOrientation(fix_orientation);
2361 if (NumOfBdrElements == 0)
2363 GetElementToFaceTable();
2365 GenerateBoundaryElements();
2370 DSTable v_to_v(NumOfVertices);
2371 GetVertexToVertexTable(v_to_v);
2372 MarkTetMeshForRefinement(v_to_v);
2375 GetElementToFaceTable();
2378 CheckBdrElementOrientation();
2380 if (generate_edges == 1)
2382 el_to_edge =
new Table;
2383 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2397 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2398 bool fix_orientation)
2401 CheckElementOrientation(fix_orientation);
2403 if (NumOfBdrElements == 0)
2405 GetElementToFaceTable();
2407 GenerateBoundaryElements();
2410 GetElementToFaceTable();
2413 CheckBdrElementOrientation();
2415 if (generate_edges == 1)
2417 el_to_edge =
new Table;
2418 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2432 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2435 CheckElementOrientation(fix_orientation);
2437 GetElementToFaceTable();
2440 if (NumOfBdrElements == 0)
2442 GenerateBoundaryElements();
2445 CheckBdrElementOrientation();
2449 el_to_edge =
new Table;
2450 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2462 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2466 Finalize(refine, fix_orientation);
2469 void Mesh::FinalizeTopology(
bool generate_bdr)
2481 bool generate_edges =
true;
2483 if (spaceDim == 0) { spaceDim = Dim; }
2484 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2492 GetElementToFaceTable();
2494 if (NumOfBdrElements == 0 && generate_bdr)
2496 GenerateBoundaryElements();
2497 GetElementToFaceTable();
2506 if (Dim > 1 && generate_edges)
2509 if (!el_to_edge) { el_to_edge =
new Table; }
2510 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2514 if (NumOfBdrElements == 0 && generate_bdr)
2516 GenerateBoundaryElements();
2533 ncmesh->OnMeshUpdated(
this);
2536 GenerateNCFaceInfo();
2544 if (tmp_vertex_parents.Size())
2546 MFEM_VERIFY(ncmesh == NULL,
"");
2548 tmp_vertex_parents.DeleteAll();
2552 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2554 if (NURBSext || ncmesh)
2556 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
2557 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
2566 const bool check_orientation =
true;
2567 const bool curved = (Nodes != NULL);
2568 const bool may_change_topology =
2569 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
2570 ( check_orientation && fix_orientation &&
2571 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
2574 Table *old_elem_vert = NULL;
2576 if (curved && may_change_topology)
2578 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
2581 if (check_orientation)
2584 CheckElementOrientation(fix_orientation);
2588 MarkForRefinement();
2591 if (may_change_topology)
2595 DoNodeReorder(old_v_to_v, old_elem_vert);
2596 delete old_elem_vert;
2609 CheckBdrElementOrientation();
2614 if (Dim >= 2 && Dim == spaceDim)
2616 const int num_faces = GetNumFaces();
2617 for (
int i = 0; i < num_faces; i++)
2619 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
2620 faces_info[i].Elem2Inf%2 != 0,
"invalid mesh topology");
2627 double sx,
double sy,
double sz,
bool sfc_ordering)
2631 int NVert, NElem, NBdrElem;
2633 NVert = (nx+1) * (ny+1) * (nz+1);
2634 NElem = nx * ny * nz;
2635 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
2636 if (type == Element::TETRAHEDRON)
2641 else if (type == Element::WEDGE)
2644 NBdrElem += 2*nx*ny;
2647 InitMesh(3, 3, NVert, NElem, NBdrElem);
2653 for (z = 0; z <= nz; z++)
2655 coord[2] = ((double) z / nz) * sz;
2656 for (y = 0; y <= ny; y++)
2658 coord[1] = ((double) y / ny) * sy;
2659 for (x = 0; x <= nx; x++)
2661 coord[0] = ((double) x / nx) * sx;
2667 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
2670 if (sfc_ordering && type == Element::HEXAHEDRON)
2673 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
2674 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
2676 for (
int k = 0; k < nx*ny*nz; k++)
2682 ind[0] = VTX(x , y , z );
2683 ind[1] = VTX(x+1, y , z );
2684 ind[2] = VTX(x+1, y+1, z );
2685 ind[3] = VTX(x , y+1, z );
2686 ind[4] = VTX(x , y , z+1);
2687 ind[5] = VTX(x+1, y , z+1);
2688 ind[6] = VTX(x+1, y+1, z+1);
2689 ind[7] = VTX(x , y+1, z+1);
2696 for (z = 0; z < nz; z++)
2698 for (y = 0; y < ny; y++)
2700 for (x = 0; x < nx; x++)
2702 ind[0] = VTX(x , y , z );
2703 ind[1] = VTX(x+1, y , z );
2704 ind[2] = VTX(x+1, y+1, z );
2705 ind[3] = VTX(x , y+1, z );
2706 ind[4] = VTX(x , y , z+1);
2707 ind[5] = VTX(x+1, y , z+1);
2708 ind[6] = VTX(x+1, y+1, z+1);
2709 ind[7] = VTX( x, y+1, z+1);
2710 if (type == Element::TETRAHEDRON)
2712 AddHexAsTets(ind, 1);
2714 else if (type == Element::WEDGE)
2716 AddHexAsWedges(ind, 1);
2729 for (y = 0; y < ny; y++)
2731 for (x = 0; x < nx; x++)
2733 ind[0] = VTX(x , y , 0);
2734 ind[1] = VTX(x , y+1, 0);
2735 ind[2] = VTX(x+1, y+1, 0);
2736 ind[3] = VTX(x+1, y , 0);
2737 if (type == Element::TETRAHEDRON)
2739 AddBdrQuadAsTriangles(ind, 1);
2741 else if (type == Element::WEDGE)
2743 AddBdrQuadAsTriangles(ind, 1);
2752 for (y = 0; y < ny; y++)
2754 for (x = 0; x < nx; x++)
2756 ind[0] = VTX(x , y , nz);
2757 ind[1] = VTX(x+1, y , nz);
2758 ind[2] = VTX(x+1, y+1, nz);
2759 ind[3] = VTX(x , y+1, nz);
2760 if (type == Element::TETRAHEDRON)
2762 AddBdrQuadAsTriangles(ind, 6);
2764 else if (type == Element::WEDGE)
2766 AddBdrQuadAsTriangles(ind, 1);
2775 for (z = 0; z < nz; z++)
2777 for (y = 0; y < ny; y++)
2779 ind[0] = VTX(0 , y , z );
2780 ind[1] = VTX(0 , y , z+1);
2781 ind[2] = VTX(0 , y+1, z+1);
2782 ind[3] = VTX(0 , y+1, z );
2783 if (type == Element::TETRAHEDRON)
2785 AddBdrQuadAsTriangles(ind, 5);
2794 for (z = 0; z < nz; z++)
2796 for (y = 0; y < ny; y++)
2798 ind[0] = VTX(nx, y , z );
2799 ind[1] = VTX(nx, y+1, z );
2800 ind[2] = VTX(nx, y+1, z+1);
2801 ind[3] = VTX(nx, y , z+1);
2802 if (type == Element::TETRAHEDRON)
2804 AddBdrQuadAsTriangles(ind, 3);
2813 for (x = 0; x < nx; x++)
2815 for (z = 0; z < nz; z++)
2817 ind[0] = VTX(x , 0, z );
2818 ind[1] = VTX(x+1, 0, z );
2819 ind[2] = VTX(x+1, 0, z+1);
2820 ind[3] = VTX(x , 0, z+1);
2821 if (type == Element::TETRAHEDRON)
2823 AddBdrQuadAsTriangles(ind, 2);
2832 for (x = 0; x < nx; x++)
2834 for (z = 0; z < nz; z++)
2836 ind[0] = VTX(x , ny, z );
2837 ind[1] = VTX(x , ny, z+1);
2838 ind[2] = VTX(x+1, ny, z+1);
2839 ind[3] = VTX(x+1, ny, z );
2840 if (type == Element::TETRAHEDRON)
2842 AddBdrQuadAsTriangles(ind, 4);
2854 ofstream test_stream(
"debug.mesh");
2856 test_stream.close();
2865 double sx,
double sy,
2866 bool generate_edges,
bool sfc_ordering)
2875 if (type == Element::QUADRILATERAL)
2877 NumOfVertices = (nx+1) * (ny+1);
2878 NumOfElements = nx * ny;
2879 NumOfBdrElements = 2 * nx + 2 * ny;
2881 vertices.SetSize(NumOfVertices);
2882 elements.SetSize(NumOfElements);
2883 boundary.SetSize(NumOfBdrElements);
2890 for (j = 0; j < ny+1; j++)
2892 cy = ((double) j / ny) * sy;
2893 for (i = 0; i < nx+1; i++)
2895 cx = ((double) i / nx) * sx;
2896 vertices[k](0) = cx;
2897 vertices[k](1) = cy;
2906 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
2907 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
2909 for (k = 0; k < nx*ny; k++)
2913 ind[0] = i + j*(nx+1);
2914 ind[1] = i + 1 +j*(nx+1);
2915 ind[2] = i + 1 + (j+1)*(nx+1);
2916 ind[3] = i + (j+1)*(nx+1);
2923 for (j = 0; j < ny; j++)
2925 for (i = 0; i < nx; i++)
2927 ind[0] = i + j*(nx+1);
2928 ind[1] = i + 1 +j*(nx+1);
2929 ind[2] = i + 1 + (j+1)*(nx+1);
2930 ind[3] = i + (j+1)*(nx+1);
2939 for (i = 0; i < nx; i++)
2941 boundary[i] =
new Segment(i, i+1, 1);
2942 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
2945 for (j = 0; j < ny; j++)
2947 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
2948 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
2952 else if (type == Element::TRIANGLE)
2954 NumOfVertices = (nx+1) * (ny+1);
2955 NumOfElements = 2 * nx * ny;
2956 NumOfBdrElements = 2 * nx + 2 * ny;
2958 vertices.SetSize(NumOfVertices);
2959 elements.SetSize(NumOfElements);
2960 boundary.SetSize(NumOfBdrElements);
2967 for (j = 0; j < ny+1; j++)
2969 cy = ((double) j / ny) * sy;
2970 for (i = 0; i < nx+1; i++)
2972 cx = ((double) i / nx) * sx;
2973 vertices[k](0) = cx;
2974 vertices[k](1) = cy;
2981 for (j = 0; j < ny; j++)
2983 for (i = 0; i < nx; i++)
2985 ind[0] = i + j*(nx+1);
2986 ind[1] = i + 1 + (j+1)*(nx+1);
2987 ind[2] = i + (j+1)*(nx+1);
2990 ind[1] = i + 1 + j*(nx+1);
2991 ind[2] = i + 1 + (j+1)*(nx+1);
2999 for (i = 0; i < nx; i++)
3001 boundary[i] =
new Segment(i, i+1, 1);
3002 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
3005 for (j = 0; j < ny; j++)
3007 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
3008 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
3015 MFEM_ABORT(
"Unsupported element type.");
3019 CheckElementOrientation();
3021 if (generate_edges == 1)
3023 el_to_edge =
new Table;
3024 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3026 CheckBdrElementOrientation();
3035 attributes.Append(1);
3036 bdr_attributes.Append(1); bdr_attributes.Append(2);
3037 bdr_attributes.Append(3); bdr_attributes.Append(4);
3042 void Mesh::Make1D(
int n,
double sx)
3051 NumOfVertices = n + 1;
3053 NumOfBdrElements = 2;
3054 vertices.SetSize(NumOfVertices);
3055 elements.SetSize(NumOfElements);
3056 boundary.SetSize(NumOfBdrElements);
3059 for (j = 0; j < n+1; j++)
3061 vertices[j](0) = ((
double) j / n) * sx;
3065 for (j = 0; j < n; j++)
3067 elements[j] =
new Segment(j, j+1, 1);
3072 boundary[0] =
new Point(ind, 1);
3074 boundary[1] =
new Point(ind, 2);
3082 attributes.Append(1);
3083 bdr_attributes.Append(1); bdr_attributes.Append(2);
3086 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
3104 last_operation = Mesh::NONE;
3107 elements.SetSize(NumOfElements);
3108 for (
int i = 0; i < NumOfElements; i++)
3110 elements[i] = mesh.
elements[i]->Duplicate(
this);
3117 boundary.SetSize(NumOfBdrElements);
3118 for (
int i = 0; i < NumOfBdrElements; i++)
3120 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
3139 faces.SetSize(mesh.
faces.Size());
3140 for (
int i = 0; i < faces.Size(); i++)
3143 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
3177 if (dynamic_cast<const ParMesh*>(&mesh))
3189 if (mesh.
Nodes && copy_nodes)
3194 FiniteElementCollection::New(fec->
Name());
3198 Nodes->MakeOwner(fec_copy);
3199 *Nodes = *mesh.
Nodes;
3209 Mesh::Mesh(
const char *filename,
int generate_edges,
int refine,
3210 bool fix_orientation)
3219 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
3223 Load(imesh, generate_edges, refine, fix_orientation);
3227 Mesh::Mesh(std::istream &input,
int generate_edges,
int refine,
3228 bool fix_orientation)
3231 Load(input, generate_edges, refine, fix_orientation);
3234 void Mesh::ChangeVertexDataOwnership(
double *vertex_data,
int len_vertex_data,
3239 MFEM_VERIFY(len_vertex_data >= NumOfVertices * 3,
3240 "Not enough vertices in external array : "
3241 "len_vertex_data = "<< len_vertex_data <<
", "
3242 "NumOfVertices * 3 = " << NumOfVertices * 3);
3244 if (vertex_data == (
double *)(vertices.GetData()))
3246 MFEM_ASSERT(!vertices.OwnsData(),
"invalid ownership");
3251 memcpy(vertex_data, vertices.GetData(),
3252 NumOfVertices * 3 *
sizeof(double));
3255 vertices.MakeRef(reinterpret_cast<Vertex*>(vertex_data), NumOfVertices);
3258 Mesh::Mesh(
double *_vertices,
int num_vertices,
3260 int *element_attributes,
int num_elements,
3262 int *boundary_attributes,
int num_boundary_elements,
3263 int dimension,
int space_dimension)
3265 if (space_dimension == -1)
3267 space_dimension = dimension;
3270 InitMesh(dimension, space_dimension, 0, num_elements,
3271 num_boundary_elements);
3273 int element_index_stride = Geometry::NumVerts[element_type];
3274 int boundary_index_stride = num_boundary_elements > 0 ?
3275 Geometry::NumVerts[boundary_type] : 0;
3278 vertices.MakeRef(reinterpret_cast<Vertex*>(_vertices), num_vertices);
3279 NumOfVertices = num_vertices;
3281 for (
int i = 0; i < num_elements; i++)
3283 elements[i] = NewElement(element_type);
3284 elements[i]->SetVertices(element_indices + i * element_index_stride);
3285 elements[i]->SetAttribute(element_attributes[i]);
3287 NumOfElements = num_elements;
3289 for (
int i = 0; i < num_boundary_elements; i++)
3291 boundary[i] = NewElement(boundary_type);
3292 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3293 boundary[i]->SetAttribute(boundary_attributes[i]);
3295 NumOfBdrElements = num_boundary_elements;
3304 case Geometry::POINT:
return (
new Point);
3305 case Geometry::SEGMENT:
return (
new Segment);
3306 case Geometry::TRIANGLE:
return (
new Triangle);
3308 case Geometry::TETRAHEDRON:
3309 #ifdef MFEM_USE_MEMALLOC
3310 return TetMemory.Alloc();
3314 case Geometry::CUBE:
return (
new Hexahedron);
3315 case Geometry::PRISM:
return (
new Wedge);
3317 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3323 Element *Mesh::ReadElementWithoutAttr(std::istream &input)
3329 el = NewElement(geom);
3330 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3333 for (
int i = 0; i < nv; i++)
3341 void Mesh::PrintElementWithoutAttr(
const Element *el, std::ostream &
out)
3346 for (
int j = 0; j < nv; j++)
3359 el = ReadElementWithoutAttr(input);
3368 PrintElementWithoutAttr(el, out);
3371 void Mesh::SetMeshGen()
3373 meshgen = mesh_geoms = 0;
3374 for (
int i = 0; i < NumOfElements; i++)
3379 case Element::TETRAHEDRON:
3380 mesh_geoms |= (1 << Geometry::TETRAHEDRON);
3381 case Element::TRIANGLE:
3382 mesh_geoms |= (1 << Geometry::TRIANGLE);
3383 case Element::SEGMENT:
3384 mesh_geoms |= (1 << Geometry::SEGMENT);
3385 case Element::POINT:
3386 mesh_geoms |= (1 << Geometry::POINT);
3390 case Element::HEXAHEDRON:
3391 mesh_geoms |= (1 << Geometry::CUBE);
3392 case Element::QUADRILATERAL:
3393 mesh_geoms |= (1 << Geometry::SQUARE);
3394 mesh_geoms |= (1 << Geometry::SEGMENT);
3395 mesh_geoms |= (1 << Geometry::POINT);
3399 case Element::WEDGE:
3400 mesh_geoms |= (1 << Geometry::PRISM);
3401 mesh_geoms |= (1 << Geometry::SQUARE);
3402 mesh_geoms |= (1 << Geometry::TRIANGLE);
3403 mesh_geoms |= (1 << Geometry::SEGMENT);
3404 mesh_geoms |= (1 << Geometry::POINT);
3409 MFEM_ABORT(
"invalid element type: " << type);
3415 void Mesh::Loader(std::istream &input,
int generate_edges,
3416 std::string parse_tag)
3418 int curved = 0, read_gf = 1;
3419 bool finalize_topo =
true;
3423 MFEM_ABORT(
"Input stream is not open");
3430 getline(input, mesh_type);
3434 bool mfem_v10 = (mesh_type ==
"MFEM mesh v1.0");
3435 bool mfem_v11 = (mesh_type ==
"MFEM mesh v1.1");
3436 bool mfem_v12 = (mesh_type ==
"MFEM mesh v1.2");
3437 if (mfem_v10 || mfem_v11 || mfem_v12)
3443 if ( mfem_v12 && parse_tag.empty() )
3445 parse_tag =
"mfem_mesh_end";
3447 ReadMFEMMesh(input, mfem_v11, curved);
3449 else if (mesh_type ==
"linemesh")
3451 ReadLineMesh(input);
3453 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
3455 if (mesh_type ==
"curved_areamesh2")
3459 ReadNetgen2DMesh(input, curved);
3461 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
3463 ReadNetgen3DMesh(input);
3465 else if (mesh_type ==
"TrueGrid")
3467 ReadTrueGridMesh(input);
3469 else if (mesh_type ==
"# vtk DataFile Version 3.0" ||
3470 mesh_type ==
"# vtk DataFile Version 2.0")
3472 ReadVTKMesh(input, curved, read_gf, finalize_topo);
3474 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
3476 ReadNURBSMesh(input, curved, read_gf);
3478 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
3480 ReadInlineMesh(input, generate_edges);
3483 else if (mesh_type ==
"$MeshFormat")
3485 ReadGmshMesh(input, curved, read_gf);
3488 ((mesh_type.size() > 2 &&
3489 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
3490 (mesh_type.size() > 3 &&
3491 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
3496 #ifdef MFEM_USE_NETCDF
3497 ReadCubit(mesh_input->
filename.c_str(), curved, read_gf);
3499 MFEM_ABORT(
"NetCDF support requires configuration with"
3500 " MFEM_USE_NETCDF=YES");
3506 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
3507 " Use mfem::named_ifgzstream for input.");
3513 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
3541 if (curved && read_gf)
3545 spaceDim = Nodes->VectorDim();
3546 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
3548 for (
int i = 0; i < spaceDim; i++)
3551 Nodes->GetNodalValues(vert_val, i+1);
3552 for (
int j = 0; j < NumOfVertices; j++)
3554 vertices[j](i) = vert_val(j);
3567 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
3568 getline(input, line);
3574 if (line ==
"mfem_mesh_end") {
break; }
3576 while (line != parse_tag);
3582 Mesh::Mesh(
Mesh *mesh_array[],
int num_pieces)
3584 int i, j, ie, ib, iv, *v, nv;
3593 if (mesh_array[0]->NURBSext)
3598 NumOfVertices = NURBSext->GetNV();
3599 NumOfElements = NURBSext->GetNE();
3601 NURBSext->GetElementTopo(elements);
3614 NumOfBdrElements = 0;
3615 for (i = 0; i < num_pieces; i++)
3617 NumOfBdrElements += mesh_array[i]->
GetNBE();
3619 boundary.SetSize(NumOfBdrElements);
3620 vertices.SetSize(NumOfVertices);
3622 for (i = 0; i < num_pieces; i++)
3628 for (j = 0; j < m->
GetNE(); j++)
3630 elements[lelem_elem[j]]->SetAttribute(m->
GetAttribute(j));
3633 for (j = 0; j < m->
GetNBE(); j++)
3638 for (
int k = 0; k < nv; k++)
3640 v[k] = lvert_vert[v[k]];
3642 boundary[ib++] = el;
3645 for (j = 0; j < m->
GetNV(); j++)
3655 NumOfBdrElements = 0;
3657 for (i = 0; i < num_pieces; i++)
3660 NumOfElements += m->
GetNE();
3661 NumOfBdrElements += m->
GetNBE();
3662 NumOfVertices += m->
GetNV();
3664 elements.SetSize(NumOfElements);
3665 boundary.SetSize(NumOfBdrElements);
3666 vertices.SetSize(NumOfVertices);
3668 for (i = 0; i < num_pieces; i++)
3672 for (j = 0; j < m->
GetNE(); j++)
3677 for (
int k = 0; k < nv; k++)
3681 elements[ie++] = el;
3684 for (j = 0; j < m->
GetNBE(); j++)
3689 for (
int k = 0; k < nv; k++)
3693 boundary[ib++] = el;
3696 for (j = 0; j < m->
GetNV(); j++)
3710 for (i = 0; i < num_pieces; i++)
3712 gf_array[i] = mesh_array[i]->
GetNodes();
3719 CheckElementOrientation(
false);
3720 CheckBdrElementOrientation(
false);
3724 Mesh::Mesh(
Mesh *orig_mesh,
int ref_factor,
int ref_type)
3727 MFEM_VERIFY(ref_factor >= 1,
"the refinement factor must be >= 1");
3728 MFEM_VERIFY(ref_type == BasisType::ClosedUniform ||
3729 ref_type == BasisType::GaussLobatto,
"invalid refinement type");
3730 MFEM_VERIFY(Dim == 1 || Dim == 2 || Dim == 3,
3731 "only implemented for Segment, Quadrilateral and Hexahedron "
3732 "elements in 1D/2D/3D");
3734 "meshes with mixed elements are not supported");
3741 int r_bndr_factor = pow(ref_factor, Dim - 1);
3742 int r_elem_factor = ref_factor * r_bndr_factor;
3745 int r_num_elem = orig_mesh->
GetNE() * r_elem_factor;
3746 int r_num_bndr = orig_mesh->
GetNBE() * r_bndr_factor;
3748 InitMesh(Dim, orig_mesh->
SpaceDimension(), r_num_vert, r_num_elem,
3752 NumOfVertices = r_num_vert;
3758 DenseMatrix node_coordinates(spaceDim*pow(2, Dim), r_num_elem);
3761 for (
int el = 0; el < orig_mesh->
GetNE(); el++)
3765 int nvert = Geometry::NumVerts[
geom];
3768 max_nv = std::max(max_nv, nvert);
3774 const int *c2h_map = rfec.GetDofMap(geom);
3775 const int *vertex_map = vertex_fec.
GetDofMap(geom);
3776 for (
int i = 0; i < phys_pts.
Width(); i++)
3778 vertices[rdofs[i]].SetCoords(spaceDim, phys_pts.
GetColumn(i));
3782 Element *elem = NewElement(geom);
3785 for (
int k = 0; k < nvert; k++)
3788 v[k] = rdofs[c2h_map[cid]];
3790 for (
int k = 0; k < nvert; k++)
3792 for (
int j = 0; j < spaceDim; ++j)
3794 node_coordinates(k*spaceDim + j, NumOfElements)
3795 = vertices[v[vertex_map[k]]](j);
3809 SetCurvature(1, discont, spaceDim, dof_ordering);
3810 Nodes->ProjectGridFunction(nodes_dg);
3814 for (
int el = 0; el < orig_mesh->
GetNBE(); el++)
3818 int nvert = Geometry::NumVerts[
geom];
3829 Element *elem = NewElement(geom);
3832 v[0] = rdofs[RG.
RefGeoms[nvert*j]];
3833 AddBdrElement(elem);
3838 const int *c2h_map = rfec.GetDofMap(geom);
3841 Element *elem = NewElement(geom);
3844 for (
int k = 0; k < nvert; k++)
3847 v[k] = rdofs[c2h_map[cid]];
3849 AddBdrElement(elem);
3854 FinalizeTopology(
false);
3856 last_operation = Mesh::REFINE;
3859 CoarseFineTr.embeddings.SetSize(GetNE());
3860 if (orig_mesh->
GetNE() > 0)
3864 CoarseFineTr.point_matrices[
geom].SetSize(Dim, max_nv, r_elem_factor);
3865 int nvert = Geometry::NumVerts[
geom];
3867 const int *c2h_map = rfec.GetDofMap(geom);
3872 for (
int k = 0; k < nvert; k++)
3880 for (
int el = 0; el < GetNE(); el++)
3882 Embedding &emb = CoarseFineTr.embeddings[el];
3883 emb.
parent = el / r_elem_factor;
3884 emb.
matrix = el % r_elem_factor;
3887 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
3888 MFEM_ASSERT(CheckBdrElementOrientation(
false) == 0,
"");
3893 if (NURBSext == NULL)
3895 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
3898 if (kv.
Size() != NURBSext->GetNKV())
3900 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
3903 NURBSext->ConvertToPatches(*Nodes);
3905 NURBSext->KnotInsert(kv);
3907 last_operation = Mesh::NONE;
3915 if (NURBSext == NULL)
3917 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
3920 if (kv.
Size() != NURBSext->GetNKV())
3922 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
3925 NURBSext->ConvertToPatches(*Nodes);
3927 NURBSext->KnotInsert(kv);
3929 last_operation = Mesh::NONE;
3935 void Mesh::NURBSUniformRefinement()
3938 NURBSext->ConvertToPatches(*Nodes);
3940 NURBSext->UniformRefinement();
3942 last_operation = Mesh::NONE;
3948 void Mesh::DegreeElevate(
int rel_degree,
int degree)
3950 if (NURBSext == NULL)
3952 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
3955 NURBSext->ConvertToPatches(*Nodes);
3957 NURBSext->DegreeElevate(rel_degree, degree);
3959 last_operation = Mesh::NONE;
3965 void Mesh::UpdateNURBS()
3969 NURBSext->SetKnotsFromPatches();
3971 Dim = NURBSext->Dimension();
3974 if (NumOfElements != NURBSext->GetNE())
3976 for (
int i = 0; i < elements.Size(); i++)
3978 FreeElement(elements[i]);
3980 NumOfElements = NURBSext->GetNE();
3981 NURBSext->GetElementTopo(elements);
3984 if (NumOfBdrElements != NURBSext->GetNBE())
3986 for (
int i = 0; i < boundary.Size(); i++)
3988 FreeElement(boundary[i]);
3990 NumOfBdrElements = NURBSext->GetNBE();
3991 NURBSext->GetBdrElementTopo(boundary);
3994 Nodes->FESpace()->Update();
3996 NURBSext->SetCoordsFromPatches(*Nodes);
3998 if (NumOfVertices != NURBSext->GetNV())
4000 NumOfVertices = NURBSext->GetNV();
4001 vertices.SetSize(NumOfVertices);
4002 int vd = Nodes->VectorDim();
4003 for (
int i = 0; i < vd; i++)
4006 Nodes->GetNodalValues(vert_val, i+1);
4007 for (
int j = 0; j < NumOfVertices; j++)
4009 vertices[j](i) = vert_val(j);
4016 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
4025 GetElementToFaceTable();
4030 void Mesh::LoadPatchTopo(std::istream &input,
Array<int> &edge_to_knot)
4046 input >> NumOfElements;
4047 elements.SetSize(NumOfElements);
4048 for (
int j = 0; j < NumOfElements; j++)
4050 elements[j] = ReadElement(input);
4056 input >> NumOfBdrElements;
4057 boundary.SetSize(NumOfBdrElements);
4058 for (
int j = 0; j < NumOfBdrElements; j++)
4060 boundary[j] = ReadElement(input);
4066 input >> NumOfEdges;
4067 edge_vertex =
new Table(NumOfEdges, 2);
4068 edge_to_knot.
SetSize(NumOfEdges);
4069 for (
int j = 0; j < NumOfEdges; j++)
4071 int *v = edge_vertex->GetRow(j);
4072 input >> edge_to_knot[j] >> v[0] >> v[1];
4075 edge_to_knot[j] = -1 - edge_to_knot[j];
4082 input >> NumOfVertices;
4083 vertices.SetSize(0);
4086 CheckBdrElementOrientation();
4093 for (
int d = 0; d < v.
Size(); d++)
4101 for (d = 0; d < p.
Size(); d++)
4105 for ( ; d < v.
Size(); d++)
4114 if (Nodes == NULL || Nodes->FESpace() != nodes.
FESpace())
4129 SetNodalGridFunction(nodes,
true);
4132 void Mesh::EnsureNodes()
4137 if (dynamic_cast<const H1_FECollection*>(fec)
4138 || dynamic_cast<const L2_FECollection*>(fec))
4144 const int order = GetNodalFESpace()->GetOrder(0);
4145 SetCurvature(order,
false, -1, Ordering::byVDIM);
4150 SetCurvature(1,
false, -1, Ordering::byVDIM);
4157 NewNodes(*nodes, make_owner);
4162 return ((Nodes) ? Nodes->FESpace() : NULL);
4165 void Mesh::SetCurvature(
int order,
bool discont,
int space_dim,
int ordering)
4167 space_dim = (space_dim == -1) ? spaceDim : space_dim;
4180 SetNodalFESpace(nfes);
4181 Nodes->MakeOwner(nfec);
4184 int Mesh::GetNumFaces()
const
4188 case 1:
return GetNV();
4189 case 2:
return GetNEdges();
4190 case 3:
return GetNFaces();
4195 static int CountFacesByType(
const Mesh &mesh,
const FaceType type)
4204 if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
4205 (type==FaceType::Boundary && e2<0 && inf2<0) ) { nf++; }
4212 const bool isInt = type==FaceType::Interior;
4213 int &nf = isInt ? nbInteriorFaces : nbBoundaryFaces;
4214 if (nf<0) { nf = CountFacesByType(*
this, type); }
4218 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
4219 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
4222 int Mesh::CheckElementOrientation(
bool fix_it)
4224 int i, j, k, wo = 0, fo = 0, *vi = 0;
4227 if (Dim == 2 && spaceDim == 2)
4231 for (i = 0; i < NumOfElements; i++)
4235 vi = elements[i]->GetVertices();
4236 for (j = 0; j < 3; j++)
4238 v[j] = vertices[vi[j]]();
4240 for (j = 0; j < 2; j++)
4241 for (k = 0; k < 2; k++)
4243 J(j, k) = v[j+1][k] - v[0][k];
4249 GetElementJacobian(i, J);
4255 switch (GetElementType(i))
4257 case Element::TRIANGLE:
4260 case Element::QUADRILATERAL:
4264 MFEM_ABORT(
"Invalid 2D element type \""
4265 << GetElementType(i) <<
"\"");
4279 for (i = 0; i < NumOfElements; i++)
4281 vi = elements[i]->GetVertices();
4282 switch (GetElementType(i))
4284 case Element::TETRAHEDRON:
4287 for (j = 0; j < 4; j++)
4289 v[j] = vertices[vi[j]]();
4291 for (j = 0; j < 3; j++)
4292 for (k = 0; k < 3; k++)
4294 J(j, k) = v[j+1][k] - v[0][k];
4300 GetElementJacobian(i, J);
4313 case Element::WEDGE:
4315 GetElementJacobian(i, J);
4326 case Element::HEXAHEDRON:
4328 GetElementJacobian(i, J);
4340 MFEM_ABORT(
"Invalid 3D element type \""
4341 << GetElementType(i) <<
"\"");
4346 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
4349 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
4350 << NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
4357 int Mesh::GetTriOrientation(
const int *base,
const int *test)
4365 if (test[0] == base[0])
4366 if (test[1] == base[1])
4374 else if (test[0] == base[1])
4375 if (test[1] == base[0])
4384 if (test[1] == base[0])
4394 const int *aor = tri_t::Orient[orient];
4395 for (
int j = 0; j < 3; j++)
4396 if (test[aor[j]] != base[j])
4405 int Mesh::GetQuadOrientation(
const int *base,
const int *test)
4409 for (i = 0; i < 4; i++)
4410 if (test[i] == base[0])
4417 if (test[(i+1)%4] == base[1])
4425 const int *aor = quad_t::Orient[orient];
4426 for (
int j = 0; j < 4; j++)
4427 if (test[aor[j]] != base[j])
4429 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
4431 for (
int k = 0; k < 4; k++)
4436 for (
int k = 0; k < 4; k++)
4445 if (test[(i+1)%4] == base[1])
4453 int Mesh::GetTetOrientation(
const int *base,
const int *test)
4461 if (test[0] == base[0])
4462 if (test[1] == base[1])
4463 if (test[2] == base[2])
4471 else if (test[2] == base[1])
4472 if (test[3] == base[2])
4481 if (test[1] == base[2])
4489 else if (test[1] == base[0])
4490 if (test[2] == base[1])
4491 if (test[0] == base[2])
4499 else if (test[3] == base[1])
4500 if (test[2] == base[2])
4509 if (test[3] == base[2])
4517 else if (test[2] == base[0])
4518 if (test[3] == base[1])
4519 if (test[0] == base[2])
4527 else if (test[0] == base[1])
4528 if (test[1] == base[2])
4537 if (test[3] == base[2])
4546 if (test[0] == base[1])
4547 if (test[2] == base[2])
4555 else if (test[1] == base[1])
4556 if (test[0] == base[2])
4565 if (test[1] == base[2])
4575 const int *aor = tet_t::Orient[orient];
4576 for (
int j = 0; j < 4; j++)
4577 if (test[aor[j]] != base[j])
4586 int Mesh::CheckBdrElementOrientation(
bool fix_it)
4592 if (el_to_edge == NULL)
4594 el_to_edge =
new Table;
4595 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
4598 for (
int i = 0; i < NumOfBdrElements; i++)
4600 if (faces_info[be_to_edge[i]].Elem2No < 0)
4602 int *bv = boundary[i]->GetVertices();
4603 int *fv = faces[be_to_edge[i]]->GetVertices();
4608 mfem::Swap<int>(bv[0], bv[1]);
4618 for (
int i = 0; i < NumOfBdrElements; i++)
4620 const int fi = be_to_face[i];
4622 if (faces_info[fi].Elem2No >= 0) {
continue; }
4625 int *bv = boundary[i]->GetVertices();
4627 MFEM_ASSERT(fi < faces.Size(),
"internal error");
4628 const int *fv = faces[fi]->GetVertices();
4634 case Element::TRIANGLE:
4636 orientation = GetTriOrientation(fv, bv);
4639 case Element::QUADRILATERAL:
4641 orientation = GetQuadOrientation(fv, bv);
4645 MFEM_ABORT(
"Invalid 2D boundary element type \""
4646 << bdr_type <<
"\"");
4651 if (orientation % 2 == 0) {
continue; }
4653 if (!fix_it) {
continue; }
4657 case Element::TRIANGLE:
4661 mfem::Swap<int>(bv[0], bv[1]);
4664 int *be = bel_to_edge->GetRow(i);
4665 mfem::Swap<int>(be[1], be[2]);
4669 case Element::QUADRILATERAL:
4671 mfem::Swap<int>(bv[0], bv[2]);
4674 int *be = bel_to_edge->GetRow(i);
4675 mfem::Swap<int>(be[0], be[1]);
4676 mfem::Swap<int>(be[2], be[3]);
4689 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
4690 << NumOfBdrElements <<
" (" << fixed_or_not[fix_it ? 0 : 1]
4697 int Mesh::GetNumGeometries(
int dim)
const
4699 MFEM_ASSERT(0 <= dim && dim <= Dim,
"invalid dim: " << dim);
4701 for (
int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
4710 MFEM_ASSERT(0 <= dim && dim <= Dim,
"invalid dim: " << dim);
4712 for (
int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
4725 el_to_edge->GetRow(i, edges);
4729 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
4730 "is not generated.");
4733 const int *v = elements[i]->GetVertices();
4734 const int ne = elements[i]->GetNEdges();
4736 for (
int j = 0; j < ne; j++)
4738 const int *e = elements[i]->GetEdgeVertices(j);
4739 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4749 edges[0] = be_to_edge[i];
4750 const int *v = boundary[i]->GetVertices();
4751 cor[0] = (v[0] < v[1]) ? (1) : (-1);
4757 bel_to_edge->GetRow(i, edges);
4764 const int *v = boundary[i]->GetVertices();
4765 const int ne = boundary[i]->GetNEdges();
4767 for (
int j = 0; j < ne; j++)
4769 const int *e = boundary[i]->GetEdgeVertices(j);
4770 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4782 const int *v = faces[i]->GetVertices();
4783 o[0] = (v[0] < v[1]) ? (1) : (-1);
4793 face_edge->GetRow(i, edges);
4795 const int *v = faces[i]->GetVertices();
4796 const int ne = faces[i]->GetNEdges();
4798 for (
int j = 0; j < ne; j++)
4800 const int *e = faces[i]->GetEdgeVertices(j);
4801 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4810 if (!edge_vertex) { GetEdgeVertexTable(); }
4811 edge_vertex->GetRow(i, vert);
4827 if (faces.Size() != NumOfFaces)
4829 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
4833 DSTable v_to_v(NumOfVertices);
4834 GetVertexToVertexTable(v_to_v);
4836 face_edge =
new Table;
4837 GetElementArrayEdgeTable(faces, v_to_v, *face_edge);
4849 DSTable v_to_v(NumOfVertices);
4850 GetVertexToVertexTable(v_to_v);
4853 edge_vertex =
new Table(nedges, 2);
4854 for (
int i = 0; i < NumOfVertices; i++)
4859 edge_vertex->Push(j, i);
4860 edge_vertex->Push(j, it.Column());
4863 edge_vertex->Finalize();
4874 vert_elem->
MakeI(NumOfVertices);
4876 for (i = 0; i < NumOfElements; i++)
4878 nv = elements[i]->GetNVertices();
4879 v = elements[i]->GetVertices();
4880 for (j = 0; j < nv; j++)
4888 for (i = 0; i < NumOfElements; i++)
4890 nv = elements[i]->GetNVertices();
4891 v = elements[i]->GetVertices();
4892 for (j = 0; j < nv; j++)
4903 Table *Mesh::GetFaceToElementTable()
const
4907 face_elem->
MakeI(faces_info.Size());
4909 for (
int i = 0; i < faces_info.Size(); i++)
4911 if (faces_info[i].Elem2No >= 0)
4923 for (
int i = 0; i < faces_info.Size(); i++)
4926 if (faces_info[i].Elem2No >= 0)
4944 el_to_face->
GetRow(i, fcs);
4948 mfem_error(
"Mesh::GetElementFaces(...) : el_to_face not generated.");
4953 for (j = 0; j < n; j++)
4954 if (faces_info[fcs[j]].Elem1No == i)
4956 cor[j] = faces_info[fcs[j]].Elem1Inf % 64;
4959 else if (faces_info[fcs[j]].Elem2No == i)
4961 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
4965 mfem_error(
"Mesh::GetElementFaces(...) : 2");
4970 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
4975 void Mesh::GetBdrElementFace(
int i,
int *
f,
int *o)
const
4980 bv = boundary[i]->GetVertices();
4981 fv = faces[be_to_face[i]]->GetVertices();
4985 switch (GetBdrElementType(i))
4987 case Element::TRIANGLE:
4988 *o = GetTriOrientation(fv, bv);
4990 case Element::QUADRILATERAL:
4991 *o = GetQuadOrientation(fv, bv);
4994 mfem_error(
"Mesh::GetBdrElementFace(...) 2");
4998 int Mesh::GetBdrElementEdgeIndex(
int i)
const
5002 case 1:
return boundary[i]->GetVertices()[0];
5003 case 2:
return be_to_edge[i];
5004 case 3:
return be_to_face[i];
5005 default:
mfem_error(
"Mesh::GetBdrElementEdgeIndex: invalid dimension!");
5010 void Mesh::GetBdrElementAdjacentElement(
int bdr_el,
int &el,
int &info)
const
5012 int fid = GetBdrElementEdgeIndex(bdr_el);
5013 const FaceInfo &fi = faces_info[fid];
5014 MFEM_ASSERT(fi.
Elem1Inf%64 == 0,
"internal error");
5015 const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
5016 const int *bv = boundary[bdr_el]->GetVertices();
5018 switch (GetBdrElementBaseGeometry(bdr_el))
5020 case Geometry::POINT: ori = 0;
break;
5021 case Geometry::SEGMENT: ori = (fv[0] == bv[0]) ? 0 : 1;
break;
5022 case Geometry::TRIANGLE: ori = GetTriOrientation(fv, bv);
break;
5023 case Geometry::SQUARE: ori = GetQuadOrientation(fv, bv);
break;
5024 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
5032 return elements[i]->GetType();
5037 return boundary[i]->GetType();
5045 v = elements[i]->GetVertices();
5046 nv = elements[i]->GetNVertices();
5048 pointmat.
SetSize(spaceDim, nv);
5049 for (k = 0; k < spaceDim; k++)
5051 for (j = 0; j < nv; j++)
5053 pointmat(k, j) = vertices[v[j]](k);
5063 v = boundary[i]->GetVertices();
5064 nv = boundary[i]->GetNVertices();
5066 pointmat.
SetSize(spaceDim, nv);
5067 for (k = 0; k < spaceDim; k++)
5068 for (j = 0; j < nv; j++)
5070 pointmat(k, j) = vertices[v[j]](k);
5074 double Mesh::GetLength(
int i,
int j)
const
5076 const double *vi = vertices[i]();
5077 const double *vj = vertices[j]();
5080 for (
int k = 0; k < spaceDim; k++)
5082 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
5085 return sqrt(length);
5093 for (
int i = 0; i < elem_array.
Size(); i++)
5098 for (
int i = 0; i < elem_array.
Size(); i++)
5100 const int *v = elem_array[i]->GetVertices();
5101 const int ne = elem_array[i]->GetNEdges();
5102 for (
int j = 0; j < ne; j++)
5104 const int *e = elem_array[i]->GetEdgeVertices(j);
5111 void Mesh::GetVertexToVertexTable(
DSTable &v_to_v)
const
5115 for (
int i = 0; i < edge_vertex->Size(); i++)
5117 const int *v = edge_vertex->GetRow(i);
5118 v_to_v.
Push(v[0], v[1]);
5123 for (
int i = 0; i < NumOfElements; i++)
5125 const int *v = elements[i]->GetVertices();
5126 const int ne = elements[i]->GetNEdges();
5127 for (
int j = 0; j < ne; j++)
5129 const int *e = elements[i]->GetEdgeVertices(j);
5130 v_to_v.
Push(v[e[0]], v[e[1]]);
5138 int i, NumberOfEdges;
5140 DSTable v_to_v(NumOfVertices);
5141 GetVertexToVertexTable(v_to_v);
5146 GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
5151 be_to_f.
SetSize(NumOfBdrElements);
5152 for (i = 0; i < NumOfBdrElements; i++)
5154 const int *v = boundary[i]->GetVertices();
5155 be_to_f[i] = v_to_v(v[0], v[1]);
5160 if (bel_to_edge == NULL)
5162 bel_to_edge =
new Table;
5164 GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
5168 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
5172 return NumberOfEdges;
5175 const Table & Mesh::ElementToElementTable()
5183 MFEM_ASSERT(faces_info.Size() >= GetNumFaces(),
"faces were not generated!");
5186 conn.
Reserve(2*faces_info.Size());
5188 for (
int i = 0; i < faces_info.Size(); i++)
5190 const FaceInfo &fi = faces_info[i];
5198 int nbr_elem_idx = NumOfElements - 1 - fi.
Elem2No;
5206 el_to_el =
new Table(NumOfElements, conn);
5211 const Table & Mesh::ElementToFaceTable()
const
5213 if (el_to_face == NULL)
5220 const Table & Mesh::ElementToEdgeTable()
const
5222 if (el_to_edge == NULL)
5229 void Mesh::AddPointFaceElement(
int lf,
int gf,
int el)
5231 if (faces_info[gf].Elem1No == -1)
5234 faces_info[gf].Elem1No = el;
5235 faces_info[gf].Elem1Inf = 64 * lf;
5236 faces_info[gf].Elem2No = -1;
5237 faces_info[gf].Elem2Inf = -1;
5241 faces_info[gf].Elem2No = el;
5242 faces_info[gf].Elem2Inf = 64 * lf + 1;
5246 void Mesh::AddSegmentFaceElement(
int lf,
int gf,
int el,
int v0,
int v1)
5248 if (faces[gf] == NULL)
5250 faces[gf] =
new Segment(v0, v1);
5251 faces_info[gf].Elem1No = el;
5252 faces_info[gf].Elem1Inf = 64 * lf;
5253 faces_info[gf].Elem2No = -1;
5254 faces_info[gf].Elem2Inf = -1;
5258 int *v = faces[gf]->GetVertices();
5259 faces_info[gf].Elem2No = el;
5260 if ( v[1] == v0 && v[0] == v1 )
5262 faces_info[gf].Elem2Inf = 64 * lf + 1;
5264 else if ( v[0] == v0 && v[1] == v1 )
5270 faces_info[gf].Elem2Inf = 64 * lf;
5274 MFEM_ABORT(
"internal error");
5279 void Mesh::AddTriangleFaceElement(
int lf,
int gf,
int el,
5280 int v0,
int v1,
int v2)
5282 if (faces[gf] == NULL)
5284 faces[gf] =
new Triangle(v0, v1, v2);
5285 faces_info[gf].Elem1No = el;
5286 faces_info[gf].Elem1Inf = 64 * lf;
5287 faces_info[gf].Elem2No = -1;
5288 faces_info[gf].Elem2Inf = -1;
5292 int orientation, vv[3] = { v0, v1, v2 };
5293 orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
5298 faces_info[gf].Elem2No = el;
5299 faces_info[gf].Elem2Inf = 64 * lf + orientation;
5303 void Mesh::AddQuadFaceElement(
int lf,
int gf,
int el,
5304 int v0,
int v1,
int v2,
int v3)
5306 if (faces_info[gf].Elem1No < 0)
5309 faces_info[gf].Elem1No = el;
5310 faces_info[gf].Elem1Inf = 64 * lf;
5311 faces_info[gf].Elem2No = -1;
5312 faces_info[gf].Elem2Inf = -1;
5316 int vv[4] = { v0, v1, v2, v3 };
5317 int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
5321 faces_info[gf].Elem2No = el;
5322 faces_info[gf].Elem2Inf = 64 * lf + oo;
5326 void Mesh::GenerateFaces()
5328 int i, nfaces = GetNumFaces();
5330 for (i = 0; i < faces.Size(); i++)
5332 FreeElement(faces[i]);
5336 faces.SetSize(nfaces);
5337 faces_info.SetSize(nfaces);
5338 for (i = 0; i < nfaces; i++)
5341 faces_info[i].Elem1No = -1;
5342 faces_info[i].NCFace = -1;
5344 for (i = 0; i < NumOfElements; i++)
5346 const int *v = elements[i]->GetVertices();
5350 AddPointFaceElement(0, v[0], i);
5351 AddPointFaceElement(1, v[1], i);
5355 ef = el_to_edge->GetRow(i);
5356 const int ne = elements[i]->GetNEdges();
5357 for (
int j = 0; j < ne; j++)
5359 const int *e = elements[i]->GetEdgeVertices(j);
5360 AddSegmentFaceElement(j, ef[j], i, v[e[0]], v[e[1]]);
5365 ef = el_to_face->GetRow(i);
5366 switch (GetElementType(i))
5368 case Element::TETRAHEDRON:
5370 for (
int j = 0; j < 4; j++)
5372 const int *fv = tet_t::FaceVert[j];
5373 AddTriangleFaceElement(j, ef[j], i,
5374 v[fv[0]], v[fv[1]], v[fv[2]]);
5378 case Element::WEDGE:
5380 for (
int j = 0; j < 2; j++)
5382 const int *fv = pri_t::FaceVert[j];
5383 AddTriangleFaceElement(j, ef[j], i,
5384 v[fv[0]], v[fv[1]], v[fv[2]]);
5386 for (
int j = 2; j < 5; j++)
5388 const int *fv = pri_t::FaceVert[j];
5389 AddQuadFaceElement(j, ef[j], i,
5390 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5394 case Element::HEXAHEDRON:
5396 for (
int j = 0; j < 6; j++)
5398 const int *fv = hex_t::FaceVert[j];
5399 AddQuadFaceElement(j, ef[j], i,
5400 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5405 MFEM_ABORT(
"Unexpected type of Element.");
5411 void Mesh::GenerateNCFaceInfo()
5413 MFEM_VERIFY(ncmesh,
"missing NCMesh.");
5415 for (
int i = 0; i < faces_info.Size(); i++)
5417 faces_info[i].NCFace = -1;
5421 (Dim == 2) ? ncmesh->GetEdgeList() : ncmesh->GetFaceList();
5423 nc_faces_info.SetSize(0);
5424 nc_faces_info.Reserve(list.
masters.Size() + list.
slaves.Size());
5426 int nfaces = GetNumFaces();
5429 for (
int i = 0; i < list.
masters.Size(); i++)
5432 if (master.
index >= nfaces) {
continue; }
5434 faces_info[master.
index].NCFace = nc_faces_info.Size();
5440 for (
int i = 0; i < list.
slaves.Size(); i++)
5444 if (slave.
index < 0 ||
5445 slave.
index >= nfaces ||
5455 slave_fi.
NCFace = nc_faces_info.Size();
5463 nc_faces_info.Append(
5472 for (
int i = 0; i < NumOfElements; i++)
5474 const int *v = elements[i]->GetVertices();
5475 switch (GetElementType(i))
5477 case Element::TETRAHEDRON:
5479 for (
int j = 0; j < 4; j++)
5481 const int *fv = tet_t::FaceVert[j];
5482 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
5486 case Element::WEDGE:
5488 for (
int j = 0; j < 2; j++)
5490 const int *fv = pri_t::FaceVert[j];
5491 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
5493 for (
int j = 2; j < 5; j++)
5495 const int *fv = pri_t::FaceVert[j];
5496 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5500 case Element::HEXAHEDRON:
5504 for (
int j = 0; j < 6; j++)
5506 const int *fv = hex_t::FaceVert[j];
5507 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5512 MFEM_ABORT(
"Unexpected type of Element.");
5523 if (el_to_face != NULL)
5527 el_to_face =
new Table(NumOfElements, 6);
5528 faces_tbl =
new STable3D(NumOfVertices);
5529 for (i = 0; i < NumOfElements; i++)
5531 v = elements[i]->GetVertices();
5532 switch (GetElementType(i))
5534 case Element::TETRAHEDRON:
5536 for (
int j = 0; j < 4; j++)
5538 const int *fv = tet_t::FaceVert[j];
5540 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
5544 case Element::WEDGE:
5546 for (
int j = 0; j < 2; j++)
5548 const int *fv = pri_t::FaceVert[j];
5550 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
5552 for (
int j = 2; j < 5; j++)
5554 const int *fv = pri_t::FaceVert[j];
5556 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
5560 case Element::HEXAHEDRON:
5564 for (
int j = 0; j < 6; j++)
5566 const int *fv = hex_t::FaceVert[j];
5568 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
5573 MFEM_ABORT(
"Unexpected type of Element.");
5576 el_to_face->Finalize();
5578 be_to_face.SetSize(NumOfBdrElements);
5579 for (i = 0; i < NumOfBdrElements; i++)
5581 v = boundary[i]->GetVertices();
5582 switch (GetBdrElementType(i))
5584 case Element::TRIANGLE:
5586 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
5589 case Element::QUADRILATERAL:
5591 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
5595 MFEM_ABORT(
"Unexpected type of boundary Element.");
5609 void Rotate3(
int &
a,
int &
b,
int &c)
5631 void Mesh::ReorientTetMesh()
5633 if (Dim != 3 || !(meshgen & 1))
5641 Table *old_elem_vert = NULL;
5645 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
5648 for (
int i = 0; i < NumOfElements; i++)
5650 if (GetElementType(i) == Element::TETRAHEDRON)
5652 int *v = elements[i]->GetVertices();
5654 Rotate3(v[0], v[1], v[2]);
5657 Rotate3(v[1], v[2], v[3]);
5666 for (
int i = 0; i < NumOfBdrElements; i++)
5668 if (GetBdrElementType(i) == Element::TRIANGLE)
5670 int *v = boundary[i]->GetVertices();
5672 Rotate3(v[0], v[1], v[2]);
5678 GetElementToFaceTable();
5682 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5687 DoNodeReorder(old_v_to_v, old_elem_vert);
5688 delete old_elem_vert;
5693 int *Mesh::CartesianPartitioning(
int nxyz[])
5699 for (
int vi = 0; vi < NumOfVertices; vi++)
5701 const double *
p = vertices[vi]();
5702 for (
int i = 0; i < spaceDim; i++)
5704 if (p[i] < pmin[i]) { pmin[i] = p[i]; }
5705 if (p[i] > pmax[i]) { pmax[i] = p[i]; }
5709 partitioning =
new int[NumOfElements];
5713 Vector pt(ppt, spaceDim);
5714 for (
int el = 0; el < NumOfElements; el++)
5716 GetElementTransformation(el)->Transform(
5719 for (
int i = spaceDim-1; i >= 0; i--)
5721 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
5722 if (idx < 0) { idx = 0; }
5723 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
5724 part = part * nxyz[i] + idx;
5726 partitioning[el] = part;
5729 return partitioning;
5732 int *Mesh::GeneratePartitioning(
int nparts,
int part_method)
5734 #ifdef MFEM_USE_METIS
5736 int print_messages = 1;
5739 int init_flag, fin_flag;
5740 MPI_Initialized(&init_flag);
5741 MPI_Finalized(&fin_flag);
5742 if (init_flag && !fin_flag)
5746 if (rank != 0) { print_messages = 0; }
5750 int i, *partitioning;
5752 ElementToElementTable();
5754 partitioning =
new int[NumOfElements];
5758 for (i = 0; i < NumOfElements; i++)
5760 partitioning[i] = 0;
5763 else if (NumOfElements <= nparts)
5765 for (i = 0; i < NumOfElements; i++)
5767 partitioning[i] = i;
5773 #ifndef MFEM_USE_METIS_5
5785 bool freedata =
false;
5787 idx_t *mpartitioning;
5790 if (
sizeof(
idx_t) ==
sizeof(int))
5792 I = (
idx_t*) el_to_el->GetI();
5793 J = (
idx_t*) el_to_el->GetJ();
5794 mpartitioning = (
idx_t*) partitioning;
5798 int *iI = el_to_el->GetI();
5799 int *iJ = el_to_el->GetJ();
5803 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
5804 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
5805 mpartitioning =
new idx_t[n];
5808 #ifndef MFEM_USE_METIS_5
5811 METIS_SetDefaultOptions(options);
5812 options[METIS_OPTION_CONTIG] = 1;
5816 if (part_method >= 0 && part_method <= 2)
5818 for (i = 0; i < n; i++)
5824 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
5830 if (part_method == 0 || part_method == 3)
5832 #ifndef MFEM_USE_METIS_5
5861 " error in METIS_PartGraphRecursive!");
5868 if (part_method == 1 || part_method == 4)
5870 #ifndef MFEM_USE_METIS_5
5899 " error in METIS_PartGraphKway!");
5906 if (part_method == 2 || part_method == 5)
5908 #ifndef MFEM_USE_METIS_5
5921 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
5938 " error in METIS_PartGraphKway!");
5946 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
5950 nparts = (int) mparts;
5951 if (mpartitioning != (
idx_t*)partitioning)
5953 for (
int k = 0; k<NumOfElements; k++)
5955 partitioning[k] = mpartitioning[k];
5962 delete[] mpartitioning;
5973 for (i = 0; i < nparts; i++)
5979 for (i = 0; i < NumOfElements; i++)
5981 psize[partitioning[i]].one++;
5985 for (i = 0; i < nparts; i++)
5987 if (psize[i].one == 0) { empty_parts++; }
5996 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
5997 << empty_parts <<
" empty parts!"
5998 <<
" Applying a simple fix ..." << endl;
6001 SortPairs<int,int>(psize, nparts);
6003 for (i = nparts-1; i > nparts-1-empty_parts; i--)
6008 for (
int j = 0; j < NumOfElements; j++)
6010 for (i = nparts-1; i > nparts-1-empty_parts; i--)
6012 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
6018 partitioning[j] = psize[nparts-1-i].two;
6025 for (i = 0; i < nparts; i++)
6030 for (i = 0; i < NumOfElements; i++)
6032 psize[partitioning[i]].one++;
6036 for (i = 0; i < nparts; i++)
6038 if (psize[i].one == 0) { empty_parts++; }
6044 return partitioning;
6048 mfem_error(
"Mesh::GeneratePartitioning(...): "
6049 "MFEM was compiled without Metis.");
6063 int num_elem, *i_elem_elem, *j_elem_elem;
6065 num_elem = elem_elem.
Size();
6066 i_elem_elem = elem_elem.
GetI();
6067 j_elem_elem = elem_elem.
GetJ();
6072 int stack_p, stack_top_p, elem;
6076 for (i = 0; i < num_elem; i++)
6078 if (partitioning[i] > num_part)
6080 num_part = partitioning[i];
6087 for (i = 0; i < num_part; i++)
6094 for (elem = 0; elem < num_elem; elem++)
6096 if (component[elem] >= 0)
6101 component[elem] = num_comp[partitioning[elem]]++;
6103 elem_stack[stack_top_p++] = elem;
6105 for ( ; stack_p < stack_top_p; stack_p++)
6107 i = elem_stack[stack_p];
6108 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
6111 if (partitioning[k] == partitioning[i])
6113 if (component[k] < 0)
6115 component[k] = component[i];
6116 elem_stack[stack_top_p++] = k;
6118 else if (component[k] != component[i])
6128 void Mesh::CheckPartitioning(
int *partitioning)
6130 int i, n_empty, n_mcomp;
6132 const Array<int> _partitioning(partitioning, GetNE());
6134 ElementToElementTable();
6138 n_empty = n_mcomp = 0;
6139 for (i = 0; i < num_comp.
Size(); i++)
6140 if (num_comp[i] == 0)
6144 else if (num_comp[i] > 1)
6151 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
6152 <<
"The following subdomains are empty :\n";
6153 for (i = 0; i < num_comp.
Size(); i++)
6154 if (num_comp[i] == 0)
6162 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
6163 <<
"The following subdomains are NOT connected :\n";
6164 for (i = 0; i < num_comp.
Size(); i++)
6165 if (num_comp[i] > 1)
6171 if (n_empty == 0 && n_mcomp == 0)
6172 mfem::out <<
"Mesh::CheckPartitioning(...) : "
6173 "All subdomains are connected." << endl;
6187 const double *a = A.
Data();
6188 const double *b = B.
Data();
6197 c(0) = a[0]*a[3]-a[1]*a[2];
6198 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
6199 c(2) = b[0]*b[3]-b[1]*b[2];
6220 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
6221 a[1] * (a[5] * a[6] - a[3] * a[8]) +
6222 a[2] * (a[3] * a[7] - a[4] * a[6]));
6224 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
6225 b[1] * (a[5] * a[6] - a[3] * a[8]) +
6226 b[2] * (a[3] * a[7] - a[4] * a[6]) +
6228 a[0] * (b[4] * a[8] - b[5] * a[7]) +
6229 a[1] * (b[5] * a[6] - b[3] * a[8]) +
6230 a[2] * (b[3] * a[7] - b[4] * a[6]) +
6232 a[0] * (a[4] * b[8] - a[5] * b[7]) +
6233 a[1] * (a[5] * b[6] - a[3] * b[8]) +
6234 a[2] * (a[3] * b[7] - a[4] * b[6]));
6236 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
6237 a[1] * (b[5] * b[6] - b[3] * b[8]) +
6238 a[2] * (b[3] * b[7] - b[4] * b[6]) +
6240 b[0] * (a[4] * b[8] - a[5] * b[7]) +
6241 b[1] * (a[5] * b[6] - a[3] * b[8]) +
6242 b[2] * (a[3] * b[7] - a[4] * b[6]) +
6244 b[0] * (b[4] * a[8] - b[5] * a[7]) +
6245 b[1] * (b[5] * a[6] - b[3] * a[8]) +
6246 b[2] * (b[3] * a[7] - b[4] * a[6]));
6248 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
6249 b[1] * (b[5] * b[6] - b[3] * b[8]) +
6250 b[2] * (b[3] * b[7] - b[4] * b[6]));
6296 double a = z(2), b = z(1), c = z(0);
6297 double D = b*b-4*a*c;
6304 x(0) = x(1) = -0.5 * b /
a;
6309 x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
6317 t = -0.5 * (b + sqrt(D));
6321 t = -0.5 * (b - sqrt(D));
6327 Swap<double>(x(0), x(1));
6335 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
6338 double Q = (a * a - 3 *
b) / 9;
6339 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
6340 double Q3 = Q * Q * Q;
6347 x(0) = x(1) = x(2) = - a / 3;
6351 double sqrtQ = sqrt(Q);
6355 x(0) = -2 * sqrtQ - a / 3;
6356 x(1) = x(2) = sqrtQ - a / 3;
6360 x(0) = x(1) = - sqrtQ - a / 3;
6361 x(2) = 2 * sqrtQ - a / 3;
6368 double theta = acos(R / sqrt(Q3));
6369 double A = -2 * sqrt(Q);
6371 x0 = A * cos(theta / 3) - a / 3;
6372 x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
6373 x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
6378 Swap<double>(x0, x1);
6382 Swap<double>(x1, x2);
6385 Swap<double>(x0, x1);
6398 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
6402 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
6404 x(0) = A + Q / A - a / 3;
6413 const double factor,
const int Dim)
6415 const double c0 = c(0);
6416 c(0) = c0 * (1.0 - pow(factor, -Dim));
6418 for (
int j = 0; j < nr; j++)
6430 c(0) = c0 * (1.0 - pow(factor, Dim));
6432 for (
int j = 0; j < nr; j++)
6446 void Mesh::CheckDisplacements(
const Vector &displacements,
double &tmax)
6448 int nvs = vertices.Size();
6449 DenseMatrix P, V, DS, PDS(spaceDim), VDS(spaceDim);
6450 Vector c(spaceDim+1), x(spaceDim);
6451 const double factor = 2.0;
6458 for (
int i = 0; i < NumOfElements; i++)
6465 for (
int j = 0; j < spaceDim; j++)
6466 for (
int k = 0; k < nv; k++)
6468 P(j, k) = vertices[v[k]](j);
6469 V(j, k) = displacements(v[k]+j*nvs);
6473 GetTransformationFEforElementType(el->
GetType());
6477 case Element::TRIANGLE:
6478 case Element::TETRAHEDRON:
6496 case Element::QUADRILATERAL:
6499 for (
int j = 0; j < nv; j++)
6523 void Mesh::MoveVertices(
const Vector &displacements)
6525 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
6526 for (
int j = 0; j < spaceDim; j++)
6528 vertices[i](j) += displacements(j*nv+i);
6532 void Mesh::GetVertices(
Vector &vert_coord)
const
6534 int nv = vertices.Size();
6535 vert_coord.
SetSize(nv*spaceDim);
6536 for (
int i = 0; i < nv; i++)
6537 for (
int j = 0; j < spaceDim; j++)
6539 vert_coord(j*nv+i) = vertices[i](j);
6543 void Mesh::SetVertices(
const Vector &vert_coord)
6545 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
6546 for (
int j = 0; j < spaceDim; j++)
6548 vertices[i](j) = vert_coord(j*nv+i);
6552 void Mesh::GetNode(
int i,
double *coord)
const
6557 for (
int j = 0; j < spaceDim; j++)
6559 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
6564 for (
int j = 0; j < spaceDim; j++)
6566 coord[j] = vertices[i](j);
6571 void Mesh::SetNode(
int i,
const double *coord)
6576 for (
int j = 0; j < spaceDim; j++)
6578 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
6583 for (
int j = 0; j < spaceDim; j++)
6585 vertices[i](j) = coord[j];
6591 void Mesh::MoveNodes(
const Vector &displacements)
6595 (*Nodes) += displacements;
6599 MoveVertices(displacements);
6603 void Mesh::GetNodes(
Vector &node_coord)
const
6607 node_coord = (*Nodes);
6611 GetVertices(node_coord);
6615 void Mesh::SetNodes(
const Vector &node_coord)
6619 (*Nodes) = node_coord;
6623 SetVertices(node_coord);
6629 if (own_nodes) {
delete Nodes; }
6632 own_nodes = (int)make_owner;
6643 mfem::Swap<GridFunction*>(Nodes, nodes);
6644 mfem::Swap<int>(own_nodes, own_nodes_);
6651 void Mesh::AverageVertices(
const int *indexes,
int n,
int result)
6655 for (k = 0; k < spaceDim; k++)
6657 vertices[result](k) = vertices[indexes[0]](k);
6660 for (j = 1; j < n; j++)
6661 for (k = 0; k < spaceDim; k++)
6663 vertices[result](k) += vertices[indexes[j]](k);
6666 for (k = 0; k < spaceDim; k++)
6668 vertices[result](k) *= (1.0 / n);
6672 void Mesh::UpdateNodes()
6676 Nodes->FESpace()->Update();
6681 void Mesh::UniformRefinement2D_base(
bool update_nodes)
6685 if (el_to_edge == NULL)
6687 el_to_edge =
new Table;
6688 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6691 int quad_counter = 0;
6692 for (
int i = 0; i < NumOfElements; i++)
6694 if (elements[i]->GetType() == Element::QUADRILATERAL)
6700 const int oedge = NumOfVertices;
6701 const int oelem = oedge + NumOfEdges;
6706 vertices.
SetSize(oelem + quad_counter);
6707 new_elements.
SetSize(4 * NumOfElements);
6710 for (
int i = 0, j = 0; i < NumOfElements; i++)
6713 const int attr = elements[i]->GetAttribute();
6714 int *v = elements[i]->GetVertices();
6715 const int *e = el_to_edge->GetRow(i);
6718 if (el_type == Element::TRIANGLE)
6720 for (
int ei = 0; ei < 3; ei++)
6722 for (
int k = 0; k < 2; k++)
6724 vv[k] = v[tri_t::Edges[ei][k]];
6726 AverageVertices(vv, 2, oedge+e[ei]);
6730 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
6732 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
6734 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
6736 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
6738 else if (el_type == Element::QUADRILATERAL)
6740 const int qe = quad_counter;
6742 AverageVertices(v, 4, oelem+qe);
6744 for (
int ei = 0; ei < 4; ei++)
6746 for (
int k = 0; k < 2; k++)
6748 vv[k] = v[quad_t::Edges[ei][k]];
6750 AverageVertices(vv, 2, oedge+e[ei]);
6754 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
6756 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
6758 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
6760 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
6764 MFEM_ABORT(
"unknown element type: " << el_type);
6766 FreeElement(elements[i]);
6771 new_boundary.
SetSize(2 * NumOfBdrElements);
6772 for (
int i = 0, j = 0; i < NumOfBdrElements; i++)
6774 const int attr = boundary[i]->GetAttribute();
6775 int *v = boundary[i]->GetVertices();
6777 new_boundary[j++] =
new Segment(v[0], oedge+be_to_edge[i], attr);
6778 new_boundary[j++] =
new Segment(oedge+be_to_edge[i], v[1], attr);
6780 FreeElement(boundary[i]);
6784 static const double A = 0.0, B = 0.5, C = 1.0;
6785 static double tri_children[2*3*4] =
6792 static double quad_children[2*4*4] =
6800 CoarseFineTr.point_matrices[Geometry::TRIANGLE]
6801 .UseExternalData(tri_children, 2, 3, 4);
6802 CoarseFineTr.point_matrices[Geometry::SQUARE]
6803 .UseExternalData(quad_children, 2, 4, 4);
6804 CoarseFineTr.embeddings.SetSize(elements.Size());
6806 for (
int i = 0; i < elements.Size(); i++)
6808 Embedding &emb = CoarseFineTr.embeddings[i];
6813 NumOfVertices = vertices.Size();
6814 NumOfElements = 4 * NumOfElements;
6815 NumOfBdrElements = 2 * NumOfBdrElements;
6818 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6821 last_operation = Mesh::REFINE;
6824 if (update_nodes) { UpdateNodes(); }
6827 if (!Nodes || update_nodes)
6829 CheckElementOrientation(
false);
6831 CheckBdrElementOrientation(
false);
6835 static inline double sqr(
const double &x)
6845 if (el_to_edge == NULL)
6847 el_to_edge =
new Table;
6848 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6851 if (el_to_face == NULL)
6853 GetElementToFaceTable();
6857 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
6860 int NumOfQuadFaces = 0;
6861 if (HasGeometry(Geometry::SQUARE))
6863 if (HasGeometry(Geometry::TRIANGLE))
6866 for (
int i = 0; i < faces.Size(); i++)
6868 if (faces[i]->GetType() == Element::QUADRILATERAL)
6870 f2qf[i] = NumOfQuadFaces;
6877 NumOfQuadFaces = faces.
Size();
6881 int hex_counter = 0;
6882 if (HasGeometry(Geometry::CUBE))
6884 for (
int i = 0; i < elements.Size(); i++)
6886 if (elements[i]->GetType() == Element::HEXAHEDRON)
6896 if (HasGeometry(Geometry::TETRAHEDRON))
6900 DSTable *v_to_v_ptr = v_to_v_p;
6903 v_to_v_ptr =
new DSTable(NumOfVertices);
6904 GetVertexToVertexTable(*v_to_v_ptr);
6909 for (
int i = 0; i < NumOfVertices; i++)
6916 std::sort(row_start, J_v2v.
end());
6919 for (
int i = 0; i < J_v2v.
Size(); i++)
6921 e2v[J_v2v[i].two] = i;
6930 for (
int i = 0; i < NumOfVertices; i++)
6934 it.SetIndex(e2v[it.Index()]);
6942 const int oedge = NumOfVertices;
6943 const int oface = oedge + NumOfEdges;
6944 const int oelem = oface + NumOfQuadFaces;
6949 vertices.
SetSize(oelem + hex_counter);
6950 new_elements.
SetSize(8 * NumOfElements);
6951 CoarseFineTr.embeddings.SetSize(new_elements.
Size());
6954 for (
int i = 0, j = 0; i < NumOfElements; i++)
6957 const int attr = elements[i]->GetAttribute();
6958 int *v = elements[i]->GetVertices();
6959 const int *e = el_to_edge->GetRow(i);
6964 const int ne = el_to_edge->RowSize(i);
6965 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
6971 case Element::TETRAHEDRON:
6973 for (
int ei = 0; ei < 6; ei++)
6975 for (
int k = 0; k < 2; k++)
6977 vv[k] = v[tet_t::Edges[ei][k]];
6979 AverageVertices(vv, 2, oedge+e[ei]);
6985 const int rt_algo = 1;
6996 double len_sqr, min_len;
6998 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
6999 sqr(J(1,0)-J(1,1)-J(1,2)) +
7000 sqr(J(2,0)-J(2,1)-J(2,2));
7003 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
7004 sqr(J(1,1)-J(1,0)-J(1,2)) +
7005 sqr(J(2,1)-J(2,0)-J(2,2));
7006 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
7008 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
7009 sqr(J(1,2)-J(1,0)-J(1,1)) +
7010 sqr(J(2,2)-J(2,0)-J(2,1));
7011 if (len_sqr < min_len) { rt = 2; }
7016 double Em_data[18], Js_data[9], Jp_data[9];
7019 double ar1, ar2,
kappa, kappa_min;
7021 for (
int s = 0; s < 3; s++)
7023 for (
int t = 0; t < 3; t++)
7025 Em(t,s) = 0.5*J(t,s);
7028 for (
int t = 0; t < 3; t++)
7030 Em(t,3) = 0.5*(J(t,0)+J(t,1));
7031 Em(t,4) = 0.5*(J(t,0)+J(t,2));
7032 Em(t,5) = 0.5*(J(t,1)+J(t,2));
7036 for (
int t = 0; t < 3; t++)
7038 Js(t,0) = Em(t,5)-Em(t,0);
7039 Js(t,1) = Em(t,1)-Em(t,0);
7040 Js(t,2) = Em(t,2)-Em(t,0);
7044 for (
int t = 0; t < 3; t++)
7046 Js(t,0) = Em(t,5)-Em(t,0);
7047 Js(t,1) = Em(t,2)-Em(t,0);
7048 Js(t,2) = Em(t,4)-Em(t,0);
7052 kappa_min = std::max(ar1, ar2);
7056 for (
int t = 0; t < 3; t++)
7058 Js(t,0) = Em(t,0)-Em(t,1);
7059 Js(t,1) = Em(t,4)-Em(t,1);
7060 Js(t,2) = Em(t,2)-Em(t,1);
7064 for (
int t = 0; t < 3; t++)
7066 Js(t,0) = Em(t,2)-Em(t,1);
7067 Js(t,1) = Em(t,4)-Em(t,1);
7068 Js(t,2) = Em(t,5)-Em(t,1);
7072 kappa = std::max(ar1, ar2);
7073 if (kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
7076 for (
int t = 0; t < 3; t++)
7078 Js(t,0) = Em(t,0)-Em(t,2);
7079 Js(t,1) = Em(t,1)-Em(t,2);
7080 Js(t,2) = Em(t,3)-Em(t,2);
7084 for (
int t = 0; t < 3; t++)
7086 Js(t,0) = Em(t,1)-Em(t,2);
7087 Js(t,1) = Em(t,5)-Em(t,2);
7088 Js(t,2) = Em(t,3)-Em(t,2);
7092 kappa = std::max(ar1, ar2);
7093 if (kappa < kappa_min) { rt = 2; }
7096 static const int mv_all[3][4][4] =
7098 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
7099 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
7100 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
7102 const int (&mv)[4][4] = mv_all[rt];
7104 #ifndef MFEM_USE_MEMALLOC
7106 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
7108 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
7110 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
7112 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
7114 for (
int k = 0; k < 4; k++)
7116 new_elements[j+4+k] =
7117 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
7118 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
7122 new_elements[j+0] = tet = TetMemory.Alloc();
7123 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
7125 new_elements[j+1] = tet = TetMemory.Alloc();
7126 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
7128 new_elements[j+2] = tet = TetMemory.Alloc();
7129 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
7131 new_elements[j+3] = tet = TetMemory.Alloc();
7132 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
7134 for (
int k = 0; k < 4; k++)
7136 new_elements[j+4+k] = tet = TetMemory.Alloc();
7137 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
7138 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
7141 for (
int k = 0; k < 4; k++)
7143 CoarseFineTr.embeddings[j+k].parent = i;
7144 CoarseFineTr.embeddings[j+k].matrix = k;
7146 for (
int k = 0; k < 4; k++)
7148 CoarseFineTr.embeddings[j+4+k].parent = i;
7149 CoarseFineTr.embeddings[j+4+k].matrix = 4*(rt+1)+k;
7156 case Element::WEDGE:
7158 const int *
f = el_to_face->GetRow(i);
7160 for (
int fi = 2; fi < 5; fi++)
7162 for (
int k = 0; k < 4; k++)
7164 vv[k] = v[pri_t::FaceVert[fi][k]];
7166 AverageVertices(vv, 4, oface + f2qf[f[fi]]);
7169 for (
int ei = 0; ei < 9; ei++)
7171 for (
int k = 0; k < 2; k++)
7173 vv[k] = v[pri_t::Edges[ei][k]];
7175 AverageVertices(vv, 2, oedge+e[ei]);
7178 const int qf2 = f2qf[f[2]];
7179 const int qf3 = f2qf[f[3]];
7180 const int qf4 = f2qf[f[4]];
7183 new Wedge(v[0], oedge+e[0], oedge+e[2],
7184 oedge+e[6], oface+qf2, oface+qf4, attr);
7187 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
7188 oface+qf3, oface+qf4, oface+qf2, attr);
7191 new Wedge(oedge+e[0], v[1], oedge+e[1],
7192 oface+qf2, oedge+e[7], oface+qf3, attr);
7195 new Wedge(oedge+e[2], oedge+e[1], v[2],
7196 oface+qf4, oface+qf3, oedge+e[8], attr);
7199 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
7200 v[3], oedge+e[3], oedge+e[5], attr);
7203 new Wedge(oface+qf3, oface+qf4, oface+qf2,
7204 oedge+e[4], oedge+e[5], oedge+e[3], attr);
7207 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
7208 oedge+e[3], v[4], oedge+e[4], attr);
7211 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
7212 oedge+e[5], oedge+e[4], v[5], attr);
7216 case Element::HEXAHEDRON:
7218 const int *
f = el_to_face->GetRow(i);
7219 const int he = hex_counter;
7224 if (f2qf.
Size() == 0)
7230 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[f[k]]; }
7234 AverageVertices(v, 8, oelem+he);
7236 for (
int fi = 0; fi < 6; fi++)
7238 for (
int k = 0; k < 4; k++)
7240 vv[k] = v[hex_t::FaceVert[fi][k]];
7242 AverageVertices(vv, 4, oface + qf[fi]);
7245 for (
int ei = 0; ei < 12; ei++)
7247 for (
int k = 0; k < 2; k++)
7249 vv[k] = v[hex_t::Edges[ei][k]];
7251 AverageVertices(vv, 2, oedge+e[ei]);
7255 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
7256 oedge+e[3], oedge+e[8], oface+qf[1],
7257 oelem+he, oface+qf[4], attr);
7260 oface+qf[0], oface+qf[1], oedge+e[9],
7261 oface+qf[2], oelem+he, attr);
7263 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
7264 oedge+e[2], oelem+he, oface+qf[2],
7265 oedge+e[10], oface+qf[3], attr);
7267 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
7268 v[3], oface+qf[4], oelem+he,
7269 oface+qf[3], oedge+e[11], attr);
7271 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
7272 oface+qf[4], v[4], oedge+e[4],
7273 oface+qf[5], oedge+e[7], attr);
7275 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
7276 oelem+he, oedge+e[4], v[5],
7277 oedge+e[5], oface+qf[5], attr);
7279 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
7280 oface+qf[3], oface+qf[5], oedge+e[5],
7281 v[6], oedge+e[6], attr);
7283 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
7284 oedge+e[11], oedge+e[7], oface+qf[5],
7285 oedge+e[6], v[7], attr);
7290 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
7293 FreeElement(elements[i]);
7298 new_boundary.
SetSize(4 * NumOfBdrElements);
7299 for (
int i = 0, j = 0; i < NumOfBdrElements; i++)
7302 const int attr = boundary[i]->GetAttribute();
7303 int *v = boundary[i]->GetVertices();
7304 const int *e = bel_to_edge->GetRow(i);
7309 const int ne = bel_to_edge->RowSize(i);
7310 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
7314 if (bdr_el_type == Element::TRIANGLE)
7317 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
7319 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
7321 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
7323 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
7325 else if (bdr_el_type == Element::QUADRILATERAL)
7328 (f2qf.
Size() == 0) ? be_to_face[i] : f2qf[be_to_face[i]];
7331 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
7333 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
7335 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
7337 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
7341 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
7343 FreeElement(boundary[i]);
7347 static const double A = 0.0, B = 0.5, C = 1.0;
7348 static double tet_children[3*4*16] =
7350 A,A,A, B,A,A, A,B,A, A,A,B,
7351 B,A,A, C,A,A, B,B,A, B,A,B,
7352 A,B,A, B,B,A, A,C,A, A,B,B,
7353 A,A,B, B,A,B, A,B,B, A,A,C,
7358 B,A,A, A,B,B, A,B,A, A,A,B,
7359 B,A,A, A,B,B, A,A,B, B,A,B,
7360 B,A,A, A,B,B, B,A,B, B,B,A,
7361 B,A,A, A,B,B, B,B,A, A,B,A,
7363 A,B,A, B,A,A, B,A,B, A,A,B,
7364 A,B,A, A,A,B, B,A,B, A,B,B,
7365 A,B,A, A,B,B, B,A,B, B,B,A,
7366 A,B,A, B,B,A, B,A,B, B,A,A,
7368 A,A,B, B,A,A, A,B,A, B,B,A,
7369 A,A,B, A,B,A, A,B,B, B,B,A,
7370 A,A,B, A,B,B, B,A,B, B,B,A,
7371 A,A,B, B,A,B, B,A,A, B,B,A
7373 static double pri_children[3*6*8] =
7375 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
7376 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
7377 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
7378 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
7379 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
7380 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
7381 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
7382 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
7384 static double hex_children[3*8*8] =
7386 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
7387 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
7388 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
7389 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
7390 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
7391 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
7392 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
7393 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
7396 CoarseFineTr.point_matrices[Geometry::TETRAHEDRON]
7397 .UseExternalData(tet_children, 3, 4, 16);
7398 CoarseFineTr.point_matrices[Geometry::PRISM]
7399 .UseExternalData(pri_children, 3, 6, 8);
7400 CoarseFineTr.point_matrices[Geometry::CUBE]
7401 .UseExternalData(hex_children, 3, 8, 8);
7403 for (
int i = 0; i < elements.Size(); i++)
7406 if (elements[i]->GetType() == Element::TETRAHEDRON) {
continue; }
7408 Embedding &emb = CoarseFineTr.embeddings[i];
7413 NumOfVertices = vertices.Size();
7414 NumOfElements = 8 * NumOfElements;
7415 NumOfBdrElements = 4 * NumOfBdrElements;
7417 GetElementToFaceTable();
7421 CheckBdrElementOrientation(
false);
7424 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7426 last_operation = Mesh::REFINE;
7429 if (update_nodes) { UpdateNodes(); }
7432 void Mesh::LocalRefinement(
const Array<int> &marked_el,
int type)
7434 int i, j, ind, nedges;
7441 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
7444 InitRefinementTransforms();
7448 int cne = NumOfElements, cnv = NumOfVertices;
7449 NumOfVertices += marked_el.
Size();
7450 NumOfElements += marked_el.
Size();
7451 vertices.SetSize(NumOfVertices);
7452 elements.SetSize(NumOfElements);
7453 CoarseFineTr.embeddings.SetSize(NumOfElements);
7455 for (j = 0; j < marked_el.
Size(); j++)
7460 int new_v = cnv + j, new_e = cne + j;
7461 AverageVertices(vert, 2, new_v);
7462 elements[new_e] =
new Segment(new_v, vert[1], attr);
7465 CoarseFineTr.embeddings[i] =
Embedding(i, 1);
7466 CoarseFineTr.embeddings[new_e] =
Embedding(i, 2);
7469 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
7470 CoarseFineTr.point_matrices[Geometry::SEGMENT].
7471 UseExternalData(seg_children, 1, 2, 3);
7479 DSTable v_to_v(NumOfVertices);
7480 GetVertexToVertexTable(v_to_v);
7484 int *edge1 =
new int[nedges];
7485 int *edge2 =
new int[nedges];
7486 int *middle =
new int[nedges];
7488 for (i = 0; i < nedges; i++)
7490 edge1[i] = edge2[i] = middle[i] = -1;
7493 for (i = 0; i < NumOfElements; i++)
7495 elements[i]->GetVertices(v);
7496 for (j = 1; j < v.
Size(); j++)
7498 ind = v_to_v(v[j-1], v[j]);
7499 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
7501 ind = v_to_v(v[0], v[v.
Size()-1]);
7502 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
7506 for (i = 0; i < marked_el.
Size(); i++)
7508 RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
7512 int need_refinement;
7515 need_refinement = 0;
7516 for (i = 0; i < nedges; i++)
7518 if (middle[i] != -1 && edge1[i] != -1)
7520 need_refinement = 1;
7521 GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
7525 while (need_refinement == 1);
7528 int v1[2], v2[2], bisect, temp;
7529 temp = NumOfBdrElements;
7530 for (i = 0; i < temp; i++)
7532 boundary[i]->GetVertices(v);
7533 bisect = v_to_v(v[0], v[1]);
7534 if (middle[bisect] != -1)
7536 if (boundary[i]->GetType() == Element::SEGMENT)
7538 v1[0] = v[0]; v1[1] = middle[bisect];
7539 v2[0] = middle[bisect]; v2[1] = v[1];
7541 boundary[i]->SetVertices(v1);
7542 boundary.
Append(
new Segment(v2, boundary[i]->GetAttribute()));
7545 mfem_error(
"Only bisection of segment is implemented"
7549 NumOfBdrElements = boundary.Size();
7556 if (el_to_edge != NULL)
7558 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7569 MFEM_VERIFY(GetNE() == 0 ||
7570 ((
Tetrahedron*)elements[0])->GetRefinementFlag() != 0,
7571 "tetrahedral mesh is not marked for refinement:"
7572 " call Finalize(true)");
7579 for (i = 0; i < marked_el.
Size(); i++)
7581 Bisection(marked_el[i], v_to_v);
7585 for (i = 0; i < marked_el.
Size(); i++)
7587 Bisection(marked_el[i], v_to_v);
7589 Bisection(NumOfElements - 1, v_to_v);
7590 Bisection(marked_el[i], v_to_v);
7594 for (i = 0; i < marked_el.
Size(); i++)
7596 Bisection(marked_el[i], v_to_v);
7598 ii = NumOfElements - 1;
7599 Bisection(ii, v_to_v);
7600 Bisection(NumOfElements - 1, v_to_v);
7601 Bisection(ii, v_to_v);
7603 Bisection(marked_el[i], v_to_v);
7604 Bisection(NumOfElements-1, v_to_v);
7605 Bisection(marked_el[i], v_to_v);
7611 int need_refinement;
7616 need_refinement = 0;
7619 for (i = 0; i < NumOfElements; i++)
7624 if (elements[i]->NeedRefinement(v_to_v))
7626 need_refinement = 1;
7627 Bisection(i, v_to_v);
7631 while (need_refinement == 1);
7638 need_refinement = 0;
7639 for (i = 0; i < NumOfBdrElements; i++)
7640 if (boundary[i]->NeedRefinement(v_to_v))
7642 need_refinement = 1;
7643 BdrBisection(i, v_to_v);
7646 while (need_refinement == 1);
7648 NumOfVertices = vertices.Size();
7649 NumOfBdrElements = boundary.Size();
7652 if (el_to_edge != NULL)
7654 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7656 if (el_to_face != NULL)
7658 GetElementToFaceTable();
7664 last_operation = Mesh::REFINE;
7670 CheckElementOrientation(
false);
7677 MFEM_VERIFY(!NURBSext,
"Nonconforming refinement of NURBS meshes is "
7678 "not supported. Project the NURBS to Nodes first.");
7685 ncmesh =
new NCMesh(
this);
7688 if (!refinements.
Size())
7690 last_operation = Mesh::NONE;
7695 ncmesh->MarkCoarseLevel();
7696 ncmesh->Refine(refinements);
7700 ncmesh->LimitNCLevel(nc_limit);
7705 ncmesh->OnMeshUpdated(mesh2);
7709 Swap(*mesh2,
false);
7712 GenerateNCFaceInfo();
7714 last_operation = Mesh::REFINE;
7719 Nodes->FESpace()->Update();
7725 const int *fine,
int nfine,
int op)
7728 for (
int i = 0; i < nfine; i++)
7730 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
7732 double err_fine = elem_error[fine[i]];
7735 case 0: error = std::min(error, err_fine);
break;
7736 case 1: error += err_fine;
break;
7737 case 2: error = std::max(error, err_fine);
break;
7744 double threshold,
int nc_limit,
int op)
7746 MFEM_VERIFY(ncmesh,
"Only supported for non-conforming meshes.");
7747 MFEM_VERIFY(!NURBSext,
"Derefinement of NURBS meshes is not supported. "
7748 "Project the NURBS to Nodes first.");
7752 const Table &dt = ncmesh->GetDerefinementTable();
7757 ncmesh->CheckDerefinementNCLevel(dt, level_ok, nc_limit);
7761 for (
int i = 0; i < dt.
Size(); i++)
7763 if (nc_limit > 0 && !level_ok[i]) {
continue; }
7766 AggregateError(elem_error, dt.
GetRow(i), dt.
RowSize(i), op);
7768 if (error < threshold) { derefs.
Append(i); }
7771 if (!derefs.
Size()) {
return false; }
7773 ncmesh->Derefine(derefs);
7776 ncmesh->OnMeshUpdated(mesh2);
7778 Swap(*mesh2,
false);
7781 GenerateNCFaceInfo();
7783 last_operation = Mesh::DEREFINE;
7792 int nc_limit,
int op)
7796 if (Nonconforming())
7798 return NonconformingDerefinement(elem_error, threshold, nc_limit, op);
7802 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
7808 bool Mesh::DerefineByError(
const Vector &elem_error,
double threshold,
7809 int nc_limit,
int op)
7812 for (
int i = 0; i < tmp.Size(); i++)
7814 tmp[i] = elem_error(i);
7816 return DerefineByError(tmp, threshold, nc_limit, op);
7820 void Mesh::InitFromNCMesh(
const NCMesh &ncmesh)
7829 NumOfVertices = vertices.Size();
7830 NumOfElements = elements.Size();
7831 NumOfBdrElements = boundary.Size();
7835 NumOfEdges = NumOfFaces = 0;
7836 nbInteriorFaces = nbBoundaryFaces = -1;
7840 el_to_edge =
new Table;
7841 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7845 GetElementToFaceTable();
7849 CheckBdrElementOrientation(
false);
7860 InitFromNCMesh(ncmesh);
7899 #ifdef MFEM_USE_MEMALLOC
7917 const int nv = Geometry::NumVerts[
geom];
7919 for (
int i = 0; i < elem_array.
Size(); i++)
7921 if (elem_array[i]->GetGeometryType() ==
geom)
7926 elem_vtx.
SetSize(nv*num_elems);
7930 for (
int i = 0; i < elem_array.
Size(); i++)
7936 elem_vtx.
Append(loc_vtx);
7944 for (
int i = 0; i < nelem; i++) { list[i] = i; }
7948 void Mesh::UniformRefinement(
int ref_algo)
7954 NURBSUniformRefinement();
7958 GeneralRefinement(AllElements(list, GetNE()));
7960 else if (ref_algo == 1 && meshgen == 1 && Dim == 3)
7963 LocalRefinement(AllElements(list, GetNE()));
7969 case 1: LocalRefinement(AllElements(list, GetNE()));
break;
7970 case 2: UniformRefinement2D();
break;
7971 case 3: UniformRefinement3D();
break;
7972 default: MFEM_ABORT(
"internal error");
7978 int nonconforming,
int nc_limit)
7984 else if (Dim == 1 || (Dim == 3 && (meshgen & 1)))
7988 else if (nonconforming < 0)
7991 if ((meshgen & 2) || (meshgen & 4))
8004 NonconformingRefinement(refinements, nc_limit);
8009 for (
int i = 0; i < refinements.
Size(); i++)
8011 el_to_refine[i] = refinements[i].index;
8015 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
8016 if (rt == 1 || rt == 2 || rt == 4)
8020 else if (rt == 3 || rt == 5 || rt == 6)
8030 LocalRefinement(el_to_refine, type);
8034 void Mesh::GeneralRefinement(
const Array<int> &el_to_refine,
int nonconforming,
8038 for (
int i = 0; i < el_to_refine.
Size(); i++)
8040 refinements[i] =
Refinement(el_to_refine[i]);
8042 GeneralRefinement(refinements, nonconforming, nc_limit);
8045 void Mesh::EnsureNCMesh(
bool simplices_nonconforming)
8047 MFEM_VERIFY(!NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
8048 "Project the NURBS to Nodes first.");
8052 if ((meshgen & 0x2) ||
8054 (simplices_nonconforming && (meshgen & 0x1)) )
8056 ncmesh =
new NCMesh(
this);
8057 ncmesh->OnMeshUpdated(
this);
8058 GenerateNCFaceInfo();
8063 void Mesh::RandomRefinement(
double prob,
bool aniso,
int nonconforming,
8067 for (
int i = 0; i < GetNE(); i++)
8069 if ((
double) rand() / RAND_MAX <
prob)
8074 type = (Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
8079 GeneralRefinement(refs, nonconforming, nc_limit);
8082 void Mesh::RefineAtVertex(
const Vertex& vert,
double eps,
int nonconforming)
8086 for (
int i = 0; i < GetNE(); i++)
8088 GetElementVertices(i, v);
8089 bool refine =
false;
8090 for (
int j = 0; j < v.
Size(); j++)
8093 for (
int l = 0; l < spaceDim; l++)
8095 double d = vert(l) - vertices[v[j]](l);
8098 if (dist <= eps*eps) { refine =
true;
break; }
8105 GeneralRefinement(refs, nonconforming);
8109 int nonconforming,
int nc_limit)
8111 MFEM_VERIFY(elem_error.
Size() == GetNE(),
"");
8113 for (
int i = 0; i < GetNE(); i++)
8115 if (elem_error[i] > threshold)
8120 if (ReduceInt(refs.Size()))
8122 GeneralRefinement(refs, nonconforming, nc_limit);
8128 bool Mesh::RefineByError(
const Vector &elem_error,
double threshold,
8129 int nonconforming,
int nc_limit)
8133 return RefineByError(tmp, threshold, nonconforming, nc_limit);
8137 void Mesh::Bisection(
int i,
const DSTable &v_to_v,
8138 int *edge1,
int *edge2,
int *middle)
8141 int v[2][4], v_new, bisect, t;
8146 if (t == Element::TRIANGLE)
8153 bisect = v_to_v(vert[0], vert[1]);
8154 MFEM_ASSERT(bisect >= 0,
"");
8156 if (middle[bisect] == -1)
8158 v_new = NumOfVertices++;
8159 for (
int d = 0; d < spaceDim; d++)
8161 V(d) = 0.5 * (vertices[vert[0]](d) + vertices[vert[1]](d));
8167 if (edge1[bisect] == i)
8169 edge1[bisect] = edge2[bisect];
8172 middle[bisect] = v_new;
8176 v_new = middle[bisect];
8184 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
8185 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
8190 elements.Append(tri_new);
8199 int coarse = FindCoarseElement(i);
8200 CoarseFineTr.embeddings[i].parent = coarse;
8201 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8206 bisect = v_to_v(v[1][0], v[1][1]);
8207 MFEM_ASSERT(bisect >= 0,
"");
8209 if (edge1[bisect] == i)
8211 edge1[bisect] = NumOfElements;
8213 else if (edge2[bisect] == i)
8215 edge2[bisect] = NumOfElements;
8222 MFEM_ABORT(
"Bisection for now works only for triangles.");
8229 int v[2][4], v_new, bisect, t;
8234 if (t == Element::TETRAHEDRON)
8236 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
8240 "TETRAHEDRON element is not marked for refinement.");
8245 bisect = v_to_v.
FindId(vert[0], vert[1]);
8248 v_new = NumOfVertices + v_to_v.
GetId(vert[0],vert[1]);
8249 for (j = 0; j < 3; j++)
8251 V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
8257 v_new = NumOfVertices + bisect;
8266 new_redges[0][0] = 2;
8267 new_redges[0][1] = 1;
8268 new_redges[1][0] = 2;
8269 new_redges[1][1] = 1;
8270 int tr1 = -1, tr2 = -1;
8271 switch (old_redges[0])
8274 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
8275 if (type == Tetrahedron::TYPE_PF) { new_redges[0][1] = 4; }
8279 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
8283 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
8286 switch (old_redges[1])
8289 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
8290 if (type == Tetrahedron::TYPE_PF) { new_redges[1][0] = 3; }
8294 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
8298 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
8305 #ifdef MFEM_USE_MEMALLOC
8313 elements.Append(tet2);
8319 int coarse = FindCoarseElement(i);
8320 CoarseFineTr.embeddings[i].parent = coarse;
8321 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8326 case Tetrahedron::TYPE_PU:
8327 new_type = Tetrahedron::TYPE_PF;
break;
8328 case Tetrahedron::TYPE_PF:
8329 new_type = Tetrahedron::TYPE_A;
break;
8331 new_type = Tetrahedron::TYPE_PU;
8341 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
8348 int v[2][3], v_new, bisect, t;
8349 Element *bdr_el = boundary[i];
8352 if (t == Element::TRIANGLE)
8359 bisect = v_to_v.
FindId(vert[0], vert[1]);
8360 MFEM_ASSERT(bisect >= 0,
"");
8361 v_new = NumOfVertices + bisect;
8362 MFEM_ASSERT(v_new != -1,
"");
8366 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
8367 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
8377 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
8382 void Mesh::UniformRefinement(
int i,
const DSTable &v_to_v,
8383 int *edge1,
int *edge2,
int *middle)
8386 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
8389 if (elements[i]->GetType() == Element::TRIANGLE)
8395 bisect[0] = v_to_v(v[0],v[1]);
8396 bisect[1] = v_to_v(v[1],v[2]);
8397 bisect[2] = v_to_v(v[0],v[2]);
8398 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
8400 for (j = 0; j < 3; j++)
8402 if (middle[bisect[j]] == -1)
8404 v_new[j] = NumOfVertices++;
8405 for (
int d = 0; d < spaceDim; d++)
8407 V(d) = (vertices[v[j]](d) + vertices[v[(j+1)%3]](d))/2.;
8413 if (edge1[bisect[j]] == i)
8415 edge1[bisect[j]] = edge2[bisect[j]];
8418 middle[bisect[j]] = v_new[j];
8422 v_new[j] = middle[bisect[j]];
8425 edge1[bisect[j]] = -1;
8431 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
8432 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
8433 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
8434 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
8440 elements.Append(tri1);
8441 elements.Append(tri2);
8442 elements.Append(tri3);
8449 tri2->ResetTransform(code);
8450 tri3->ResetTransform(code);
8454 tri2->PushTransform(1);
8455 tri3->PushTransform(2);
8458 int coarse = FindCoarseElement(i);
8459 CoarseFineTr.embeddings[i] =
Embedding(coarse);
8460 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8461 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8462 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8468 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
8472 void Mesh::InitRefinementTransforms()
8475 CoarseFineTr.Clear();
8476 CoarseFineTr.embeddings.SetSize(NumOfElements);
8477 for (
int i = 0; i < NumOfElements; i++)
8479 elements[i]->ResetTransform(0);
8480 CoarseFineTr.embeddings[i] =
Embedding(i);
8484 int Mesh::FindCoarseElement(
int i)
8487 while ((coarse = CoarseFineTr.embeddings[i].parent) != i)
8496 MFEM_VERIFY(GetLastOperation() == Mesh::REFINE,
"");
8500 return ncmesh->GetRefinementTransforms();
8504 for (
int i = 0; i < elem_geoms.
Size(); i++)
8507 if (CoarseFineTr.point_matrices[geom].SizeK()) {
continue; }
8509 if (geom == Geometry::TRIANGLE ||
8510 geom == Geometry::TETRAHEDRON)
8512 std::map<unsigned, int> mat_no;
8516 for (
int i = 0; i < elements.Size(); i++)
8519 unsigned code = elements[i]->GetTransform();
8522 int &matrix = mat_no[code];
8523 if (!matrix) { matrix = mat_no.size(); }
8526 CoarseFineTr.embeddings[i].matrix =
index;
8530 pmats.
SetSize(Dim, Dim+1, mat_no.size());
8533 std::map<unsigned, int>::iterator it;
8534 for (it = mat_no.begin(); it != mat_no.end(); ++it)
8536 if (geom == Geometry::TRIANGLE)
8538 Triangle::GetPointMatrix(it->first, pmats(it->second-1));
8542 Tetrahedron::GetPointMatrix(it->first, pmats(it->second-1));
8548 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
8549 " geom = " << geom);
8554 return CoarseFineTr;
8557 void Mesh::PrintXG(std::ostream &
out)
const
8559 MFEM_ASSERT(Dim==spaceDim,
"2D Manifold meshes not supported");
8568 out <<
"areamesh2\n\n";
8572 out <<
"curved_areamesh2\n\n";
8576 out << NumOfBdrElements <<
'\n';
8577 for (i = 0; i < NumOfBdrElements; i++)
8579 boundary[i]->GetVertices(v);
8581 out << boundary[i]->GetAttribute();
8582 for (j = 0; j < v.
Size(); j++)
8584 out <<
' ' << v[j] + 1;
8590 out << NumOfElements <<
'\n';
8591 for (i = 0; i < NumOfElements; i++)
8593 elements[i]->GetVertices(v);
8595 out << elements[i]->GetAttribute() <<
' ' << v.
Size();
8596 for (j = 0; j < v.
Size(); j++)
8598 out <<
' ' << v[j] + 1;
8606 out << NumOfVertices <<
'\n';
8607 for (i = 0; i < NumOfVertices; i++)
8609 out << vertices[i](0);
8610 for (j = 1; j < Dim; j++)
8612 out <<
' ' << vertices[i](j);
8619 out << NumOfVertices <<
'\n';
8627 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
8635 out <<
"NETGEN_Neutral_Format\n";
8637 out << NumOfVertices <<
'\n';
8638 for (i = 0; i < NumOfVertices; i++)
8640 for (j = 0; j < Dim; j++)
8642 out <<
' ' << vertices[i](j);
8648 out << NumOfElements <<
'\n';
8649 for (i = 0; i < NumOfElements; i++)
8651 nv = elements[i]->GetNVertices();
8652 ind = elements[i]->GetVertices();
8653 out << elements[i]->GetAttribute();
8654 for (j = 0; j < nv; j++)
8656 out <<
' ' << ind[j]+1;
8662 out << NumOfBdrElements <<
'\n';
8663 for (i = 0; i < NumOfBdrElements; i++)
8665 nv = boundary[i]->GetNVertices();
8666 ind = boundary[i]->GetVertices();
8667 out << boundary[i]->GetAttribute();
8668 for (j = 0; j < nv; j++)
8670 out <<
' ' << ind[j]+1;
8675 else if (meshgen == 2)
8681 <<
"1 " << NumOfVertices <<
" " << NumOfElements
8682 <<
" 0 0 0 0 0 0 0\n"
8683 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
8684 <<
"0 0 " << NumOfBdrElements <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
8685 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
8686 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
8688 for (i = 0; i < NumOfVertices; i++)
8689 out << i+1 <<
" 0.0 " << vertices[i](0) <<
' ' << vertices[i](1)
8690 <<
' ' << vertices[i](2) <<
" 0.0\n";
8692 for (i = 0; i < NumOfElements; i++)
8694 nv = elements[i]->GetNVertices();
8695 ind = elements[i]->GetVertices();
8696 out << i+1 <<
' ' << elements[i]->GetAttribute();
8697 for (j = 0; j < nv; j++)
8699 out <<
' ' << ind[j]+1;
8704 for (i = 0; i < NumOfBdrElements; i++)
8706 nv = boundary[i]->GetNVertices();
8707 ind = boundary[i]->GetVertices();
8708 out << boundary[i]->GetAttribute();
8709 for (j = 0; j < nv; j++)
8711 out <<
' ' << ind[j]+1;
8713 out <<
" 1.0 1.0 1.0 1.0\n";
8721 void Mesh::Printer(std::ostream &
out, std::string section_delimiter)
const
8728 NURBSext->Print(out);
8739 out << (ncmesh ?
"MFEM mesh v1.1\n" :
8740 section_delimiter.empty() ?
"MFEM mesh v1.0\n" :
8741 "MFEM mesh v1.2\n");
8745 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8750 "# TETRAHEDRON = 4\n"
8755 out <<
"\ndimension\n" << Dim
8756 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8757 for (i = 0; i < NumOfElements; i++)
8759 PrintElement(elements[i], out);
8762 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
8763 for (i = 0; i < NumOfBdrElements; i++)
8765 PrintElement(boundary[i], out);
8770 out <<
"\nvertex_parents\n";
8771 ncmesh->PrintVertexParents(out);
8773 out <<
"\ncoarse_elements\n";
8774 ncmesh->PrintCoarseElements(out);
8777 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8780 out << spaceDim <<
'\n';
8781 for (i = 0; i < NumOfVertices; i++)
8783 out << vertices[i](0);
8784 for (j = 1; j < spaceDim; j++)
8786 out <<
' ' << vertices[i](j);
8798 if (!ncmesh && !section_delimiter.empty())
8800 out << section_delimiter << endl;
8809 out <<
"MFEM NURBS mesh v1.0\n";
8813 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8819 out <<
"\ndimension\n" << Dim
8820 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8821 for (i = 0; i < NumOfElements; i++)
8823 PrintElement(elements[i], out);
8826 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
8827 for (i = 0; i < NumOfBdrElements; i++)
8829 PrintElement(boundary[i], out);
8832 out <<
"\nedges\n" << NumOfEdges <<
'\n';
8833 for (i = 0; i < NumOfEdges; i++)
8835 edge_vertex->GetRow(i, vert);
8841 out << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
8843 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8846 #ifdef MFEM_USE_ADIOS2
8853 void Mesh::PrintVTK(std::ostream &
out)
8856 "# vtk DataFile Version 3.0\n"
8857 "Generated by MFEM\n"
8859 "DATASET UNSTRUCTURED_GRID\n";
8863 out <<
"POINTS " << NumOfVertices <<
" double\n";
8864 for (
int i = 0; i < NumOfVertices; i++)
8866 out << vertices[i](0);
8868 for (j = 1; j < spaceDim; j++)
8870 out <<
' ' << vertices[i](j);
8882 out <<
"POINTS " << Nodes->FESpace()->GetNDofs() <<
" double\n";
8883 for (
int i = 0; i < Nodes->FESpace()->GetNDofs(); i++)
8887 Nodes->FESpace()->DofsToVDofs(vdofs);
8888 out << (*Nodes)(vdofs[0]);
8890 for (j = 1; j < spaceDim; j++)
8892 out <<
' ' << (*Nodes)(vdofs[j]);
8906 for (
int i = 0; i < NumOfElements; i++)
8908 size += elements[i]->GetNVertices() + 1;
8910 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
8911 for (
int i = 0; i < NumOfElements; i++)
8913 const int *v = elements[i]->GetVertices();
8914 const int nv = elements[i]->GetNVertices();
8916 for (
int j = 0; j < nv; j++)
8928 for (
int i = 0; i < NumOfElements; i++)
8930 Nodes->FESpace()->GetElementDofs(i, dofs);
8931 MFEM_ASSERT(Dim != 0 || dofs.
Size() == 1,
8932 "Point meshes should have a single dof per element");
8933 size += dofs.
Size() + 1;
8935 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
8936 const char *fec_name = Nodes->FESpace()->FEColl()->Name();
8938 if (!strcmp(fec_name,
"Linear") ||
8939 !strcmp(fec_name,
"H1_0D_P1") ||
8940 !strcmp(fec_name,
"H1_1D_P1") ||
8941 !strcmp(fec_name,
"H1_2D_P1") ||
8942 !strcmp(fec_name,
"H1_3D_P1"))
8946 else if (!strcmp(fec_name,
"Quadratic") ||
8947 !strcmp(fec_name,
"H1_1D_P2") ||
8948 !strcmp(fec_name,
"H1_2D_P2") ||
8949 !strcmp(fec_name,
"H1_3D_P2"))
8955 mfem::err <<
"Mesh::PrintVTK : can not save '"
8956 << fec_name <<
"' elements!" << endl;
8959 for (
int i = 0; i < NumOfElements; i++)
8961 Nodes->FESpace()->GetElementDofs(i, dofs);
8965 for (
int j = 0; j < dofs.
Size(); j++)
8967 out <<
' ' << dofs[j];
8970 else if (order == 2)
8972 const int *vtk_mfem;
8973 switch (elements[i]->GetGeometryType())
8975 case Geometry::SEGMENT:
8976 case Geometry::TRIANGLE:
8977 case Geometry::SQUARE:
8978 vtk_mfem = vtk_quadratic_hex;
break;
8979 case Geometry::TETRAHEDRON:
8980 vtk_mfem = vtk_quadratic_tet;
break;
8981 case Geometry::PRISM:
8982 vtk_mfem = vtk_quadratic_wedge;
break;
8983 case Geometry::CUBE:
8985 vtk_mfem = vtk_quadratic_hex;
break;
8987 for (
int j = 0; j < dofs.
Size(); j++)
8989 out <<
' ' << dofs[vtk_mfem[j]];
8996 out <<
"CELL_TYPES " << NumOfElements <<
'\n';
8997 for (
int i = 0; i < NumOfElements; i++)
8999 int vtk_cell_type = 5;
9005 case Geometry::POINT: vtk_cell_type = 1;
break;
9006 case Geometry::SEGMENT: vtk_cell_type = 3;
break;
9007 case Geometry::TRIANGLE: vtk_cell_type = 5;
break;
9008 case Geometry::SQUARE: vtk_cell_type = 9;
break;
9009 case Geometry::TETRAHEDRON: vtk_cell_type = 10;
break;
9010 case Geometry::CUBE: vtk_cell_type = 12;
break;
9011 case Geometry::PRISM: vtk_cell_type = 13;
break;
9015 else if (order == 2)
9019 case Geometry::SEGMENT: vtk_cell_type = 21;
break;
9020 case Geometry::TRIANGLE: vtk_cell_type = 22;
break;
9021 case Geometry::SQUARE: vtk_cell_type = 28;
break;
9022 case Geometry::TETRAHEDRON: vtk_cell_type = 24;
break;
9023 case Geometry::CUBE: vtk_cell_type = 29;
break;
9024 case Geometry::PRISM: vtk_cell_type = 32;
break;
9029 out << vtk_cell_type <<
'\n';
9033 out <<
"CELL_DATA " << NumOfElements <<
'\n'
9034 <<
"SCALARS material int\n"
9035 <<
"LOOKUP_TABLE default\n";
9036 for (
int i = 0; i < NumOfElements; i++)
9038 out << elements[i]->GetAttribute() <<
'\n';
9043 void Mesh::PrintVTU(std::string fname,
9045 bool high_order_output,
9046 int compression_level,
9049 int ref = (high_order_output && Nodes) ? Nodes->FESpace()->GetOrder(0) : 1;
9050 fname = fname +
".vtu";
9052 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
9053 if (compression_level != 0)
9055 out <<
" compressor=\"vtkZLibDataCompressor\"";
9058 out <<
"<UnstructuredGrid>\n";
9059 PrintVTU(out, ref, format, high_order_output, compression_level, bdr);
9060 out <<
"</Piece>\n";
9061 out <<
"</UnstructuredGrid>\n";
9062 out <<
"</VTKFile>" << std::endl;
9067 void Mesh::PrintBdrVTU(std::string fname,
9069 bool high_order_output,
9070 int compression_level)
9072 PrintVTU(fname, format, high_order_output, compression_level,
true);
9075 template <
typename T>
9079 if (format == VTKFormat::ASCII) { out << val << suffix; }
9086 const uint8_t &val,
const char *suffix,
9089 if (format == VTKFormat::ASCII) { out << static_cast<int>(val) << suffix; }
9095 const double &val,
const char *suffix,
9098 if (format == VTKFormat::BINARY32)
9100 bin_io::AppendBytes<float>(buf, float(val));
9102 else if (format == VTKFormat::BINARY)
9108 out << val << suffix;
9114 const float &val,
const char *suffix,
9117 if (format == VTKFormat::BINARY) { bin_io::AppendBytes<double>(buf, val); }
9119 else {
out << val << suffix; }
9123 int compression_level)
9131 bool high_order_output,
int compression_level,
9137 const char *fmt_str = (format == VTKFormat::ASCII) ?
"ascii" :
"binary";
9138 const char *type_str = (format != VTKFormat::BINARY32) ?
"Float64" :
"Float32";
9139 std::vector<char> buf;
9141 auto get_geom = [&](
int i)
9143 if (bdr_elements) {
return GetBdrElementBaseGeometry(i); }
9144 else {
return GetElementBaseGeometry(i); }
9147 int ne = bdr_elements ? GetNBE() : GetNE();
9149 int np = 0, nc_ref = 0, size = 0;
9150 for (
int i = 0; i < ne; i++)
9160 out <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
9161 << (high_order_output ? ne : nc_ref) <<
"\">\n";
9164 out <<
"<Points>\n";
9165 out <<
"<DataArray type=\"" << type_str
9166 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
9167 for (
int i = 0; i < ne; i++)
9173 GetBdrElementTransformation(i)->Transform(RefG->
RefPts, pmat);
9177 GetElementTransformation(i)->Transform(RefG->
RefPts, pmat);
9180 for (
int j = 0; j < pmat.
Width(); j++)
9199 if (format == VTKFormat::ASCII) { out <<
'\n'; }
9202 if (format != VTKFormat::ASCII)
9206 out <<
"</DataArray>" << std::endl;
9207 out <<
"</Points>" << std::endl;
9209 out <<
"<Cells>" << std::endl;
9210 out <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
9211 << fmt_str <<
"\">" << std::endl;
9213 std::vector<int> offset;
9216 if (high_order_output)
9219 for (
int iel = 0; iel < ne; iel++)
9223 int nnodes = local_connectivity.
Size();
9224 for (
int i=0; i<nnodes; ++i)
9228 if (format == VTKFormat::ASCII) { out <<
'\n'; }
9230 offset.push_back(np);
9236 for (
int i = 0; i < ne; i++)
9242 for (
int j = 0; j < RG.
Size(); )
9245 offset.push_back(coff);
9246 for (
int k = 0; k < nv; k++, j++)
9250 if (format == VTKFormat::ASCII) { out <<
'\n'; }
9255 if (format != VTKFormat::ASCII)
9259 out <<
"</DataArray>" << std::endl;
9261 out <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
9262 << fmt_str <<
"\">" << std::endl;
9264 for (
size_t ii=0; ii<offset.size(); ii++)
9268 if (format != VTKFormat::ASCII)
9272 out <<
"</DataArray>" << std::endl;
9273 out <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
9274 << fmt_str <<
"\">" << std::endl;
9276 for (
int i = 0; i < ne; i++)
9279 uint8_t vtk_cell_type = 5;
9284 case Geometry::POINT:
9287 case Geometry::SEGMENT:
9288 vtk_cell_type = high_order_output ? 68 : 3;
9290 case Geometry::TRIANGLE:
9291 vtk_cell_type = high_order_output ? 69 : 5;
9293 case Geometry::SQUARE:
9294 vtk_cell_type = high_order_output ? 70 : 9;
9296 case Geometry::TETRAHEDRON:
9297 vtk_cell_type = high_order_output ? 71 : 10;
9299 case Geometry::CUBE:
9300 vtk_cell_type = high_order_output ? 72 : 12;
9302 case Geometry::PRISM:
9303 vtk_cell_type = high_order_output ? 73 : 13;
9306 MFEM_ABORT(
"Unrecognized VTK element type \"" << geom <<
"\"");
9310 if (high_order_output)
9319 for (
int j = 0; j < RG.
Size(); j += nv)
9325 if (format != VTKFormat::ASCII)
9329 out <<
"</DataArray>" << std::endl;
9330 out <<
"</Cells>" << std::endl;
9332 out <<
"<CellData Scalars=\"attribute\">" << std::endl;
9333 out <<
"<DataArray type=\"Int32\" Name=\"attribute\" format=\""
9334 << fmt_str <<
"\">" << std::endl;
9335 for (
int i = 0; i < ne; i++)
9337 int attr = bdr_elements ? GetBdrAttribute(i) : GetAttribute(i);
9338 if (high_order_output)
9353 if (format != VTKFormat::ASCII)
9357 out <<
"</DataArray>" << std::endl;
9358 out <<
"</CellData>" << std::endl;
9362 void Mesh::PrintVTK(std::ostream &
out,
int ref,
int field_data)
9369 "# vtk DataFile Version 3.0\n"
9370 "Generated by MFEM\n"
9372 "DATASET UNSTRUCTURED_GRID\n";
9377 out <<
"FIELD FieldData 1\n"
9378 <<
"MaterialIds " << 1 <<
" " << attributes.
Size() <<
" int\n";
9379 for (
int i = 0; i < attributes.Size(); i++)
9381 out <<
' ' << attributes[i];
9388 for (
int i = 0; i < GetNE(); i++)
9397 out <<
"POINTS " << np <<
" double\n";
9399 for (
int i = 0; i < GetNE(); i++)
9402 GetElementBaseGeometry(i), ref, 1);
9404 GetElementTransformation(i)->Transform(RefG->
RefPts, pmat);
9406 for (
int j = 0; j < pmat.
Width(); j++)
9408 out << pmat(0, j) <<
' ';
9411 out << pmat(1, j) <<
' ';
9423 out << 0.0 <<
' ' << 0.0;
9430 out <<
"CELLS " << nc <<
' ' << size <<
'\n';
9432 for (
int i = 0; i < GetNE(); i++)
9439 for (
int j = 0; j < RG.
Size(); )
9442 for (
int k = 0; k < nv; k++, j++)
9444 out <<
' ' << np + RG[j];
9450 out <<
"CELL_TYPES " << nc <<
'\n';
9451 for (
int i = 0; i < GetNE(); i++)
9457 int vtk_cell_type = 5;
9461 case Geometry::POINT: vtk_cell_type = 1;
break;
9462 case Geometry::SEGMENT: vtk_cell_type = 3;
break;
9463 case Geometry::TRIANGLE: vtk_cell_type = 5;
break;
9464 case Geometry::SQUARE: vtk_cell_type = 9;
break;
9465 case Geometry::TETRAHEDRON: vtk_cell_type = 10;
break;
9466 case Geometry::CUBE: vtk_cell_type = 12;
break;
9467 case Geometry::PRISM: vtk_cell_type = 13;
break;
9469 MFEM_ABORT(
"Unrecognized VTK element type \"" << geom <<
"\"");
9473 for (
int j = 0; j < RG.
Size(); j += nv)
9475 out << vtk_cell_type <<
'\n';
9479 out <<
"CELL_DATA " << nc <<
'\n'
9480 <<
"SCALARS material int\n"
9481 <<
"LOOKUP_TABLE default\n";
9482 for (
int i = 0; i < GetNE(); i++)
9487 int attr = GetAttribute(i);
9490 out << attr <<
'\n';
9497 srand((
unsigned)time(0));
9498 double a = double(rand()) / (double(RAND_MAX) + 1.);
9499 int el0 = (int)floor(a * GetNE());
9500 GetElementColoring(coloring, el0);
9501 out <<
"SCALARS element_coloring int\n"
9502 <<
"LOOKUP_TABLE default\n";
9503 for (
int i = 0; i < GetNE(); i++)
9510 out << coloring[i] + 1 <<
'\n';
9516 out <<
"POINT_DATA " << np <<
'\n' << flush;
9521 int delete_el_to_el = (el_to_el) ? (0) : (1);
9522 const Table &el_el = ElementToElementTable();
9523 int num_el = GetNE(), stack_p, stack_top_p, max_num_col;
9526 const int *i_el_el = el_el.
GetI();
9527 const int *j_el_el = el_el.
GetJ();
9532 stack_p = stack_top_p = 0;
9533 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
9535 if (colors[el] != -2)
9541 el_stack[stack_top_p++] = el;
9543 for ( ; stack_p < stack_top_p; stack_p++)
9545 int i = el_stack[stack_p];
9546 int num_nb = i_el_el[i+1] - i_el_el[i];
9547 if (max_num_col < num_nb + 1)
9549 max_num_col = num_nb + 1;
9551 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
9554 if (colors[k] == -2)
9557 el_stack[stack_top_p++] = k;
9565 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
9567 int i = el_stack[stack_p], col;
9569 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
9571 col = colors[j_el_el[j]];
9574 col_marker[col] = 1;
9578 for (col = 0; col < max_num_col; col++)
9579 if (col_marker[col] == 0)
9587 if (delete_el_to_el)
9594 void Mesh::PrintWithPartitioning(
int *partitioning, std::ostream &
out,
9595 int elem_attr)
const
9597 if (Dim != 3 && Dim != 2) {
return; }
9599 int i, j, k, l, nv, nbe, *v;
9601 out <<
"MFEM mesh v1.0\n";
9605 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
9610 "# TETRAHEDRON = 4\n"
9615 out <<
"\ndimension\n" << Dim
9616 <<
"\n\nelements\n" << NumOfElements <<
'\n';
9617 for (i = 0; i < NumOfElements; i++)
9619 out << int((elem_attr) ? partitioning[i]+1 : elements[i]->GetAttribute())
9620 <<
' ' << elements[i]->GetGeometryType();
9621 nv = elements[i]->GetNVertices();
9622 v = elements[i]->GetVertices();
9623 for (j = 0; j < nv; j++)
9630 for (i = 0; i < faces_info.Size(); i++)
9632 if ((l = faces_info[i].Elem2No) >= 0)
9634 k = partitioning[faces_info[i].Elem1No];
9635 l = partitioning[l];
9639 if (!Nonconforming() || !IsSlaveFace(faces_info[i]))
9650 out <<
"\nboundary\n" << nbe <<
'\n';
9651 for (i = 0; i < faces_info.Size(); i++)
9653 if ((l = faces_info[i].Elem2No) >= 0)
9655 k = partitioning[faces_info[i].Elem1No];
9656 l = partitioning[l];
9659 nv = faces[i]->GetNVertices();
9660 v = faces[i]->GetVertices();
9661 out << k+1 <<
' ' << faces[i]->GetGeometryType();
9662 for (j = 0; j < nv; j++)
9667 if (!Nonconforming() || !IsSlaveFace(faces_info[i]))
9669 out << l+1 <<
' ' << faces[i]->GetGeometryType();
9670 for (j = nv-1; j >= 0; j--)
9680 k = partitioning[faces_info[i].Elem1No];
9681 nv = faces[i]->GetNVertices();
9682 v = faces[i]->GetVertices();
9683 out << k+1 <<
' ' << faces[i]->GetGeometryType();
9684 for (j = 0; j < nv; j++)
9691 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
9694 out << spaceDim <<
'\n';
9695 for (i = 0; i < NumOfVertices; i++)
9697 out << vertices[i](0);
9698 for (j = 1; j < spaceDim; j++)
9700 out <<
' ' << vertices[i](j);
9713 void Mesh::PrintElementsWithPartitioning(
int *partitioning,
9717 MFEM_ASSERT(Dim == spaceDim,
"2D Manifolds not supported\n");
9718 if (Dim != 3 && Dim != 2) {
return; }
9725 int *vcount =
new int[NumOfVertices];
9726 for (i = 0; i < NumOfVertices; i++)
9730 for (i = 0; i < NumOfElements; i++)
9732 nv = elements[i]->GetNVertices();
9733 ind = elements[i]->GetVertices();
9734 for (j = 0; j < nv; j++)
9740 int *voff =
new int[NumOfVertices+1];
9742 for (i = 1; i <= NumOfVertices; i++)
9744 voff[i] = vcount[i-1] + voff[i-1];
9747 int **vown =
new int*[NumOfVertices];
9748 for (i = 0; i < NumOfVertices; i++)
9750 vown[i] =
new int[vcount[i]];
9760 Transpose(ElementToEdgeTable(), edge_el);
9763 for (i = 0; i < NumOfElements; i++)
9765 nv = elements[i]->GetNVertices();
9766 ind = elements[i]->GetVertices();
9767 for (j = 0; j < nv; j++)
9770 vown[ind[j]][vcount[ind[j]]] = i;
9774 for (i = 0; i < NumOfVertices; i++)
9776 vcount[i] = voff[i+1] - voff[i];
9780 for (i = 0; i < edge_el.
Size(); i++)
9782 const int *el = edge_el.
GetRow(i);
9785 k = partitioning[el[0]];
9786 l = partitioning[el[1]];
9787 if (interior_faces || k != l)
9799 out <<
"areamesh2\n\n" << nbe <<
'\n';
9801 for (i = 0; i < edge_el.
Size(); i++)
9803 const int *el = edge_el.
GetRow(i);
9806 k = partitioning[el[0]];
9807 l = partitioning[el[1]];
9808 if (interior_faces || k != l)
9811 GetEdgeVertices(i,ev);
9813 for (j = 0; j < 2; j++)
9814 for (s = 0; s < vcount[ev[j]]; s++)
9815 if (vown[ev[j]][s] == el[0])
9817 out <<
' ' << voff[ev[j]]+s+1;
9821 for (j = 1; j >= 0; j--)
9822 for (s = 0; s < vcount[ev[j]]; s++)
9823 if (vown[ev[j]][s] == el[1])
9825 out <<
' ' << voff[ev[j]]+s+1;
9832 k = partitioning[el[0]];
9834 GetEdgeVertices(i,ev);
9836 for (j = 0; j < 2; j++)
9837 for (s = 0; s < vcount[ev[j]]; s++)
9838 if (vown[ev[j]][s] == el[0])
9840 out <<
' ' << voff[ev[j]]+s+1;
9847 out << NumOfElements <<
'\n';
9848 for (i = 0; i < NumOfElements; i++)
9850 nv = elements[i]->GetNVertices();
9851 ind = elements[i]->GetVertices();
9852 out << partitioning[i]+1 <<
' ';
9854 for (j = 0; j < nv; j++)
9856 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
9857 vown[ind[j]][vcount[ind[j]]] = i;
9862 for (i = 0; i < NumOfVertices; i++)
9864 vcount[i] = voff[i+1] - voff[i];
9868 out << voff[NumOfVertices] <<
'\n';
9869 for (i = 0; i < NumOfVertices; i++)
9870 for (k = 0; k < vcount[i]; k++)
9872 for (j = 0; j < Dim; j++)
9874 out << vertices[i](j) <<
' ';
9880 else if (meshgen == 1)
9882 out <<
"NETGEN_Neutral_Format\n";
9884 out << voff[NumOfVertices] <<
'\n';
9885 for (i = 0; i < NumOfVertices; i++)
9886 for (k = 0; k < vcount[i]; k++)
9888 for (j = 0; j < Dim; j++)
9890 out <<
' ' << vertices[i](j);
9896 out << NumOfElements <<
'\n';
9897 for (i = 0; i < NumOfElements; i++)
9899 nv = elements[i]->GetNVertices();
9900 ind = elements[i]->GetVertices();
9901 out << partitioning[i]+1;
9902 for (j = 0; j < nv; j++)
9904 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
9905 vown[ind[j]][vcount[ind[j]]] = i;
9910 for (i = 0; i < NumOfVertices; i++)
9912 vcount[i] = voff[i+1] - voff[i];
9918 for (i = 0; i < NumOfFaces; i++)
9919 if ((l = faces_info[i].Elem2No) >= 0)
9921 k = partitioning[faces_info[i].Elem1No];
9922 l = partitioning[l];
9923 if (interior_faces || k != l)
9934 for (i = 0; i < NumOfFaces; i++)
9935 if ((l = faces_info[i].Elem2No) >= 0)
9937 k = partitioning[faces_info[i].Elem1No];
9938 l = partitioning[l];
9939 if (interior_faces || k != l)
9941 nv = faces[i]->GetNVertices();
9942 ind = faces[i]->GetVertices();
9944 for (j = 0; j < nv; j++)
9945 for (s = 0; s < vcount[ind[j]]; s++)
9946 if (vown[ind[j]][s] == faces_info[i].Elem1No)
9948 out <<
' ' << voff[ind[j]]+s+1;
9952 for (j = nv-1; j >= 0; j--)
9953 for (s = 0; s < vcount[ind[j]]; s++)
9954 if (vown[ind[j]][s] == faces_info[i].Elem2No)
9956 out <<
' ' << voff[ind[j]]+s+1;
9963 k = partitioning[faces_info[i].Elem1No];
9964 nv = faces[i]->GetNVertices();
9965 ind = faces[i]->GetVertices();
9967 for (j = 0; j < nv; j++)
9968 for (s = 0; s < vcount[ind[j]]; s++)
9969 if (vown[ind[j]][s] == faces_info[i].Elem1No)
9971 out <<
' ' << voff[ind[j]]+s+1;
9977 else if (meshgen == 2)
9982 for (i = 0; i < NumOfFaces; i++)
9983 if ((l = faces_info[i].Elem2No) >= 0)
9985 k = partitioning[faces_info[i].Elem1No];
9986 l = partitioning[l];
9987 if (interior_faces || k != l)
9999 <<
"1 " << voff[NumOfVertices] <<
" " << NumOfElements
10000 <<
" 0 0 0 0 0 0 0\n"
10001 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
10002 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
10003 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
10004 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
10006 for (i = 0; i < NumOfVertices; i++)
10007 for (k = 0; k < vcount[i]; k++)
10008 out << voff[i]+k <<
" 0.0 " << vertices[i](0) <<
' '
10009 << vertices[i](1) <<
' ' << vertices[i](2) <<
" 0.0\n";
10011 for (i = 0; i < NumOfElements; i++)
10013 nv = elements[i]->GetNVertices();
10014 ind = elements[i]->GetVertices();
10015 out << i+1 <<
' ' << partitioning[i]+1;
10016 for (j = 0; j < nv; j++)
10018 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
10019 vown[ind[j]][vcount[ind[j]]] = i;
10024 for (i = 0; i < NumOfVertices; i++)
10026 vcount[i] = voff[i+1] - voff[i];
10030 for (i = 0; i < NumOfFaces; i++)
10031 if ((l = faces_info[i].Elem2No) >= 0)
10033 k = partitioning[faces_info[i].Elem1No];
10034 l = partitioning[l];
10035 if (interior_faces || k != l)
10037 nv = faces[i]->GetNVertices();
10038 ind = faces[i]->GetVertices();
10040 for (j = 0; j < nv; j++)
10041 for (s = 0; s < vcount[ind[j]]; s++)
10042 if (vown[ind[j]][s] == faces_info[i].Elem1No)
10044 out <<
' ' << voff[ind[j]]+s+1;
10046 out <<
" 1.0 1.0 1.0 1.0\n";
10048 for (j = nv-1; j >= 0; j--)
10049 for (s = 0; s < vcount[ind[j]]; s++)
10050 if (vown[ind[j]][s] == faces_info[i].Elem2No)
10052 out <<
' ' << voff[ind[j]]+s+1;
10054 out <<
" 1.0 1.0 1.0 1.0\n";
10059 k = partitioning[faces_info[i].Elem1No];
10060 nv = faces[i]->GetNVertices();
10061 ind = faces[i]->GetVertices();
10063 for (j = 0; j < nv; j++)
10064 for (s = 0; s < vcount[ind[j]]; s++)
10065 if (vown[ind[j]][s] == faces_info[i].Elem1No)
10067 out <<
' ' << voff[ind[j]]+s+1;
10069 out <<
" 1.0 1.0 1.0 1.0\n";
10075 for (i = 0; i < NumOfVertices; i++)
10085 void Mesh::PrintSurfaces(
const Table & Aface_face, std::ostream &
out)
const
10092 " NURBS mesh is not supported!");
10096 out <<
"MFEM mesh v1.0\n";
10100 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
10105 "# TETRAHEDRON = 4\n"
10110 out <<
"\ndimension\n" << Dim
10111 <<
"\n\nelements\n" << NumOfElements <<
'\n';
10112 for (i = 0; i < NumOfElements; i++)
10114 PrintElement(elements[i], out);
10118 const int *
const i_AF_f = Aface_face.
GetI();
10119 const int *
const j_AF_f = Aface_face.
GetJ();
10121 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
10122 for (
const int * iface = j_AF_f + i_AF_f[iAF];
10123 iface < j_AF_f + i_AF_f[iAF+1];
10126 out << iAF+1 <<
' ';
10127 PrintElementWithoutAttr(faces[*iface],out);
10130 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
10133 out << spaceDim <<
'\n';
10134 for (i = 0; i < NumOfVertices; i++)
10136 out << vertices[i](0);
10137 for (j = 1; j < spaceDim; j++)
10139 out <<
' ' << vertices[i](j);
10147 out <<
"\nnodes\n";
10152 void Mesh::ScaleSubdomains(
double sf)
10157 int na = attributes.
Size();
10158 double *cg =
new double[na*spaceDim];
10159 int *nbea =
new int[na];
10161 int *vn =
new int[NumOfVertices];
10162 for (i = 0; i < NumOfVertices; i++)
10166 for (i = 0; i < na; i++)
10168 for (j = 0; j < spaceDim; j++)
10170 cg[i*spaceDim+j] = 0.0;
10175 for (i = 0; i < NumOfElements; i++)
10177 GetElementVertices(i, vert);
10178 for (k = 0; k < vert.
Size(); k++)
10184 for (i = 0; i < NumOfElements; i++)
10186 int bea = GetAttribute(i)-1;
10187 GetPointMatrix(i, pointmat);
10188 GetElementVertices(i, vert);
10190 for (k = 0; k < vert.
Size(); k++)
10191 if (vn[vert[k]] == 1)
10194 for (j = 0; j < spaceDim; j++)
10196 cg[bea*spaceDim+j] += pointmat(j,k);
10202 for (i = 0; i < NumOfElements; i++)
10204 int bea = GetAttribute(i)-1;
10205 GetElementVertices (i, vert);
10207 for (k = 0; k < vert.
Size(); k++)
10210 for (j = 0; j < spaceDim; j++)
10211 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
10212 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
10222 void Mesh::ScaleElements(
double sf)
10227 int na = NumOfElements;
10228 double *cg =
new double[na*spaceDim];
10229 int *nbea =
new int[na];
10231 int *vn =
new int[NumOfVertices];
10232 for (i = 0; i < NumOfVertices; i++)
10236 for (i = 0; i < na; i++)
10238 for (j = 0; j < spaceDim; j++)
10240 cg[i*spaceDim+j] = 0.0;
10245 for (i = 0; i < NumOfElements; i++)
10247 GetElementVertices(i, vert);
10248 for (k = 0; k < vert.
Size(); k++)
10254 for (i = 0; i < NumOfElements; i++)
10257 GetPointMatrix(i, pointmat);
10258 GetElementVertices(i, vert);
10260 for (k = 0; k < vert.
Size(); k++)
10261 if (vn[vert[k]] == 1)
10264 for (j = 0; j < spaceDim; j++)
10266 cg[bea*spaceDim+j] += pointmat(j,k);
10272 for (i = 0; i < NumOfElements; i++)
10275 GetElementVertices(i, vert);
10277 for (k = 0; k < vert.
Size(); k++)
10280 for (j = 0; j < spaceDim; j++)
10281 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
10282 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
10297 Vector vold(spaceDim), vnew(NULL, spaceDim);
10298 for (
int i = 0; i < vertices.Size(); i++)
10300 for (
int j = 0; j < spaceDim; j++)
10302 vold(j) = vertices[i](j);
10312 xnew.ProjectCoefficient(f_pert);
10319 MFEM_VERIFY(spaceDim == deformation.
GetVDim(),
10320 "incompatible vector dimensions");
10327 for (
int i = 0; i < NumOfVertices; i++)
10328 for (
int d = 0; d < spaceDim; d++)
10330 vertices[i](d) = xnew(d + spaceDim*i);
10341 void Mesh::RemoveUnusedVertices()
10343 if (NURBSext || ncmesh) {
return; }
10347 for (
int i = 0; i < GetNE(); i++)
10352 for (
int j = 0; j < nv; j++)
10357 for (
int i = 0; i < GetNBE(); i++)
10359 Element *el = GetBdrElement(i);
10362 for (
int j = 0; j < nv; j++)
10368 for (
int i = 0; i < v2v.
Size(); i++)
10372 vertices[num_vert] = vertices[i];
10373 v2v[i] = num_vert++;
10377 if (num_vert == v2v.
Size()) {
return; }
10379 Vector nodes_by_element;
10384 for (
int i = 0; i < GetNE(); i++)
10386 Nodes->FESpace()->GetElementVDofs(i, vdofs);
10391 for (
int i = 0; i < GetNE(); i++)
10393 Nodes->FESpace()->GetElementVDofs(i, vdofs);
10394 Nodes->GetSubVector(vdofs, &nodes_by_element(s));
10398 vertices.SetSize(num_vert);
10399 NumOfVertices = num_vert;
10400 for (
int i = 0; i < GetNE(); i++)
10405 for (
int j = 0; j < nv; j++)
10410 for (
int i = 0; i < GetNBE(); i++)
10412 Element *el = GetBdrElement(i);
10415 for (
int j = 0; j < nv; j++)
10424 el_to_edge =
new Table;
10425 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
10430 GetElementToFaceTable();
10436 Nodes->FESpace()->Update();
10439 for (
int i = 0; i < GetNE(); i++)
10441 Nodes->FESpace()->GetElementVDofs(i, vdofs);
10442 Nodes->SetSubVector(vdofs, &nodes_by_element(s));
10448 void Mesh::RemoveInternalBoundaries()
10450 if (NURBSext || ncmesh) {
return; }
10452 int num_bdr_elem = 0;
10453 int new_bel_to_edge_nnz = 0;
10454 for (
int i = 0; i < GetNBE(); i++)
10456 if (FaceIsInterior(GetBdrElementEdgeIndex(i)))
10458 FreeElement(boundary[i]);
10465 new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
10470 if (num_bdr_elem == GetNBE()) {
return; }
10474 Table *new_bel_to_edge = NULL;
10478 new_be_to_edge.
Reserve(num_bdr_elem);
10482 new_be_to_face.
Reserve(num_bdr_elem);
10483 new_bel_to_edge =
new Table;
10484 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
10486 for (
int i = 0; i < GetNBE(); i++)
10488 if (!FaceIsInterior(GetBdrElementEdgeIndex(i)))
10490 new_boundary.
Append(boundary[i]);
10493 new_be_to_edge.
Append(be_to_edge[i]);
10497 int row = new_be_to_face.
Size();
10498 new_be_to_face.
Append(be_to_face[i]);
10499 int *e = bel_to_edge->GetRow(i);
10500 int ne = bel_to_edge->RowSize(i);
10501 int *new_e = new_bel_to_edge->
GetRow(row);
10502 for (
int j = 0; j < ne; j++)
10506 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
10511 NumOfBdrElements = new_boundary.
Size();
10521 delete bel_to_edge;
10522 bel_to_edge = new_bel_to_edge;
10526 for (
int i = 0; i < attribs.
Size(); i++)
10528 attribs[i] = GetBdrAttribute(i);
10532 bdr_attributes.DeleteAll();
10533 attribs.
Copy(bdr_attributes);
10538 #ifdef MFEM_USE_MEMALLOC
10541 if (E->
GetType() == Element::TETRAHEDRON)
10565 const int npts = point_mat.
Width();
10566 if (!npts) {
return 0; }
10567 MFEM_VERIFY(point_mat.
Height() == spaceDim,
"Invalid points matrix");
10571 if (!GetNE()) {
return 0; }
10573 double *data = point_mat.
GetData();
10580 min_dist = std::numeric_limits<double>::max();
10584 for (
int i = 0; i < GetNE(); i++)
10586 GetElementTransformation(i)->Transform(
10588 for (
int k = 0; k < npts; k++)
10590 double dist = pt.
DistanceTo(data+k*spaceDim);
10591 if (dist < min_dist(k))
10593 min_dist(k) = dist;
10602 for (
int k = 0; k < npts; k++)
10606 int res = inv_tr->
Transform(pt, ips[k]);
10607 if (res == InverseElementTransformation::Inside)
10609 elem_ids[k] = e_idx[k];
10613 if (pts_found != npts)
10616 Table *vtoel = GetVertexToElementTable();
10617 for (
int k = 0; k < npts; k++)
10619 if (elem_ids[k] != -1) {
continue; }
10622 GetElementVertices(e_idx[k], vertices);
10623 for (
int v = 0; v < vertices.
Size(); v++)
10625 int vv = vertices[v];
10627 const int* els = vtoel->
GetRow(vv);
10628 for (
int e = 0; e < ne; e++)
10630 if (els[e] == e_idx[k]) {
continue; }
10632 int res = inv_tr->
Transform(pt, ips[k]);
10633 if (res == InverseElementTransformation::Inside)
10635 elem_ids[k] = els[e];
10645 int le = ncmesh->leaf_elements[e_idx[k]];
10646 ncmesh->FindNeighbors(le,neigh);
10647 for (
int e = 0; e < neigh.
Size(); e++)
10650 if (ncmesh->IsGhost(ncmesh->elements[nn])) {
continue; }
10651 int el = ncmesh->elements[nn].index;
10653 int res = inv_tr->
Transform(pt, ips[k]);
10654 if (res == InverseElementTransformation::Inside)
10666 if (inv_trans == NULL) {
delete inv_tr; }
10668 if (warn && pts_found != npts)
10670 MFEM_WARNING((npts-pts_found) <<
" points were not found");
10681 computed_factors = flags;
10687 const int vdim = fespace->
GetVDim();
10688 const int NE = fespace->
GetNE();
10689 const int ND = fe->
GetDof();
10694 ElementDofOrdering::NATIVE);
10696 unsigned eval_flags = 0;
10697 if (flags & GeometricFactors::COORDINATES)
10699 X.SetSize(vdim*NQ*NE);
10700 eval_flags |= QuadratureInterpolator::VALUES;
10702 if (flags & GeometricFactors::JACOBIANS)
10704 J.SetSize(dim*vdim*NQ*NE);
10705 eval_flags |= QuadratureInterpolator::DERIVATIVES;
10707 if (flags & GeometricFactors::DETERMINANTS)
10709 detJ.SetSize(NQ*NE);
10710 eval_flags |= QuadratureInterpolator::DETERMINANTS;
10719 Vector Enodes(vdim*ND*NE);
10720 elem_restr->
Mult(*nodes, Enodes);
10721 qi->
Mult(Enodes, eval_flags, X, J, detJ);
10725 qi->
Mult(*nodes, eval_flags, X, J, detJ);
10729 FaceGeometricFactors::FaceGeometricFactors(
const Mesh *mesh,
10740 const int vdim = fespace->
GetVDim();
10749 face_restr->
Mult(*nodes, Fnodes);
10751 unsigned eval_flags = 0;
10792 V(1) = s * ((ip.
y + layer) / n);
10797 V(2) = s * ((ip.
z + layer) / n);
10806 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
10810 int nvy = (closed) ? (ny) : (ny + 1);
10811 int nvt = mesh->
GetNV() * nvy;
10820 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
10825 for (
int i = 0; i < mesh->
GetNV(); i++)
10828 for (
int j = 0; j < nvy; j++)
10830 vc[1] = sy * (double(j) / ny);
10836 for (
int i = 0; i < mesh->
GetNE(); i++)
10841 for (
int j = 0; j < ny; j++)
10844 qv[0] = vert[0] * nvy + j;
10845 qv[1] = vert[1] * nvy + j;
10846 qv[2] = vert[1] * nvy + (j + 1) % nvy;
10847 qv[3] = vert[0] * nvy + (j + 1) % nvy;
10853 for (
int i = 0; i < mesh->
GetNBE(); i++)
10858 for (
int j = 0; j < ny; j++)
10861 sv[0] = vert[0] * nvy + j;
10862 sv[1] = vert[0] * nvy + (j + 1) % nvy;
10866 Swap<int>(sv[0], sv[1]);
10878 for (
int i = 0; i < mesh->
GetNE(); i++)
10884 sv[0] = vert[0] * nvy;
10885 sv[1] = vert[1] * nvy;
10889 sv[0] = vert[1] * nvy + ny;
10890 sv[1] = vert[0] * nvy + ny;
10906 string cname = name;
10907 if (cname ==
"Linear")
10911 else if (cname ==
"Quadratic")
10915 else if (cname ==
"Cubic")
10919 else if (!strncmp(name,
"H1_", 3))
10923 else if (!strncmp(name,
"L2_T", 4))
10927 else if (!strncmp(name,
"L2_", 3))
10934 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
10946 for (
int i = 0; i < mesh->
GetNE(); i++)
10949 for (
int j = ny-1; j >= 0; j--)
10966 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
10971 int nvt = mesh->
GetNV() * nvz;
10976 bool wdgMesh =
false;
10977 bool hexMesh =
false;
10981 for (
int i = 0; i < mesh->
GetNV(); i++)
10985 for (
int j = 0; j < nvz; j++)
10987 vc[2] = sz * (double(j) / nz);
10993 for (
int i = 0; i < mesh->
GetNE(); i++)
11003 for (
int j = 0; j < nz; j++)
11006 pv[0] = vert[0] * nvz + j;
11007 pv[1] = vert[1] * nvz + j;
11008 pv[2] = vert[2] * nvz + j;
11009 pv[3] = vert[0] * nvz + (j + 1) % nvz;
11010 pv[4] = vert[1] * nvz + (j + 1) % nvz;
11011 pv[5] = vert[2] * nvz + (j + 1) % nvz;
11018 for (
int j = 0; j < nz; j++)
11021 hv[0] = vert[0] * nvz + j;
11022 hv[1] = vert[1] * nvz + j;
11023 hv[2] = vert[2] * nvz + j;
11024 hv[3] = vert[3] * nvz + j;
11025 hv[4] = vert[0] * nvz + (j + 1) % nvz;
11026 hv[5] = vert[1] * nvz + (j + 1) % nvz;
11027 hv[6] = vert[2] * nvz + (j + 1) % nvz;
11028 hv[7] = vert[3] * nvz + (j + 1) % nvz;
11030 mesh3d->
AddHex(hv, attr);
11034 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
11035 << geom <<
"\'" << endl;
11041 for (
int i = 0; i < mesh->
GetNBE(); i++)
11046 for (
int j = 0; j < nz; j++)
11049 qv[0] = vert[0] * nvz + j;
11050 qv[1] = vert[1] * nvz + j;
11051 qv[2] = vert[1] * nvz + (j + 1) % nvz;
11052 qv[3] = vert[0] * nvz + (j + 1) % nvz;
11061 for (
int i = 0; i < mesh->
GetNE(); i++)
11072 tv[0] = vert[0] * nvz;
11073 tv[1] = vert[2] * nvz;
11074 tv[2] = vert[1] * nvz;
11078 tv[0] = vert[0] * nvz + nz;
11079 tv[1] = vert[1] * nvz + nz;
11080 tv[2] = vert[2] * nvz + nz;
11088 qv[0] = vert[0] * nvz;
11089 qv[1] = vert[3] * nvz;
11090 qv[2] = vert[2] * nvz;
11091 qv[3] = vert[1] * nvz;
11095 qv[0] = vert[0] * nvz + nz;
11096 qv[1] = vert[1] * nvz + nz;
11097 qv[2] = vert[2] * nvz + nz;
11098 qv[3] = vert[3] * nvz + nz;
11104 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
11105 << geom <<
"\'" << endl;
11111 if ( hexMesh && wdgMesh )
11115 else if ( hexMesh )
11119 else if ( wdgMesh )
11132 string cname = name;
11133 if (cname ==
"Linear")
11137 else if (cname ==
"Quadratic")
11141 else if (cname ==
"Cubic")
11145 else if (!strncmp(name,
"H1_", 3))
11149 else if (!strncmp(name,
"L2_T", 4))
11153 else if (!strncmp(name,
"L2_", 3))
11160 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
11172 for (
int i = 0; i < mesh->
GetNE(); i++)
11175 for (
int j = nz-1; j >= 0; j--)
11196 out << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
11197 <<
" 0 0 " << i <<
" -1 0\n";
11204 double mid[3] = {0, 0, 0};
11205 for (
int j = 0; j < 2; j++)
11207 for (
int k = 0; k <
spaceDim; k++)
11212 out << NumOfVertices+i <<
" "
11213 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
11214 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
int GetNPoints() const
Returns the number of the points in the integration rule.
Abstract class for all finite elements.
const IntegrationRule * IntRule
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
int Size() const
For backward compatibility define Size to be synonym of Width()
Geometry::Type GetGeometryType() const
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 Get(double *p, const int dim) const
virtual void Print(std::ostream &out=mfem::out) const
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.
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...
T * end()
STL-like end. Returns pointer after the last element of the array.
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)
virtual Element * Duplicate(Mesh *m) const =0
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Vector J
Jacobians of the element transformations at all quadrature points.
void AddColumnsInRow(int r, int ncol)
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 ...
Base class for vector Coefficients that optionally depend on time and space.
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 GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
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 MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of fec and fes.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Lists all edges/faces in the nonconforming mesh.
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().
void ShiftRight(int &a, int &b, int &c)
void SetDims(int rows, int nnz)
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
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.
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.
void GetMeshComponents(Mesh &mesh) const
Fill Mesh::{vertices,elements,boundary} for the current finest level.
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 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.
void Save(std::ostream &out, int fmt=0) const
Save the Array to the stream out using the format fmt. The format fmt can be:
void SetSize(int i, int j, int k)
int GetNE() const
Returns number of elements.
void FinalizeWedgeMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a wedge Mesh.
void order(Functional *functional, uint iterations=1, uint window=2, uint period=2, uint seed=0, Progress *progress=0)
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
void GetElementLocalToGlobal(Array< int > &lelem_elem)
void SetOutputLayout(QVectorLayout out_layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
Data type quadrilateral element.
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)
double * GetData() const
Return a pointer to the beginning of the Vector data.
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
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.
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
void DeleteAll()
Delete the whole array.
void AddConnections(int r, const int *c, int nc)
int master
master number (in Mesh numbering)
Array< NCFaceInfo > nc_faces_info
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Structure for storing face geometric factors: coordinates, Jacobians, determinants of the Jacobians...
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)
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 DisableTensorProducts(bool disable=true) const
Disable the use of tensor product evaluations, for tensor-product elements, e.g. quads and hexes...
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.
const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Data type hexahedron element.
void WriteBase64WithSizeAndClear(std::ostream &out, std::vector< char > &buf, int compression_level)
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)
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.
Symmetric 3D Table stored as an array of rows each of which has a stack of column, floor, number nodes. The number of the node is assigned by counting the nodes from zero as they are pushed into the table. Diagonals of any kind are not allowed so the row, column and floor must all be different for each node. Only one node is stored for all 6 symmetric entries that are indexable by unique triplets of row, column, and floor.
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Append(const T &el)
Append element 'el' to array, resize if necessary.
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
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 METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void AddConnection(int r, int c)
void WriteBinaryOrASCII< uint8_t >(std::ostream &out, std::vector< char > &buf, const uint8_t &val, const char *suffix, VTKFormat format)
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.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
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)
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
GeometryRefiner GlobGeometryRefiner
void SetLayer(const int l)
int GetAttribute() const
Return element's attribute.
T * begin()
STL-like begin. Returns pointer to the first element of the array.
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
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.
signed char local
local number within 'element'
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
int GetVDim()
Returns dimension of the vector.
void SetNodalFESpace(FiniteElementSpace *nfes)
virtual void ResetTransform(int tr)
Set current coarse-fine transformation number.
int Size() const
Returns the number of TYPE I elements.
int GetVDim() const
Returns vector dimension.
Vector normal
Normals at all quadrature points.
A class for non-conforming AMR on higher-order hexahedral, prismatic, quadrilateral or triangular mes...
FiniteElementSpace * FESpace()
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...
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
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 SpaceDimension() const
virtual const Operator * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
virtual unsigned GetTransform() const
Return current coarse-fine transformation.
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
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 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)
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...
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
Vector X
Mapped (physical) coordinates of all quadrature points.
void AddAColumnInRow(int r)
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
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.
int FindRoots(const Vector &z, Vector &x)
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.
Helper struct for defining a connectivity table, see Table::MakeFromList.
class Linear3DFiniteElement TetrahedronFE
A class to initialize the size of a Tensor.
Array< Element * > elements
Linear2DFiniteElement TriangleFE
Evaluate the values at quadrature points.
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i'th boundary element.
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...
virtual const char * Name() const
Class for integration point with weight.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
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.
NodeExtrudeCoefficient(const int dim, const int _n, const double _s)
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
void Mult(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_der, Vector &q_det) const
Interpolate the E-vector e_vec to quadrature points.
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
const char * VTKByteOrder()
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)
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)
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.
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)
Piecewise-(bi)quadratic continuous finite elements.
void AppendBytes(std::vector< char > &vec, const T &val)
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
NCMesh * ncmesh
Optional non-conforming mesh extension.
Geometry::Type GetBdrElementBaseGeometry(int i) const
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
BiLinear2DFiniteElement QuadrilateralFE
Evaluate the derivatives at quadrature points.
const FiniteElementCollection * FEColl() const
void GetNodes(Vector &node_coord) const
int NumberOfEntries() const
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)
const std::string filename
MemAlloc< Tetrahedron, 1024 > TetMemory
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
Node::Index insert_node(Float length=1)
TriLinear3DFiniteElement HexahedronFE
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
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...
Rank 3 tensor (array of matrices)
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.
Data type line segment element.
Linear1DFiniteElement SegmentFE
Array< GeometricFactors * > geom_factors
Optional geometric factors.
int GetNFDofs() const
Number of all scalar face-interior dofs.
MPI_Comm GetGlobalMPI_Comm()
Get MFEM's "global" MPI communicator.
int NumberOfElements()
Return the number of elements added to the table.
virtual Type GetType() const =0
Returns element's type.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
const Element * GetBdrElement(int i) const
Defines the position of a fine element within a coarse element.
int matrix
Index into the DenseTensor corresponding to the parent Geometry::Type stored in CoarseFineTransformat...
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
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...