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 "../fem/quadinterpolator.hpp"
33 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
38 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
43 int*,
int*,
int*,
int*,
int*,
idxtype*);
45 int*,
int*,
int*,
int*,
int*,
idxtype*);
47 int*,
int*,
int*,
int*,
int*,
idxtype*);
68 void Mesh::GetElementCenter(
int i,
Vector ¢er)
71 int geom = GetElementBaseGeometry(i);
76 double Mesh::GetElementSize(
int i,
int type)
79 GetElementJacobian(i, J);
82 return pow(fabs(J.
Det()), 1./Dim);
94 double Mesh::GetElementSize(
int i,
const Vector &dir)
98 GetElementJacobian(i, J);
100 return sqrt((d_hat * d_hat) / (dir * dir));
103 double Mesh::GetElementVolume(
int i)
125 for (
int d = 0; d < spaceDim; d++)
134 for (
int i = 0; i < NumOfVertices; i++)
136 coord = GetVertex(i);
137 for (
int d = 0; d < spaceDim; d++)
139 if (coord[d] < min(d)) { min(d) = coord[d]; }
140 if (coord[d] > max(d)) { max(d) = coord[d]; }
146 const bool use_boundary =
false;
147 int ne = use_boundary ? GetNBE() : GetNE();
155 for (
int i = 0; i < ne; i++)
159 GetBdrElementFace(i, &fn, &fo);
161 Tr = GetFaceElementTransformations(fn, 5);
168 T = GetElementTransformation(i);
172 for (
int j = 0; j < pointmat.
Width(); j++)
174 for (
int d = 0; d < pointmat.
Height(); d++)
176 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
177 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
184 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
185 double &kappa_min,
double &kappa_max,
193 sdim = SpaceDimension();
195 if (Vh) { Vh->
SetSize(NumOfElements); }
196 if (Vk) { Vk->
SetSize(NumOfElements); }
199 h_max = kappa_max = -h_min;
200 if (dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
202 for (i = 0; i < NumOfElements; i++)
204 GetElementJacobian(i, J);
205 h = pow(fabs(J.
Weight()), 1.0/
double(dim));
206 kappa = (dim == sdim) ?
208 if (Vh) { (*Vh)(i) = h; }
209 if (Vk) { (*Vk)(i) = kappa; }
211 if (h < h_min) { h_min = h; }
212 if (h > h_max) { h_max = h; }
213 if (kappa < kappa_min) { kappa_min =
kappa; }
214 if (kappa > kappa_max) { kappa_max =
kappa; }
219 void Mesh::PrintElementsByGeometry(
int dim,
223 for (
int g = Geometry::DimStart[dim], first = 1;
224 g < Geometry::DimStart[dim+1]; g++)
226 if (!num_elems_by_geom[g]) {
continue; }
227 if (!first) { out <<
" + "; }
229 out << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
235 double h_min, h_max, kappa_min, kappa_max;
237 out <<
"Mesh Characteristics:";
239 this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
241 Array<int> num_elems_by_geom(Geometry::NumGeom);
242 num_elems_by_geom = 0;
243 for (
int i = 0; i < GetNE(); i++)
245 num_elems_by_geom[GetElementBaseGeometry(i)]++;
249 <<
"Dimension : " << Dimension() <<
'\n'
250 <<
"Space dimension : " << SpaceDimension();
254 <<
"Number of vertices : " << GetNV() <<
'\n'
255 <<
"Number of elements : " << GetNE() <<
'\n'
256 <<
"Number of bdr elem : " << GetNBE() <<
'\n';
261 <<
"Number of vertices : " << GetNV() <<
'\n'
262 <<
"Number of elements : " << GetNE() <<
'\n'
263 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
264 <<
"h_min : " << h_min <<
'\n'
265 <<
"h_max : " << h_max <<
'\n';
270 <<
"Number of vertices : " << GetNV() <<
'\n'
271 <<
"Number of edges : " << GetNEdges() <<
'\n'
272 <<
"Number of elements : " << GetNE() <<
" -- ";
273 PrintElementsByGeometry(2, num_elems_by_geom, out);
275 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
276 <<
"Euler Number : " << EulerNumber2D() <<
'\n'
277 <<
"h_min : " << h_min <<
'\n'
278 <<
"h_max : " << h_max <<
'\n'
279 <<
"kappa_min : " << kappa_min <<
'\n'
280 <<
"kappa_max : " << kappa_max <<
'\n';
284 Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
285 num_bdr_elems_by_geom = 0;
286 for (
int i = 0; i < GetNBE(); i++)
288 num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
290 Array<int> num_faces_by_geom(Geometry::NumGeom);
291 num_faces_by_geom = 0;
292 for (
int i = 0; i < GetNFaces(); i++)
294 num_faces_by_geom[GetFaceBaseGeometry(i)]++;
298 <<
"Number of vertices : " << GetNV() <<
'\n'
299 <<
"Number of edges : " << GetNEdges() <<
'\n'
300 <<
"Number of faces : " << GetNFaces() <<
" -- ";
301 PrintElementsByGeometry(Dim-1, num_faces_by_geom, out);
303 <<
"Number of elements : " << GetNE() <<
" -- ";
304 PrintElementsByGeometry(Dim, num_elems_by_geom, out);
306 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
307 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, out);
309 <<
"Euler Number : " << EulerNumber() <<
'\n'
310 <<
"h_min : " << h_min <<
'\n'
311 <<
"h_max : " << h_max <<
'\n'
312 <<
"kappa_min : " << kappa_min <<
'\n'
313 <<
"kappa_max : " << kappa_max <<
'\n';
315 out <<
'\n' << std::flush;
322 case Element::POINT :
return &
PointFE;
323 case Element::SEGMENT :
return &
SegmentFE;
328 case Element::WEDGE :
return &
WedgeFE;
330 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
333 MFEM_ABORT(
"Unknown element type");
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));
368 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
377 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
378 int nv = elements[i]->GetNVertices();
379 const int *v = elements[i]->GetVertices();
380 int n = vertices.Size();
382 for (
int k = 0; k < spaceDim; k++)
384 for (
int j = 0; j < nv; j++)
386 pm(k, j) = nodes(k*n+v[j]);
389 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
393 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
395 Nodes->FESpace()->GetElementVDofs(i, vdofs);
396 int n = vdofs.
Size()/spaceDim;
398 for (
int k = 0; k < spaceDim; k++)
400 for (
int j = 0; j < n; j++)
402 pm(k,j) = nodes(vdofs[n*k+j]);
405 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
412 GetElementTransformation(i, &Transformation);
414 return &Transformation;
419 GetBdrElementTransformation(i, &BdrTransformation);
420 return &BdrTransformation;
430 GetBdrPointMatrix(i, pm);
431 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
441 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
442 int n = vdofs.
Size()/spaceDim;
444 for (
int k = 0; k < spaceDim; k++)
446 for (
int j = 0; j < n; j++)
448 pm(k,j) = nodes(vdofs[n*k+j]);
455 int elem_id, face_info;
456 GetBdrElementAdjacentElement(i, elem_id, face_info);
458 GetLocalFaceTransformation(GetBdrElementType(i),
459 GetElementType(elem_id),
460 FaceElemTr.Loc1.Transf, face_info);
464 Nodes->FESpace()->GetTraceElement(elem_id,
465 GetBdrElementBaseGeometry(i));
468 FaceElemTr.Loc1.Transf.ElementNo = elem_id;
469 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
470 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
472 ElTr->
SetFE(face_el);
480 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
485 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
486 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
488 for (
int i = 0; i < spaceDim; i++)
490 for (
int j = 0; j < nv; j++)
492 pm(i, j) = vertices[v[j]](i);
495 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
499 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
505 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
506 int n = vdofs.
Size()/spaceDim;
508 for (
int i = 0; i < spaceDim; i++)
510 for (
int j = 0; j < n; j++)
512 pm(i, j) = nodes(vdofs[n*i+j]);
519 FaceInfo &face_info = faces_info[FaceNo];
524 GetLocalFaceTransformation(face_type,
525 GetElementType(face_info.
Elem1No),
526 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
529 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
533 FaceElemTr.Loc1.Transf.ElementNo = face_info.
Elem1No;
534 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
535 Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
545 GetFaceTransformation(FaceNo, &FaceTransformation);
546 return &FaceTransformation;
553 GetFaceTransformation(EdgeNo, EdTr);
558 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
567 GetEdgeVertices(EdgeNo, v);
570 for (
int i = 0; i < spaceDim; i++)
572 for (
int j = 0; j < nv; j++)
574 pm(i, j) = vertices[v[j]](i);
577 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
581 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
585 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
586 int n = vdofs.
Size()/spaceDim;
588 for (
int i = 0; i < spaceDim; i++)
590 for (
int j = 0; j < n; j++)
592 pm(i, j) = (*Nodes)(vdofs[n*i+j]);
595 EdTr->
SetFE(edge_el);
599 MFEM_ABORT(
"Not implemented.");
607 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
608 return &EdgeTransformation;
612 void Mesh::GetLocalPtToSegTransformation(
627 void Mesh::GetLocalSegToTriTransformation(
635 tv = tri_t::Edges[i/64];
636 so = seg_t::Orient[i%64];
639 for (
int j = 0; j < 2; j++)
641 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
642 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
647 void Mesh::GetLocalSegToQuadTransformation(
655 qv = quad_t::Edges[i/64];
656 so = seg_t::Orient[i%64];
659 for (
int j = 0; j < 2; j++)
661 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
662 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
667 void Mesh::GetLocalTriToTetTransformation(
674 const int *tv = tet_t::FaceVert[i/64];
677 const int *to = tri_t::Orient[i%64];
681 for (
int j = 0; j < 3; j++)
684 locpm(0, j) = vert.
x;
685 locpm(1, j) = vert.
y;
686 locpm(2, j) = vert.
z;
691 void Mesh::GetLocalTriToWdgTransformation(
698 MFEM_VERIFY(i < 128,
"Local face index " << i/64
699 <<
" is not a triangular face of a wedge.");
700 const int *pv = pri_t::FaceVert[i/64];
703 const int *to = tri_t::Orient[i%64];
707 for (
int j = 0; j < 3; j++)
710 locpm(0, j) = vert.
x;
711 locpm(1, j) = vert.
y;
712 locpm(2, j) = vert.
z;
717 void Mesh::GetLocalQuadToHexTransformation(
724 const int *hv = hex_t::FaceVert[i/64];
726 const int *qo = quad_t::Orient[i%64];
729 for (
int j = 0; j < 4; j++)
732 locpm(0, j) = vert.
x;
733 locpm(1, j) = vert.
y;
734 locpm(2, j) = vert.
z;
739 void Mesh::GetLocalQuadToWdgTransformation(
746 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
747 <<
" is not a quadrilateral face of a wedge.");
748 const int *pv = pri_t::FaceVert[i/64];
750 const int *qo = quad_t::Orient[i%64];
753 for (
int j = 0; j < 4; j++)
756 locpm(0, j) = vert.
x;
757 locpm(1, j) = vert.
y;
758 locpm(2, j) = vert.
z;
766 for (
int i = 0; i < geom_factors.Size(); i++)
778 geom_factors.Append(gf);
786 for (
int i = 0; i < face_geom_factors.Size(); i++)
799 face_geom_factors.Append(gf);
803 void Mesh::DeleteGeometricFactors()
805 for (
int i = 0; i < geom_factors.Size(); i++)
807 delete geom_factors[i];
809 geom_factors.SetSize(0);
810 for (
int i = 0; i < face_geom_factors.Size(); i++)
812 delete face_geom_factors[i];
814 face_geom_factors.SetSize(0);
817 void Mesh::GetLocalFaceTransformation(
823 GetLocalPtToSegTransformation(Transf, info);
826 case Element::SEGMENT:
827 if (elem_type == Element::TRIANGLE)
829 GetLocalSegToTriTransformation(Transf, info);
833 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
834 GetLocalSegToQuadTransformation(Transf, info);
838 case Element::TRIANGLE:
839 if (elem_type == Element::TETRAHEDRON)
841 GetLocalTriToTetTransformation(Transf, info);
845 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
846 GetLocalTriToWdgTransformation(Transf, info);
850 case Element::QUADRILATERAL:
851 if (elem_type == Element::HEXAHEDRON)
853 GetLocalQuadToHexTransformation(Transf, info);
857 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
858 GetLocalQuadToWdgTransformation(Transf, info);
867 FaceInfo &face_info = faces_info[FaceNo];
869 FaceElemTr.Elem1 = NULL;
870 FaceElemTr.Elem2 = NULL;
876 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
877 FaceElemTr.Elem1 = &Transformation;
883 FaceElemTr.Elem2No = face_info.
Elem2No;
884 if ((mask & 2) && FaceElemTr.Elem2No >= 0)
887 if (NURBSext && (mask & 1)) { MFEM_ABORT(
"NURBS mesh not supported!"); }
889 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
890 FaceElemTr.Elem2 = &Transformation2;
894 FaceElemTr.FaceGeom = GetFaceGeometryType(FaceNo);
895 FaceElemTr.Face = (mask & 16) ? GetFaceTransformation(FaceNo) : NULL;
898 int face_type = GetFaceElementType(FaceNo);
901 int elem_type = GetElementType(face_info.
Elem1No);
902 GetLocalFaceTransformation(face_type, elem_type,
903 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
905 if ((mask & 8) && FaceElemTr.Elem2No >= 0)
907 int elem_type = GetElementType(face_info.
Elem2No);
908 GetLocalFaceTransformation(face_type, elem_type,
909 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
912 if (Nonconforming() && IsSlaveFace(face_info))
914 ApplyLocalSlaveTransformation(FaceElemTr.Loc2.Transf, face_info);
916 if (face_type == Element::SEGMENT)
919 DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
920 std::swap(pm(0,0), pm(0,1));
921 std::swap(pm(1,0), pm(1,1));
931 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
937 #ifdef MFEM_THREAD_SAFE
942 MFEM_ASSERT(fi.
NCFace >= 0,
"");
954 fn = be_to_face[BdrElemNo];
958 fn = be_to_edge[BdrElemNo];
962 fn = boundary[BdrElemNo]->GetVertices()[0];
965 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
969 tr = GetFaceElementTransformations(fn);
974 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const
976 *Elem1 = faces_info[Face].
Elem1No;
977 *Elem2 = faces_info[Face].Elem2No;
980 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const
982 *Inf1 = faces_info[Face].Elem1Inf;
983 *Inf2 = faces_info[Face].Elem2Inf;
988 return (Dim == 1) ? Geometry::POINT : faces[Face]->GetGeometryType();
993 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
1001 NumOfElements = NumOfBdrElements = 0;
1002 NumOfEdges = NumOfFaces = 0;
1003 nbInteriorFaces = -1;
1004 nbBoundaryFaces = -1;
1005 meshgen = mesh_geoms = 0;
1011 last_operation = Mesh::NONE;
1014 void Mesh::InitTables()
1017 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
1020 void Mesh::SetEmpty()
1026 void Mesh::DestroyTables()
1031 DeleteGeometricFactors();
1042 void Mesh::DestroyPointers()
1044 if (own_nodes) {
delete Nodes; }
1050 for (
int i = 0; i < NumOfElements; i++)
1052 FreeElement(elements[i]);
1055 for (
int i = 0; i < NumOfBdrElements; i++)
1057 FreeElement(boundary[i]);
1060 for (
int i = 0; i < faces.Size(); i++)
1062 FreeElement(faces[i]);
1068 void Mesh::Destroy()
1072 elements.DeleteAll();
1073 vertices.DeleteAll();
1074 boundary.DeleteAll();
1076 faces_info.DeleteAll();
1077 nc_faces_info.DeleteAll();
1078 be_to_edge.DeleteAll();
1079 be_to_face.DeleteAll();
1087 CoarseFineTr.Clear();
1089 #ifdef MFEM_USE_MEMALLOC
1093 attributes.DeleteAll();
1094 bdr_attributes.DeleteAll();
1097 void Mesh::ResetLazyData()
1099 delete el_to_el; el_to_el = NULL;
1100 delete face_edge; face_edge = NULL;
1101 delete edge_vertex; edge_vertex = NULL;
1102 DeleteGeometricFactors();
1103 nbInteriorFaces = -1;
1104 nbBoundaryFaces = -1;
1107 void Mesh::SetAttributes()
1112 for (
int i = 0; i < attribs.
Size(); i++)
1114 attribs[i] = GetBdrAttribute(i);
1118 attribs.
Copy(bdr_attributes);
1119 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1121 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1125 for (
int i = 0; i < attribs.
Size(); i++)
1127 attribs[i] = GetAttribute(i);
1131 attribs.
Copy(attributes);
1132 if (attributes.Size() > 0 && attributes[0] <= 0)
1134 MFEM_WARNING(
"Non-positive attributes in the domain!");
1138 void Mesh::InitMesh(
int _Dim,
int _spaceDim,
int NVert,
int NElem,
int NBdrElem)
1143 spaceDim = _spaceDim;
1146 vertices.SetSize(NVert);
1149 elements.SetSize(NElem);
1151 NumOfBdrElements = 0;
1152 boundary.SetSize(NBdrElem);
1155 void Mesh::AddVertex(
const double *x)
1157 double *y = vertices[NumOfVertices]();
1159 for (
int i = 0; i < spaceDim; i++)
1166 void Mesh::AddSegment(
const int *vi,
int attr)
1168 elements[NumOfElements++] =
new Segment(vi, attr);
1171 void Mesh::AddTri(
const int *vi,
int attr)
1173 elements[NumOfElements++] =
new Triangle(vi, attr);
1176 void Mesh::AddTriangle(
const int *vi,
int attr)
1178 elements[NumOfElements++] =
new Triangle(vi, attr);
1181 void Mesh::AddQuad(
const int *vi,
int attr)
1186 void Mesh::AddTet(
const int *vi,
int attr)
1188 #ifdef MFEM_USE_MEMALLOC
1190 tet = TetMemory.Alloc();
1193 elements[NumOfElements++] = tet;
1195 elements[NumOfElements++] =
new Tetrahedron(vi, attr);
1199 void Mesh::AddWedge(
const int *vi,
int attr)
1201 elements[NumOfElements++] =
new Wedge(vi, attr);
1204 void Mesh::AddHex(
const int *vi,
int attr)
1206 elements[NumOfElements++] =
new Hexahedron(vi, attr);
1209 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1211 static const int hex_to_tet[6][4] =
1213 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1214 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1218 for (
int i = 0; i < 6; i++)
1220 for (
int j = 0; j < 4; j++)
1222 ti[j] = vi[hex_to_tet[i][j]];
1228 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1230 static const int hex_to_wdg[2][6] =
1232 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1236 for (
int i = 0; i < 2; i++)
1238 for (
int j = 0; j < 6; j++)
1240 ti[j] = vi[hex_to_wdg[i][j]];
1246 void Mesh::AddBdrSegment(
const int *vi,
int attr)
1248 boundary[NumOfBdrElements++] =
new Segment(vi, attr);
1251 void Mesh::AddBdrTriangle(
const int *vi,
int attr)
1253 boundary[NumOfBdrElements++] =
new Triangle(vi, attr);
1256 void Mesh::AddBdrQuad(
const int *vi,
int attr)
1261 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1263 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1266 for (
int i = 0; i < 2; i++)
1268 for (
int j = 0; j < 3; j++)
1270 ti[j] = vi[quad_to_tri[i][j]];
1272 AddBdrTriangle(ti, attr);
1276 void Mesh::GenerateBoundaryElements()
1279 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1283 for (i = 0; i < boundary.Size(); i++)
1285 FreeElement(boundary[i]);
1295 NumOfBdrElements = 0;
1296 for (i = 0; i < faces_info.Size(); i++)
1298 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1301 boundary.SetSize(NumOfBdrElements);
1302 be2face.
SetSize(NumOfBdrElements);
1303 for (j = i = 0; i < faces_info.Size(); i++)
1305 if (faces_info[i].Elem2No < 0)
1307 boundary[j] = faces[i]->Duplicate(
this);
1314 void Mesh::FinalizeCheck()
1316 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1317 vertices.Size() == 0,
1318 "incorrect number of vertices: preallocated: " << vertices.Size()
1319 <<
", actually added: " << NumOfVertices);
1320 MFEM_VERIFY(elements.Size() == NumOfElements,
1321 "incorrect number of elements: preallocated: " << elements.Size()
1322 <<
", actually added: " << NumOfElements);
1323 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1324 "incorrect number of boundary elements: preallocated: "
1325 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1328 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1331 CheckElementOrientation(fix_orientation);
1335 MarkTriMeshForRefinement();
1340 el_to_edge =
new Table;
1341 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1343 CheckBdrElementOrientation();
1357 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1358 bool fix_orientation)
1361 if (fix_orientation)
1363 CheckElementOrientation(fix_orientation);
1368 el_to_edge =
new Table;
1369 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1371 CheckBdrElementOrientation();
1386 #ifdef MFEM_USE_GECKO
1388 int iterations,
int window,
1389 int period,
int seed)
1393 Gecko::Functional *functional =
1394 new Gecko::FunctionalGeometric();
1397 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1404 const Table &my_el_to_el = ElementToElementTable();
1405 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1407 const int *neighid = my_el_to_el.
GetRow(elemid);
1408 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
1410 graph.insert(elemid + 1, neighid[i] + 1);
1415 graph.order(functional, iterations, window, period, seed);
1418 Gecko::Node::Index NE = GetNE();
1419 for (Gecko::Node::Index gnodeid = 1; gnodeid <= NE; ++gnodeid)
1421 ordering[gnodeid - 1] = graph.rank(gnodeid);
1436 HilbertCmp(
int coord,
bool dir,
const Array<double> &points,
double mid)
1437 : coord(coord), dir(dir), points(points), mid(mid) {}
1439 bool operator()(
int i)
const
1441 return (points[3*i + coord] < mid) != dir;
1445 static void HilbertSort2D(
int coord1,
1448 const Array<double> &points,
int *beg,
int *end,
1449 double xmin,
double ymin,
double xmax,
double ymax)
1451 if (end - beg <= 1) {
return; }
1453 double xmid = (xmin + xmax)*0.5;
1454 double ymid = (ymin + ymax)*0.5;
1456 int coord2 = (coord1 + 1) % 2;
1459 int *p0 = beg, *p4 = end;
1460 int *p2 = std::partition(p0, p4, HilbertCmp(coord1, dir1, points, xmid));
1461 int *p1 = std::partition(p0, p2, HilbertCmp(coord2, dir2, points, ymid));
1462 int *p3 = std::partition(p2, p4, HilbertCmp(coord2, !dir2, points, ymid));
1466 HilbertSort2D(coord2, dir2, dir1, points, p0, p1,
1467 ymin, xmin, ymid, xmid);
1469 if (p1 != p0 || p2 != p4)
1471 HilbertSort2D(coord1, dir1, dir2, points, p1, p2,
1472 xmin, ymid, xmid, ymax);
1474 if (p2 != p0 || p3 != p4)
1476 HilbertSort2D(coord1, dir1, dir2, points, p2, p3,
1477 xmid, ymid, xmax, ymax);
1481 HilbertSort2D(coord2, !dir2, !dir1, points, p3, p4,
1482 ymid, xmax, ymin, xmid);
1486 static void HilbertSort3D(
int coord1,
bool dir1,
bool dir2,
bool dir3,
1487 const Array<double> &points,
int *beg,
int *end,
1488 double xmin,
double ymin,
double zmin,
1489 double xmax,
double ymax,
double zmax)
1491 if (end - beg <= 1) {
return; }
1493 double xmid = (xmin + xmax)*0.5;
1494 double ymid = (ymin + ymax)*0.5;
1495 double zmid = (zmin + zmax)*0.5;
1497 int coord2 = (coord1 + 1) % 3;
1498 int coord3 = (coord1 + 2) % 3;
1501 int *p0 = beg, *p8 = end;
1502 int *p4 = std::partition(p0, p8, HilbertCmp(coord1, dir1, points, xmid));
1503 int *p2 = std::partition(p0, p4, HilbertCmp(coord2, dir2, points, ymid));
1504 int *p6 = std::partition(p4, p8, HilbertCmp(coord2, !dir2, points, ymid));
1505 int *p1 = std::partition(p0, p2, HilbertCmp(coord3, dir3, points, zmid));
1506 int *p3 = std::partition(p2, p4, HilbertCmp(coord3, !dir3, points, zmid));
1507 int *p5 = std::partition(p4, p6, HilbertCmp(coord3, dir3, points, zmid));
1508 int *p7 = std::partition(p6, p8, HilbertCmp(coord3, !dir3, points, zmid));
1512 HilbertSort3D(coord3, dir3, dir1, dir2, points, p0, p1,
1513 zmin, xmin, ymin, zmid, xmid, ymid);
1515 if (p1 != p0 || p2 != p8)
1517 HilbertSort3D(coord2, dir2, dir3, dir1, points, p1, p2,
1518 ymin, zmid, xmin, ymid, zmax, xmid);
1520 if (p2 != p0 || p3 != p8)
1522 HilbertSort3D(coord2, dir2, dir3, dir1, points, p2, p3,
1523 ymid, zmid, xmin, ymax, zmax, xmid);
1525 if (p3 != p0 || p4 != p8)
1527 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p3, p4,
1528 xmin, ymax, zmid, xmid, ymid, zmin);
1530 if (p4 != p0 || p5 != p8)
1532 HilbertSort3D(coord1, dir1, !dir2, !dir3, points, p4, p5,
1533 xmid, ymax, zmid, xmax, ymid, zmin);
1535 if (p5 != p0 || p6 != p8)
1537 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p5, p6,
1538 ymax, zmid, xmax, ymid, zmax, xmid);
1540 if (p6 != p0 || p7 != p8)
1542 HilbertSort3D(coord2, !dir2, dir3, !dir1, points, p6, p7,
1543 ymid, zmid, xmax, ymin, zmax, xmid);
1547 HilbertSort3D(coord3, !dir3, !dir1, dir2, points, p7, p8,
1548 zmid, xmax, ymin, zmin, xmid, ymid);
1554 MFEM_VERIFY(spaceDim <= 3,
"");
1557 GetBoundingBox(min, max);
1562 if (spaceDim < 3) { points = 0.0; }
1565 for (
int i = 0; i < GetNE(); i++)
1567 GetElementCenter(i, center);
1568 for (
int j = 0; j < spaceDim; j++)
1570 points[3*i + j] = center(j);
1577 indices.
Sort([&](
int a,
int b)
1578 {
return points[3*
a] < points[3*
b]; });
1580 else if (spaceDim == 2)
1583 HilbertSort2D(0,
false,
false,
1584 points, indices.
begin(), indices.
end(),
1585 min(0), min(1), max(0), max(1));
1590 HilbertSort3D(0,
false,
false,
false,
1591 points, indices.
begin(), indices.
end(),
1592 min(0), min(1), min(2), max(0), max(1), max(2));
1597 for (
int i = 0; i < GetNE(); i++)
1599 ordering[indices[i]] = i;
1604 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
1608 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
1613 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
1617 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
1647 old_elem_node_vals.SetSize(GetNE());
1648 nodes_fes = Nodes->FESpace();
1651 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1653 nodes_fes->GetElementVDofs(old_elid, old_dofs);
1655 old_elem_node_vals[old_elid] =
new Vector(vals);
1661 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
1663 int new_elid = ordering[old_elid];
1664 new_elements[new_elid] = elements[old_elid];
1669 if (reorder_vertices)
1674 vertex_ordering = -1;
1676 int new_vertex_ind = 0;
1677 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1679 int *elem_vert = elements[new_elid]->GetVertices();
1680 int nv = elements[new_elid]->GetNVertices();
1681 for (
int vi = 0; vi < nv; ++vi)
1683 int old_vertex_ind = elem_vert[vi];
1684 if (vertex_ordering[old_vertex_ind] == -1)
1686 vertex_ordering[old_vertex_ind] = new_vertex_ind;
1687 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
1697 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1699 int *elem_vert = elements[new_elid]->GetVertices();
1700 int nv = elements[new_elid]->GetNVertices();
1701 for (
int vi = 0; vi < nv; ++vi)
1703 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
1708 for (
int belid = 0; belid < GetNBE(); ++belid)
1710 int *be_vert = boundary[belid]->GetVertices();
1711 int nv = boundary[belid]->GetNVertices();
1712 for (
int vi = 0; vi < nv; ++vi)
1714 be_vert[vi] = vertex_ordering[be_vert[vi]];
1725 el_to_edge =
new Table;
1726 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1731 GetElementToFaceTable();
1741 last_operation = Mesh::NONE;
1742 nodes_fes->Update(
false);
1745 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1747 int new_elid = ordering[old_elid];
1748 nodes_fes->GetElementVDofs(new_elid, new_dofs);
1749 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
1750 delete old_elem_node_vals[old_elid];
1756 void Mesh::MarkForRefinement()
1762 MarkTriMeshForRefinement();
1766 DSTable v_to_v(NumOfVertices);
1767 GetVertexToVertexTable(v_to_v);
1768 MarkTetMeshForRefinement(v_to_v);
1773 void Mesh::MarkTriMeshForRefinement()
1778 for (
int i = 0; i < NumOfElements; i++)
1780 if (elements[i]->GetType() == Element::TRIANGLE)
1782 GetPointMatrix(i, pmat);
1783 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
1794 for (
int i = 0; i < NumOfVertices; i++)
1799 length_idx[j].one = GetLength(i, it.Column());
1800 length_idx[j].two = j;
1807 for (
int i = 0; i < NumOfEdges; i++)
1809 order[length_idx[i].two] = i;
1813 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
1818 GetEdgeOrdering(v_to_v, order);
1820 for (
int i = 0; i < NumOfElements; i++)
1822 if (elements[i]->GetType() == Element::TETRAHEDRON)
1824 elements[i]->MarkEdge(v_to_v, order);
1827 for (
int i = 0; i < NumOfBdrElements; i++)
1829 if (boundary[i]->GetType() == Element::TRIANGLE)
1831 boundary[i]->MarkEdge(v_to_v, order);
1838 if (*old_v_to_v && *old_elem_vert)
1845 if (*old_v_to_v == NULL)
1847 bool need_v_to_v =
false;
1849 for (
int i = 0; i < GetNEdges(); i++)
1855 if (dofs.
Size() > 0)
1863 *old_v_to_v =
new DSTable(NumOfVertices);
1864 GetVertexToVertexTable(*(*old_v_to_v));
1867 if (*old_elem_vert == NULL)
1869 bool need_elem_vert =
false;
1871 for (
int i = 0; i < GetNE(); i++)
1877 if (dofs.
Size() > 1)
1879 need_elem_vert =
true;
1885 *old_elem_vert =
new Table;
1886 (*old_elem_vert)->
MakeI(GetNE());
1887 for (
int i = 0; i < GetNE(); i++)
1889 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
1891 (*old_elem_vert)->MakeJ();
1892 for (
int i = 0; i < GetNE(); i++)
1894 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
1895 elements[i]->GetNVertices());
1897 (*old_elem_vert)->ShiftUpI();
1910 const int num_edge_dofs = old_dofs.
Size();
1913 const Vector onodes = *Nodes;
1917 int offset = NumOfVertices * old_dofs.
Size();
1921 if (num_edge_dofs > 0)
1923 DSTable new_v_to_v(NumOfVertices);
1924 GetVertexToVertexTable(new_v_to_v);
1926 for (
int i = 0; i < NumOfVertices; i++)
1930 const int old_i = (*old_v_to_v)(i, it.Column());
1931 const int new_i = it.Index();
1932 if (new_i == old_i) {
continue; }
1934 old_dofs.
SetSize(num_edge_dofs);
1935 new_dofs.
SetSize(num_edge_dofs);
1936 for (
int j = 0; j < num_edge_dofs; j++)
1938 old_dofs[j] = offset + old_i * num_edge_dofs + j;
1939 new_dofs[j] = offset + new_i * num_edge_dofs + j;
1943 for (
int j = 0; j < old_dofs.
Size(); j++)
1945 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1949 offset += NumOfEdges * num_edge_dofs;
1957 Table old_face_vertex;
1958 old_face_vertex.
MakeI(NumOfFaces);
1959 for (
int i = 0; i < NumOfFaces; i++)
1963 old_face_vertex.
MakeJ();
1964 for (
int i = 0; i < NumOfFaces; i++)
1966 faces[i]->GetNVertices());
1970 STable3D *faces_tbl = GetElementToFaceTable(1);
1976 for (
int i = 0; i < NumOfFaces; i++)
1978 const int *old_v = old_face_vertex.
GetRow(i);
1980 switch (old_face_vertex.
RowSize(i))
1983 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
1987 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
1991 new_fdofs[new_i+1] = old_dofs.
Size();
1996 for (
int i = 0; i < NumOfFaces; i++)
1998 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
2001 switch (old_face_vertex.
RowSize(i))
2004 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
2005 new_v = faces[new_i]->GetVertices();
2006 new_or = GetTriOrientation(old_v, new_v);
2011 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
2012 new_v = faces[new_i]->GetVertices();
2013 new_or = GetQuadOrientation(old_v, new_v);
2020 for (
int j = 0; j < old_dofs.
Size(); j++)
2023 const int old_j = dof_ord[j];
2024 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
2028 for (
int j = 0; j < old_dofs.
Size(); j++)
2030 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2052 for (
int i = 0; i < GetNE(); i++)
2054 const int *old_v = old_elem_vert->
GetRow(i);
2055 const int *new_v = elements[i]->GetVertices();
2061 case Geometry::SEGMENT:
2062 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
2064 case Geometry::TRIANGLE:
2065 new_or = GetTriOrientation(old_v, new_v);
2067 case Geometry::SQUARE:
2068 new_or = GetQuadOrientation(old_v, new_v);
2072 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
2073 <<
" FE collection) are not supported yet!");
2077 MFEM_VERIFY(dof_ord != NULL,
2078 "FE collection '" << fec->
Name()
2079 <<
"' does not define reordering for "
2080 << Geometry::Name[
geom] <<
" elements!");
2083 for (
int j = 0; j < new_dofs.
Size(); j++)
2086 const int old_j = dof_ord[j];
2087 new_dofs[old_j] = offset + j;
2089 offset += new_dofs.
Size();
2092 for (
int j = 0; j < old_dofs.
Size(); j++)
2094 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
2106 GetElementToFaceTable();
2109 CheckBdrElementOrientation();
2114 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2119 CheckBdrElementOrientation();
2124 last_operation = Mesh::NONE;
2129 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
2132 CheckElementOrientation(fix_orientation);
2134 if (NumOfBdrElements == 0)
2136 GetElementToFaceTable();
2138 GenerateBoundaryElements();
2143 DSTable v_to_v(NumOfVertices);
2144 GetVertexToVertexTable(v_to_v);
2145 MarkTetMeshForRefinement(v_to_v);
2148 GetElementToFaceTable();
2151 CheckBdrElementOrientation();
2153 if (generate_edges == 1)
2155 el_to_edge =
new Table;
2156 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2170 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
2171 bool fix_orientation)
2174 CheckElementOrientation(fix_orientation);
2176 if (NumOfBdrElements == 0)
2178 GetElementToFaceTable();
2180 GenerateBoundaryElements();
2183 GetElementToFaceTable();
2186 CheckBdrElementOrientation();
2188 if (generate_edges == 1)
2190 el_to_edge =
new Table;
2191 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2205 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
2208 CheckElementOrientation(fix_orientation);
2210 GetElementToFaceTable();
2213 if (NumOfBdrElements == 0)
2215 GenerateBoundaryElements();
2218 CheckBdrElementOrientation();
2222 el_to_edge =
new Table;
2223 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2235 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2239 Finalize(refine, fix_orientation);
2242 void Mesh::FinalizeTopology(
bool generate_bdr)
2254 bool generate_edges =
true;
2256 if (spaceDim == 0) { spaceDim = Dim; }
2257 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2265 GetElementToFaceTable();
2267 if (NumOfBdrElements == 0 && generate_bdr)
2269 GenerateBoundaryElements();
2270 GetElementToFaceTable();
2279 if (Dim > 1 && generate_edges)
2282 if (!el_to_edge) { el_to_edge =
new Table; }
2283 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2287 if (NumOfBdrElements == 0 && generate_bdr)
2289 GenerateBoundaryElements();
2306 ncmesh->OnMeshUpdated(
this);
2309 GenerateNCFaceInfo();
2316 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2318 if (NURBSext || ncmesh)
2320 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
2321 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
2330 const bool check_orientation =
true;
2331 const bool curved = (Nodes != NULL);
2332 const bool may_change_topology =
2333 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
2334 ( check_orientation && fix_orientation &&
2335 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
2338 Table *old_elem_vert = NULL;
2340 if (curved && may_change_topology)
2342 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
2345 if (check_orientation)
2348 CheckElementOrientation(fix_orientation);
2352 MarkForRefinement();
2355 if (may_change_topology)
2359 DoNodeReorder(old_v_to_v, old_elem_vert);
2360 delete old_elem_vert;
2373 CheckBdrElementOrientation();
2378 if (Dim >= 2 && Dim == spaceDim)
2380 const int num_faces = GetNumFaces();
2381 for (
int i = 0; i < num_faces; i++)
2383 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
2384 faces_info[i].Elem2Inf%2 != 0,
"invalid mesh topology");
2391 double sx,
double sy,
double sz,
bool sfc_ordering)
2395 int NVert, NElem, NBdrElem;
2397 NVert = (nx+1) * (ny+1) * (nz+1);
2398 NElem = nx * ny * nz;
2399 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
2400 if (type == Element::TETRAHEDRON)
2405 else if (type == Element::WEDGE)
2408 NBdrElem += 2*nx*ny;
2411 InitMesh(3, 3, NVert, NElem, NBdrElem);
2417 for (z = 0; z <= nz; z++)
2419 coord[2] = ((double) z / nz) * sz;
2420 for (y = 0; y <= ny; y++)
2422 coord[1] = ((double) y / ny) * sy;
2423 for (x = 0; x <= nx; x++)
2425 coord[0] = ((double) x / nx) * sx;
2431 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
2434 if (sfc_ordering && type == Element::HEXAHEDRON)
2437 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
2438 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
2440 for (
int k = 0; k < nx*ny*nz; k++)
2446 ind[0] = VTX(x , y , z );
2447 ind[1] = VTX(x+1, y , z );
2448 ind[2] = VTX(x+1, y+1, z );
2449 ind[3] = VTX(x , y+1, z );
2450 ind[4] = VTX(x , y , z+1);
2451 ind[5] = VTX(x+1, y , z+1);
2452 ind[6] = VTX(x+1, y+1, z+1);
2453 ind[7] = VTX(x , y+1, z+1);
2460 for (z = 0; z < nz; z++)
2462 for (y = 0; y < ny; y++)
2464 for (x = 0; x < nx; x++)
2466 ind[0] = VTX(x , y , z );
2467 ind[1] = VTX(x+1, y , z );
2468 ind[2] = VTX(x+1, y+1, z );
2469 ind[3] = VTX(x , y+1, z );
2470 ind[4] = VTX(x , y , z+1);
2471 ind[5] = VTX(x+1, y , z+1);
2472 ind[6] = VTX(x+1, y+1, z+1);
2473 ind[7] = VTX( x, y+1, z+1);
2474 if (type == Element::TETRAHEDRON)
2476 AddHexAsTets(ind, 1);
2478 else if (type == Element::WEDGE)
2480 AddHexAsWedges(ind, 1);
2493 for (y = 0; y < ny; y++)
2495 for (x = 0; x < nx; x++)
2497 ind[0] = VTX(x , y , 0);
2498 ind[1] = VTX(x , y+1, 0);
2499 ind[2] = VTX(x+1, y+1, 0);
2500 ind[3] = VTX(x+1, y , 0);
2501 if (type == Element::TETRAHEDRON)
2503 AddBdrQuadAsTriangles(ind, 1);
2505 else if (type == Element::WEDGE)
2507 AddBdrQuadAsTriangles(ind, 1);
2516 for (y = 0; y < ny; y++)
2518 for (x = 0; x < nx; x++)
2520 ind[0] = VTX(x , y , nz);
2521 ind[1] = VTX(x+1, y , nz);
2522 ind[2] = VTX(x+1, y+1, nz);
2523 ind[3] = VTX(x , y+1, nz);
2524 if (type == Element::TETRAHEDRON)
2526 AddBdrQuadAsTriangles(ind, 6);
2528 else if (type == Element::WEDGE)
2530 AddBdrQuadAsTriangles(ind, 1);
2539 for (z = 0; z < nz; z++)
2541 for (y = 0; y < ny; y++)
2543 ind[0] = VTX(0 , y , z );
2544 ind[1] = VTX(0 , y , z+1);
2545 ind[2] = VTX(0 , y+1, z+1);
2546 ind[3] = VTX(0 , y+1, z );
2547 if (type == Element::TETRAHEDRON)
2549 AddBdrQuadAsTriangles(ind, 5);
2558 for (z = 0; z < nz; z++)
2560 for (y = 0; y < ny; y++)
2562 ind[0] = VTX(nx, y , z );
2563 ind[1] = VTX(nx, y+1, z );
2564 ind[2] = VTX(nx, y+1, z+1);
2565 ind[3] = VTX(nx, y , z+1);
2566 if (type == Element::TETRAHEDRON)
2568 AddBdrQuadAsTriangles(ind, 3);
2577 for (x = 0; x < nx; x++)
2579 for (z = 0; z < nz; z++)
2581 ind[0] = VTX(x , 0, z );
2582 ind[1] = VTX(x+1, 0, z );
2583 ind[2] = VTX(x+1, 0, z+1);
2584 ind[3] = VTX(x , 0, z+1);
2585 if (type == Element::TETRAHEDRON)
2587 AddBdrQuadAsTriangles(ind, 2);
2596 for (x = 0; x < nx; x++)
2598 for (z = 0; z < nz; z++)
2600 ind[0] = VTX(x , ny, z );
2601 ind[1] = VTX(x , ny, z+1);
2602 ind[2] = VTX(x+1, ny, z+1);
2603 ind[3] = VTX(x+1, ny, z );
2604 if (type == Element::TETRAHEDRON)
2606 AddBdrQuadAsTriangles(ind, 4);
2618 ofstream test_stream(
"debug.mesh");
2620 test_stream.close();
2629 double sx,
double sy,
2630 bool generate_edges,
bool sfc_ordering)
2639 if (type == Element::QUADRILATERAL)
2641 NumOfVertices = (nx+1) * (ny+1);
2642 NumOfElements = nx * ny;
2643 NumOfBdrElements = 2 * nx + 2 * ny;
2645 vertices.SetSize(NumOfVertices);
2646 elements.SetSize(NumOfElements);
2647 boundary.SetSize(NumOfBdrElements);
2654 for (j = 0; j < ny+1; j++)
2656 cy = ((double) j / ny) * sy;
2657 for (i = 0; i < nx+1; i++)
2659 cx = ((double) i / nx) * sx;
2660 vertices[k](0) = cx;
2661 vertices[k](1) = cy;
2670 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
2671 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
2673 for (k = 0; k < nx*ny; k++)
2677 ind[0] = i + j*(nx+1);
2678 ind[1] = i + 1 +j*(nx+1);
2679 ind[2] = i + 1 + (j+1)*(nx+1);
2680 ind[3] = i + (j+1)*(nx+1);
2687 for (j = 0; j < ny; j++)
2689 for (i = 0; i < nx; i++)
2691 ind[0] = i + j*(nx+1);
2692 ind[1] = i + 1 +j*(nx+1);
2693 ind[2] = i + 1 + (j+1)*(nx+1);
2694 ind[3] = i + (j+1)*(nx+1);
2703 for (i = 0; i < nx; i++)
2705 boundary[i] =
new Segment(i, i+1, 1);
2706 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
2709 for (j = 0; j < ny; j++)
2711 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
2712 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
2716 else if (type == Element::TRIANGLE)
2718 NumOfVertices = (nx+1) * (ny+1);
2719 NumOfElements = 2 * nx * ny;
2720 NumOfBdrElements = 2 * nx + 2 * ny;
2722 vertices.SetSize(NumOfVertices);
2723 elements.SetSize(NumOfElements);
2724 boundary.SetSize(NumOfBdrElements);
2731 for (j = 0; j < ny+1; j++)
2733 cy = ((double) j / ny) * sy;
2734 for (i = 0; i < nx+1; i++)
2736 cx = ((double) i / nx) * sx;
2737 vertices[k](0) = cx;
2738 vertices[k](1) = cy;
2745 for (j = 0; j < ny; j++)
2747 for (i = 0; i < nx; i++)
2749 ind[0] = i + j*(nx+1);
2750 ind[1] = i + 1 + (j+1)*(nx+1);
2751 ind[2] = i + (j+1)*(nx+1);
2754 ind[1] = i + 1 + j*(nx+1);
2755 ind[2] = i + 1 + (j+1)*(nx+1);
2763 for (i = 0; i < nx; i++)
2765 boundary[i] =
new Segment(i, i+1, 1);
2766 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
2769 for (j = 0; j < ny; j++)
2771 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
2772 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
2779 MFEM_ABORT(
"Unsupported element type.");
2783 CheckElementOrientation();
2785 if (generate_edges == 1)
2787 el_to_edge =
new Table;
2788 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2790 CheckBdrElementOrientation();
2799 attributes.Append(1);
2800 bdr_attributes.Append(1); bdr_attributes.Append(2);
2801 bdr_attributes.Append(3); bdr_attributes.Append(4);
2806 void Mesh::Make1D(
int n,
double sx)
2815 NumOfVertices = n + 1;
2817 NumOfBdrElements = 2;
2818 vertices.SetSize(NumOfVertices);
2819 elements.SetSize(NumOfElements);
2820 boundary.SetSize(NumOfBdrElements);
2823 for (j = 0; j < n+1; j++)
2825 vertices[j](0) = ((
double) j / n) * sx;
2829 for (j = 0; j < n; j++)
2831 elements[j] =
new Segment(j, j+1, 1);
2836 boundary[0] =
new Point(ind, 1);
2838 boundary[1] =
new Point(ind, 2);
2846 attributes.Append(1);
2847 bdr_attributes.Append(1); bdr_attributes.Append(2);
2850 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
2868 last_operation = Mesh::NONE;
2871 elements.SetSize(NumOfElements);
2872 for (
int i = 0; i < NumOfElements; i++)
2874 elements[i] = mesh.
elements[i]->Duplicate(
this);
2881 boundary.SetSize(NumOfBdrElements);
2882 for (
int i = 0; i < NumOfBdrElements; i++)
2884 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
2903 faces.SetSize(mesh.
faces.Size());
2904 for (
int i = 0; i < faces.Size(); i++)
2907 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
2941 if (dynamic_cast<const ParMesh*>(&mesh))
2953 if (mesh.
Nodes && copy_nodes)
2958 FiniteElementCollection::New(fec->
Name());
2962 Nodes->MakeOwner(fec_copy);
2963 *Nodes = *mesh.
Nodes;
2973 Mesh::Mesh(
const char *filename,
int generate_edges,
int refine,
2974 bool fix_orientation)
2983 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
2987 Load(imesh, generate_edges, refine, fix_orientation);
2991 Mesh::Mesh(std::istream &input,
int generate_edges,
int refine,
2992 bool fix_orientation)
2995 Load(input, generate_edges, refine, fix_orientation);
2998 void Mesh::ChangeVertexDataOwnership(
double *vertex_data,
int len_vertex_data,
3003 MFEM_VERIFY(len_vertex_data >= NumOfVertices * 3,
3004 "Not enough vertices in external array : "
3005 "len_vertex_data = "<< len_vertex_data <<
", "
3006 "NumOfVertices * 3 = " << NumOfVertices * 3);
3008 if (vertex_data == (
double *)(vertices.GetData()))
3010 MFEM_ASSERT(!vertices.OwnsData(),
"invalid ownership");
3015 memcpy(vertex_data, vertices.GetData(),
3016 NumOfVertices * 3 *
sizeof(double));
3019 vertices.MakeRef(reinterpret_cast<Vertex*>(vertex_data), NumOfVertices);
3022 Mesh::Mesh(
double *_vertices,
int num_vertices,
3024 int *element_attributes,
int num_elements,
3026 int *boundary_attributes,
int num_boundary_elements,
3027 int dimension,
int space_dimension)
3029 if (space_dimension == -1)
3031 space_dimension = dimension;
3034 InitMesh(dimension, space_dimension, 0, num_elements,
3035 num_boundary_elements);
3037 int element_index_stride = Geometry::NumVerts[element_type];
3038 int boundary_index_stride = num_boundary_elements > 0 ?
3039 Geometry::NumVerts[boundary_type] : 0;
3042 vertices.MakeRef(reinterpret_cast<Vertex*>(_vertices), num_vertices);
3043 NumOfVertices = num_vertices;
3045 for (
int i = 0; i < num_elements; i++)
3047 elements[i] = NewElement(element_type);
3048 elements[i]->SetVertices(element_indices + i * element_index_stride);
3049 elements[i]->SetAttribute(element_attributes[i]);
3051 NumOfElements = num_elements;
3053 for (
int i = 0; i < num_boundary_elements; i++)
3055 boundary[i] = NewElement(boundary_type);
3056 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
3057 boundary[i]->SetAttribute(boundary_attributes[i]);
3059 NumOfBdrElements = num_boundary_elements;
3068 case Geometry::POINT:
return (
new Point);
3069 case Geometry::SEGMENT:
return (
new Segment);
3070 case Geometry::TRIANGLE:
return (
new Triangle);
3072 case Geometry::TETRAHEDRON:
3073 #ifdef MFEM_USE_MEMALLOC
3074 return TetMemory.Alloc();
3078 case Geometry::CUBE:
return (
new Hexahedron);
3079 case Geometry::PRISM:
return (
new Wedge);
3081 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
3087 Element *Mesh::ReadElementWithoutAttr(std::istream &input)
3093 el = NewElement(geom);
3094 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
3097 for (
int i = 0; i < nv; i++)
3105 void Mesh::PrintElementWithoutAttr(
const Element *el, std::ostream &
out)
3110 for (
int j = 0; j < nv; j++)
3123 el = ReadElementWithoutAttr(input);
3132 PrintElementWithoutAttr(el, out);
3135 void Mesh::SetMeshGen()
3137 meshgen = mesh_geoms = 0;
3138 for (
int i = 0; i < NumOfElements; i++)
3143 case Element::TETRAHEDRON:
3144 mesh_geoms |= (1 << Geometry::TETRAHEDRON);
3145 case Element::TRIANGLE:
3146 mesh_geoms |= (1 << Geometry::TRIANGLE);
3147 case Element::SEGMENT:
3148 mesh_geoms |= (1 << Geometry::SEGMENT);
3149 case Element::POINT:
3150 mesh_geoms |= (1 << Geometry::POINT);
3154 case Element::HEXAHEDRON:
3155 mesh_geoms |= (1 << Geometry::CUBE);
3156 case Element::QUADRILATERAL:
3157 mesh_geoms |= (1 << Geometry::SQUARE);
3158 mesh_geoms |= (1 << Geometry::SEGMENT);
3159 mesh_geoms |= (1 << Geometry::POINT);
3163 case Element::WEDGE:
3164 mesh_geoms |= (1 << Geometry::PRISM);
3165 mesh_geoms |= (1 << Geometry::SQUARE);
3166 mesh_geoms |= (1 << Geometry::TRIANGLE);
3167 mesh_geoms |= (1 << Geometry::SEGMENT);
3168 mesh_geoms |= (1 << Geometry::POINT);
3173 MFEM_ABORT(
"invalid element type: " << type);
3179 void Mesh::Loader(std::istream &input,
int generate_edges,
3180 std::string parse_tag)
3182 int curved = 0, read_gf = 1;
3183 bool finalize_topo =
true;
3187 MFEM_ABORT(
"Input stream is not open");
3194 getline(input, mesh_type);
3198 bool mfem_v10 = (mesh_type ==
"MFEM mesh v1.0");
3199 bool mfem_v11 = (mesh_type ==
"MFEM mesh v1.1");
3200 bool mfem_v12 = (mesh_type ==
"MFEM mesh v1.2");
3201 if (mfem_v10 || mfem_v11 || mfem_v12)
3207 if ( mfem_v12 && parse_tag.empty() )
3209 parse_tag =
"mfem_mesh_end";
3211 ReadMFEMMesh(input, mfem_v11, curved);
3213 else if (mesh_type ==
"linemesh")
3215 ReadLineMesh(input);
3217 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
3219 if (mesh_type ==
"curved_areamesh2")
3223 ReadNetgen2DMesh(input, curved);
3225 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
3227 ReadNetgen3DMesh(input);
3229 else if (mesh_type ==
"TrueGrid")
3231 ReadTrueGridMesh(input);
3233 else if (mesh_type ==
"# vtk DataFile Version 3.0" ||
3234 mesh_type ==
"# vtk DataFile Version 2.0")
3236 ReadVTKMesh(input, curved, read_gf, finalize_topo);
3238 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
3240 ReadNURBSMesh(input, curved, read_gf);
3242 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
3244 ReadInlineMesh(input, generate_edges);
3247 else if (mesh_type ==
"$MeshFormat")
3249 ReadGmshMesh(input);
3252 ((mesh_type.size() > 2 &&
3253 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
3254 (mesh_type.size() > 3 &&
3255 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
3260 #ifdef MFEM_USE_NETCDF
3261 ReadCubit(mesh_input->
filename.c_str(), curved, read_gf);
3263 MFEM_ABORT(
"NetCDF support requires configuration with"
3264 " MFEM_USE_NETCDF=YES");
3270 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
3271 " Use mfem::named_ifgzstream for input.");
3277 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
3305 if (curved && read_gf)
3309 spaceDim = Nodes->VectorDim();
3310 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
3312 for (
int i = 0; i < spaceDim; i++)
3315 Nodes->GetNodalValues(vert_val, i+1);
3316 for (
int j = 0; j < NumOfVertices; j++)
3318 vertices[j](i) = vert_val(j);
3331 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
3332 getline(input, line);
3338 if (line ==
"mfem_mesh_end") {
break; }
3340 while (line != parse_tag);
3346 Mesh::Mesh(
Mesh *mesh_array[],
int num_pieces)
3348 int i, j, ie, ib, iv, *v, nv;
3357 if (mesh_array[0]->NURBSext)
3362 NumOfVertices = NURBSext->GetNV();
3363 NumOfElements = NURBSext->GetNE();
3365 NURBSext->GetElementTopo(elements);
3378 NumOfBdrElements = 0;
3379 for (i = 0; i < num_pieces; i++)
3381 NumOfBdrElements += mesh_array[i]->
GetNBE();
3383 boundary.SetSize(NumOfBdrElements);
3384 vertices.SetSize(NumOfVertices);
3386 for (i = 0; i < num_pieces; i++)
3392 for (j = 0; j < m->
GetNE(); j++)
3394 elements[lelem_elem[j]]->SetAttribute(m->
GetAttribute(j));
3397 for (j = 0; j < m->
GetNBE(); j++)
3402 for (
int k = 0; k < nv; k++)
3404 v[k] = lvert_vert[v[k]];
3406 boundary[ib++] = el;
3409 for (j = 0; j < m->
GetNV(); j++)
3419 NumOfBdrElements = 0;
3421 for (i = 0; i < num_pieces; i++)
3424 NumOfElements += m->
GetNE();
3425 NumOfBdrElements += m->
GetNBE();
3426 NumOfVertices += m->
GetNV();
3428 elements.SetSize(NumOfElements);
3429 boundary.SetSize(NumOfBdrElements);
3430 vertices.SetSize(NumOfVertices);
3432 for (i = 0; i < num_pieces; i++)
3436 for (j = 0; j < m->
GetNE(); j++)
3441 for (
int k = 0; k < nv; k++)
3445 elements[ie++] = el;
3448 for (j = 0; j < m->
GetNBE(); j++)
3453 for (
int k = 0; k < nv; k++)
3457 boundary[ib++] = el;
3460 for (j = 0; j < m->
GetNV(); j++)
3474 for (i = 0; i < num_pieces; i++)
3476 gf_array[i] = mesh_array[i]->
GetNodes();
3483 CheckElementOrientation(
false);
3484 CheckBdrElementOrientation(
false);
3488 Mesh::Mesh(
Mesh *orig_mesh,
int ref_factor,
int ref_type)
3491 MFEM_VERIFY(ref_factor >= 1,
"the refinement factor must be >= 1");
3492 MFEM_VERIFY(ref_type == BasisType::ClosedUniform ||
3493 ref_type == BasisType::GaussLobatto,
"invalid refinement type");
3494 MFEM_VERIFY(Dim == 1 || Dim == 2 || Dim == 3,
3495 "only implemented for Segment, Quadrilateral and Hexahedron "
3496 "elements in 1D/2D/3D");
3498 "meshes with mixed elements are not supported");
3505 int r_bndr_factor = pow(ref_factor, Dim - 1);
3506 int r_elem_factor = ref_factor * r_bndr_factor;
3509 int r_num_elem = orig_mesh->
GetNE() * r_elem_factor;
3510 int r_num_bndr = orig_mesh->
GetNBE() * r_bndr_factor;
3512 InitMesh(Dim, orig_mesh->
SpaceDimension(), r_num_vert, r_num_elem,
3516 NumOfVertices = r_num_vert;
3522 DenseMatrix node_coordinates(spaceDim*pow(2, Dim), r_num_elem);
3525 for (
int el = 0; el < orig_mesh->
GetNE(); el++)
3529 int nvert = Geometry::NumVerts[
geom];
3532 max_nv = std::max(max_nv, nvert);
3538 const int *c2h_map = rfec.GetDofMap(geom);
3539 const int *vertex_map = vertex_fec.
GetDofMap(geom);
3540 for (
int i = 0; i < phys_pts.
Width(); i++)
3542 vertices[rdofs[i]].SetCoords(spaceDim, phys_pts.
GetColumn(i));
3546 Element *elem = NewElement(geom);
3549 for (
int k = 0; k < nvert; k++)
3552 v[k] = rdofs[c2h_map[cid]];
3554 for (
int k = 0; k < nvert; k++)
3556 for (
int j = 0; j < spaceDim; ++j)
3558 node_coordinates(k*spaceDim + j, NumOfElements)
3559 = vertices[v[vertex_map[k]]](j);
3566 SetCurvature(1,
true, spaceDim);
3567 Vector node_coordinates_vec(
3568 node_coordinates.
Data(),
3569 node_coordinates.
Width()*node_coordinates.
Height());
3570 SetNodes(node_coordinates_vec);
3573 for (
int el = 0; el < orig_mesh->
GetNBE(); el++)
3577 int nvert = Geometry::NumVerts[
geom];
3588 Element *elem = NewElement(geom);
3591 v[0] = rdofs[RG.
RefGeoms[nvert*j]];
3592 AddBdrElement(elem);
3597 const int *c2h_map = rfec.GetDofMap(geom);
3600 Element *elem = NewElement(geom);
3603 for (
int k = 0; k < nvert; k++)
3606 v[k] = rdofs[c2h_map[cid]];
3608 AddBdrElement(elem);
3613 FinalizeTopology(
false);
3615 last_operation = Mesh::REFINE;
3618 CoarseFineTr.embeddings.SetSize(GetNE());
3619 if (orig_mesh->
GetNE() > 0)
3623 CoarseFineTr.point_matrices[
geom].SetSize(Dim, max_nv, r_elem_factor);
3624 int nvert = Geometry::NumVerts[
geom];
3626 const int *c2h_map = rfec.GetDofMap(geom);
3631 for (
int k = 0; k < nvert; k++)
3639 for (
int el = 0; el < GetNE(); el++)
3641 Embedding &emb = CoarseFineTr.embeddings[el];
3642 emb.
parent = el / r_elem_factor;
3643 emb.
matrix = el % r_elem_factor;
3646 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
3647 MFEM_ASSERT(CheckBdrElementOrientation(
false) == 0,
"");
3652 if (NURBSext == NULL)
3654 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
3657 if (kv.
Size() != NURBSext->GetNKV())
3659 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
3662 NURBSext->ConvertToPatches(*Nodes);
3664 NURBSext->KnotInsert(kv);
3666 last_operation = Mesh::NONE;
3674 if (NURBSext == NULL)
3676 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
3679 if (kv.
Size() != NURBSext->GetNKV())
3681 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
3684 NURBSext->ConvertToPatches(*Nodes);
3686 NURBSext->KnotInsert(kv);
3688 last_operation = Mesh::NONE;
3694 void Mesh::NURBSUniformRefinement()
3697 NURBSext->ConvertToPatches(*Nodes);
3699 NURBSext->UniformRefinement();
3701 last_operation = Mesh::NONE;
3707 void Mesh::DegreeElevate(
int rel_degree,
int degree)
3709 if (NURBSext == NULL)
3711 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
3714 NURBSext->ConvertToPatches(*Nodes);
3716 NURBSext->DegreeElevate(rel_degree, degree);
3718 last_operation = Mesh::NONE;
3724 void Mesh::UpdateNURBS()
3728 NURBSext->SetKnotsFromPatches();
3730 Dim = NURBSext->Dimension();
3733 if (NumOfElements != NURBSext->GetNE())
3735 for (
int i = 0; i < elements.Size(); i++)
3737 FreeElement(elements[i]);
3739 NumOfElements = NURBSext->GetNE();
3740 NURBSext->GetElementTopo(elements);
3743 if (NumOfBdrElements != NURBSext->GetNBE())
3745 for (
int i = 0; i < boundary.Size(); i++)
3747 FreeElement(boundary[i]);
3749 NumOfBdrElements = NURBSext->GetNBE();
3750 NURBSext->GetBdrElementTopo(boundary);
3753 Nodes->FESpace()->Update();
3755 NURBSext->SetCoordsFromPatches(*Nodes);
3757 if (NumOfVertices != NURBSext->GetNV())
3759 NumOfVertices = NURBSext->GetNV();
3760 vertices.SetSize(NumOfVertices);
3761 int vd = Nodes->VectorDim();
3762 for (
int i = 0; i < vd; i++)
3765 Nodes->GetNodalValues(vert_val, i+1);
3766 for (
int j = 0; j < NumOfVertices; j++)
3768 vertices[j](i) = vert_val(j);
3775 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3784 GetElementToFaceTable();
3789 void Mesh::LoadPatchTopo(std::istream &input,
Array<int> &edge_to_knot)
3805 input >> NumOfElements;
3806 elements.SetSize(NumOfElements);
3807 for (
int j = 0; j < NumOfElements; j++)
3809 elements[j] = ReadElement(input);
3815 input >> NumOfBdrElements;
3816 boundary.SetSize(NumOfBdrElements);
3817 for (
int j = 0; j < NumOfBdrElements; j++)
3819 boundary[j] = ReadElement(input);
3825 input >> NumOfEdges;
3826 edge_vertex =
new Table(NumOfEdges, 2);
3827 edge_to_knot.
SetSize(NumOfEdges);
3828 for (
int j = 0; j < NumOfEdges; j++)
3830 int *v = edge_vertex->GetRow(j);
3831 input >> edge_to_knot[j] >> v[0] >> v[1];
3834 edge_to_knot[j] = -1 - edge_to_knot[j];
3841 input >> NumOfVertices;
3842 vertices.SetSize(0);
3845 CheckBdrElementOrientation();
3852 for (
int d = 0; d < v.
Size(); d++)
3860 for (d = 0; d < p.
Size(); d++)
3864 for ( ; d < v.
Size(); d++)
3873 if (Nodes == NULL || Nodes->FESpace() != nodes.
FESpace())
3888 SetNodalGridFunction(nodes,
true);
3891 void Mesh::EnsureNodes()
3896 if (dynamic_cast<const H1_FECollection*>(fec)
3897 || dynamic_cast<const L2_FECollection*>(fec))
3903 const int order = GetNodalFESpace()->GetOrder(0);
3904 SetCurvature(order,
false, -1, Ordering::byVDIM);
3909 SetCurvature(1,
false, -1, Ordering::byVDIM);
3916 NewNodes(*nodes, make_owner);
3921 return ((Nodes) ? Nodes->FESpace() : NULL);
3924 void Mesh::SetCurvature(
int order,
bool discont,
int space_dim,
int ordering)
3926 space_dim = (space_dim == -1) ? spaceDim : space_dim;
3939 SetNodalFESpace(nfes);
3940 Nodes->MakeOwner(nfec);
3943 int Mesh::GetNumFaces()
const
3947 case 1:
return GetNV();
3948 case 2:
return GetNEdges();
3949 case 3:
return GetNFaces();
3954 static int CountFacesByType(
const Mesh &mesh,
const FaceType type)
3963 if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
3964 (type==FaceType::Boundary && e2<0 && inf2<0) ) { nf++; }
3971 const bool isInt = type==FaceType::Interior;
3972 int &nf = isInt ? nbInteriorFaces : nbBoundaryFaces;
3973 if (nf<0) { nf = CountFacesByType(*
this, type); }
3977 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3978 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
3981 int Mesh::CheckElementOrientation(
bool fix_it)
3983 int i, j, k, wo = 0, fo = 0, *vi = 0;
3986 if (Dim == 2 && spaceDim == 2)
3990 for (i = 0; i < NumOfElements; i++)
3994 vi = elements[i]->GetVertices();
3995 for (j = 0; j < 3; j++)
3997 v[j] = vertices[vi[j]]();
3999 for (j = 0; j < 2; j++)
4000 for (k = 0; k < 2; k++)
4002 J(j, k) = v[j+1][k] - v[0][k];
4008 GetElementJacobian(i, J);
4014 switch (GetElementType(i))
4016 case Element::TRIANGLE:
4019 case Element::QUADRILATERAL:
4023 MFEM_ABORT(
"Invalid 2D element type \""
4024 << GetElementType(i) <<
"\"");
4038 for (i = 0; i < NumOfElements; i++)
4040 vi = elements[i]->GetVertices();
4041 switch (GetElementType(i))
4043 case Element::TETRAHEDRON:
4046 for (j = 0; j < 4; j++)
4048 v[j] = vertices[vi[j]]();
4050 for (j = 0; j < 3; j++)
4051 for (k = 0; k < 3; k++)
4053 J(j, k) = v[j+1][k] - v[0][k];
4059 GetElementJacobian(i, J);
4072 case Element::WEDGE:
4074 GetElementJacobian(i, J);
4085 case Element::HEXAHEDRON:
4087 GetElementJacobian(i, J);
4099 MFEM_ABORT(
"Invalid 3D element type \""
4100 << GetElementType(i) <<
"\"");
4105 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
4108 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
4109 << NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
4116 int Mesh::GetTriOrientation(
const int *base,
const int *test)
4124 if (test[0] == base[0])
4125 if (test[1] == base[1])
4133 else if (test[0] == base[1])
4134 if (test[1] == base[0])
4143 if (test[1] == base[0])
4153 const int *aor = tri_t::Orient[orient];
4154 for (
int j = 0; j < 3; j++)
4155 if (test[aor[j]] != base[j])
4164 int Mesh::GetQuadOrientation(
const int *base,
const int *test)
4168 for (i = 0; i < 4; i++)
4169 if (test[i] == base[0])
4176 if (test[(i+1)%4] == base[1])
4184 const int *aor = quad_t::Orient[orient];
4185 for (
int j = 0; j < 4; j++)
4186 if (test[aor[j]] != base[j])
4188 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
4190 for (
int k = 0; k < 4; k++)
4195 for (
int k = 0; k < 4; k++)
4204 if (test[(i+1)%4] == base[1])
4212 int Mesh::CheckBdrElementOrientation(
bool fix_it)
4218 if (el_to_edge == NULL)
4220 el_to_edge =
new Table;
4221 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
4224 for (
int i = 0; i < NumOfBdrElements; i++)
4226 if (faces_info[be_to_edge[i]].Elem2No < 0)
4228 int *bv = boundary[i]->GetVertices();
4229 int *fv = faces[be_to_edge[i]]->GetVertices();
4234 mfem::Swap<int>(bv[0], bv[1]);
4244 for (
int i = 0; i < NumOfBdrElements; i++)
4246 const int fi = be_to_face[i];
4248 if (faces_info[fi].Elem2No >= 0) {
continue; }
4251 int *bv = boundary[i]->GetVertices();
4253 MFEM_ASSERT(fi < faces.Size(),
"internal error");
4254 const int *fv = faces[fi]->GetVertices();
4260 case Element::TRIANGLE:
4262 orientation = GetTriOrientation(fv, bv);
4265 case Element::QUADRILATERAL:
4267 orientation = GetQuadOrientation(fv, bv);
4271 MFEM_ABORT(
"Invalid 2D boundary element type \""
4272 << bdr_type <<
"\"");
4277 if (orientation % 2 == 0) {
continue; }
4279 if (!fix_it) {
continue; }
4283 case Element::TRIANGLE:
4287 mfem::Swap<int>(bv[0], bv[1]);
4290 int *be = bel_to_edge->GetRow(i);
4291 mfem::Swap<int>(be[1], be[2]);
4295 case Element::QUADRILATERAL:
4297 mfem::Swap<int>(bv[0], bv[2]);
4300 int *be = bel_to_edge->GetRow(i);
4301 mfem::Swap<int>(be[0], be[1]);
4302 mfem::Swap<int>(be[2], be[3]);
4315 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
4316 << NumOfBdrElements <<
" (" << fixed_or_not[fix_it ? 0 : 1]
4323 int Mesh::GetNumGeometries(
int dim)
const
4325 MFEM_ASSERT(0 <= dim && dim <= Dim,
"invalid dim: " << dim);
4327 for (
int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
4336 MFEM_ASSERT(0 <= dim && dim <= Dim,
"invalid dim: " << dim);
4338 for (
int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
4351 el_to_edge->GetRow(i, edges);
4355 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
4356 "is not generated.");
4359 const int *v = elements[i]->GetVertices();
4360 const int ne = elements[i]->GetNEdges();
4362 for (
int j = 0; j < ne; j++)
4364 const int *e = elements[i]->GetEdgeVertices(j);
4365 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4375 edges[0] = be_to_edge[i];
4376 const int *v = boundary[i]->GetVertices();
4377 cor[0] = (v[0] < v[1]) ? (1) : (-1);
4383 bel_to_edge->GetRow(i, edges);
4390 const int *v = boundary[i]->GetVertices();
4391 const int ne = boundary[i]->GetNEdges();
4393 for (
int j = 0; j < ne; j++)
4395 const int *e = boundary[i]->GetEdgeVertices(j);
4396 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4408 const int *v = faces[i]->GetVertices();
4409 o[0] = (v[0] < v[1]) ? (1) : (-1);
4419 face_edge->GetRow(i, edges);
4421 const int *v = faces[i]->GetVertices();
4422 const int ne = faces[i]->GetNEdges();
4424 for (
int j = 0; j < ne; j++)
4426 const int *e = faces[i]->GetEdgeVertices(j);
4427 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4436 if (!edge_vertex) { GetEdgeVertexTable(); }
4437 edge_vertex->GetRow(i, vert);
4453 if (faces.Size() != NumOfFaces)
4455 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
4459 DSTable v_to_v(NumOfVertices);
4460 GetVertexToVertexTable(v_to_v);
4462 face_edge =
new Table;
4463 GetElementArrayEdgeTable(faces, v_to_v, *face_edge);
4475 DSTable v_to_v(NumOfVertices);
4476 GetVertexToVertexTable(v_to_v);
4479 edge_vertex =
new Table(nedges, 2);
4480 for (
int i = 0; i < NumOfVertices; i++)
4485 edge_vertex->Push(j, i);
4486 edge_vertex->Push(j, it.Column());
4489 edge_vertex->Finalize();
4500 vert_elem->
MakeI(NumOfVertices);
4502 for (i = 0; i < NumOfElements; i++)
4504 nv = elements[i]->GetNVertices();
4505 v = elements[i]->GetVertices();
4506 for (j = 0; j < nv; j++)
4514 for (i = 0; i < NumOfElements; i++)
4516 nv = elements[i]->GetNVertices();
4517 v = elements[i]->GetVertices();
4518 for (j = 0; j < nv; j++)
4529 Table *Mesh::GetFaceToElementTable()
const
4533 face_elem->
MakeI(faces_info.Size());
4535 for (
int i = 0; i < faces_info.Size(); i++)
4537 if (faces_info[i].Elem2No >= 0)
4549 for (
int i = 0; i < faces_info.Size(); i++)
4552 if (faces_info[i].Elem2No >= 0)
4570 el_to_face->
GetRow(i, fcs);
4574 mfem_error(
"Mesh::GetElementFaces(...) : el_to_face not generated.");
4579 for (j = 0; j < n; j++)
4580 if (faces_info[fcs[j]].Elem1No == i)
4582 cor[j] = faces_info[fcs[j]].Elem1Inf % 64;
4585 else if (faces_info[fcs[j]].Elem2No == i)
4587 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
4591 mfem_error(
"Mesh::GetElementFaces(...) : 2");
4596 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
4601 void Mesh::GetBdrElementFace(
int i,
int *f,
int *o)
const
4606 bv = boundary[i]->GetVertices();
4607 fv = faces[be_to_face[i]]->GetVertices();
4611 switch (GetBdrElementType(i))
4613 case Element::TRIANGLE:
4614 *o = GetTriOrientation(fv, bv);
4616 case Element::QUADRILATERAL:
4617 *o = GetQuadOrientation(fv, bv);
4620 mfem_error(
"Mesh::GetBdrElementFace(...) 2");
4624 int Mesh::GetBdrElementEdgeIndex(
int i)
const
4628 case 1:
return boundary[i]->GetVertices()[0];
4629 case 2:
return be_to_edge[i];
4630 case 3:
return be_to_face[i];
4631 default:
mfem_error(
"Mesh::GetBdrElementEdgeIndex: invalid dimension!");
4636 void Mesh::GetBdrElementAdjacentElement(
int bdr_el,
int &el,
int &info)
const
4638 int fid = GetBdrElementEdgeIndex(bdr_el);
4639 const FaceInfo &fi = faces_info[fid];
4640 MFEM_ASSERT(fi.
Elem1Inf%64 == 0,
"internal error");
4641 const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
4642 const int *bv = boundary[bdr_el]->GetVertices();
4644 switch (GetBdrElementBaseGeometry(bdr_el))
4646 case Geometry::POINT: ori = 0;
break;
4647 case Geometry::SEGMENT: ori = (fv[0] == bv[0]) ? 0 : 1;
break;
4648 case Geometry::TRIANGLE: ori = GetTriOrientation(fv, bv);
break;
4649 case Geometry::SQUARE: ori = GetQuadOrientation(fv, bv);
break;
4650 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
4658 return elements[i]->GetType();
4663 return boundary[i]->GetType();
4671 v = elements[i]->GetVertices();
4672 nv = elements[i]->GetNVertices();
4674 pointmat.
SetSize(spaceDim, nv);
4675 for (k = 0; k < spaceDim; k++)
4676 for (j = 0; j < nv; j++)
4678 pointmat(k, j) = vertices[v[j]](k);
4687 v = boundary[i]->GetVertices();
4688 nv = boundary[i]->GetNVertices();
4690 pointmat.
SetSize(spaceDim, nv);
4691 for (k = 0; k < spaceDim; k++)
4692 for (j = 0; j < nv; j++)
4694 pointmat(k, j) = vertices[v[j]](k);
4698 double Mesh::GetLength(
int i,
int j)
const
4700 const double *vi = vertices[i]();
4701 const double *vj = vertices[j]();
4704 for (
int k = 0; k < spaceDim; k++)
4706 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
4709 return sqrt(length);
4717 for (
int i = 0; i < elem_array.
Size(); i++)
4722 for (
int i = 0; i < elem_array.
Size(); i++)
4724 const int *v = elem_array[i]->GetVertices();
4725 const int ne = elem_array[i]->GetNEdges();
4726 for (
int j = 0; j < ne; j++)
4728 const int *e = elem_array[i]->GetEdgeVertices(j);
4735 void Mesh::GetVertexToVertexTable(
DSTable &v_to_v)
const
4739 for (
int i = 0; i < edge_vertex->Size(); i++)
4741 const int *v = edge_vertex->GetRow(i);
4742 v_to_v.
Push(v[0], v[1]);
4747 for (
int i = 0; i < NumOfElements; i++)
4749 const int *v = elements[i]->GetVertices();
4750 const int ne = elements[i]->GetNEdges();
4751 for (
int j = 0; j < ne; j++)
4753 const int *e = elements[i]->GetEdgeVertices(j);
4754 v_to_v.
Push(v[e[0]], v[e[1]]);
4762 int i, NumberOfEdges;
4764 DSTable v_to_v(NumOfVertices);
4765 GetVertexToVertexTable(v_to_v);
4770 GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
4775 be_to_f.
SetSize(NumOfBdrElements);
4776 for (i = 0; i < NumOfBdrElements; i++)
4778 const int *v = boundary[i]->GetVertices();
4779 be_to_f[i] = v_to_v(v[0], v[1]);
4784 if (bel_to_edge == NULL)
4786 bel_to_edge =
new Table;
4788 GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
4792 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
4796 return NumberOfEdges;
4799 const Table & Mesh::ElementToElementTable()
4807 MFEM_ASSERT(faces_info.Size() >= GetNumFaces(),
"faces were not generated!");
4810 conn.
Reserve(2*faces_info.Size());
4812 for (
int i = 0; i < faces_info.Size(); i++)
4814 const FaceInfo &fi = faces_info[i];
4822 int nbr_elem_idx = NumOfElements - 1 - fi.
Elem2No;
4830 el_to_el =
new Table(NumOfElements, conn);
4835 const Table & Mesh::ElementToFaceTable()
const
4837 if (el_to_face == NULL)
4844 const Table & Mesh::ElementToEdgeTable()
const
4846 if (el_to_edge == NULL)
4853 void Mesh::AddPointFaceElement(
int lf,
int gf,
int el)
4855 if (faces_info[gf].Elem1No == -1)
4858 faces_info[gf].Elem1No = el;
4859 faces_info[gf].Elem1Inf = 64 * lf;
4860 faces_info[gf].Elem2No = -1;
4861 faces_info[gf].Elem2Inf = -1;
4865 faces_info[gf].Elem2No = el;
4866 faces_info[gf].Elem2Inf = 64 * lf + 1;
4870 void Mesh::AddSegmentFaceElement(
int lf,
int gf,
int el,
int v0,
int v1)
4872 if (faces[gf] == NULL)
4874 faces[gf] =
new Segment(v0, v1);
4875 faces_info[gf].Elem1No = el;
4876 faces_info[gf].Elem1Inf = 64 * lf;
4877 faces_info[gf].Elem2No = -1;
4878 faces_info[gf].Elem2Inf = -1;
4882 int *v = faces[gf]->GetVertices();
4883 faces_info[gf].Elem2No = el;
4884 if ( v[1] == v0 && v[0] == v1 )
4886 faces_info[gf].Elem2Inf = 64 * lf + 1;
4888 else if ( v[0] == v0 && v[1] == v1 )
4894 faces_info[gf].Elem2Inf = 64 * lf;
4898 MFEM_ABORT(
"internal error");
4903 void Mesh::AddTriangleFaceElement(
int lf,
int gf,
int el,
4904 int v0,
int v1,
int v2)
4906 if (faces[gf] == NULL)
4908 faces[gf] =
new Triangle(v0, v1, v2);
4909 faces_info[gf].Elem1No = el;
4910 faces_info[gf].Elem1Inf = 64 * lf;
4911 faces_info[gf].Elem2No = -1;
4912 faces_info[gf].Elem2Inf = -1;
4916 int orientation, vv[3] = { v0, v1, v2 };
4917 orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
4922 faces_info[gf].Elem2No = el;
4923 faces_info[gf].Elem2Inf = 64 * lf + orientation;
4927 void Mesh::AddQuadFaceElement(
int lf,
int gf,
int el,
4928 int v0,
int v1,
int v2,
int v3)
4930 if (faces_info[gf].Elem1No < 0)
4933 faces_info[gf].Elem1No = el;
4934 faces_info[gf].Elem1Inf = 64 * lf;
4935 faces_info[gf].Elem2No = -1;
4936 faces_info[gf].Elem2Inf = -1;
4940 int vv[4] = { v0, v1, v2, v3 };
4941 int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
4945 faces_info[gf].Elem2No = el;
4946 faces_info[gf].Elem2Inf = 64 * lf + oo;
4950 void Mesh::GenerateFaces()
4952 int i, nfaces = GetNumFaces();
4954 for (i = 0; i < faces.Size(); i++)
4956 FreeElement(faces[i]);
4960 faces.SetSize(nfaces);
4961 faces_info.SetSize(nfaces);
4962 for (i = 0; i < nfaces; i++)
4965 faces_info[i].Elem1No = -1;
4966 faces_info[i].NCFace = -1;
4968 for (i = 0; i < NumOfElements; i++)
4970 const int *v = elements[i]->GetVertices();
4974 AddPointFaceElement(0, v[0], i);
4975 AddPointFaceElement(1, v[1], i);
4979 ef = el_to_edge->GetRow(i);
4980 const int ne = elements[i]->GetNEdges();
4981 for (
int j = 0; j < ne; j++)
4983 const int *e = elements[i]->GetEdgeVertices(j);
4984 AddSegmentFaceElement(j, ef[j], i, v[e[0]], v[e[1]]);
4989 ef = el_to_face->GetRow(i);
4990 switch (GetElementType(i))
4992 case Element::TETRAHEDRON:
4994 for (
int j = 0; j < 4; j++)
4996 const int *fv = tet_t::FaceVert[j];
4997 AddTriangleFaceElement(j, ef[j], i,
4998 v[fv[0]], v[fv[1]], v[fv[2]]);
5002 case Element::WEDGE:
5004 for (
int j = 0; j < 2; j++)
5006 const int *fv = pri_t::FaceVert[j];
5007 AddTriangleFaceElement(j, ef[j], i,
5008 v[fv[0]], v[fv[1]], v[fv[2]]);
5010 for (
int j = 2; j < 5; j++)
5012 const int *fv = pri_t::FaceVert[j];
5013 AddQuadFaceElement(j, ef[j], i,
5014 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5018 case Element::HEXAHEDRON:
5020 for (
int j = 0; j < 6; j++)
5022 const int *fv = hex_t::FaceVert[j];
5023 AddQuadFaceElement(j, ef[j], i,
5024 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5029 MFEM_ABORT(
"Unexpected type of Element.");
5035 void Mesh::GenerateNCFaceInfo()
5037 MFEM_VERIFY(ncmesh,
"missing NCMesh.");
5039 for (
int i = 0; i < faces_info.Size(); i++)
5041 faces_info[i].NCFace = -1;
5045 (Dim == 2) ? ncmesh->GetEdgeList() : ncmesh->GetFaceList();
5047 nc_faces_info.SetSize(0);
5048 nc_faces_info.Reserve(list.
masters.size() + list.
slaves.size());
5050 int nfaces = GetNumFaces();
5053 for (
unsigned i = 0; i < list.
masters.size(); i++)
5056 if (master.
index >= nfaces) {
continue; }
5058 faces_info[master.
index].NCFace = nc_faces_info.Size();
5064 for (
unsigned i = 0; i < list.
slaves.size(); i++)
5068 if (slave.
index < 0 ||
5069 slave.
index >= nfaces ||
5079 slave_fi.
NCFace = nc_faces_info.Size();
5091 for (
int i = 0; i < NumOfElements; i++)
5093 const int *v = elements[i]->GetVertices();
5094 switch (GetElementType(i))
5096 case Element::TETRAHEDRON:
5098 for (
int j = 0; j < 4; j++)
5100 const int *fv = tet_t::FaceVert[j];
5101 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
5105 case Element::WEDGE:
5107 for (
int j = 0; j < 2; j++)
5109 const int *fv = pri_t::FaceVert[j];
5110 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
5112 for (
int j = 2; j < 5; j++)
5114 const int *fv = pri_t::FaceVert[j];
5115 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5119 case Element::HEXAHEDRON:
5123 for (
int j = 0; j < 6; j++)
5125 const int *fv = hex_t::FaceVert[j];
5126 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
5131 MFEM_ABORT(
"Unexpected type of Element.");
5142 if (el_to_face != NULL)
5146 el_to_face =
new Table(NumOfElements, 6);
5147 faces_tbl =
new STable3D(NumOfVertices);
5148 for (i = 0; i < NumOfElements; i++)
5150 v = elements[i]->GetVertices();
5151 switch (GetElementType(i))
5153 case Element::TETRAHEDRON:
5155 for (
int j = 0; j < 4; j++)
5157 const int *fv = tet_t::FaceVert[j];
5159 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
5163 case Element::WEDGE:
5165 for (
int j = 0; j < 2; j++)
5167 const int *fv = pri_t::FaceVert[j];
5169 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
5171 for (
int j = 2; j < 5; j++)
5173 const int *fv = pri_t::FaceVert[j];
5175 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
5179 case Element::HEXAHEDRON:
5183 for (
int j = 0; j < 6; j++)
5185 const int *fv = hex_t::FaceVert[j];
5187 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
5192 MFEM_ABORT(
"Unexpected type of Element.");
5195 el_to_face->Finalize();
5197 be_to_face.SetSize(NumOfBdrElements);
5198 for (i = 0; i < NumOfBdrElements; i++)
5200 v = boundary[i]->GetVertices();
5201 switch (GetBdrElementType(i))
5203 case Element::TRIANGLE:
5205 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
5208 case Element::QUADRILATERAL:
5210 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
5214 MFEM_ABORT(
"Unexpected type of boundary Element.");
5228 void Rotate3(
int &
a,
int &
b,
int &c)
5250 void Mesh::ReorientTetMesh()
5252 if (Dim != 3 || !(meshgen & 1))
5260 Table *old_elem_vert = NULL;
5264 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
5267 for (
int i = 0; i < NumOfElements; i++)
5269 if (GetElementType(i) == Element::TETRAHEDRON)
5271 int *v = elements[i]->GetVertices();
5273 Rotate3(v[0], v[1], v[2]);
5276 Rotate3(v[1], v[2], v[3]);
5285 for (
int i = 0; i < NumOfBdrElements; i++)
5287 if (GetBdrElementType(i) == Element::TRIANGLE)
5289 int *v = boundary[i]->GetVertices();
5291 Rotate3(v[0], v[1], v[2]);
5297 GetElementToFaceTable();
5301 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5306 DoNodeReorder(old_v_to_v, old_elem_vert);
5307 delete old_elem_vert;
5312 int *Mesh::CartesianPartitioning(
int nxyz[])
5318 for (
int vi = 0; vi < NumOfVertices; vi++)
5320 const double *p = vertices[vi]();
5321 for (
int i = 0; i < spaceDim; i++)
5323 if (p[i] < pmin[i]) { pmin[i] = p[i]; }
5324 if (p[i] > pmax[i]) { pmax[i] = p[i]; }
5328 partitioning =
new int[NumOfElements];
5332 Vector pt(ppt, spaceDim);
5333 for (
int el = 0; el < NumOfElements; el++)
5335 GetElementTransformation(el)->Transform(
5338 for (
int i = spaceDim-1; i >= 0; i--)
5340 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
5341 if (idx < 0) { idx = 0; }
5342 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
5343 part = part * nxyz[i] + idx;
5345 partitioning[el] = part;
5348 return partitioning;
5351 int *Mesh::GeneratePartitioning(
int nparts,
int part_method)
5353 #ifdef MFEM_USE_METIS
5355 int print_messages = 1;
5358 int init_flag, fin_flag;
5359 MPI_Initialized(&init_flag);
5360 MPI_Finalized(&fin_flag);
5361 if (init_flag && !fin_flag)
5365 if (rank != 0) { print_messages = 0; }
5369 int i, *partitioning;
5371 ElementToElementTable();
5373 partitioning =
new int[NumOfElements];
5377 for (i = 0; i < NumOfElements; i++)
5379 partitioning[i] = 0;
5382 else if (NumOfElements <= nparts)
5384 for (i = 0; i < NumOfElements; i++)
5386 partitioning[i] = i;
5392 #ifndef MFEM_USE_METIS_5
5404 bool freedata =
false;
5406 idx_t *mpartitioning;
5409 if (
sizeof(
idx_t) ==
sizeof(int))
5411 I = (
idx_t*) el_to_el->GetI();
5412 J = (
idx_t*) el_to_el->GetJ();
5413 mpartitioning = (
idx_t*) partitioning;
5417 int *iI = el_to_el->GetI();
5418 int *iJ = el_to_el->GetJ();
5422 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
5423 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
5424 mpartitioning =
new idx_t[n];
5427 #ifndef MFEM_USE_METIS_5
5430 METIS_SetDefaultOptions(options);
5431 options[METIS_OPTION_CONTIG] = 1;
5435 if (part_method >= 0 && part_method <= 2)
5437 for (i = 0; i < n; i++)
5443 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
5449 if (part_method == 0 || part_method == 3)
5451 #ifndef MFEM_USE_METIS_5
5480 " error in METIS_PartGraphRecursive!");
5487 if (part_method == 1 || part_method == 4)
5489 #ifndef MFEM_USE_METIS_5
5518 " error in METIS_PartGraphKway!");
5525 if (part_method == 2 || part_method == 5)
5527 #ifndef MFEM_USE_METIS_5
5540 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
5557 " error in METIS_PartGraphKway!");
5565 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
5569 nparts = (int) mparts;
5570 if (mpartitioning != (
idx_t*)partitioning)
5572 for (
int k = 0; k<NumOfElements; k++)
5574 partitioning[k] = mpartitioning[k];
5581 delete[] mpartitioning;
5591 for (i = 0; i < nparts; i++)
5597 for (i = 0; i < NumOfElements; i++)
5599 psize[partitioning[i]].one++;
5602 int empty_parts = 0;
5603 for (i = 0; i < nparts; i++)
5605 if (psize[i].one == 0) { empty_parts++; }
5614 mfem::err <<
"Mesh::GeneratePartitioning(...): METIS returned "
5615 << empty_parts <<
" empty parts!"
5616 <<
" Applying a simple fix ..." << endl;
5619 SortPairs<int,int>(psize, nparts);
5621 for (i = nparts-1; i > nparts-1-empty_parts; i--)
5626 for (
int j = 0; j < NumOfElements; j++)
5628 for (i = nparts-1; i > nparts-1-empty_parts; i--)
5630 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
5636 partitioning[j] = psize[nparts-1-i].two;
5644 return partitioning;
5648 mfem_error(
"Mesh::GeneratePartitioning(...): "
5649 "MFEM was compiled without Metis.");
5663 int num_elem, *i_elem_elem, *j_elem_elem;
5665 num_elem = elem_elem.
Size();
5666 i_elem_elem = elem_elem.
GetI();
5667 j_elem_elem = elem_elem.
GetJ();
5672 int stack_p, stack_top_p, elem;
5676 for (i = 0; i < num_elem; i++)
5678 if (partitioning[i] > num_part)
5680 num_part = partitioning[i];
5687 for (i = 0; i < num_part; i++)
5694 for (elem = 0; elem < num_elem; elem++)
5696 if (component[elem] >= 0)
5701 component[elem] = num_comp[partitioning[elem]]++;
5703 elem_stack[stack_top_p++] = elem;
5705 for ( ; stack_p < stack_top_p; stack_p++)
5707 i = elem_stack[stack_p];
5708 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
5711 if (partitioning[k] == partitioning[i])
5713 if (component[k] < 0)
5715 component[k] = component[i];
5716 elem_stack[stack_top_p++] = k;
5718 else if (component[k] != component[i])
5728 void Mesh::CheckPartitioning(
int *partitioning)
5730 int i, n_empty, n_mcomp;
5732 const Array<int> _partitioning(partitioning, GetNE());
5734 ElementToElementTable();
5738 n_empty = n_mcomp = 0;
5739 for (i = 0; i < num_comp.
Size(); i++)
5740 if (num_comp[i] == 0)
5744 else if (num_comp[i] > 1)
5751 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
5752 <<
"The following subdomains are empty :\n";
5753 for (i = 0; i < num_comp.
Size(); i++)
5754 if (num_comp[i] == 0)
5762 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
5763 <<
"The following subdomains are NOT connected :\n";
5764 for (i = 0; i < num_comp.
Size(); i++)
5765 if (num_comp[i] > 1)
5771 if (n_empty == 0 && n_mcomp == 0)
5772 mfem::out <<
"Mesh::CheckPartitioning(...) : "
5773 "All subdomains are connected." << endl;
5787 const double *a = A.
Data();
5788 const double *b = B.
Data();
5797 c(0) = a[0]*a[3]-a[1]*a[2];
5798 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
5799 c(2) = b[0]*b[3]-b[1]*b[2];
5820 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
5821 a[1] * (a[5] * a[6] - a[3] * a[8]) +
5822 a[2] * (a[3] * a[7] - a[4] * a[6]));
5824 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
5825 b[1] * (a[5] * a[6] - a[3] * a[8]) +
5826 b[2] * (a[3] * a[7] - a[4] * a[6]) +
5828 a[0] * (b[4] * a[8] - b[5] * a[7]) +
5829 a[1] * (b[5] * a[6] - b[3] * a[8]) +
5830 a[2] * (b[3] * a[7] - b[4] * a[6]) +
5832 a[0] * (a[4] * b[8] - a[5] * b[7]) +
5833 a[1] * (a[5] * b[6] - a[3] * b[8]) +
5834 a[2] * (a[3] * b[7] - a[4] * b[6]));
5836 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
5837 a[1] * (b[5] * b[6] - b[3] * b[8]) +
5838 a[2] * (b[3] * b[7] - b[4] * b[6]) +
5840 b[0] * (a[4] * b[8] - a[5] * b[7]) +
5841 b[1] * (a[5] * b[6] - a[3] * b[8]) +
5842 b[2] * (a[3] * b[7] - a[4] * b[6]) +
5844 b[0] * (b[4] * a[8] - b[5] * a[7]) +
5845 b[1] * (b[5] * a[6] - b[3] * a[8]) +
5846 b[2] * (b[3] * a[7] - b[4] * a[6]));
5848 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
5849 b[1] * (b[5] * b[6] - b[3] * b[8]) +
5850 b[2] * (b[3] * b[7] - b[4] * b[6]));
5896 double a = z(2), b = z(1), c = z(0);
5897 double D = b*b-4*a*c;
5904 x(0) = x(1) = -0.5 * b /
a;
5909 x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
5917 t = -0.5 * (b + sqrt(D));
5921 t = -0.5 * (b - sqrt(D));
5927 Swap<double>(x(0), x(1));
5935 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
5938 double Q = (a * a - 3 *
b) / 9;
5939 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
5940 double Q3 = Q * Q * Q;
5947 x(0) = x(1) = x(2) = - a / 3;
5951 double sqrtQ = sqrt(Q);
5955 x(0) = -2 * sqrtQ - a / 3;
5956 x(1) = x(2) = sqrtQ - a / 3;
5960 x(0) = x(1) = - sqrtQ - a / 3;
5961 x(2) = 2 * sqrtQ - a / 3;
5968 double theta = acos(R / sqrt(Q3));
5969 double A = -2 * sqrt(Q);
5971 x0 = A * cos(theta / 3) - a / 3;
5972 x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
5973 x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
5978 Swap<double>(x0, x1);
5982 Swap<double>(x1, x2);
5985 Swap<double>(x0, x1);
5998 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
6002 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
6004 x(0) = A + Q / A - a / 3;
6013 const double factor,
const int Dim)
6015 const double c0 = c(0);
6016 c(0) = c0 * (1.0 - pow(factor, -Dim));
6018 for (
int j = 0; j < nr; j++)
6030 c(0) = c0 * (1.0 - pow(factor, Dim));
6032 for (
int j = 0; j < nr; j++)
6046 void Mesh::CheckDisplacements(
const Vector &displacements,
double &tmax)
6048 int nvs = vertices.Size();
6049 DenseMatrix P, V, DS, PDS(spaceDim), VDS(spaceDim);
6050 Vector c(spaceDim+1), x(spaceDim);
6051 const double factor = 2.0;
6058 for (
int i = 0; i < NumOfElements; i++)
6065 for (
int j = 0; j < spaceDim; j++)
6066 for (
int k = 0; k < nv; k++)
6068 P(j, k) = vertices[v[k]](j);
6069 V(j, k) = displacements(v[k]+j*nvs);
6073 GetTransformationFEforElementType(el->
GetType());
6077 case Element::TRIANGLE:
6078 case Element::TETRAHEDRON:
6096 case Element::QUADRILATERAL:
6099 for (
int j = 0; j < nv; j++)
6123 void Mesh::MoveVertices(
const Vector &displacements)
6125 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
6126 for (
int j = 0; j < spaceDim; j++)
6128 vertices[i](j) += displacements(j*nv+i);
6132 void Mesh::GetVertices(
Vector &vert_coord)
const
6134 int nv = vertices.Size();
6135 vert_coord.
SetSize(nv*spaceDim);
6136 for (
int i = 0; i < nv; i++)
6137 for (
int j = 0; j < spaceDim; j++)
6139 vert_coord(j*nv+i) = vertices[i](j);
6143 void Mesh::SetVertices(
const Vector &vert_coord)
6145 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
6146 for (
int j = 0; j < spaceDim; j++)
6148 vertices[i](j) = vert_coord(j*nv+i);
6152 void Mesh::GetNode(
int i,
double *coord)
const
6157 for (
int j = 0; j < spaceDim; j++)
6159 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
6164 for (
int j = 0; j < spaceDim; j++)
6166 coord[j] = vertices[i](j);
6171 void Mesh::SetNode(
int i,
const double *coord)
6176 for (
int j = 0; j < spaceDim; j++)
6178 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
6183 for (
int j = 0; j < spaceDim; j++)
6185 vertices[i](j) = coord[j];
6191 void Mesh::MoveNodes(
const Vector &displacements)
6195 (*Nodes) += displacements;
6199 MoveVertices(displacements);
6203 void Mesh::GetNodes(
Vector &node_coord)
const
6207 node_coord = (*Nodes);
6211 GetVertices(node_coord);
6215 void Mesh::SetNodes(
const Vector &node_coord)
6219 (*Nodes) = node_coord;
6223 SetVertices(node_coord);
6229 if (own_nodes) {
delete Nodes; }
6232 own_nodes = (int)make_owner;
6243 mfem::Swap<GridFunction*>(Nodes, nodes);
6244 mfem::Swap<int>(own_nodes, own_nodes_);
6251 void Mesh::AverageVertices(
const int *indexes,
int n,
int result)
6255 for (k = 0; k < spaceDim; k++)
6257 vertices[result](k) = vertices[indexes[0]](k);
6260 for (j = 1; j < n; j++)
6261 for (k = 0; k < spaceDim; k++)
6263 vertices[result](k) += vertices[indexes[j]](k);
6266 for (k = 0; k < spaceDim; k++)
6268 vertices[result](k) *= (1.0 / n);
6272 void Mesh::UpdateNodes()
6276 Nodes->FESpace()->Update();
6281 void Mesh::UniformRefinement2D_base(
bool update_nodes)
6285 if (el_to_edge == NULL)
6287 el_to_edge =
new Table;
6288 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6291 int quad_counter = 0;
6292 for (
int i = 0; i < NumOfElements; i++)
6294 if (elements[i]->GetType() == Element::QUADRILATERAL)
6300 const int oedge = NumOfVertices;
6301 const int oelem = oedge + NumOfEdges;
6306 vertices.
SetSize(oelem + quad_counter);
6307 new_elements.
SetSize(4 * NumOfElements);
6310 for (
int i = 0, j = 0; i < NumOfElements; i++)
6313 const int attr = elements[i]->GetAttribute();
6314 int *v = elements[i]->GetVertices();
6315 const int *e = el_to_edge->GetRow(i);
6318 if (el_type == Element::TRIANGLE)
6320 for (
int ei = 0; ei < 3; ei++)
6322 for (
int k = 0; k < 2; k++)
6324 vv[k] = v[tri_t::Edges[ei][k]];
6326 AverageVertices(vv, 2, oedge+e[ei]);
6330 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
6332 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
6334 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
6336 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
6338 else if (el_type == Element::QUADRILATERAL)
6340 const int qe = quad_counter;
6342 AverageVertices(v, 4, oelem+qe);
6344 for (
int ei = 0; ei < 4; ei++)
6346 for (
int k = 0; k < 2; k++)
6348 vv[k] = v[quad_t::Edges[ei][k]];
6350 AverageVertices(vv, 2, oedge+e[ei]);
6354 new Quadrilateral(v[0], oedge+e[0], oelem+qe, oedge+e[3], attr);
6356 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oelem+qe, attr);
6358 new Quadrilateral(oelem+qe, oedge+e[1], v[2], oedge+e[2], attr);
6360 new Quadrilateral(oedge+e[3], oelem+qe, oedge+e[2], v[3], attr);
6364 MFEM_ABORT(
"unknown element type: " << el_type);
6366 FreeElement(elements[i]);
6371 new_boundary.
SetSize(2 * NumOfBdrElements);
6372 for (
int i = 0, j = 0; i < NumOfBdrElements; i++)
6374 const int attr = boundary[i]->GetAttribute();
6375 int *v = boundary[i]->GetVertices();
6377 new_boundary[j++] =
new Segment(v[0], oedge+be_to_edge[i], attr);
6378 new_boundary[j++] =
new Segment(oedge+be_to_edge[i], v[1], attr);
6380 FreeElement(boundary[i]);
6384 static const double A = 0.0, B = 0.5, C = 1.0;
6385 static double tri_children[2*3*4] =
6392 static double quad_children[2*4*4] =
6400 CoarseFineTr.point_matrices[Geometry::TRIANGLE]
6401 .UseExternalData(tri_children, 2, 3, 4);
6402 CoarseFineTr.point_matrices[Geometry::SQUARE]
6403 .UseExternalData(quad_children, 2, 4, 4);
6404 CoarseFineTr.embeddings.SetSize(elements.Size());
6406 for (
int i = 0; i < elements.Size(); i++)
6408 Embedding &emb = CoarseFineTr.embeddings[i];
6413 NumOfVertices = vertices.Size();
6414 NumOfElements = 4 * NumOfElements;
6415 NumOfBdrElements = 2 * NumOfBdrElements;
6418 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6421 last_operation = Mesh::REFINE;
6424 if (update_nodes) { UpdateNodes(); }
6427 if (!Nodes || update_nodes)
6429 CheckElementOrientation(
false);
6431 CheckBdrElementOrientation(
false);
6435 static inline double sqr(
const double &x)
6445 if (el_to_edge == NULL)
6447 el_to_edge =
new Table;
6448 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6451 if (el_to_face == NULL)
6453 GetElementToFaceTable();
6457 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
6460 int NumOfQuadFaces = 0;
6461 if (HasGeometry(Geometry::SQUARE))
6463 if (HasGeometry(Geometry::TRIANGLE))
6466 for (
int i = 0; i < faces.Size(); i++)
6468 if (faces[i]->GetType() == Element::QUADRILATERAL)
6470 f2qf[i] = NumOfQuadFaces;
6477 NumOfQuadFaces = faces.
Size();
6481 int hex_counter = 0;
6482 if (HasGeometry(Geometry::CUBE))
6484 for (
int i = 0; i < elements.Size(); i++)
6486 if (elements[i]->GetType() == Element::HEXAHEDRON)
6496 if (HasGeometry(Geometry::TETRAHEDRON))
6500 DSTable *v_to_v_ptr = v_to_v_p;
6503 v_to_v_ptr =
new DSTable(NumOfVertices);
6504 GetVertexToVertexTable(*v_to_v_ptr);
6509 for (
int i = 0; i < NumOfVertices; i++)
6516 std::sort(row_start, J_v2v.
end());
6519 for (
int i = 0; i < J_v2v.
Size(); i++)
6521 e2v[J_v2v[i].two] = i;
6530 for (
int i = 0; i < NumOfVertices; i++)
6534 it.SetIndex(e2v[it.Index()]);
6542 const int oedge = NumOfVertices;
6543 const int oface = oedge + NumOfEdges;
6544 const int oelem = oface + NumOfQuadFaces;
6549 vertices.
SetSize(oelem + hex_counter);
6550 new_elements.
SetSize(8 * NumOfElements);
6551 CoarseFineTr.embeddings.SetSize(new_elements.
Size());
6554 for (
int i = 0, j = 0; i < NumOfElements; i++)
6557 const int attr = elements[i]->GetAttribute();
6558 int *v = elements[i]->GetVertices();
6559 const int *e = el_to_edge->GetRow(i);
6564 const int ne = el_to_edge->RowSize(i);
6565 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
6571 case Element::TETRAHEDRON:
6573 for (
int ei = 0; ei < 6; ei++)
6575 for (
int k = 0; k < 2; k++)
6577 vv[k] = v[tet_t::Edges[ei][k]];
6579 AverageVertices(vv, 2, oedge+e[ei]);
6585 const int rt_algo = 1;
6596 double len_sqr, min_len;
6598 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
6599 sqr(J(1,0)-J(1,1)-J(1,2)) +
6600 sqr(J(2,0)-J(2,1)-J(2,2));
6603 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
6604 sqr(J(1,1)-J(1,0)-J(1,2)) +
6605 sqr(J(2,1)-J(2,0)-J(2,2));
6606 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
6608 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
6609 sqr(J(1,2)-J(1,0)-J(1,1)) +
6610 sqr(J(2,2)-J(2,0)-J(2,1));
6611 if (len_sqr < min_len) { rt = 2; }
6616 double Em_data[18], Js_data[9], Jp_data[9];
6619 double ar1, ar2,
kappa, kappa_min;
6621 for (
int s = 0; s < 3; s++)
6623 for (
int t = 0; t < 3; t++)
6625 Em(t,s) = 0.5*J(t,s);
6628 for (
int t = 0; t < 3; t++)
6630 Em(t,3) = 0.5*(J(t,0)+J(t,1));
6631 Em(t,4) = 0.5*(J(t,0)+J(t,2));
6632 Em(t,5) = 0.5*(J(t,1)+J(t,2));
6636 for (
int t = 0; t < 3; t++)
6638 Js(t,0) = Em(t,5)-Em(t,0);
6639 Js(t,1) = Em(t,1)-Em(t,0);
6640 Js(t,2) = Em(t,2)-Em(t,0);
6644 for (
int t = 0; t < 3; t++)
6646 Js(t,0) = Em(t,5)-Em(t,0);
6647 Js(t,1) = Em(t,2)-Em(t,0);
6648 Js(t,2) = Em(t,4)-Em(t,0);
6652 kappa_min = std::max(ar1, ar2);
6656 for (
int t = 0; t < 3; t++)
6658 Js(t,0) = Em(t,0)-Em(t,1);
6659 Js(t,1) = Em(t,4)-Em(t,1);
6660 Js(t,2) = Em(t,2)-Em(t,1);
6664 for (
int t = 0; t < 3; t++)
6666 Js(t,0) = Em(t,2)-Em(t,1);
6667 Js(t,1) = Em(t,4)-Em(t,1);
6668 Js(t,2) = Em(t,5)-Em(t,1);
6672 kappa = std::max(ar1, ar2);
6673 if (kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
6676 for (
int t = 0; t < 3; t++)
6678 Js(t,0) = Em(t,0)-Em(t,2);
6679 Js(t,1) = Em(t,1)-Em(t,2);
6680 Js(t,2) = Em(t,3)-Em(t,2);
6684 for (
int t = 0; t < 3; t++)
6686 Js(t,0) = Em(t,1)-Em(t,2);
6687 Js(t,1) = Em(t,5)-Em(t,2);
6688 Js(t,2) = Em(t,3)-Em(t,2);
6692 kappa = std::max(ar1, ar2);
6693 if (kappa < kappa_min) { rt = 2; }
6696 static const int mv_all[3][4][4] =
6698 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
6699 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
6700 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
6702 const int (&mv)[4][4] = mv_all[rt];
6704 #ifndef MFEM_USE_MEMALLOC
6706 new Tetrahedron(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
6708 new Tetrahedron(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
6710 new Tetrahedron(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
6712 new Tetrahedron(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
6714 for (
int k = 0; k < 4; k++)
6716 new_elements[j+4+k] =
6717 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
6718 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
6722 new_elements[j+0] = tet = TetMemory.Alloc();
6723 tet->
Init(v[0], oedge+e[0], oedge+e[1], oedge+e[2], attr);
6725 new_elements[j+1] = tet = TetMemory.Alloc();
6726 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
6728 new_elements[j+2] = tet = TetMemory.Alloc();
6729 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
6731 new_elements[j+3] = tet = TetMemory.Alloc();
6732 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
6734 for (
int k = 0; k < 4; k++)
6736 new_elements[j+4+k] = tet = TetMemory.Alloc();
6737 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
6738 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
6741 for (
int k = 0; k < 4; k++)
6743 CoarseFineTr.embeddings[j+k].parent = i;
6744 CoarseFineTr.embeddings[j+k].matrix = k;
6746 for (
int k = 0; k < 4; k++)
6748 CoarseFineTr.embeddings[j+4+k].parent = i;
6749 CoarseFineTr.embeddings[j+4+k].matrix = 4*(rt+1)+k;
6756 case Element::WEDGE:
6758 const int *f = el_to_face->GetRow(i);
6760 for (
int fi = 2; fi < 5; fi++)
6762 for (
int k = 0; k < 4; k++)
6764 vv[k] = v[pri_t::FaceVert[fi][k]];
6766 AverageVertices(vv, 4, oface + f2qf[f[fi]]);
6769 for (
int ei = 0; ei < 9; ei++)
6771 for (
int k = 0; k < 2; k++)
6773 vv[k] = v[pri_t::Edges[ei][k]];
6775 AverageVertices(vv, 2, oedge+e[ei]);
6778 const int qf2 = f2qf[f[2]];
6779 const int qf3 = f2qf[f[3]];
6780 const int qf4 = f2qf[f[4]];
6783 new Wedge(v[0], oedge+e[0], oedge+e[2],
6784 oedge+e[6], oface+qf2, oface+qf4, attr);
6787 new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
6788 oface+qf3, oface+qf4, oface+qf2, attr);
6791 new Wedge(oedge+e[0], v[1], oedge+e[1],
6792 oface+qf2, oedge+e[7], oface+qf3, attr);
6795 new Wedge(oedge+e[2], oedge+e[1], v[2],
6796 oface+qf4, oface+qf3, oedge+e[8], attr);
6799 new Wedge(oedge+e[6], oface+qf2, oface+qf4,
6800 v[3], oedge+e[3], oedge+e[5], attr);
6803 new Wedge(oface+qf3, oface+qf4, oface+qf2,
6804 oedge+e[4], oedge+e[5], oedge+e[3], attr);
6807 new Wedge(oface+qf2, oedge+e[7], oface+qf3,
6808 oedge+e[3], v[4], oedge+e[4], attr);
6811 new Wedge(oface+qf4, oface+qf3, oedge+e[8],
6812 oedge+e[5], oedge+e[4], v[5], attr);
6816 case Element::HEXAHEDRON:
6818 const int *f = el_to_face->GetRow(i);
6819 const int he = hex_counter;
6824 if (f2qf.
Size() == 0)
6830 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[f[k]]; }
6834 AverageVertices(v, 8, oelem+he);
6836 for (
int fi = 0; fi < 6; fi++)
6838 for (
int k = 0; k < 4; k++)
6840 vv[k] = v[hex_t::FaceVert[fi][k]];
6842 AverageVertices(vv, 4, oface + qf[fi]);
6845 for (
int ei = 0; ei < 12; ei++)
6847 for (
int k = 0; k < 2; k++)
6849 vv[k] = v[hex_t::Edges[ei][k]];
6851 AverageVertices(vv, 2, oedge+e[ei]);
6855 new Hexahedron(v[0], oedge+e[0], oface+qf[0],
6856 oedge+e[3], oedge+e[8], oface+qf[1],
6857 oelem+he, oface+qf[4], attr);
6860 oface+qf[0], oface+qf[1], oedge+e[9],
6861 oface+qf[2], oelem+he, attr);
6863 new Hexahedron(oface+qf[0], oedge+e[1], v[2],
6864 oedge+e[2], oelem+he, oface+qf[2],
6865 oedge+e[10], oface+qf[3], attr);
6867 new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
6868 v[3], oface+qf[4], oelem+he,
6869 oface+qf[3], oedge+e[11], attr);
6871 new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
6872 oface+qf[4], v[4], oedge+e[4],
6873 oface+qf[5], oedge+e[7], attr);
6875 new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
6876 oelem+he, oedge+e[4], v[5],
6877 oedge+e[5], oface+qf[5], attr);
6879 new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
6880 oface+qf[3], oface+qf[5], oedge+e[5],
6881 v[6], oedge+e[6], attr);
6883 new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
6884 oedge+e[11], oedge+e[7], oface+qf[5],
6885 oedge+e[6], v[7], attr);
6890 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
6893 FreeElement(elements[i]);
6898 new_boundary.
SetSize(4 * NumOfBdrElements);
6899 for (
int i = 0, j = 0; i < NumOfBdrElements; i++)
6902 const int attr = boundary[i]->GetAttribute();
6903 int *v = boundary[i]->GetVertices();
6904 const int *e = bel_to_edge->GetRow(i);
6909 const int ne = bel_to_edge->RowSize(i);
6910 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
6914 if (bdr_el_type == Element::TRIANGLE)
6917 new Triangle(v[0], oedge+e[0], oedge+e[2], attr);
6919 new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
6921 new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
6923 new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
6925 else if (bdr_el_type == Element::QUADRILATERAL)
6928 (f2qf.
Size() == 0) ? be_to_face[i] : f2qf[be_to_face[i]];
6931 new Quadrilateral(v[0], oedge+e[0], oface+qf, oedge+e[3], attr);
6933 new Quadrilateral(oedge+e[0], v[1], oedge+e[1], oface+qf, attr);
6935 new Quadrilateral(oface+qf, oedge+e[1], v[2], oedge+e[2], attr);
6937 new Quadrilateral(oedge+e[3], oface+qf, oedge+e[2], v[3], attr);
6941 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
6943 FreeElement(boundary[i]);
6947 static const double A = 0.0, B = 0.5, C = 1.0;
6948 static double tet_children[3*4*16] =
6950 A,A,A, B,A,A, A,B,A, A,A,B,
6951 B,A,A, C,A,A, B,B,A, B,A,B,
6952 A,B,A, B,B,A, A,C,A, A,B,B,
6953 A,A,B, B,A,B, A,B,B, A,A,C,
6958 B,A,A, A,B,B, A,B,A, A,A,B,
6959 B,A,A, A,B,B, A,A,B, B,A,B,
6960 B,A,A, A,B,B, B,A,B, B,B,A,
6961 B,A,A, A,B,B, B,B,A, A,B,A,
6963 A,B,A, B,A,A, B,A,B, A,A,B,
6964 A,B,A, A,A,B, B,A,B, A,B,B,
6965 A,B,A, A,B,B, B,A,B, B,B,A,
6966 A,B,A, B,B,A, B,A,B, B,A,A,
6968 A,A,B, B,A,A, A,B,A, B,B,A,
6969 A,A,B, A,B,A, A,B,B, B,B,A,
6970 A,A,B, A,B,B, B,A,B, B,B,A,
6971 A,A,B, B,A,B, B,A,A, B,B,A
6973 static double pri_children[3*6*8] =
6975 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
6976 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
6977 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
6978 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
6979 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
6980 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
6981 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
6982 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
6984 static double hex_children[3*8*8] =
6986 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
6987 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
6988 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
6989 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
6990 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
6991 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
6992 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
6993 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
6996 CoarseFineTr.point_matrices[Geometry::TETRAHEDRON]
6997 .UseExternalData(tet_children, 3, 4, 16);
6998 CoarseFineTr.point_matrices[Geometry::PRISM]
6999 .UseExternalData(pri_children, 3, 6, 8);
7000 CoarseFineTr.point_matrices[Geometry::CUBE]
7001 .UseExternalData(hex_children, 3, 8, 8);
7003 for (
int i = 0; i < elements.Size(); i++)
7006 if (elements[i]->GetType() == Element::TETRAHEDRON) {
continue; }
7008 Embedding &emb = CoarseFineTr.embeddings[i];
7013 NumOfVertices = vertices.Size();
7014 NumOfElements = 8 * NumOfElements;
7015 NumOfBdrElements = 4 * NumOfBdrElements;
7017 GetElementToFaceTable();
7021 CheckBdrElementOrientation(
false);
7024 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7026 last_operation = Mesh::REFINE;
7029 if (update_nodes) { UpdateNodes(); }
7032 void Mesh::LocalRefinement(
const Array<int> &marked_el,
int type)
7034 int i, j, ind, nedges;
7041 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
7044 InitRefinementTransforms();
7048 int cne = NumOfElements, cnv = NumOfVertices;
7049 NumOfVertices += marked_el.
Size();
7050 NumOfElements += marked_el.
Size();
7051 vertices.SetSize(NumOfVertices);
7052 elements.SetSize(NumOfElements);
7053 CoarseFineTr.embeddings.SetSize(NumOfElements);
7055 for (j = 0; j < marked_el.
Size(); j++)
7060 int new_v = cnv + j, new_e = cne + j;
7061 AverageVertices(vert, 2, new_v);
7062 elements[new_e] =
new Segment(new_v, vert[1], attr);
7065 CoarseFineTr.embeddings[i] =
Embedding(i, 1);
7066 CoarseFineTr.embeddings[new_e] =
Embedding(i, 2);
7069 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
7070 CoarseFineTr.point_matrices[Geometry::SEGMENT].
7071 UseExternalData(seg_children, 1, 2, 3);
7079 DSTable v_to_v(NumOfVertices);
7080 GetVertexToVertexTable(v_to_v);
7084 int *edge1 =
new int[nedges];
7085 int *edge2 =
new int[nedges];
7086 int *middle =
new int[nedges];
7088 for (i = 0; i < nedges; i++)
7090 edge1[i] = edge2[i] = middle[i] = -1;
7093 for (i = 0; i < NumOfElements; i++)
7095 elements[i]->GetVertices(v);
7096 for (j = 1; j < v.
Size(); j++)
7098 ind = v_to_v(v[j-1], v[j]);
7099 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
7101 ind = v_to_v(v[0], v[v.
Size()-1]);
7102 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
7106 for (i = 0; i < marked_el.
Size(); i++)
7108 RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
7112 int need_refinement;
7115 need_refinement = 0;
7116 for (i = 0; i < nedges; i++)
7118 if (middle[i] != -1 && edge1[i] != -1)
7120 need_refinement = 1;
7121 GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
7125 while (need_refinement == 1);
7128 int v1[2], v2[2], bisect, temp;
7129 temp = NumOfBdrElements;
7130 for (i = 0; i < temp; i++)
7132 boundary[i]->GetVertices(v);
7133 bisect = v_to_v(v[0], v[1]);
7134 if (middle[bisect] != -1)
7136 if (boundary[i]->GetType() == Element::SEGMENT)
7138 v1[0] = v[0]; v1[1] = middle[bisect];
7139 v2[0] = middle[bisect]; v2[1] = v[1];
7141 boundary[i]->SetVertices(v1);
7142 boundary.
Append(
new Segment(v2, boundary[i]->GetAttribute()));
7145 mfem_error(
"Only bisection of segment is implemented"
7149 NumOfBdrElements = boundary.Size();
7156 if (el_to_edge != NULL)
7158 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7169 MFEM_VERIFY(GetNE() == 0 ||
7170 ((
Tetrahedron*)elements[0])->GetRefinementFlag() != 0,
7171 "tetrahedral mesh is not marked for refinement:"
7172 " call Finalize(true)");
7179 for (i = 0; i < marked_el.
Size(); i++)
7181 Bisection(marked_el[i], v_to_v);
7185 for (i = 0; i < marked_el.
Size(); i++)
7187 Bisection(marked_el[i], v_to_v);
7189 Bisection(NumOfElements - 1, v_to_v);
7190 Bisection(marked_el[i], v_to_v);
7194 for (i = 0; i < marked_el.
Size(); i++)
7196 Bisection(marked_el[i], v_to_v);
7198 ii = NumOfElements - 1;
7199 Bisection(ii, v_to_v);
7200 Bisection(NumOfElements - 1, v_to_v);
7201 Bisection(ii, v_to_v);
7203 Bisection(marked_el[i], v_to_v);
7204 Bisection(NumOfElements-1, v_to_v);
7205 Bisection(marked_el[i], v_to_v);
7211 int need_refinement;
7216 need_refinement = 0;
7219 for (i = 0; i < NumOfElements; i++)
7224 if (elements[i]->NeedRefinement(v_to_v))
7226 need_refinement = 1;
7227 Bisection(i, v_to_v);
7231 while (need_refinement == 1);
7238 need_refinement = 0;
7239 for (i = 0; i < NumOfBdrElements; i++)
7240 if (boundary[i]->NeedRefinement(v_to_v))
7242 need_refinement = 1;
7243 BdrBisection(i, v_to_v);
7246 while (need_refinement == 1);
7248 NumOfVertices = vertices.Size();
7249 NumOfBdrElements = boundary.Size();
7252 if (el_to_edge != NULL)
7254 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7256 if (el_to_face != NULL)
7258 GetElementToFaceTable();
7264 last_operation = Mesh::REFINE;
7270 CheckElementOrientation(
false);
7277 MFEM_VERIFY(!NURBSext,
"Nonconforming refinement of NURBS meshes is "
7278 "not supported. Project the NURBS to Nodes first.");
7285 ncmesh =
new NCMesh(
this);
7288 if (!refinements.
Size())
7290 last_operation = Mesh::NONE;
7295 ncmesh->MarkCoarseLevel();
7296 ncmesh->Refine(refinements);
7300 ncmesh->LimitNCLevel(nc_limit);
7305 ncmesh->OnMeshUpdated(mesh2);
7309 Swap(*mesh2,
false);
7312 GenerateNCFaceInfo();
7314 last_operation = Mesh::REFINE;
7319 Nodes->FESpace()->Update();
7325 const int *fine,
int nfine,
int op)
7328 for (
int i = 0; i < nfine; i++)
7330 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
7332 double err_fine = elem_error[fine[i]];
7335 case 0: error = std::min(error, err_fine);
break;
7336 case 1: error += err_fine;
break;
7337 case 2: error = std::max(error, err_fine);
break;
7344 double threshold,
int nc_limit,
int op)
7346 MFEM_VERIFY(ncmesh,
"Only supported for non-conforming meshes.");
7347 MFEM_VERIFY(!NURBSext,
"Derefinement of NURBS meshes is not supported. "
7348 "Project the NURBS to Nodes first.");
7352 const Table &dt = ncmesh->GetDerefinementTable();
7357 ncmesh->CheckDerefinementNCLevel(dt, level_ok, nc_limit);
7361 for (
int i = 0; i < dt.
Size(); i++)
7363 if (nc_limit > 0 && !level_ok[i]) {
continue; }
7366 AggregateError(elem_error, dt.
GetRow(i), dt.
RowSize(i), op);
7368 if (error < threshold) { derefs.
Append(i); }
7371 if (!derefs.
Size()) {
return false; }
7373 ncmesh->Derefine(derefs);
7376 ncmesh->OnMeshUpdated(mesh2);
7378 Swap(*mesh2,
false);
7381 GenerateNCFaceInfo();
7383 last_operation = Mesh::DEREFINE;
7392 int nc_limit,
int op)
7396 if (Nonconforming())
7398 return NonconformingDerefinement(elem_error, threshold, nc_limit, op);
7402 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
7408 bool Mesh::DerefineByError(
const Vector &elem_error,
double threshold,
7409 int nc_limit,
int op)
7412 for (
int i = 0; i < tmp.Size(); i++)
7414 tmp[i] = elem_error(i);
7416 return DerefineByError(tmp, threshold, nc_limit, op);
7420 void Mesh::InitFromNCMesh(
const NCMesh &ncmesh)
7429 NumOfVertices = vertices.Size();
7430 NumOfElements = elements.Size();
7431 NumOfBdrElements = boundary.Size();
7435 NumOfEdges = NumOfFaces = 0;
7436 nbInteriorFaces = nbBoundaryFaces = -1;
7440 el_to_edge =
new Table;
7441 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7445 GetElementToFaceTable();
7449 CheckBdrElementOrientation(
false);
7460 InitFromNCMesh(ncmesh);
7499 #ifdef MFEM_USE_MEMALLOC
7517 const int nv = Geometry::NumVerts[
geom];
7519 for (
int i = 0; i < elem_array.
Size(); i++)
7521 if (elem_array[i]->GetGeometryType() ==
geom)
7526 elem_vtx.
SetSize(nv*num_elems);
7530 for (
int i = 0; i < elem_array.
Size(); i++)
7536 elem_vtx.
Append(loc_vtx);
7544 for (
int i = 0; i < nelem; i++) { list[i] = i; }
7548 void Mesh::UniformRefinement(
int ref_algo)
7554 NURBSUniformRefinement();
7558 GeneralRefinement(AllElements(list, GetNE()));
7560 else if (ref_algo == 1 && meshgen == 1 && Dim == 3)
7563 LocalRefinement(AllElements(list, GetNE()));
7569 case 1: LocalRefinement(AllElements(list, GetNE()));
break;
7570 case 2: UniformRefinement2D();
break;
7571 case 3: UniformRefinement3D();
break;
7572 default: MFEM_ABORT(
"internal error");
7578 int nonconforming,
int nc_limit)
7584 else if (Dim == 1 || (Dim == 3 && (meshgen & 1)))
7588 else if (nonconforming < 0)
7591 if ((meshgen & 2) || (meshgen & 4))
7604 NonconformingRefinement(refinements, nc_limit);
7609 for (
int i = 0; i < refinements.
Size(); i++)
7611 el_to_refine[i] = refinements[i].index;
7615 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
7616 if (rt == 1 || rt == 2 || rt == 4)
7620 else if (rt == 3 || rt == 5 || rt == 6)
7630 LocalRefinement(el_to_refine, type);
7634 void Mesh::GeneralRefinement(
const Array<int> &el_to_refine,
int nonconforming,
7638 for (
int i = 0; i < el_to_refine.
Size(); i++)
7640 refinements[i] =
Refinement(el_to_refine[i]);
7642 GeneralRefinement(refinements, nonconforming, nc_limit);
7645 void Mesh::EnsureNCMesh(
bool simplices_nonconforming)
7647 MFEM_VERIFY(!NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
7648 "Project the NURBS to Nodes first.");
7652 if ((meshgen & 0x2) ||
7654 (simplices_nonconforming && (meshgen & 0x1)) )
7656 ncmesh =
new NCMesh(
this);
7657 ncmesh->OnMeshUpdated(
this);
7658 GenerateNCFaceInfo();
7663 void Mesh::RandomRefinement(
double prob,
bool aniso,
int nonconforming,
7667 for (
int i = 0; i < GetNE(); i++)
7669 if ((
double) rand() / RAND_MAX < prob)
7674 type = (Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
7679 GeneralRefinement(refs, nonconforming, nc_limit);
7682 void Mesh::RefineAtVertex(
const Vertex& vert,
double eps,
int nonconforming)
7686 for (
int i = 0; i < GetNE(); i++)
7688 GetElementVertices(i, v);
7689 bool refine =
false;
7690 for (
int j = 0; j < v.
Size(); j++)
7693 for (
int l = 0; l < spaceDim; l++)
7695 double d = vert(l) - vertices[v[j]](l);
7698 if (dist <= eps*eps) { refine =
true;
break; }
7705 GeneralRefinement(refs, nonconforming);
7709 int nonconforming,
int nc_limit)
7711 MFEM_VERIFY(elem_error.
Size() == GetNE(),
"");
7713 for (
int i = 0; i < GetNE(); i++)
7715 if (elem_error[i] > threshold)
7720 if (ReduceInt(refs.Size()))
7722 GeneralRefinement(refs, nonconforming, nc_limit);
7728 bool Mesh::RefineByError(
const Vector &elem_error,
double threshold,
7729 int nonconforming,
int nc_limit)
7733 return RefineByError(tmp, threshold, nonconforming, nc_limit);
7737 void Mesh::Bisection(
int i,
const DSTable &v_to_v,
7738 int *edge1,
int *edge2,
int *middle)
7741 int v[2][4], v_new, bisect, t;
7746 if (t == Element::TRIANGLE)
7753 bisect = v_to_v(vert[0], vert[1]);
7754 MFEM_ASSERT(bisect >= 0,
"");
7756 if (middle[bisect] == -1)
7758 v_new = NumOfVertices++;
7759 for (
int d = 0; d < spaceDim; d++)
7761 V(d) = 0.5 * (vertices[vert[0]](d) + vertices[vert[1]](d));
7767 if (edge1[bisect] == i)
7769 edge1[bisect] = edge2[bisect];
7772 middle[bisect] = v_new;
7776 v_new = middle[bisect];
7784 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
7785 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
7790 elements.Append(tri_new);
7799 int coarse = FindCoarseElement(i);
7800 CoarseFineTr.embeddings[i].parent = coarse;
7801 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7806 bisect = v_to_v(v[1][0], v[1][1]);
7807 MFEM_ASSERT(bisect >= 0,
"");
7809 if (edge1[bisect] == i)
7811 edge1[bisect] = NumOfElements;
7813 else if (edge2[bisect] == i)
7815 edge2[bisect] = NumOfElements;
7822 MFEM_ABORT(
"Bisection for now works only for triangles.");
7829 int v[2][4], v_new, bisect, t;
7834 if (t == Element::TETRAHEDRON)
7836 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
7840 "TETRAHEDRON element is not marked for refinement.");
7845 bisect = v_to_v.
FindId(vert[0], vert[1]);
7848 v_new = NumOfVertices + v_to_v.
GetId(vert[0],vert[1]);
7849 for (j = 0; j < 3; j++)
7851 V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
7857 v_new = NumOfVertices + bisect;
7866 new_redges[0][0] = 2;
7867 new_redges[0][1] = 1;
7868 new_redges[1][0] = 2;
7869 new_redges[1][1] = 1;
7870 int tr1 = -1, tr2 = -1;
7871 switch (old_redges[0])
7874 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
7875 if (type == Tetrahedron::TYPE_PF) { new_redges[0][1] = 4; }
7879 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
7883 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
7886 switch (old_redges[1])
7889 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
7890 if (type == Tetrahedron::TYPE_PF) { new_redges[1][0] = 3; }
7894 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
7898 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
7905 #ifdef MFEM_USE_MEMALLOC
7913 elements.Append(tet2);
7919 int coarse = FindCoarseElement(i);
7920 CoarseFineTr.embeddings[i].parent = coarse;
7921 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7926 case Tetrahedron::TYPE_PU:
7927 new_type = Tetrahedron::TYPE_PF;
break;
7928 case Tetrahedron::TYPE_PF:
7929 new_type = Tetrahedron::TYPE_A;
break;
7931 new_type = Tetrahedron::TYPE_PU;
7941 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
7948 int v[2][3], v_new, bisect, t;
7949 Element *bdr_el = boundary[i];
7952 if (t == Element::TRIANGLE)
7959 bisect = v_to_v.
FindId(vert[0], vert[1]);
7960 MFEM_ASSERT(bisect >= 0,
"");
7961 v_new = NumOfVertices + bisect;
7962 MFEM_ASSERT(v_new != -1,
"");
7966 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
7967 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
7977 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
7982 void Mesh::UniformRefinement(
int i,
const DSTable &v_to_v,
7983 int *edge1,
int *edge2,
int *middle)
7986 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
7989 if (elements[i]->GetType() == Element::TRIANGLE)
7995 bisect[0] = v_to_v(v[0],v[1]);
7996 bisect[1] = v_to_v(v[1],v[2]);
7997 bisect[2] = v_to_v(v[0],v[2]);
7998 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
8000 for (j = 0; j < 3; j++)
8002 if (middle[bisect[j]] == -1)
8004 v_new[j] = NumOfVertices++;
8005 for (
int d = 0; d < spaceDim; d++)
8007 V(d) = (vertices[v[j]](d) + vertices[v[(j+1)%3]](d))/2.;
8013 if (edge1[bisect[j]] == i)
8015 edge1[bisect[j]] = edge2[bisect[j]];
8018 middle[bisect[j]] = v_new[j];
8022 v_new[j] = middle[bisect[j]];
8025 edge1[bisect[j]] = -1;
8031 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
8032 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
8033 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
8034 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
8040 elements.Append(tri1);
8041 elements.Append(tri2);
8042 elements.Append(tri3);
8049 tri2->ResetTransform(code);
8050 tri3->ResetTransform(code);
8054 tri2->PushTransform(1);
8055 tri3->PushTransform(2);
8058 int coarse = FindCoarseElement(i);
8059 CoarseFineTr.embeddings[i] =
Embedding(coarse);
8060 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8061 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8062 CoarseFineTr.embeddings.Append(
Embedding(coarse));
8068 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
8072 void Mesh::InitRefinementTransforms()
8075 CoarseFineTr.Clear();
8076 CoarseFineTr.embeddings.SetSize(NumOfElements);
8077 for (
int i = 0; i < NumOfElements; i++)
8079 elements[i]->ResetTransform(0);
8080 CoarseFineTr.embeddings[i] =
Embedding(i);
8084 int Mesh::FindCoarseElement(
int i)
8087 while ((coarse = CoarseFineTr.embeddings[i].parent) != i)
8096 MFEM_VERIFY(GetLastOperation() == Mesh::REFINE,
"");
8100 return ncmesh->GetRefinementTransforms();
8104 for (
int i = 0; i < elem_geoms.
Size(); i++)
8107 if (CoarseFineTr.point_matrices[geom].SizeK()) {
continue; }
8109 if (geom == Geometry::TRIANGLE ||
8110 geom == Geometry::TETRAHEDRON)
8112 std::map<unsigned, int> mat_no;
8116 for (
int i = 0; i < elements.Size(); i++)
8119 unsigned code = elements[i]->GetTransform();
8122 int &matrix = mat_no[code];
8123 if (!matrix) { matrix = mat_no.size(); }
8126 CoarseFineTr.embeddings[i].matrix =
index;
8130 pmats.
SetSize(Dim, Dim+1, mat_no.size());
8133 std::map<unsigned, int>::iterator it;
8134 for (it = mat_no.begin(); it != mat_no.end(); ++it)
8136 if (geom == Geometry::TRIANGLE)
8138 Triangle::GetPointMatrix(it->first, pmats(it->second-1));
8142 Tetrahedron::GetPointMatrix(it->first, pmats(it->second-1));
8148 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
8149 " geom = " << geom);
8154 return CoarseFineTr;
8157 void Mesh::PrintXG(std::ostream &
out)
const
8159 MFEM_ASSERT(Dim==spaceDim,
"2D Manifold meshes not supported");
8168 out <<
"areamesh2\n\n";
8172 out <<
"curved_areamesh2\n\n";
8176 out << NumOfBdrElements <<
'\n';
8177 for (i = 0; i < NumOfBdrElements; i++)
8179 boundary[i]->GetVertices(v);
8181 out << boundary[i]->GetAttribute();
8182 for (j = 0; j < v.
Size(); j++)
8184 out <<
' ' << v[j] + 1;
8190 out << NumOfElements <<
'\n';
8191 for (i = 0; i < NumOfElements; i++)
8193 elements[i]->GetVertices(v);
8195 out << elements[i]->GetAttribute() <<
' ' << v.
Size();
8196 for (j = 0; j < v.
Size(); j++)
8198 out <<
' ' << v[j] + 1;
8206 out << NumOfVertices <<
'\n';
8207 for (i = 0; i < NumOfVertices; i++)
8209 out << vertices[i](0);
8210 for (j = 1; j < Dim; j++)
8212 out <<
' ' << vertices[i](j);
8219 out << NumOfVertices <<
'\n';
8227 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
8235 out <<
"NETGEN_Neutral_Format\n";
8237 out << NumOfVertices <<
'\n';
8238 for (i = 0; i < NumOfVertices; i++)
8240 for (j = 0; j < Dim; j++)
8242 out <<
' ' << vertices[i](j);
8248 out << NumOfElements <<
'\n';
8249 for (i = 0; i < NumOfElements; i++)
8251 nv = elements[i]->GetNVertices();
8252 ind = elements[i]->GetVertices();
8253 out << elements[i]->GetAttribute();
8254 for (j = 0; j < nv; j++)
8256 out <<
' ' << ind[j]+1;
8262 out << NumOfBdrElements <<
'\n';
8263 for (i = 0; i < NumOfBdrElements; i++)
8265 nv = boundary[i]->GetNVertices();
8266 ind = boundary[i]->GetVertices();
8267 out << boundary[i]->GetAttribute();
8268 for (j = 0; j < nv; j++)
8270 out <<
' ' << ind[j]+1;
8275 else if (meshgen == 2)
8281 <<
"1 " << NumOfVertices <<
" " << NumOfElements
8282 <<
" 0 0 0 0 0 0 0\n"
8283 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
8284 <<
"0 0 " << NumOfBdrElements <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
8285 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
8286 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
8288 for (i = 0; i < NumOfVertices; i++)
8289 out << i+1 <<
" 0.0 " << vertices[i](0) <<
' ' << vertices[i](1)
8290 <<
' ' << vertices[i](2) <<
" 0.0\n";
8292 for (i = 0; i < NumOfElements; i++)
8294 nv = elements[i]->GetNVertices();
8295 ind = elements[i]->GetVertices();
8296 out << i+1 <<
' ' << elements[i]->GetAttribute();
8297 for (j = 0; j < nv; j++)
8299 out <<
' ' << ind[j]+1;
8304 for (i = 0; i < NumOfBdrElements; i++)
8306 nv = boundary[i]->GetNVertices();
8307 ind = boundary[i]->GetVertices();
8308 out << boundary[i]->GetAttribute();
8309 for (j = 0; j < nv; j++)
8311 out <<
' ' << ind[j]+1;
8313 out <<
" 1.0 1.0 1.0 1.0\n";
8321 void Mesh::Printer(std::ostream &
out, std::string section_delimiter)
const
8328 NURBSext->Print(out);
8339 out << (ncmesh ?
"MFEM mesh v1.1\n" :
8340 section_delimiter.empty() ?
"MFEM mesh v1.0\n" :
8341 "MFEM mesh v1.2\n");
8345 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8350 "# TETRAHEDRON = 4\n"
8355 out <<
"\ndimension\n" << Dim
8356 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8357 for (i = 0; i < NumOfElements; i++)
8359 PrintElement(elements[i], out);
8362 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
8363 for (i = 0; i < NumOfBdrElements; i++)
8365 PrintElement(boundary[i], out);
8370 out <<
"\nvertex_parents\n";
8371 ncmesh->PrintVertexParents(out);
8373 out <<
"\ncoarse_elements\n";
8374 ncmesh->PrintCoarseElements(out);
8377 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8380 out << spaceDim <<
'\n';
8381 for (i = 0; i < NumOfVertices; i++)
8383 out << vertices[i](0);
8384 for (j = 1; j < spaceDim; j++)
8386 out <<
' ' << vertices[i](j);
8398 if (!ncmesh && !section_delimiter.empty())
8400 out << section_delimiter << endl;
8409 out <<
"MFEM NURBS mesh v1.0\n";
8413 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8419 out <<
"\ndimension\n" << Dim
8420 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8421 for (i = 0; i < NumOfElements; i++)
8423 PrintElement(elements[i], out);
8426 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
8427 for (i = 0; i < NumOfBdrElements; i++)
8429 PrintElement(boundary[i], out);
8432 out <<
"\nedges\n" << NumOfEdges <<
'\n';
8433 for (i = 0; i < NumOfEdges; i++)
8435 edge_vertex->GetRow(i, vert);
8441 out << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
8443 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8446 void Mesh::PrintVTK(std::ostream &
out)
8449 "# vtk DataFile Version 3.0\n"
8450 "Generated by MFEM\n"
8452 "DATASET UNSTRUCTURED_GRID\n";
8456 out <<
"POINTS " << NumOfVertices <<
" double\n";
8457 for (
int i = 0; i < NumOfVertices; i++)
8459 out << vertices[i](0);
8461 for (j = 1; j < spaceDim; j++)
8463 out <<
' ' << vertices[i](j);
8475 out <<
"POINTS " << Nodes->FESpace()->GetNDofs() <<
" double\n";
8476 for (
int i = 0; i < Nodes->FESpace()->GetNDofs(); i++)
8480 Nodes->FESpace()->DofsToVDofs(vdofs);
8481 out << (*Nodes)(vdofs[0]);
8483 for (j = 1; j < spaceDim; j++)
8485 out <<
' ' << (*Nodes)(vdofs[j]);
8499 for (
int i = 0; i < NumOfElements; i++)
8501 size += elements[i]->GetNVertices() + 1;
8503 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
8504 for (
int i = 0; i < NumOfElements; i++)
8506 const int *v = elements[i]->GetVertices();
8507 const int nv = elements[i]->GetNVertices();
8509 for (
int j = 0; j < nv; j++)
8521 for (
int i = 0; i < NumOfElements; i++)
8523 Nodes->FESpace()->GetElementDofs(i, dofs);
8524 MFEM_ASSERT(Dim != 0 || dofs.
Size() == 1,
8525 "Point meshes should have a single dof per element");
8526 size += dofs.
Size() + 1;
8528 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
8529 const char *fec_name = Nodes->FESpace()->FEColl()->Name();
8531 if (!strcmp(fec_name,
"Linear") ||
8532 !strcmp(fec_name,
"H1_0D_P1") ||
8533 !strcmp(fec_name,
"H1_1D_P1") ||
8534 !strcmp(fec_name,
"H1_2D_P1") ||
8535 !strcmp(fec_name,
"H1_3D_P1"))
8539 else if (!strcmp(fec_name,
"Quadratic") ||
8540 !strcmp(fec_name,
"H1_1D_P2") ||
8541 !strcmp(fec_name,
"H1_2D_P2") ||
8542 !strcmp(fec_name,
"H1_3D_P2"))
8548 mfem::err <<
"Mesh::PrintVTK : can not save '"
8549 << fec_name <<
"' elements!" << endl;
8552 for (
int i = 0; i < NumOfElements; i++)
8554 Nodes->FESpace()->GetElementDofs(i, dofs);
8558 for (
int j = 0; j < dofs.
Size(); j++)
8560 out <<
' ' << dofs[j];
8563 else if (order == 2)
8565 const int *vtk_mfem;
8566 switch (elements[i]->GetGeometryType())
8568 case Geometry::SEGMENT:
8569 case Geometry::TRIANGLE:
8570 case Geometry::SQUARE:
8571 vtk_mfem = vtk_quadratic_hex;
break;
8572 case Geometry::TETRAHEDRON:
8573 vtk_mfem = vtk_quadratic_tet;
break;
8574 case Geometry::PRISM:
8575 vtk_mfem = vtk_quadratic_wedge;
break;
8576 case Geometry::CUBE:
8578 vtk_mfem = vtk_quadratic_hex;
break;
8580 for (
int j = 0; j < dofs.
Size(); j++)
8582 out <<
' ' << dofs[vtk_mfem[j]];
8589 out <<
"CELL_TYPES " << NumOfElements <<
'\n';
8590 for (
int i = 0; i < NumOfElements; i++)
8592 int vtk_cell_type = 5;
8598 case Geometry::POINT: vtk_cell_type = 1;
break;
8599 case Geometry::SEGMENT: vtk_cell_type = 3;
break;
8600 case Geometry::TRIANGLE: vtk_cell_type = 5;
break;
8601 case Geometry::SQUARE: vtk_cell_type = 9;
break;
8602 case Geometry::TETRAHEDRON: vtk_cell_type = 10;
break;
8603 case Geometry::CUBE: vtk_cell_type = 12;
break;
8604 case Geometry::PRISM: vtk_cell_type = 13;
break;
8608 else if (order == 2)
8612 case Geometry::SEGMENT: vtk_cell_type = 21;
break;
8613 case Geometry::TRIANGLE: vtk_cell_type = 22;
break;
8614 case Geometry::SQUARE: vtk_cell_type = 28;
break;
8615 case Geometry::TETRAHEDRON: vtk_cell_type = 24;
break;
8616 case Geometry::CUBE: vtk_cell_type = 29;
break;
8617 case Geometry::PRISM: vtk_cell_type = 32;
break;
8622 out << vtk_cell_type <<
'\n';
8626 out <<
"CELL_DATA " << NumOfElements <<
'\n'
8627 <<
"SCALARS material int\n"
8628 <<
"LOOKUP_TABLE default\n";
8629 for (
int i = 0; i < NumOfElements; i++)
8631 out << elements[i]->GetAttribute() <<
'\n';
8636 void Mesh::PrintVTU(std::string fname,
8638 bool high_order_output,
8639 int compression_level)
8641 int ref = (high_order_output && Nodes) ? Nodes->FESpace()->GetOrder(0) : 1;
8642 fname = fname +
".vtu";
8644 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
8645 if (compression_level != 0)
8647 out <<
" compressor=\"vtkZLibDataCompressor\"";
8650 out <<
"<UnstructuredGrid>\n";
8651 PrintVTU(out, ref, format, high_order_output, compression_level);
8652 out <<
"</Piece>\n";
8653 out <<
"</UnstructuredGrid>\n";
8654 out <<
"</VTKFile>" << std::endl;
8659 template <
typename T>
8663 if (format == VTKFormat::ASCII) { out << val << suffix; }
8670 const uint8_t &val,
const char *suffix,
8673 if (format == VTKFormat::ASCII) { out << static_cast<int>(val) << suffix; }
8679 const double &val,
const char *suffix,
8682 if (format == VTKFormat::BINARY32)
8684 bin_io::AppendBytes<float>(buf, float(val));
8686 else if (format == VTKFormat::BINARY)
8692 out << val << suffix;
8698 const float &val,
const char *suffix,
8701 if (format == VTKFormat::BINARY) { bin_io::AppendBytes<double>(buf, val); }
8703 else {
out << val << suffix; }
8707 int compression_level)
8715 bool high_order_output,
int compression_level)
8720 const char *fmt_str = (format == VTKFormat::ASCII) ?
"ascii" :
"binary";
8721 const char *type_str = (format != VTKFormat::BINARY32) ?
"Float64" :
"Float32";
8722 std::vector<char> buf;
8725 int np = 0, nc_ref = 0, size = 0;
8726 for (
int i = 0; i < GetNE(); i++)
8736 out <<
"<Piece NumberOfPoints=\"" << np <<
"\" NumberOfCells=\""
8737 << (high_order_output ? GetNE() : nc_ref) <<
"\">\n";
8740 out <<
"<Points>\n";
8741 out <<
"<DataArray type=\"" << type_str
8742 <<
"\" NumberOfComponents=\"3\" format=\"" << fmt_str <<
"\">\n";
8743 for (
int i = 0; i < GetNE(); i++)
8746 GetElementBaseGeometry(i), ref, 1);
8748 GetElementTransformation(i)->Transform(RefG->
RefPts, pmat);
8750 for (
int j = 0; j < pmat.
Width(); j++)
8769 if (format == VTKFormat::ASCII) { out <<
'\n'; }
8772 if (format != VTKFormat::ASCII)
8776 out <<
"</DataArray>" << std::endl;
8777 out <<
"</Points>" << std::endl;
8779 out <<
"<Cells>" << std::endl;
8780 out <<
"<DataArray type=\"Int32\" Name=\"connectivity\" format=\""
8781 << fmt_str <<
"\">" << std::endl;
8783 std::vector<int> offset;
8786 if (high_order_output)
8789 for (
int iel = 0; iel < GetNE(); iel++)
8793 int nnodes = local_connectivity.
Size();
8794 for (
int i=0; i<nnodes; ++i)
8798 if (format == VTKFormat::ASCII) { out <<
'\n'; }
8800 offset.push_back(np);
8806 for (
int i = 0; i < GetNE(); i++)
8812 for (
int j = 0; j < RG.
Size(); )
8816 offset.push_back(coff);
8818 for (
int k = 0; k < nv; k++, j++)
8822 if (format == VTKFormat::ASCII) { out <<
'\n'; }
8827 if (format != VTKFormat::ASCII)
8831 out <<
"</DataArray>" << std::endl;
8833 out <<
"<DataArray type=\"Int32\" Name=\"offsets\" format=\""
8834 << fmt_str <<
"\">" << std::endl;
8836 for (
size_t ii=0; ii<offset.size(); ii++)
8840 if (format != VTKFormat::ASCII)
8844 out <<
"</DataArray>" << std::endl;
8845 out <<
"<DataArray type=\"UInt8\" Name=\"types\" format=\""
8846 << fmt_str <<
"\">" << std::endl;
8848 for (
int i = 0; i < GetNE(); i++)
8851 uint8_t vtk_cell_type = 5;
8856 case Geometry::POINT:
8859 case Geometry::SEGMENT:
8860 vtk_cell_type = high_order_output ? 68 : 3;
8862 case Geometry::TRIANGLE:
8863 vtk_cell_type = high_order_output ? 69 : 5;
8865 case Geometry::SQUARE:
8866 vtk_cell_type = high_order_output ? 70 : 9;
8868 case Geometry::TETRAHEDRON:
8869 vtk_cell_type = high_order_output ? 71 : 10;
8871 case Geometry::CUBE:
8872 vtk_cell_type = high_order_output ? 72 : 12;
8874 case Geometry::PRISM:
8875 vtk_cell_type = high_order_output ? 73 : 13;
8878 MFEM_ABORT(
"Unrecognized VTK element type \"" << geom <<
"\"");
8882 if (high_order_output)
8891 for (
int j = 0; j < RG.
Size(); j += nv)
8897 if (format != VTKFormat::ASCII)
8901 out <<
"</DataArray>" << std::endl;
8902 out <<
"</Cells>" << std::endl;
8904 out <<
"<CellData Scalars=\"material\">" << std::endl;
8905 out <<
"<DataArray type=\"Int32\" Name=\"material\" format=\""
8906 << fmt_str <<
"\">" << std::endl;
8907 for (
int i = 0; i < GetNE(); i++)
8909 int attr = GetAttribute(i);
8910 if (high_order_output)
8925 if (format != VTKFormat::ASCII)
8929 out <<
"</DataArray>" << std::endl;
8930 out <<
"</CellData>" << std::endl;
8934 void Mesh::PrintVTK(std::ostream &
out,
int ref,
int field_data)
8941 "# vtk DataFile Version 3.0\n"
8942 "Generated by MFEM\n"
8944 "DATASET UNSTRUCTURED_GRID\n";
8949 out <<
"FIELD FieldData 1\n"
8950 <<
"MaterialIds " << 1 <<
" " << attributes.
Size() <<
" int\n";
8951 for (
int i = 0; i < attributes.Size(); i++)
8953 out <<
' ' << attributes[i];
8960 for (
int i = 0; i < GetNE(); i++)
8969 out <<
"POINTS " << np <<
" double\n";
8971 for (
int i = 0; i < GetNE(); i++)
8974 GetElementBaseGeometry(i), ref, 1);
8976 GetElementTransformation(i)->Transform(RefG->
RefPts, pmat);
8978 for (
int j = 0; j < pmat.
Width(); j++)
8980 out << pmat(0, j) <<
' ';
8983 out << pmat(1, j) <<
' ';
8995 out << 0.0 <<
' ' << 0.0;
9002 out <<
"CELLS " << nc <<
' ' << size <<
'\n';
9004 for (
int i = 0; i < GetNE(); i++)
9011 for (
int j = 0; j < RG.
Size(); )
9014 for (
int k = 0; k < nv; k++, j++)
9016 out <<
' ' << np + RG[j];
9022 out <<
"CELL_TYPES " << nc <<
'\n';
9023 for (
int i = 0; i < GetNE(); i++)
9029 int vtk_cell_type = 5;
9033 case Geometry::POINT: vtk_cell_type = 1;
break;
9034 case Geometry::SEGMENT: vtk_cell_type = 3;
break;
9035 case Geometry::TRIANGLE: vtk_cell_type = 5;
break;
9036 case Geometry::SQUARE: vtk_cell_type = 9;
break;
9037 case Geometry::TETRAHEDRON: vtk_cell_type = 10;
break;
9038 case Geometry::CUBE: vtk_cell_type = 12;
break;
9039 case Geometry::PRISM: vtk_cell_type = 13;
break;
9041 MFEM_ABORT(
"Unrecognized VTK element type \"" << geom <<
"\"");
9045 for (
int j = 0; j < RG.
Size(); j += nv)
9047 out << vtk_cell_type <<
'\n';
9051 out <<
"CELL_DATA " << nc <<
'\n'
9052 <<
"SCALARS material int\n"
9053 <<
"LOOKUP_TABLE default\n";
9054 for (
int i = 0; i < GetNE(); i++)
9059 int attr = GetAttribute(i);
9062 out << attr <<
'\n';
9069 srand((
unsigned)time(0));
9070 double a = double(rand()) / (double(RAND_MAX) + 1.);
9071 int el0 = (int)floor(a * GetNE());
9072 GetElementColoring(coloring, el0);
9073 out <<
"SCALARS element_coloring int\n"
9074 <<
"LOOKUP_TABLE default\n";
9075 for (
int i = 0; i < GetNE(); i++)
9082 out << coloring[i] + 1 <<
'\n';
9088 out <<
"POINT_DATA " << np <<
'\n' << flush;
9093 int delete_el_to_el = (el_to_el) ? (0) : (1);
9094 const Table &el_el = ElementToElementTable();
9095 int num_el = GetNE(), stack_p, stack_top_p, max_num_col;
9098 const int *i_el_el = el_el.
GetI();
9099 const int *j_el_el = el_el.
GetJ();
9104 stack_p = stack_top_p = 0;
9105 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
9107 if (colors[el] != -2)
9113 el_stack[stack_top_p++] = el;
9115 for ( ; stack_p < stack_top_p; stack_p++)
9117 int i = el_stack[stack_p];
9118 int num_nb = i_el_el[i+1] - i_el_el[i];
9119 if (max_num_col < num_nb + 1)
9121 max_num_col = num_nb + 1;
9123 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
9126 if (colors[k] == -2)
9129 el_stack[stack_top_p++] = k;
9137 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
9139 int i = el_stack[stack_p], col;
9141 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
9143 col = colors[j_el_el[j]];
9146 col_marker[col] = 1;
9150 for (col = 0; col < max_num_col; col++)
9151 if (col_marker[col] == 0)
9159 if (delete_el_to_el)
9166 void Mesh::PrintWithPartitioning(
int *partitioning, std::ostream &
out,
9167 int elem_attr)
const
9169 if (Dim != 3 && Dim != 2) {
return; }
9171 int i, j, k, l, nv, nbe, *v;
9173 out <<
"MFEM mesh v1.0\n";
9177 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
9182 "# TETRAHEDRON = 4\n"
9187 out <<
"\ndimension\n" << Dim
9188 <<
"\n\nelements\n" << NumOfElements <<
'\n';
9189 for (i = 0; i < NumOfElements; i++)
9191 out << int((elem_attr) ? partitioning[i]+1 : elements[i]->GetAttribute())
9192 <<
' ' << elements[i]->GetGeometryType();
9193 nv = elements[i]->GetNVertices();
9194 v = elements[i]->GetVertices();
9195 for (j = 0; j < nv; j++)
9202 for (i = 0; i < faces_info.Size(); i++)
9204 if ((l = faces_info[i].Elem2No) >= 0)
9206 k = partitioning[faces_info[i].Elem1No];
9207 l = partitioning[l];
9211 if (!Nonconforming() || !IsSlaveFace(faces_info[i]))
9222 out <<
"\nboundary\n" << nbe <<
'\n';
9223 for (i = 0; i < faces_info.Size(); i++)
9225 if ((l = faces_info[i].Elem2No) >= 0)
9227 k = partitioning[faces_info[i].Elem1No];
9228 l = partitioning[l];
9231 nv = faces[i]->GetNVertices();
9232 v = faces[i]->GetVertices();
9233 out << k+1 <<
' ' << faces[i]->GetGeometryType();
9234 for (j = 0; j < nv; j++)
9239 if (!Nonconforming() || !IsSlaveFace(faces_info[i]))
9241 out << l+1 <<
' ' << faces[i]->GetGeometryType();
9242 for (j = nv-1; j >= 0; j--)
9252 k = partitioning[faces_info[i].Elem1No];
9253 nv = faces[i]->GetNVertices();
9254 v = faces[i]->GetVertices();
9255 out << k+1 <<
' ' << faces[i]->GetGeometryType();
9256 for (j = 0; j < nv; j++)
9263 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
9266 out << spaceDim <<
'\n';
9267 for (i = 0; i < NumOfVertices; i++)
9269 out << vertices[i](0);
9270 for (j = 1; j < spaceDim; j++)
9272 out <<
' ' << vertices[i](j);
9285 void Mesh::PrintElementsWithPartitioning(
int *partitioning,
9289 MFEM_ASSERT(Dim == spaceDim,
"2D Manifolds not supported\n");
9290 if (Dim != 3 && Dim != 2) {
return; }
9297 int *vcount =
new int[NumOfVertices];
9298 for (i = 0; i < NumOfVertices; i++)
9302 for (i = 0; i < NumOfElements; i++)
9304 nv = elements[i]->GetNVertices();
9305 ind = elements[i]->GetVertices();
9306 for (j = 0; j < nv; j++)
9312 int *voff =
new int[NumOfVertices+1];
9314 for (i = 1; i <= NumOfVertices; i++)
9316 voff[i] = vcount[i-1] + voff[i-1];
9319 int **vown =
new int*[NumOfVertices];
9320 for (i = 0; i < NumOfVertices; i++)
9322 vown[i] =
new int[vcount[i]];
9332 Transpose(ElementToEdgeTable(), edge_el);
9335 for (i = 0; i < NumOfElements; i++)
9337 nv = elements[i]->GetNVertices();
9338 ind = elements[i]->GetVertices();
9339 for (j = 0; j < nv; j++)
9342 vown[ind[j]][vcount[ind[j]]] = i;
9346 for (i = 0; i < NumOfVertices; i++)
9348 vcount[i] = voff[i+1] - voff[i];
9352 for (i = 0; i < edge_el.
Size(); i++)
9354 const int *el = edge_el.
GetRow(i);
9357 k = partitioning[el[0]];
9358 l = partitioning[el[1]];
9359 if (interior_faces || k != l)
9371 out <<
"areamesh2\n\n" << nbe <<
'\n';
9373 for (i = 0; i < edge_el.
Size(); i++)
9375 const int *el = edge_el.
GetRow(i);
9378 k = partitioning[el[0]];
9379 l = partitioning[el[1]];
9380 if (interior_faces || k != l)
9383 GetEdgeVertices(i,ev);
9385 for (j = 0; j < 2; j++)
9386 for (s = 0; s < vcount[ev[j]]; s++)
9387 if (vown[ev[j]][s] == el[0])
9389 out <<
' ' << voff[ev[j]]+s+1;
9393 for (j = 1; j >= 0; j--)
9394 for (s = 0; s < vcount[ev[j]]; s++)
9395 if (vown[ev[j]][s] == el[1])
9397 out <<
' ' << voff[ev[j]]+s+1;
9404 k = partitioning[el[0]];
9406 GetEdgeVertices(i,ev);
9408 for (j = 0; j < 2; j++)
9409 for (s = 0; s < vcount[ev[j]]; s++)
9410 if (vown[ev[j]][s] == el[0])
9412 out <<
' ' << voff[ev[j]]+s+1;
9419 out << NumOfElements <<
'\n';
9420 for (i = 0; i < NumOfElements; i++)
9422 nv = elements[i]->GetNVertices();
9423 ind = elements[i]->GetVertices();
9424 out << partitioning[i]+1 <<
' ';
9426 for (j = 0; j < nv; j++)
9428 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
9429 vown[ind[j]][vcount[ind[j]]] = i;
9434 for (i = 0; i < NumOfVertices; i++)
9436 vcount[i] = voff[i+1] - voff[i];
9440 out << voff[NumOfVertices] <<
'\n';
9441 for (i = 0; i < NumOfVertices; i++)
9442 for (k = 0; k < vcount[i]; k++)
9444 for (j = 0; j < Dim; j++)
9446 out << vertices[i](j) <<
' ';
9452 else if (meshgen == 1)
9454 out <<
"NETGEN_Neutral_Format\n";
9456 out << voff[NumOfVertices] <<
'\n';
9457 for (i = 0; i < NumOfVertices; i++)
9458 for (k = 0; k < vcount[i]; k++)
9460 for (j = 0; j < Dim; j++)
9462 out <<
' ' << vertices[i](j);
9468 out << NumOfElements <<
'\n';
9469 for (i = 0; i < NumOfElements; i++)
9471 nv = elements[i]->GetNVertices();
9472 ind = elements[i]->GetVertices();
9473 out << partitioning[i]+1;
9474 for (j = 0; j < nv; j++)
9476 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
9477 vown[ind[j]][vcount[ind[j]]] = i;
9482 for (i = 0; i < NumOfVertices; i++)
9484 vcount[i] = voff[i+1] - voff[i];
9490 for (i = 0; i < NumOfFaces; i++)
9491 if ((l = faces_info[i].Elem2No) >= 0)
9493 k = partitioning[faces_info[i].Elem1No];
9494 l = partitioning[l];
9495 if (interior_faces || k != l)
9506 for (i = 0; i < NumOfFaces; i++)
9507 if ((l = faces_info[i].Elem2No) >= 0)
9509 k = partitioning[faces_info[i].Elem1No];
9510 l = partitioning[l];
9511 if (interior_faces || k != l)
9513 nv = faces[i]->GetNVertices();
9514 ind = faces[i]->GetVertices();
9516 for (j = 0; j < nv; j++)
9517 for (s = 0; s < vcount[ind[j]]; s++)
9518 if (vown[ind[j]][s] == faces_info[i].Elem1No)
9520 out <<
' ' << voff[ind[j]]+s+1;
9524 for (j = nv-1; j >= 0; j--)
9525 for (s = 0; s < vcount[ind[j]]; s++)
9526 if (vown[ind[j]][s] == faces_info[i].Elem2No)
9528 out <<
' ' << voff[ind[j]]+s+1;
9535 k = partitioning[faces_info[i].Elem1No];
9536 nv = faces[i]->GetNVertices();
9537 ind = faces[i]->GetVertices();
9539 for (j = 0; j < nv; j++)
9540 for (s = 0; s < vcount[ind[j]]; s++)
9541 if (vown[ind[j]][s] == faces_info[i].Elem1No)
9543 out <<
' ' << voff[ind[j]]+s+1;
9549 else if (meshgen == 2)
9554 for (i = 0; i < NumOfFaces; i++)
9555 if ((l = faces_info[i].Elem2No) >= 0)
9557 k = partitioning[faces_info[i].Elem1No];
9558 l = partitioning[l];
9559 if (interior_faces || k != l)
9571 <<
"1 " << voff[NumOfVertices] <<
" " << NumOfElements
9572 <<
" 0 0 0 0 0 0 0\n"
9573 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
9574 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
9575 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
9576 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
9578 for (i = 0; i < NumOfVertices; i++)
9579 for (k = 0; k < vcount[i]; k++)
9580 out << voff[i]+k <<
" 0.0 " << vertices[i](0) <<
' '
9581 << vertices[i](1) <<
' ' << vertices[i](2) <<
" 0.0\n";
9583 for (i = 0; i < NumOfElements; i++)
9585 nv = elements[i]->GetNVertices();
9586 ind = elements[i]->GetVertices();
9587 out << i+1 <<
' ' << partitioning[i]+1;
9588 for (j = 0; j < nv; j++)
9590 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
9591 vown[ind[j]][vcount[ind[j]]] = i;
9596 for (i = 0; i < NumOfVertices; i++)
9598 vcount[i] = voff[i+1] - voff[i];
9602 for (i = 0; i < NumOfFaces; i++)
9603 if ((l = faces_info[i].Elem2No) >= 0)
9605 k = partitioning[faces_info[i].Elem1No];
9606 l = partitioning[l];
9607 if (interior_faces || k != l)
9609 nv = faces[i]->GetNVertices();
9610 ind = faces[i]->GetVertices();
9612 for (j = 0; j < nv; j++)
9613 for (s = 0; s < vcount[ind[j]]; s++)
9614 if (vown[ind[j]][s] == faces_info[i].Elem1No)
9616 out <<
' ' << voff[ind[j]]+s+1;
9618 out <<
" 1.0 1.0 1.0 1.0\n";
9620 for (j = nv-1; j >= 0; j--)
9621 for (s = 0; s < vcount[ind[j]]; s++)
9622 if (vown[ind[j]][s] == faces_info[i].Elem2No)
9624 out <<
' ' << voff[ind[j]]+s+1;
9626 out <<
" 1.0 1.0 1.0 1.0\n";
9631 k = partitioning[faces_info[i].Elem1No];
9632 nv = faces[i]->GetNVertices();
9633 ind = faces[i]->GetVertices();
9635 for (j = 0; j < nv; j++)
9636 for (s = 0; s < vcount[ind[j]]; s++)
9637 if (vown[ind[j]][s] == faces_info[i].Elem1No)
9639 out <<
' ' << voff[ind[j]]+s+1;
9641 out <<
" 1.0 1.0 1.0 1.0\n";
9647 for (i = 0; i < NumOfVertices; i++)
9657 void Mesh::PrintSurfaces(
const Table & Aface_face, std::ostream &
out)
const
9664 " NURBS mesh is not supported!");
9668 out <<
"MFEM mesh v1.0\n";
9672 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
9677 "# TETRAHEDRON = 4\n"
9682 out <<
"\ndimension\n" << Dim
9683 <<
"\n\nelements\n" << NumOfElements <<
'\n';
9684 for (i = 0; i < NumOfElements; i++)
9686 PrintElement(elements[i], out);
9690 const int *
const i_AF_f = Aface_face.
GetI();
9691 const int *
const j_AF_f = Aface_face.
GetJ();
9693 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
9694 for (
const int * iface = j_AF_f + i_AF_f[iAF];
9695 iface < j_AF_f + i_AF_f[iAF+1];
9698 out << iAF+1 <<
' ';
9699 PrintElementWithoutAttr(faces[*iface],out);
9702 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
9705 out << spaceDim <<
'\n';
9706 for (i = 0; i < NumOfVertices; i++)
9708 out << vertices[i](0);
9709 for (j = 1; j < spaceDim; j++)
9711 out <<
' ' << vertices[i](j);
9724 void Mesh::ScaleSubdomains(
double sf)
9729 int na = attributes.
Size();
9730 double *cg =
new double[na*spaceDim];
9731 int *nbea =
new int[na];
9733 int *vn =
new int[NumOfVertices];
9734 for (i = 0; i < NumOfVertices; i++)
9738 for (i = 0; i < na; i++)
9740 for (j = 0; j < spaceDim; j++)
9742 cg[i*spaceDim+j] = 0.0;
9747 for (i = 0; i < NumOfElements; i++)
9749 GetElementVertices(i, vert);
9750 for (k = 0; k < vert.
Size(); k++)
9756 for (i = 0; i < NumOfElements; i++)
9758 int bea = GetAttribute(i)-1;
9759 GetPointMatrix(i, pointmat);
9760 GetElementVertices(i, vert);
9762 for (k = 0; k < vert.
Size(); k++)
9763 if (vn[vert[k]] == 1)
9766 for (j = 0; j < spaceDim; j++)
9768 cg[bea*spaceDim+j] += pointmat(j,k);
9774 for (i = 0; i < NumOfElements; i++)
9776 int bea = GetAttribute(i)-1;
9777 GetElementVertices (i, vert);
9779 for (k = 0; k < vert.
Size(); k++)
9782 for (j = 0; j < spaceDim; j++)
9783 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
9784 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
9794 void Mesh::ScaleElements(
double sf)
9799 int na = NumOfElements;
9800 double *cg =
new double[na*spaceDim];
9801 int *nbea =
new int[na];
9803 int *vn =
new int[NumOfVertices];
9804 for (i = 0; i < NumOfVertices; i++)
9808 for (i = 0; i < na; i++)
9810 for (j = 0; j < spaceDim; j++)
9812 cg[i*spaceDim+j] = 0.0;
9817 for (i = 0; i < NumOfElements; i++)
9819 GetElementVertices(i, vert);
9820 for (k = 0; k < vert.
Size(); k++)
9826 for (i = 0; i < NumOfElements; i++)
9829 GetPointMatrix(i, pointmat);
9830 GetElementVertices(i, vert);
9832 for (k = 0; k < vert.
Size(); k++)
9833 if (vn[vert[k]] == 1)
9836 for (j = 0; j < spaceDim; j++)
9838 cg[bea*spaceDim+j] += pointmat(j,k);
9844 for (i = 0; i < NumOfElements; i++)
9847 GetElementVertices(i, vert);
9849 for (k = 0; k < vert.
Size(); k++)
9852 for (j = 0; j < spaceDim; j++)
9853 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
9854 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
9869 Vector vold(spaceDim), vnew(NULL, spaceDim);
9870 for (
int i = 0; i < vertices.Size(); i++)
9872 for (
int j = 0; j < spaceDim; j++)
9874 vold(j) = vertices[i](j);
9884 xnew.ProjectCoefficient(f_pert);
9891 MFEM_VERIFY(spaceDim == deformation.
GetVDim(),
9892 "incompatible vector dimensions");
9899 for (
int i = 0; i < NumOfVertices; i++)
9900 for (
int d = 0; d < spaceDim; d++)
9902 vertices[i](d) = xnew(d + spaceDim*i);
9913 void Mesh::RemoveUnusedVertices()
9915 if (NURBSext || ncmesh) {
return; }
9919 for (
int i = 0; i < GetNE(); i++)
9924 for (
int j = 0; j < nv; j++)
9929 for (
int i = 0; i < GetNBE(); i++)
9931 Element *el = GetBdrElement(i);
9934 for (
int j = 0; j < nv; j++)
9940 for (
int i = 0; i < v2v.
Size(); i++)
9944 vertices[num_vert] = vertices[i];
9945 v2v[i] = num_vert++;
9949 if (num_vert == v2v.
Size()) {
return; }
9956 for (
int i = 0; i < GetNE(); i++)
9958 Nodes->FESpace()->GetElementVDofs(i, vdofs);
9963 for (
int i = 0; i < GetNE(); i++)
9965 Nodes->FESpace()->GetElementVDofs(i, vdofs);
9966 Nodes->GetSubVector(vdofs, &nodes_by_element(s));
9970 vertices.SetSize(num_vert);
9971 NumOfVertices = num_vert;
9972 for (
int i = 0; i < GetNE(); i++)
9977 for (
int j = 0; j < nv; j++)
9982 for (
int i = 0; i < GetNBE(); i++)
9984 Element *el = GetBdrElement(i);
9987 for (
int j = 0; j < nv; j++)
9996 el_to_edge =
new Table;
9997 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
10002 GetElementToFaceTable();
10008 Nodes->FESpace()->Update();
10011 for (
int i = 0; i < GetNE(); i++)
10013 Nodes->FESpace()->GetElementVDofs(i, vdofs);
10014 Nodes->SetSubVector(vdofs, &nodes_by_element(s));
10020 void Mesh::RemoveInternalBoundaries()
10022 if (NURBSext || ncmesh) {
return; }
10024 int num_bdr_elem = 0;
10025 int new_bel_to_edge_nnz = 0;
10026 for (
int i = 0; i < GetNBE(); i++)
10028 if (FaceIsInterior(GetBdrElementEdgeIndex(i)))
10030 FreeElement(boundary[i]);
10037 new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
10042 if (num_bdr_elem == GetNBE()) {
return; }
10046 Table *new_bel_to_edge = NULL;
10050 new_be_to_edge.
Reserve(num_bdr_elem);
10054 new_be_to_face.
Reserve(num_bdr_elem);
10055 new_bel_to_edge =
new Table;
10056 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
10058 for (
int i = 0; i < GetNBE(); i++)
10060 if (!FaceIsInterior(GetBdrElementEdgeIndex(i)))
10062 new_boundary.
Append(boundary[i]);
10065 new_be_to_edge.
Append(be_to_edge[i]);
10069 int row = new_be_to_face.
Size();
10070 new_be_to_face.
Append(be_to_face[i]);
10071 int *e = bel_to_edge->GetRow(i);
10072 int ne = bel_to_edge->RowSize(i);
10073 int *new_e = new_bel_to_edge->
GetRow(row);
10074 for (
int j = 0; j < ne; j++)
10078 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
10083 NumOfBdrElements = new_boundary.
Size();
10093 delete bel_to_edge;
10094 bel_to_edge = new_bel_to_edge;
10098 for (
int i = 0; i < attribs.
Size(); i++)
10100 attribs[i] = GetBdrAttribute(i);
10104 bdr_attributes.DeleteAll();
10105 attribs.
Copy(bdr_attributes);
10110 #ifdef MFEM_USE_MEMALLOC
10113 if (E->
GetType() == Element::TETRAHEDRON)
10137 const int npts = point_mat.
Width();
10138 if (!npts) {
return 0; }
10139 MFEM_VERIFY(point_mat.
Height() == spaceDim,
"Invalid points matrix");
10143 if (!GetNE()) {
return 0; }
10145 double *data = point_mat.
GetData();
10152 min_dist = std::numeric_limits<double>::max();
10156 for (
int i = 0; i < GetNE(); i++)
10158 GetElementTransformation(i)->Transform(
10160 for (
int k = 0; k < npts; k++)
10162 double dist = pt.
DistanceTo(data+k*spaceDim);
10163 if (dist < min_dist(k))
10165 min_dist(k) = dist;
10174 for (
int k = 0; k < npts; k++)
10178 int res = inv_tr->
Transform(pt, ips[k]);
10179 if (res == InverseElementTransformation::Inside)
10181 elem_ids[k] = e_idx[k];
10185 if (pts_found != npts)
10188 Table *vtoel = GetVertexToElementTable();
10189 for (
int k = 0; k < npts; k++)
10191 if (elem_ids[k] != -1) {
continue; }
10194 GetElementVertices(e_idx[k], vertices);
10195 for (
int v = 0; v < vertices.
Size(); v++)
10197 int vv = vertices[v];
10199 const int* els = vtoel->
GetRow(vv);
10200 for (
int e = 0; e < ne; e++)
10202 if (els[e] == e_idx[k]) {
continue; }
10204 int res = inv_tr->
Transform(pt, ips[k]);
10205 if (res == InverseElementTransformation::Inside)
10207 elem_ids[k] = els[e];
10217 int le = ncmesh->leaf_elements[e_idx[k]];
10218 ncmesh->FindNeighbors(le,neigh);
10219 for (
int e = 0; e < neigh.
Size(); e++)
10222 if (ncmesh->IsGhost(ncmesh->elements[nn])) {
continue; }
10223 int el = ncmesh->elements[nn].index;
10225 int res = inv_tr->
Transform(pt, ips[k]);
10226 if (res == InverseElementTransformation::Inside)
10238 if (inv_trans == NULL) {
delete inv_tr; }
10240 if (warn && pts_found != npts)
10242 MFEM_WARNING((npts-pts_found) <<
" points were not found");
10253 computed_factors = flags;
10258 const int vdim = fespace->
GetVDim();
10259 const int NE = fespace->
GetNE();
10260 const int ND = fe->
GetDof();
10265 ElementDofOrdering::NATIVE);
10267 unsigned eval_flags = 0;
10268 if (flags & GeometricFactors::COORDINATES)
10270 X.SetSize(vdim*NQ*NE);
10271 eval_flags |= QuadratureInterpolator::VALUES;
10273 if (flags & GeometricFactors::JACOBIANS)
10275 J.SetSize(vdim*vdim*NQ*NE);
10276 eval_flags |= QuadratureInterpolator::DERIVATIVES;
10278 if (flags & GeometricFactors::DETERMINANTS)
10280 detJ.SetSize(NQ*NE);
10281 eval_flags |= QuadratureInterpolator::DETERMINANTS;
10290 Vector Enodes(vdim*ND*NE);
10291 elem_restr->
Mult(*nodes, Enodes);
10292 qi->
Mult(Enodes, eval_flags, X, J, detJ);
10296 qi->
Mult(*nodes, eval_flags, X, J, detJ);
10300 FaceGeometricFactors::FaceGeometricFactors(
const Mesh *mesh,
10311 const int vdim = fespace->
GetVDim();
10320 face_restr->
Mult(*nodes, Fnodes);
10322 unsigned eval_flags = 0;
10363 V(1) = s * ((ip.
y + layer) / n);
10368 V(2) = s * ((ip.
z + layer) / n);
10377 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
10381 int nvy = (closed) ? (ny) : (ny + 1);
10382 int nvt = mesh->
GetNV() * nvy;
10391 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
10396 for (
int i = 0; i < mesh->
GetNV(); i++)
10399 for (
int j = 0; j < nvy; j++)
10401 vc[1] = sy * (double(j) / ny);
10407 for (
int i = 0; i < mesh->
GetNE(); i++)
10412 for (
int j = 0; j < ny; j++)
10415 qv[0] = vert[0] * nvy + j;
10416 qv[1] = vert[1] * nvy + j;
10417 qv[2] = vert[1] * nvy + (j + 1) % nvy;
10418 qv[3] = vert[0] * nvy + (j + 1) % nvy;
10424 for (
int i = 0; i < mesh->
GetNBE(); i++)
10429 for (
int j = 0; j < ny; j++)
10432 sv[0] = vert[0] * nvy + j;
10433 sv[1] = vert[0] * nvy + (j + 1) % nvy;
10437 Swap<int>(sv[0], sv[1]);
10449 for (
int i = 0; i < mesh->
GetNE(); i++)
10455 sv[0] = vert[0] * nvy;
10456 sv[1] = vert[1] * nvy;
10460 sv[0] = vert[1] * nvy + ny;
10461 sv[1] = vert[0] * nvy + ny;
10477 string cname = name;
10478 if (cname ==
"Linear")
10482 else if (cname ==
"Quadratic")
10486 else if (cname ==
"Cubic")
10490 else if (!strncmp(name,
"H1_", 3))
10494 else if (!strncmp(name,
"L2_T", 4))
10498 else if (!strncmp(name,
"L2_", 3))
10505 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
10517 for (
int i = 0; i < mesh->
GetNE(); i++)
10520 for (
int j = ny-1; j >= 0; j--)
10537 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
10542 int nvt = mesh->
GetNV() * nvz;
10547 bool wdgMesh =
false;
10548 bool hexMesh =
false;
10552 for (
int i = 0; i < mesh->
GetNV(); i++)
10556 for (
int j = 0; j < nvz; j++)
10558 vc[2] = sz * (double(j) / nz);
10564 for (
int i = 0; i < mesh->
GetNE(); i++)
10574 for (
int j = 0; j < nz; j++)
10577 pv[0] = vert[0] * nvz + j;
10578 pv[1] = vert[1] * nvz + j;
10579 pv[2] = vert[2] * nvz + j;
10580 pv[3] = vert[0] * nvz + (j + 1) % nvz;
10581 pv[4] = vert[1] * nvz + (j + 1) % nvz;
10582 pv[5] = vert[2] * nvz + (j + 1) % nvz;
10589 for (
int j = 0; j < nz; j++)
10592 hv[0] = vert[0] * nvz + j;
10593 hv[1] = vert[1] * nvz + j;
10594 hv[2] = vert[2] * nvz + j;
10595 hv[3] = vert[3] * nvz + j;
10596 hv[4] = vert[0] * nvz + (j + 1) % nvz;
10597 hv[5] = vert[1] * nvz + (j + 1) % nvz;
10598 hv[6] = vert[2] * nvz + (j + 1) % nvz;
10599 hv[7] = vert[3] * nvz + (j + 1) % nvz;
10601 mesh3d->
AddHex(hv, attr);
10605 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
10606 << geom <<
"\'" << endl;
10612 for (
int i = 0; i < mesh->
GetNBE(); i++)
10617 for (
int j = 0; j < nz; j++)
10620 qv[0] = vert[0] * nvz + j;
10621 qv[1] = vert[1] * nvz + j;
10622 qv[2] = vert[1] * nvz + (j + 1) % nvz;
10623 qv[3] = vert[0] * nvz + (j + 1) % nvz;
10632 for (
int i = 0; i < mesh->
GetNE(); i++)
10643 tv[0] = vert[0] * nvz;
10644 tv[1] = vert[2] * nvz;
10645 tv[2] = vert[1] * nvz;
10649 tv[0] = vert[0] * nvz + nz;
10650 tv[1] = vert[1] * nvz + nz;
10651 tv[2] = vert[2] * nvz + nz;
10659 qv[0] = vert[0] * nvz;
10660 qv[1] = vert[3] * nvz;
10661 qv[2] = vert[2] * nvz;
10662 qv[3] = vert[1] * nvz;
10666 qv[0] = vert[0] * nvz + nz;
10667 qv[1] = vert[1] * nvz + nz;
10668 qv[2] = vert[2] * nvz + nz;
10669 qv[3] = vert[3] * nvz + nz;
10675 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
10676 << geom <<
"\'" << endl;
10682 if ( hexMesh && wdgMesh )
10686 else if ( hexMesh )
10690 else if ( wdgMesh )
10703 string cname = name;
10704 if (cname ==
"Linear")
10708 else if (cname ==
"Quadratic")
10712 else if (cname ==
"Cubic")
10716 else if (!strncmp(name,
"H1_", 3))
10720 else if (!strncmp(name,
"L2_T", 4))
10724 else if (!strncmp(name,
"L2_", 3))
10731 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
10743 for (
int i = 0; i < mesh->
GetNE(); i++)
10746 for (
int j = nz-1; j >= 0; j--)
10767 out << i <<
" " << v[0] <<
" " << v[1] <<
" " << v[2]
10768 <<
" 0 0 " << i <<
" -1 0\n";
10775 double mid[3] = {0, 0, 0};
10776 for (
int j = 0; j < 2; j++)
10778 for (
int k = 0; k <
spaceDim; k++)
10783 out << NumOfVertices+i <<
" "
10784 << mid[0]/2 <<
" " << mid[1]/2 <<
" " << mid[2]/2 <<
" "
10785 << ev[0] <<
" " << ev[1] <<
" -1 " << i <<
" 0\n";
int GetNPoints() const
Returns the number of the points in the integration rule.
Abstract class for Finite Elements.
void AddHex(const int *vi, int attr=1)
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
std::ostream & operator<<(std::ostream &out, const Mesh &mesh)
int Size() const
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 GetNDofs() const
Returns number of degrees of freedom.
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
virtual void Update(bool want_transform=true)
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 ...
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
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 current array.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
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)
Piecewise-(bi)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.
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.
void AddBdrSegment(const int *vi, int attr=1)
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
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)
void DeleteAll()
Delete 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...
void AddVertex(const double *)
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.
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
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.
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 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.
std::vector< Master > masters
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
TriLinear3DFiniteElement HexahedronFE
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()
const IntegrationRule & GetNodes() const
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
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.
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. 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>.
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.
const NURBSExtension * GetNURBSext() const
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
int SpaceDimension() const
std::vector< Slave > slaves
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...
const IntegrationRule * IntRule
Vector X
Mapped (physical) coordinates of all quadrature points.
void AddAColumnInRow(int r)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void PartialSum()
Partial Sum.
int FindRoots(const Vector &z, Vector &x)
void AddQuad(const int *vi, int attr=1)
void AddWedge(const int *vi, int attr=1)
int Push4(int r, int c, int f, int t)
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.
void AddBdrTriangle(const int *vi, int attr=1)
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...
void AddBdrQuad(const int *vi, int attr=1)
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)
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()
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 associated with i'th element.
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)
NCMesh * ncmesh
Optional non-conforming mesh extension.
Geometry::Type GetBdrElementBaseGeometry(int i) const
void filter_dos(std::string &line)
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)
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...
BiLinear2DFiniteElement QuadrilateralFE
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.
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
DenseMatrix point_matrix
position within the master edge/face
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.
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...