15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
17 #include "../general/text.hpp"
18 #include "../general/device.hpp"
31 #if defined(MFEM_USE_METIS) && defined(MFEM_USE_METIS_5)
36 #if defined(MFEM_USE_METIS) && !defined(MFEM_USE_METIS_5)
41 int*,
int*,
int*,
int*,
int*,
idxtype*);
43 int*,
int*,
int*,
int*,
int*,
idxtype*);
45 int*,
int*,
int*,
int*,
int*,
idxtype*);
66 void Mesh::GetElementCenter(
int i,
Vector ¢)
69 int geom = GetElementBaseGeometry(i);
74 double Mesh::GetElementSize(
int i,
int type)
77 GetElementJacobian(i, J);
80 return pow(fabs(J.
Det()), 1./Dim);
92 double Mesh::GetElementSize(
int i,
const Vector &dir)
96 GetElementJacobian(i, J);
98 return sqrt((d_hat * d_hat) / (dir * dir));
101 double Mesh::GetElementVolume(
int i)
123 for (
int d = 0; d < spaceDim; d++)
132 for (
int i = 0; i < NumOfVertices; i++)
134 coord = GetVertex(i);
135 for (
int d = 0; d < spaceDim; d++)
137 if (coord[d] < min(d)) { min(d) = coord[d]; }
138 if (coord[d] > max(d)) { max(d) = coord[d]; }
144 const bool use_boundary =
false;
145 int ne = use_boundary ? GetNBE() : GetNE();
153 for (
int i = 0; i < ne; i++)
157 GetBdrElementFace(i, &fn, &fo);
159 Tr = GetFaceElementTransformations(fn, 5);
166 T = GetElementTransformation(i);
170 for (
int j = 0; j < pointmat.
Width(); j++)
172 for (
int d = 0; d < pointmat.
Height(); d++)
174 if (pointmat(d,j) < min(d)) { min(d) = pointmat(d,j); }
175 if (pointmat(d,j) > max(d)) { max(d) = pointmat(d,j); }
182 void Mesh::GetCharacteristics(
double &h_min,
double &h_max,
183 double &kappa_min,
double &kappa_max,
191 sdim = SpaceDimension();
193 if (Vh) { Vh->
SetSize(NumOfElements); }
194 if (Vk) { Vk->
SetSize(NumOfElements); }
197 h_max = kappa_max = -h_min;
198 if (dim == 0) {
if (Vh) { *Vh = 1.0; }
if (Vk) {*Vk = 1.0; }
return; }
200 for (i = 0; i < NumOfElements; i++)
202 GetElementJacobian(i, J);
203 h = pow(fabs(J.
Weight()), 1.0/
double(dim));
204 kappa = (dim == sdim) ?
206 if (Vh) { (*Vh)(i) = h; }
207 if (Vk) { (*Vk)(i) = kappa; }
209 if (h < h_min) { h_min = h; }
210 if (h > h_max) { h_max = h; }
211 if (kappa < kappa_min) { kappa_min =
kappa; }
212 if (kappa > kappa_max) { kappa_max =
kappa; }
217 void Mesh::PrintElementsByGeometry(
int dim,
221 for (
int g = Geometry::DimStart[dim], first = 1;
222 g < Geometry::DimStart[dim+1]; g++)
224 if (!num_elems_by_geom[g]) {
continue; }
225 if (!first) { out <<
" + "; }
227 out << num_elems_by_geom[g] <<
' ' << Geometry::Name[g] <<
"(s)";
233 double h_min, h_max, kappa_min, kappa_max;
235 out <<
"Mesh Characteristics:";
237 this->GetCharacteristics(h_min, h_max, kappa_min, kappa_max, Vh, Vk);
239 Array<int> num_elems_by_geom(Geometry::NumGeom);
240 num_elems_by_geom = 0;
241 for (
int i = 0; i < GetNE(); i++)
243 num_elems_by_geom[GetElementBaseGeometry(i)]++;
247 <<
"Dimension : " << Dimension() <<
'\n'
248 <<
"Space dimension : " << SpaceDimension();
252 <<
"Number of vertices : " << GetNV() <<
'\n'
253 <<
"Number of elements : " << GetNE() <<
'\n'
254 <<
"Number of bdr elem : " << GetNBE() <<
'\n';
259 <<
"Number of vertices : " << GetNV() <<
'\n'
260 <<
"Number of elements : " << GetNE() <<
'\n'
261 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
262 <<
"h_min : " << h_min <<
'\n'
263 <<
"h_max : " << h_max <<
'\n';
268 <<
"Number of vertices : " << GetNV() <<
'\n'
269 <<
"Number of edges : " << GetNEdges() <<
'\n'
270 <<
"Number of elements : " << GetNE() <<
" -- ";
271 PrintElementsByGeometry(2, num_elems_by_geom, out);
273 <<
"Number of bdr elem : " << GetNBE() <<
'\n'
274 <<
"Euler Number : " << EulerNumber2D() <<
'\n'
275 <<
"h_min : " << h_min <<
'\n'
276 <<
"h_max : " << h_max <<
'\n'
277 <<
"kappa_min : " << kappa_min <<
'\n'
278 <<
"kappa_max : " << kappa_max <<
'\n';
282 Array<int> num_bdr_elems_by_geom(Geometry::NumGeom);
283 num_bdr_elems_by_geom = 0;
284 for (
int i = 0; i < GetNBE(); i++)
286 num_bdr_elems_by_geom[GetBdrElementBaseGeometry(i)]++;
288 Array<int> num_faces_by_geom(Geometry::NumGeom);
289 num_faces_by_geom = 0;
290 for (
int i = 0; i < GetNFaces(); i++)
292 num_faces_by_geom[GetFaceBaseGeometry(i)]++;
296 <<
"Number of vertices : " << GetNV() <<
'\n'
297 <<
"Number of edges : " << GetNEdges() <<
'\n'
298 <<
"Number of faces : " << GetNFaces() <<
" -- ";
299 PrintElementsByGeometry(Dim-1, num_faces_by_geom, out);
301 <<
"Number of elements : " << GetNE() <<
" -- ";
302 PrintElementsByGeometry(Dim, num_elems_by_geom, out);
304 <<
"Number of bdr elem : " << GetNBE() <<
" -- ";
305 PrintElementsByGeometry(Dim-1, num_bdr_elems_by_geom, out);
307 <<
"Euler Number : " << EulerNumber() <<
'\n'
308 <<
"h_min : " << h_min <<
'\n'
309 <<
"h_max : " << h_max <<
'\n'
310 <<
"kappa_min : " << kappa_min <<
'\n'
311 <<
"kappa_max : " << kappa_max <<
'\n';
313 out <<
'\n' << std::flush;
320 case Element::POINT :
return &
PointFE;
321 case Element::SEGMENT :
return &
SegmentFE;
326 case Element::WEDGE :
return &
WedgeFE;
328 MFEM_ABORT(
"Unknown element type \"" << ElemType <<
"\"");
331 MFEM_ABORT(
"Unknown element type");
343 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
349 Nodes->FESpace()->GetElementVDofs(i, vdofs);
352 int n = vdofs.
Size()/spaceDim;
354 for (
int k = 0; k < spaceDim; k++)
356 for (
int j = 0; j < n; j++)
358 pm(k,j) = nodes(vdofs[n*k+j]);
361 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
366 void Mesh::GetElementTransformation(
int i,
const Vector &nodes,
374 MFEM_ASSERT(nodes.
Size() == spaceDim*GetNV(),
"");
375 int nv = elements[i]->GetNVertices();
376 const int *v = elements[i]->GetVertices();
377 int n = vertices.
Size();
379 for (
int k = 0; k < spaceDim; k++)
381 for (
int j = 0; j < nv; j++)
383 pm(k, j) = nodes(k*n+v[j]);
386 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
390 MFEM_ASSERT(nodes.
Size() == Nodes->Size(),
"");
392 Nodes->FESpace()->GetElementVDofs(i, vdofs);
393 int n = vdofs.
Size()/spaceDim;
395 for (
int k = 0; k < spaceDim; k++)
397 for (
int j = 0; j < n; j++)
399 pm(k,j) = nodes(vdofs[n*k+j]);
402 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
409 GetElementTransformation(i, &Transformation);
411 return &Transformation;
416 GetBdrElementTransformation(i, &BdrTransformation);
417 return &BdrTransformation;
427 GetBdrPointMatrix(i, pm);
428 ElTr->
SetFE(GetTransformationFEforElementType(GetBdrElementType(i)));
436 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
437 int n = vdofs.
Size()/spaceDim;
439 for (
int k = 0; k < spaceDim; k++)
441 for (
int j = 0; j < n; j++)
443 pm(k,j) = (*Nodes)(vdofs[n*k+j]);
450 int elem_id, face_info;
451 GetBdrElementAdjacentElement(i, elem_id, face_info);
453 GetLocalFaceTransformation(GetBdrElementType(i),
454 GetElementType(elem_id),
455 FaceElemTr.Loc1.Transf, face_info);
459 Nodes->FESpace()->GetTraceElement(elem_id,
460 GetBdrElementBaseGeometry(i));
463 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
465 Nodes->GetVectorValues(Transformation, eir, pm);
467 ElTr->
SetFE(face_el);
475 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
480 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
481 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
483 for (
int i = 0; i < spaceDim; i++)
485 for (
int j = 0; j < nv; j++)
487 pm(i, j) = vertices[v[j]](i);
490 FTr->
SetFE(GetTransformationFEforElementType(GetFaceElementType(FaceNo)));
494 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
498 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
499 int n = vdofs.
Size()/spaceDim;
501 for (
int i = 0; i < spaceDim; i++)
503 for (
int j = 0; j < n; j++)
505 pm(i, j) = (*Nodes)(vdofs[n*i+j]);
512 FaceInfo &face_info = faces_info[FaceNo];
517 GetLocalFaceTransformation(face_type,
518 GetElementType(face_info.
Elem1No),
519 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
522 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
526 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
528 Nodes->GetVectorValues(Transformation, eir, pm);
538 GetFaceTransformation(FaceNo, &FaceTransformation);
539 return &FaceTransformation;
546 GetFaceTransformation(EdgeNo, EdTr);
551 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
560 GetEdgeVertices(EdgeNo, v);
563 for (
int i = 0; i < spaceDim; i++)
565 for (
int j = 0; j < nv; j++)
567 pm(i, j) = vertices[v[j]](i);
570 EdTr->
SetFE(GetTransformationFEforElementType(Element::SEGMENT));
574 const FiniteElement *edge_el = Nodes->FESpace()->GetEdgeElement(EdgeNo);
578 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
579 int n = vdofs.
Size()/spaceDim;
581 for (
int i = 0; i < spaceDim; i++)
583 for (
int j = 0; j < n; j++)
585 pm(i, j) = (*Nodes)(vdofs[n*i+j]);
588 EdTr->
SetFE(edge_el);
592 MFEM_ABORT(
"Not implemented.");
600 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
601 return &EdgeTransformation;
605 void Mesh::GetLocalPtToSegTransformation(
620 void Mesh::GetLocalSegToTriTransformation(
628 tv = tri_t::Edges[i/64];
629 so = seg_t::Orient[i%64];
632 for (
int j = 0; j < 2; j++)
634 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
635 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
640 void Mesh::GetLocalSegToQuadTransformation(
648 qv = quad_t::Edges[i/64];
649 so = seg_t::Orient[i%64];
652 for (
int j = 0; j < 2; j++)
654 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
655 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
660 void Mesh::GetLocalTriToTetTransformation(
667 const int *tv = tet_t::FaceVert[i/64];
670 const int *to = tri_t::Orient[i%64];
674 for (
int j = 0; j < 3; j++)
677 locpm(0, j) = vert.
x;
678 locpm(1, j) = vert.
y;
679 locpm(2, j) = vert.
z;
684 void Mesh::GetLocalTriToWdgTransformation(
691 MFEM_VERIFY(i < 128,
"Local face index " << i/64
692 <<
" is not a triangular face of a wedge.");
693 const int *pv = pri_t::FaceVert[i/64];
696 const int *to = tri_t::Orient[i%64];
700 for (
int j = 0; j < 3; j++)
703 locpm(0, j) = vert.
x;
704 locpm(1, j) = vert.
y;
705 locpm(2, j) = vert.
z;
710 void Mesh::GetLocalQuadToHexTransformation(
717 const int *hv = hex_t::FaceVert[i/64];
719 const int *qo = quad_t::Orient[i%64];
722 for (
int j = 0; j < 4; j++)
725 locpm(0, j) = vert.
x;
726 locpm(1, j) = vert.
y;
727 locpm(2, j) = vert.
z;
732 void Mesh::GetLocalQuadToWdgTransformation(
739 MFEM_VERIFY(i >= 128,
"Local face index " << i/64
740 <<
" is not a quadrilateral face of a wedge.");
741 const int *pv = pri_t::FaceVert[i/64];
743 const int *qo = quad_t::Orient[i%64];
746 for (
int j = 0; j < 4; j++)
749 locpm(0, j) = vert.
x;
750 locpm(1, j) = vert.
y;
751 locpm(2, j) = vert.
z;
759 for (
int i = 0; i < geom_factors.Size(); i++)
771 geom_factors.Append(gf);
775 void Mesh::DeleteGeometricFactors()
777 for (
int i = 0; i < geom_factors.Size(); i++)
779 delete geom_factors[i];
781 geom_factors.SetSize(0);
784 void Mesh::GetLocalFaceTransformation(
790 GetLocalPtToSegTransformation(Transf, info);
793 case Element::SEGMENT:
794 if (elem_type == Element::TRIANGLE)
796 GetLocalSegToTriTransformation(Transf, info);
800 MFEM_ASSERT(elem_type == Element::QUADRILATERAL,
"");
801 GetLocalSegToQuadTransformation(Transf, info);
805 case Element::TRIANGLE:
806 if (elem_type == Element::TETRAHEDRON)
808 GetLocalTriToTetTransformation(Transf, info);
812 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
813 GetLocalTriToWdgTransformation(Transf, info);
817 case Element::QUADRILATERAL:
818 if (elem_type == Element::HEXAHEDRON)
820 GetLocalQuadToHexTransformation(Transf, info);
824 MFEM_ASSERT(elem_type == Element::WEDGE,
"");
825 GetLocalQuadToWdgTransformation(Transf, info);
834 FaceInfo &face_info = faces_info[FaceNo];
836 FaceElemTr.Elem1 = NULL;
837 FaceElemTr.Elem2 = NULL;
843 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
844 FaceElemTr.Elem1 = &Transformation;
850 FaceElemTr.Elem2No = face_info.
Elem2No;
851 if ((mask & 2) && FaceElemTr.Elem2No >= 0)
854 if (NURBSext && (mask & 1)) { MFEM_ABORT(
"NURBS mesh not supported!"); }
856 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
857 FaceElemTr.Elem2 = &Transformation2;
861 FaceElemTr.FaceGeom = GetFaceGeometryType(FaceNo);
862 FaceElemTr.Face = (mask & 16) ? GetFaceTransformation(FaceNo) : NULL;
865 int face_type = GetFaceElementType(FaceNo);
868 int elem_type = GetElementType(face_info.
Elem1No);
869 GetLocalFaceTransformation(face_type, elem_type,
870 FaceElemTr.Loc1.Transf, face_info.
Elem1Inf);
872 if ((mask & 8) && FaceElemTr.Elem2No >= 0)
874 int elem_type = GetElementType(face_info.
Elem2No);
875 GetLocalFaceTransformation(face_type, elem_type,
876 FaceElemTr.Loc2.Transf, face_info.
Elem2Inf);
879 if (Nonconforming() && IsSlaveFace(face_info))
881 ApplyLocalSlaveTransformation(FaceElemTr.Loc2.Transf, face_info);
883 if (face_type == Element::SEGMENT)
886 DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
887 std::swap(pm(0,0), pm(0,1));
888 std::swap(pm(1,0), pm(1,1));
898 return fi.
NCFace >= 0 && nc_faces_info[fi.
NCFace].Slave;
904 #ifdef MFEM_THREAD_SAFE
909 MFEM_ASSERT(fi.
NCFace >= 0,
"");
921 fn = be_to_face[BdrElemNo];
925 fn = be_to_edge[BdrElemNo];
929 fn = boundary[BdrElemNo]->GetVertices()[0];
932 if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
936 tr = GetFaceElementTransformations(fn);
941 void Mesh::GetFaceElements(
int Face,
int *Elem1,
int *Elem2)
const
943 *Elem1 = faces_info[Face].
Elem1No;
944 *Elem2 = faces_info[Face].Elem2No;
947 void Mesh::GetFaceInfos(
int Face,
int *Inf1,
int *Inf2)
const
949 *Inf1 = faces_info[Face].Elem1Inf;
950 *Inf2 = faces_info[Face].Elem2Inf;
955 return (Dim == 1) ? Geometry::POINT : faces[Face]->GetGeometryType();
960 return (Dim == 1) ? Element::POINT : faces[Face]->GetType();
968 NumOfElements = NumOfBdrElements = 0;
969 NumOfEdges = NumOfFaces = 0;
970 meshgen = mesh_geoms = 0;
976 last_operation = Mesh::NONE;
979 void Mesh::InitTables()
982 el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
985 void Mesh::SetEmpty()
991 void Mesh::DestroyTables()
996 DeleteGeometricFactors();
1007 void Mesh::DestroyPointers()
1009 if (own_nodes) {
delete Nodes; }
1015 for (
int i = 0; i < NumOfElements; i++)
1017 FreeElement(elements[i]);
1020 for (
int i = 0; i < NumOfBdrElements; i++)
1022 FreeElement(boundary[i]);
1025 for (
int i = 0; i < faces.Size(); i++)
1027 FreeElement(faces[i]);
1033 void Mesh::Destroy()
1037 elements.DeleteAll();
1038 vertices.DeleteAll();
1039 boundary.DeleteAll();
1041 faces_info.DeleteAll();
1042 nc_faces_info.DeleteAll();
1043 be_to_edge.DeleteAll();
1044 be_to_face.DeleteAll();
1051 CoarseFineTr.Clear();
1053 #ifdef MFEM_USE_MEMALLOC
1057 attributes.DeleteAll();
1058 bdr_attributes.DeleteAll();
1061 void Mesh::DeleteLazyTables()
1063 delete el_to_el; el_to_el = NULL;
1064 delete face_edge; face_edge = NULL;
1065 delete edge_vertex; edge_vertex = NULL;
1066 DeleteGeometricFactors();
1069 void Mesh::SetAttributes()
1074 for (
int i = 0; i < attribs.
Size(); i++)
1076 attribs[i] = GetBdrAttribute(i);
1080 attribs.
Copy(bdr_attributes);
1081 if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
1083 MFEM_WARNING(
"Non-positive attributes on the boundary!");
1087 for (
int i = 0; i < attribs.
Size(); i++)
1089 attribs[i] = GetAttribute(i);
1093 attribs.
Copy(attributes);
1094 if (attributes.Size() > 0 && attributes[0] <= 0)
1096 MFEM_WARNING(
"Non-positive attributes in the domain!");
1100 void Mesh::InitMesh(
int _Dim,
int _spaceDim,
int NVert,
int NElem,
int NBdrElem)
1105 spaceDim = _spaceDim;
1108 vertices.SetSize(NVert);
1111 elements.SetSize(NElem);
1113 NumOfBdrElements = 0;
1114 boundary.SetSize(NBdrElem);
1117 void Mesh::AddVertex(
const double *x)
1119 double *y = vertices[NumOfVertices]();
1121 for (
int i = 0; i < spaceDim; i++)
1128 void Mesh::AddTri(
const int *vi,
int attr)
1130 elements[NumOfElements++] =
new Triangle(vi, attr);
1133 void Mesh::AddTriangle(
const int *vi,
int attr)
1135 elements[NumOfElements++] =
new Triangle(vi, attr);
1138 void Mesh::AddQuad(
const int *vi,
int attr)
1143 void Mesh::AddTet(
const int *vi,
int attr)
1145 #ifdef MFEM_USE_MEMALLOC
1147 tet = TetMemory.Alloc();
1150 elements[NumOfElements++] = tet;
1152 elements[NumOfElements++] =
new Tetrahedron(vi, attr);
1156 void Mesh::AddWedge(
const int *vi,
int attr)
1158 elements[NumOfElements++] =
new Wedge(vi, attr);
1161 void Mesh::AddHex(
const int *vi,
int attr)
1163 elements[NumOfElements++] =
new Hexahedron(vi, attr);
1166 void Mesh::AddHexAsTets(
const int *vi,
int attr)
1168 static const int hex_to_tet[6][4] =
1170 { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
1171 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 }
1175 for (
int i = 0; i < 6; i++)
1177 for (
int j = 0; j < 4; j++)
1179 ti[j] = vi[hex_to_tet[i][j]];
1185 void Mesh::AddHexAsWedges(
const int *vi,
int attr)
1187 static const int hex_to_wdg[2][6] =
1189 { 0, 1, 2, 4, 5, 6 }, { 0, 2, 3, 4, 6, 7 }
1193 for (
int i = 0; i < 2; i++)
1195 for (
int j = 0; j < 6; j++)
1197 ti[j] = vi[hex_to_wdg[i][j]];
1203 void Mesh::AddBdrSegment(
const int *vi,
int attr)
1205 boundary[NumOfBdrElements++] =
new Segment(vi, attr);
1208 void Mesh::AddBdrTriangle(
const int *vi,
int attr)
1210 boundary[NumOfBdrElements++] =
new Triangle(vi, attr);
1213 void Mesh::AddBdrQuad(
const int *vi,
int attr)
1218 void Mesh::AddBdrQuadAsTriangles(
const int *vi,
int attr)
1220 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
1223 for (
int i = 0; i < 2; i++)
1225 for (
int j = 0; j < 3; j++)
1227 ti[j] = vi[quad_to_tri[i][j]];
1229 AddBdrTriangle(ti, attr);
1233 void Mesh::GenerateBoundaryElements()
1236 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
1240 for (i = 0; i < boundary.Size(); i++)
1242 FreeElement(boundary[i]);
1252 NumOfBdrElements = 0;
1253 for (i = 0; i < faces_info.Size(); i++)
1255 if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
1258 boundary.SetSize(NumOfBdrElements);
1259 be2face.
SetSize(NumOfBdrElements);
1260 for (j = i = 0; i < faces_info.Size(); i++)
1262 if (faces_info[i].Elem2No < 0)
1264 boundary[j] = faces[i]->Duplicate(
this);
1271 void Mesh::FinalizeCheck()
1273 MFEM_VERIFY(vertices.Size() == NumOfVertices ||
1274 vertices.Size() == 0,
1275 "incorrect number of vertices: preallocated: " << vertices.Size()
1276 <<
", actually added: " << NumOfVertices);
1277 MFEM_VERIFY(elements.Size() == NumOfElements,
1278 "incorrect number of elements: preallocated: " << elements.Size()
1279 <<
", actually added: " << NumOfElements);
1280 MFEM_VERIFY(boundary.Size() == NumOfBdrElements,
1281 "incorrect number of boundary elements: preallocated: "
1282 << boundary.Size() <<
", actually added: " << NumOfBdrElements);
1285 void Mesh::FinalizeTriMesh(
int generate_edges,
int refine,
bool fix_orientation)
1288 CheckElementOrientation(fix_orientation);
1292 MarkTriMeshForRefinement();
1297 el_to_edge =
new Table;
1298 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1300 CheckBdrElementOrientation();
1314 void Mesh::FinalizeQuadMesh(
int generate_edges,
int refine,
1315 bool fix_orientation)
1318 if (fix_orientation)
1320 CheckElementOrientation(fix_orientation);
1325 el_to_edge =
new Table;
1326 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1328 CheckBdrElementOrientation();
1343 #ifdef MFEM_USE_GECKO
1345 int iterations,
int window,
1346 int period,
int seed)
1350 Gecko::Functional *functional =
1351 new Gecko::FunctionalGeometric();
1354 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1361 const Table &my_el_to_el = ElementToElementTable();
1362 for (
int elemid = 0; elemid < GetNE(); ++elemid)
1364 const int *neighid = my_el_to_el.
GetRow(elemid);
1365 for (
int i = 0; i < my_el_to_el.
RowSize(elemid); ++i)
1367 graph.insert(elemid + 1, neighid[i] + 1);
1372 graph.order(functional, iterations, window, period, seed);
1375 Gecko::Node::Index NE = GetNE();
1376 for (Gecko::Node::Index gnodeid = 1; gnodeid <= NE; ++gnodeid)
1378 ordering[gnodeid - 1] = graph.rank(gnodeid);
1386 void Mesh::ReorderElements(
const Array<int> &ordering,
bool reorder_vertices)
1390 MFEM_WARNING(
"element reordering of NURBS meshes is not supported.");
1395 MFEM_WARNING(
"element reordering of non-conforming meshes is not"
1399 MFEM_VERIFY(ordering.
Size() == GetNE(),
"invalid reordering array.")
1429 old_elem_node_vals.SetSize(GetNE());
1430 nodes_fes = Nodes->FESpace();
1433 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1435 nodes_fes->GetElementVDofs(old_elid, old_dofs);
1437 old_elem_node_vals[old_elid] =
new Vector(vals);
1443 for (
int old_elid = 0; old_elid < ordering.
Size(); ++old_elid)
1445 int new_elid = ordering[old_elid];
1446 new_elements[new_elid] = elements[old_elid];
1451 if (reorder_vertices)
1456 vertex_ordering = -1;
1458 int new_vertex_ind = 0;
1459 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1461 int *elem_vert = elements[new_elid]->GetVertices();
1462 int nv = elements[new_elid]->GetNVertices();
1463 for (
int vi = 0; vi < nv; ++vi)
1465 int old_vertex_ind = elem_vert[vi];
1466 if (vertex_ordering[old_vertex_ind] == -1)
1468 vertex_ordering[old_vertex_ind] = new_vertex_ind;
1469 new_vertices[new_vertex_ind] = vertices[old_vertex_ind];
1479 for (
int new_elid = 0; new_elid < GetNE(); ++new_elid)
1481 int *elem_vert = elements[new_elid]->GetVertices();
1482 int nv = elements[new_elid]->GetNVertices();
1483 for (
int vi = 0; vi < nv; ++vi)
1485 elem_vert[vi] = vertex_ordering[elem_vert[vi]];
1490 for (
int belid = 0; belid < GetNBE(); ++belid)
1492 int *be_vert = boundary[belid]->GetVertices();
1493 int nv = boundary[belid]->GetNVertices();
1494 for (
int vi = 0; vi < nv; ++vi)
1496 be_vert[vi] = vertex_ordering[be_vert[vi]];
1507 el_to_edge =
new Table;
1508 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1513 GetElementToFaceTable();
1523 last_operation = Mesh::NONE;
1524 nodes_fes->Update(
false);
1527 for (
int old_elid = 0; old_elid < GetNE(); ++old_elid)
1529 int new_elid = ordering[old_elid];
1530 nodes_fes->GetElementVDofs(new_elid, new_dofs);
1531 Nodes->SetSubVector(new_dofs, *(old_elem_node_vals[old_elid]));
1532 delete old_elem_node_vals[old_elid];
1538 void Mesh::MarkForRefinement()
1544 MarkTriMeshForRefinement();
1548 DSTable v_to_v(NumOfVertices);
1549 GetVertexToVertexTable(v_to_v);
1550 MarkTetMeshForRefinement(v_to_v);
1555 void Mesh::MarkTriMeshForRefinement()
1560 for (
int i = 0; i < NumOfElements; i++)
1562 if (elements[i]->GetType() == Element::TRIANGLE)
1564 GetPointMatrix(i, pmat);
1565 static_cast<Triangle*
>(elements[i])->MarkEdge(pmat);
1576 for (
int i = 0; i < NumOfVertices; i++)
1581 length_idx[j].one = GetLength(i, it.Column());
1582 length_idx[j].two = j;
1589 for (
int i = 0; i < NumOfEdges; i++)
1591 order[length_idx[i].two] = i;
1595 void Mesh::MarkTetMeshForRefinement(
DSTable &v_to_v)
1600 GetEdgeOrdering(v_to_v, order);
1602 for (
int i = 0; i < NumOfElements; i++)
1604 if (elements[i]->GetType() == Element::TETRAHEDRON)
1606 elements[i]->MarkEdge(v_to_v, order);
1609 for (
int i = 0; i < NumOfBdrElements; i++)
1611 if (boundary[i]->GetType() == Element::TRIANGLE)
1613 boundary[i]->MarkEdge(v_to_v, order);
1620 if (*old_v_to_v && *old_elem_vert)
1627 if (*old_v_to_v == NULL)
1629 bool need_v_to_v =
false;
1631 for (
int i = 0; i < GetNEdges(); i++)
1637 if (dofs.
Size() > 0)
1645 *old_v_to_v =
new DSTable(NumOfVertices);
1646 GetVertexToVertexTable(*(*old_v_to_v));
1649 if (*old_elem_vert == NULL)
1651 bool need_elem_vert =
false;
1653 for (
int i = 0; i < GetNE(); i++)
1659 if (dofs.
Size() > 1)
1661 need_elem_vert =
true;
1667 *old_elem_vert =
new Table;
1668 (*old_elem_vert)->
MakeI(GetNE());
1669 for (
int i = 0; i < GetNE(); i++)
1671 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
1673 (*old_elem_vert)->MakeJ();
1674 for (
int i = 0; i < GetNE(); i++)
1676 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
1677 elements[i]->GetNVertices());
1679 (*old_elem_vert)->ShiftUpI();
1692 const int num_edge_dofs = old_dofs.
Size();
1695 const Vector onodes = *Nodes;
1699 int offset = NumOfVertices * old_dofs.
Size();
1703 if (num_edge_dofs > 0)
1705 DSTable new_v_to_v(NumOfVertices);
1706 GetVertexToVertexTable(new_v_to_v);
1708 for (
int i = 0; i < NumOfVertices; i++)
1712 const int old_i = (*old_v_to_v)(i, it.Column());
1713 const int new_i = it.Index();
1714 if (new_i == old_i) {
continue; }
1716 old_dofs.
SetSize(num_edge_dofs);
1717 new_dofs.
SetSize(num_edge_dofs);
1718 for (
int j = 0; j < num_edge_dofs; j++)
1720 old_dofs[j] = offset + old_i * num_edge_dofs + j;
1721 new_dofs[j] = offset + new_i * num_edge_dofs + j;
1725 for (
int j = 0; j < old_dofs.
Size(); j++)
1727 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1731 offset += NumOfEdges * num_edge_dofs;
1739 Table old_face_vertex;
1740 old_face_vertex.
MakeI(NumOfFaces);
1741 for (
int i = 0; i < NumOfFaces; i++)
1745 old_face_vertex.
MakeJ();
1746 for (
int i = 0; i < NumOfFaces; i++)
1748 faces[i]->GetNVertices());
1752 STable3D *faces_tbl = GetElementToFaceTable(1);
1758 for (
int i = 0; i < NumOfFaces; i++)
1760 const int *old_v = old_face_vertex.
GetRow(i);
1762 switch (old_face_vertex.
RowSize(i))
1765 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
1769 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
1773 new_fdofs[new_i+1] = old_dofs.
Size();
1778 for (
int i = 0; i < NumOfFaces; i++)
1780 const int *old_v = old_face_vertex.
GetRow(i), *new_v;
1783 switch (old_face_vertex.
RowSize(i))
1786 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
1787 new_v = faces[new_i]->GetVertices();
1788 new_or = GetTriOrientation(old_v, new_v);
1793 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
1794 new_v = faces[new_i]->GetVertices();
1795 new_or = GetQuadOrientation(old_v, new_v);
1802 for (
int j = 0; j < old_dofs.
Size(); j++)
1805 const int old_j = dof_ord[j];
1806 new_dofs[old_j] = offset + new_fdofs[new_i] + j;
1810 for (
int j = 0; j < old_dofs.
Size(); j++)
1812 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1834 for (
int i = 0; i < GetNE(); i++)
1836 const int *old_v = old_elem_vert->
GetRow(i);
1837 const int *new_v = elements[i]->GetVertices();
1843 case Geometry::SEGMENT:
1844 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
1846 case Geometry::TRIANGLE:
1847 new_or = GetTriOrientation(old_v, new_v);
1849 case Geometry::SQUARE:
1850 new_or = GetQuadOrientation(old_v, new_v);
1854 MFEM_ABORT(Geometry::Name[geom] <<
" elements (" << fec->
Name()
1855 <<
" FE collection) are not supported yet!");
1859 MFEM_VERIFY(dof_ord != NULL,
1860 "FE collection '" << fec->
Name()
1861 <<
"' does not define reordering for "
1862 << Geometry::Name[
geom] <<
" elements!");
1865 for (
int j = 0; j < new_dofs.
Size(); j++)
1868 const int old_j = dof_ord[j];
1869 new_dofs[old_j] = offset + j;
1871 offset += new_dofs.
Size();
1874 for (
int j = 0; j < old_dofs.
Size(); j++)
1876 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1888 GetElementToFaceTable();
1891 CheckBdrElementOrientation();
1896 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1901 CheckBdrElementOrientation();
1906 last_operation = Mesh::NONE;
1911 void Mesh::FinalizeTetMesh(
int generate_edges,
int refine,
bool fix_orientation)
1914 CheckElementOrientation(fix_orientation);
1916 if (NumOfBdrElements == 0)
1918 GetElementToFaceTable();
1920 GenerateBoundaryElements();
1925 DSTable v_to_v(NumOfVertices);
1926 GetVertexToVertexTable(v_to_v);
1927 MarkTetMeshForRefinement(v_to_v);
1930 GetElementToFaceTable();
1933 CheckBdrElementOrientation();
1935 if (generate_edges == 1)
1937 el_to_edge =
new Table;
1938 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1952 void Mesh::FinalizeWedgeMesh(
int generate_edges,
int refine,
1953 bool fix_orientation)
1956 CheckElementOrientation(fix_orientation);
1958 if (NumOfBdrElements == 0)
1960 GetElementToFaceTable();
1962 GenerateBoundaryElements();
1965 GetElementToFaceTable();
1968 CheckBdrElementOrientation();
1970 if (generate_edges == 1)
1972 el_to_edge =
new Table;
1973 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1987 void Mesh::FinalizeHexMesh(
int generate_edges,
int refine,
bool fix_orientation)
1990 CheckElementOrientation(fix_orientation);
1992 GetElementToFaceTable();
1995 if (NumOfBdrElements == 0)
1997 GenerateBoundaryElements();
2000 CheckBdrElementOrientation();
2004 el_to_edge =
new Table;
2005 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2017 void Mesh::FinalizeMesh(
int refine,
bool fix_orientation)
2021 Finalize(refine, fix_orientation);
2024 void Mesh::FinalizeTopology()
2036 bool generate_edges =
true;
2038 if (spaceDim == 0) { spaceDim = Dim; }
2039 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
2047 GetElementToFaceTable();
2049 if (NumOfBdrElements == 0)
2051 GenerateBoundaryElements();
2052 GetElementToFaceTable();
2061 if (Dim > 1 && generate_edges)
2064 if (!el_to_edge) { el_to_edge =
new Table; }
2065 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2069 if (NumOfBdrElements == 0)
2071 GenerateBoundaryElements();
2088 ncmesh->OnMeshUpdated(
this);
2091 GenerateNCFaceInfo();
2098 void Mesh::Finalize(
bool refine,
bool fix_orientation)
2100 if (NURBSext || ncmesh)
2102 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
2103 MFEM_ASSERT(CheckBdrElementOrientation() == 0,
"");
2112 const bool check_orientation =
true;
2113 const bool curved = (Nodes != NULL);
2114 const bool may_change_topology =
2115 ( refine && (Dim > 1 && (meshgen & 1)) ) ||
2116 ( check_orientation && fix_orientation &&
2117 (Dim == 2 || (Dim == 3 && (meshgen & 1))) );
2120 Table *old_elem_vert = NULL;
2122 if (curved && may_change_topology)
2124 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
2127 if (check_orientation)
2130 CheckElementOrientation(fix_orientation);
2134 MarkForRefinement();
2137 if (may_change_topology)
2141 DoNodeReorder(old_v_to_v, old_elem_vert);
2142 delete old_elem_vert;
2155 CheckBdrElementOrientation();
2160 if (Dim >= 2 && Dim == spaceDim)
2162 const int num_faces = GetNumFaces();
2163 for (
int i = 0; i < num_faces; i++)
2165 MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
2166 faces_info[i].Elem2Inf%2 != 0,
"invalid mesh topology");
2173 double sx,
double sy,
double sz,
bool sfc_ordering)
2177 int NVert, NElem, NBdrElem;
2179 NVert = (nx+1) * (ny+1) * (nz+1);
2180 NElem = nx * ny * nz;
2181 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
2182 if (type == Element::TETRAHEDRON)
2187 else if (type == Element::WEDGE)
2190 NBdrElem += 2*nx*ny;
2193 InitMesh(3, 3, NVert, NElem, NBdrElem);
2199 for (z = 0; z <= nz; z++)
2201 coord[2] = ((double) z / nz) * sz;
2202 for (y = 0; y <= ny; y++)
2204 coord[1] = ((double) y / ny) * sy;
2205 for (x = 0; x <= nx; x++)
2207 coord[0] = ((double) x / nx) * sx;
2213 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
2216 if (sfc_ordering && type == Element::HEXAHEDRON)
2219 NCMesh::GridSfcOrdering3D(nx, ny, nz, sfc);
2220 MFEM_VERIFY(sfc.
Size() == 3*nx*ny*nz,
"");
2222 for (
int k = 0; k < nx*ny*nz; k++)
2228 ind[0] = VTX(x , y , z );
2229 ind[1] = VTX(x+1, y , z );
2230 ind[2] = VTX(x+1, y+1, z );
2231 ind[3] = VTX(x , y+1, z );
2232 ind[4] = VTX(x , y , z+1);
2233 ind[5] = VTX(x+1, y , z+1);
2234 ind[6] = VTX(x+1, y+1, z+1);
2235 ind[7] = VTX(x , y+1, z+1);
2242 for (z = 0; z < nz; z++)
2244 for (y = 0; y < ny; y++)
2246 for (x = 0; x < nx; x++)
2248 ind[0] = VTX(x , y , z );
2249 ind[1] = VTX(x+1, y , z );
2250 ind[2] = VTX(x+1, y+1, z );
2251 ind[3] = VTX(x , y+1, z );
2252 ind[4] = VTX(x , y , z+1);
2253 ind[5] = VTX(x+1, y , z+1);
2254 ind[6] = VTX(x+1, y+1, z+1);
2255 ind[7] = VTX(x , y+1, z+1);
2256 if (type == Element::TETRAHEDRON)
2258 AddHexAsTets(ind, 1);
2260 else if (type == Element::WEDGE)
2262 AddHexAsWedges(ind, 1);
2275 for (y = 0; y < ny; y++)
2277 for (x = 0; x < nx; x++)
2279 ind[0] = VTX(x , y , 0);
2280 ind[1] = VTX(x , y+1, 0);
2281 ind[2] = VTX(x+1, y+1, 0);
2282 ind[3] = VTX(x+1, y , 0);
2283 if (type == Element::TETRAHEDRON)
2285 AddBdrQuadAsTriangles(ind, 1);
2287 else if (type == Element::WEDGE)
2289 AddBdrQuadAsTriangles(ind, 1);
2298 for (y = 0; y < ny; y++)
2300 for (x = 0; x < nx; x++)
2302 ind[0] = VTX(x , y , nz);
2303 ind[1] = VTX(x+1, y , nz);
2304 ind[2] = VTX(x+1, y+1, nz);
2305 ind[3] = VTX(x , y+1, nz);
2306 if (type == Element::TETRAHEDRON)
2308 AddBdrQuadAsTriangles(ind, 6);
2310 else if (type == Element::WEDGE)
2312 AddBdrQuadAsTriangles(ind, 1);
2321 for (z = 0; z < nz; z++)
2323 for (y = 0; y < ny; y++)
2325 ind[0] = VTX(0 , y , z );
2326 ind[1] = VTX(0 , y , z+1);
2327 ind[2] = VTX(0 , y+1, z+1);
2328 ind[3] = VTX(0 , y+1, z );
2329 if (type == Element::TETRAHEDRON)
2331 AddBdrQuadAsTriangles(ind, 5);
2340 for (z = 0; z < nz; z++)
2342 for (y = 0; y < ny; y++)
2344 ind[0] = VTX(nx, y , z );
2345 ind[1] = VTX(nx, y+1, z );
2346 ind[2] = VTX(nx, y+1, z+1);
2347 ind[3] = VTX(nx, y , z+1);
2348 if (type == Element::TETRAHEDRON)
2350 AddBdrQuadAsTriangles(ind, 3);
2359 for (x = 0; x < nx; x++)
2361 for (z = 0; z < nz; z++)
2363 ind[0] = VTX(x , 0, z );
2364 ind[1] = VTX(x+1, 0, z );
2365 ind[2] = VTX(x+1, 0, z+1);
2366 ind[3] = VTX(x , 0, z+1);
2367 if (type == Element::TETRAHEDRON)
2369 AddBdrQuadAsTriangles(ind, 2);
2378 for (x = 0; x < nx; x++)
2380 for (z = 0; z < nz; z++)
2382 ind[0] = VTX(x , ny, z );
2383 ind[1] = VTX(x , ny, z+1);
2384 ind[2] = VTX(x+1, ny, z+1);
2385 ind[3] = VTX(x+1, ny, z );
2386 if (type == Element::TETRAHEDRON)
2388 AddBdrQuadAsTriangles(ind, 4);
2400 ofstream test_stream(
"debug.mesh");
2402 test_stream.close();
2411 double sx,
double sy,
2412 bool generate_edges,
bool sfc_ordering)
2421 if (type == Element::QUADRILATERAL)
2423 NumOfVertices = (nx+1) * (ny+1);
2424 NumOfElements = nx * ny;
2425 NumOfBdrElements = 2 * nx + 2 * ny;
2427 vertices.SetSize(NumOfVertices);
2428 elements.SetSize(NumOfElements);
2429 boundary.SetSize(NumOfBdrElements);
2436 for (j = 0; j < ny+1; j++)
2438 cy = ((double) j / ny) * sy;
2439 for (i = 0; i < nx+1; i++)
2441 cx = ((double) i / nx) * sx;
2442 vertices[k](0) = cx;
2443 vertices[k](1) = cy;
2452 NCMesh::GridSfcOrdering2D(nx, ny, sfc);
2453 MFEM_VERIFY(sfc.
Size() == 2*nx*ny,
"");
2455 for (k = 0; k < nx*ny; k++)
2459 ind[0] = i + j*(nx+1);
2460 ind[1] = i + 1 +j*(nx+1);
2461 ind[2] = i + 1 + (j+1)*(nx+1);
2462 ind[3] = i + (j+1)*(nx+1);
2469 for (j = 0; j < ny; j++)
2471 for (i = 0; i < nx; i++)
2473 ind[0] = i + j*(nx+1);
2474 ind[1] = i + 1 +j*(nx+1);
2475 ind[2] = i + 1 + (j+1)*(nx+1);
2476 ind[3] = i + (j+1)*(nx+1);
2485 for (i = 0; i < nx; i++)
2487 boundary[i] =
new Segment(i, i+1, 1);
2488 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
2491 for (j = 0; j < ny; j++)
2493 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
2494 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
2498 else if (type == Element::TRIANGLE)
2500 NumOfVertices = (nx+1) * (ny+1);
2501 NumOfElements = 2 * nx * ny;
2502 NumOfBdrElements = 2 * nx + 2 * ny;
2504 vertices.SetSize(NumOfVertices);
2505 elements.SetSize(NumOfElements);
2506 boundary.SetSize(NumOfBdrElements);
2513 for (j = 0; j < ny+1; j++)
2515 cy = ((double) j / ny) * sy;
2516 for (i = 0; i < nx+1; i++)
2518 cx = ((double) i / nx) * sx;
2519 vertices[k](0) = cx;
2520 vertices[k](1) = cy;
2527 for (j = 0; j < ny; j++)
2529 for (i = 0; i < nx; i++)
2531 ind[0] = i + j*(nx+1);
2532 ind[1] = i + 1 + (j+1)*(nx+1);
2533 ind[2] = i + (j+1)*(nx+1);
2536 ind[1] = i + 1 + j*(nx+1);
2537 ind[2] = i + 1 + (j+1)*(nx+1);
2545 for (i = 0; i < nx; i++)
2547 boundary[i] =
new Segment(i, i+1, 1);
2548 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
2551 for (j = 0; j < ny; j++)
2553 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
2554 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
2561 MFEM_ABORT(
"Unsupported element type.");
2565 CheckElementOrientation();
2567 if (generate_edges == 1)
2569 el_to_edge =
new Table;
2570 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2572 CheckBdrElementOrientation();
2581 attributes.Append(1);
2582 bdr_attributes.Append(1); bdr_attributes.Append(2);
2583 bdr_attributes.Append(3); bdr_attributes.Append(4);
2588 void Mesh::Make1D(
int n,
double sx)
2597 NumOfVertices = n + 1;
2599 NumOfBdrElements = 2;
2600 vertices.SetSize(NumOfVertices);
2601 elements.SetSize(NumOfElements);
2602 boundary.SetSize(NumOfBdrElements);
2605 for (j = 0; j < n+1; j++)
2607 vertices[j](0) = ((
double) j / n) * sx;
2611 for (j = 0; j < n; j++)
2613 elements[j] =
new Segment(j, j+1, 1);
2618 boundary[0] =
new Point(ind, 1);
2620 boundary[1] =
new Point(ind, 2);
2628 attributes.Append(1);
2629 bdr_attributes.Append(1); bdr_attributes.Append(2);
2632 Mesh::Mesh(
const Mesh &mesh,
bool copy_nodes)
2648 last_operation = Mesh::NONE;
2651 elements.SetSize(NumOfElements);
2652 for (
int i = 0; i < NumOfElements; i++)
2654 elements[i] = mesh.
elements[i]->Duplicate(
this);
2661 boundary.SetSize(NumOfBdrElements);
2662 for (
int i = 0; i < NumOfBdrElements; i++)
2664 boundary[i] = mesh.
boundary[i]->Duplicate(
this);
2683 faces.SetSize(mesh.
faces.Size());
2684 for (
int i = 0; i < faces.Size(); i++)
2687 faces[i] = (face) ? face->
Duplicate(
this) : NULL;
2721 if (dynamic_cast<const ParMesh*>(&mesh))
2733 if (mesh.
Nodes && copy_nodes)
2738 FiniteElementCollection::New(fec->
Name());
2742 Nodes->MakeOwner(fec_copy);
2743 *Nodes = *mesh.
Nodes;
2753 Mesh::Mesh(
const char *filename,
int generate_edges,
int refine,
2754 bool fix_orientation)
2763 MFEM_ABORT(
"Mesh file not found: " << filename <<
'\n');
2767 Load(imesh, generate_edges, refine, fix_orientation);
2771 Mesh::Mesh(std::istream &input,
int generate_edges,
int refine,
2772 bool fix_orientation)
2775 Load(input, generate_edges, refine, fix_orientation);
2778 void Mesh::ChangeVertexDataOwnership(
double *vertex_data,
int len_vertex_data,
2783 MFEM_VERIFY(len_vertex_data >= NumOfVertices * 3,
2784 "Not enough vertices in external array : "
2785 "len_vertex_data = "<< len_vertex_data <<
", "
2786 "NumOfVertices * 3 = " << NumOfVertices * 3);
2788 if (vertex_data == (
double *)(vertices.GetData()))
2790 MFEM_ASSERT(!vertices.OwnsData(),
"invalid ownership");
2795 memcpy(vertex_data, vertices.GetData(),
2796 NumOfVertices * 3 *
sizeof(double));
2799 vertices.MakeRef(reinterpret_cast<Vertex*>(vertex_data), NumOfVertices);
2802 Mesh::Mesh(
double *_vertices,
int num_vertices,
2804 int *element_attributes,
int num_elements,
2806 int *boundary_attributes,
int num_boundary_elements,
2807 int dimension,
int space_dimension)
2809 if (space_dimension == -1)
2811 space_dimension = dimension;
2814 InitMesh(dimension, space_dimension, 0, num_elements,
2815 num_boundary_elements);
2817 int element_index_stride = Geometry::NumVerts[element_type];
2818 int boundary_index_stride = num_boundary_elements > 0 ?
2819 Geometry::NumVerts[boundary_type] : 0;
2822 vertices.MakeRef(reinterpret_cast<Vertex*>(_vertices), num_vertices);
2823 NumOfVertices = num_vertices;
2825 for (
int i = 0; i < num_elements; i++)
2827 elements[i] = NewElement(element_type);
2828 elements[i]->SetVertices(element_indices + i * element_index_stride);
2829 elements[i]->SetAttribute(element_attributes[i]);
2831 NumOfElements = num_elements;
2833 for (
int i = 0; i < num_boundary_elements; i++)
2835 boundary[i] = NewElement(boundary_type);
2836 boundary[i]->SetVertices(boundary_indices + i * boundary_index_stride);
2837 boundary[i]->SetAttribute(boundary_attributes[i]);
2839 NumOfBdrElements = num_boundary_elements;
2848 case Geometry::POINT:
return (
new Point);
2849 case Geometry::SEGMENT:
return (
new Segment);
2850 case Geometry::TRIANGLE:
return (
new Triangle);
2852 case Geometry::TETRAHEDRON:
2853 #ifdef MFEM_USE_MEMALLOC
2854 return TetMemory.Alloc();
2858 case Geometry::CUBE:
return (
new Hexahedron);
2859 case Geometry::PRISM:
return (
new Wedge);
2861 MFEM_ABORT(
"invalid Geometry::Type, geom = " << geom);
2867 Element *Mesh::ReadElementWithoutAttr(std::istream &input)
2873 el = NewElement(geom);
2874 MFEM_VERIFY(el,
"Unsupported element type: " << geom);
2877 for (
int i = 0; i < nv; i++)
2885 void Mesh::PrintElementWithoutAttr(
const Element *el, std::ostream &
out)
2890 for (
int j = 0; j < nv; j++)
2903 el = ReadElementWithoutAttr(input);
2912 PrintElementWithoutAttr(el, out);
2915 void Mesh::SetMeshGen()
2917 meshgen = mesh_geoms = 0;
2918 for (
int i = 0; i < NumOfElements; i++)
2923 case Element::TETRAHEDRON:
2924 mesh_geoms |= (1 << Geometry::TETRAHEDRON);
2925 case Element::TRIANGLE:
2926 mesh_geoms |= (1 << Geometry::TRIANGLE);
2927 case Element::SEGMENT:
2928 mesh_geoms |= (1 << Geometry::SEGMENT);
2929 case Element::POINT:
2930 mesh_geoms |= (1 << Geometry::POINT);
2934 case Element::HEXAHEDRON:
2935 mesh_geoms |= (1 << Geometry::CUBE);
2936 case Element::QUADRILATERAL:
2937 mesh_geoms |= (1 << Geometry::SQUARE);
2938 mesh_geoms |= (1 << Geometry::SEGMENT);
2939 mesh_geoms |= (1 << Geometry::POINT);
2943 case Element::WEDGE:
2944 mesh_geoms |= (1 << Geometry::PRISM);
2945 mesh_geoms |= (1 << Geometry::SQUARE);
2946 mesh_geoms |= (1 << Geometry::TRIANGLE);
2947 mesh_geoms |= (1 << Geometry::SEGMENT);
2948 mesh_geoms |= (1 << Geometry::POINT);
2953 MFEM_ABORT(
"invalid element type: " << type);
2959 void Mesh::Loader(std::istream &input,
int generate_edges,
2960 std::string parse_tag)
2962 int curved = 0, read_gf = 1;
2963 bool finalize_topo =
true;
2967 MFEM_ABORT(
"Input stream is not open");
2974 getline(input, mesh_type);
2978 bool mfem_v10 = (mesh_type ==
"MFEM mesh v1.0");
2979 bool mfem_v11 = (mesh_type ==
"MFEM mesh v1.1");
2980 bool mfem_v12 = (mesh_type ==
"MFEM mesh v1.2");
2981 if (mfem_v10 || mfem_v11 || mfem_v12)
2987 if ( mfem_v12 && parse_tag.empty() )
2989 parse_tag =
"mfem_mesh_end";
2991 ReadMFEMMesh(input, mfem_v11, curved);
2993 else if (mesh_type ==
"linemesh")
2995 ReadLineMesh(input);
2997 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
2999 if (mesh_type ==
"curved_areamesh2")
3003 ReadNetgen2DMesh(input, curved);
3005 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
3007 ReadNetgen3DMesh(input);
3009 else if (mesh_type ==
"TrueGrid")
3011 ReadTrueGridMesh(input);
3013 else if (mesh_type ==
"# vtk DataFile Version 3.0" ||
3014 mesh_type ==
"# vtk DataFile Version 2.0")
3016 ReadVTKMesh(input, curved, read_gf, finalize_topo);
3018 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
3020 ReadNURBSMesh(input, curved, read_gf);
3022 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
3024 ReadInlineMesh(input, generate_edges);
3027 else if (mesh_type ==
"$MeshFormat")
3029 ReadGmshMesh(input);
3032 ((mesh_type.size() > 2 &&
3033 mesh_type[0] ==
'C' && mesh_type[1] ==
'D' && mesh_type[2] ==
'F') ||
3034 (mesh_type.size() > 3 &&
3035 mesh_type[1] ==
'H' && mesh_type[2] ==
'D' && mesh_type[3] ==
'F'))
3040 #ifdef MFEM_USE_NETCDF
3041 ReadCubit(mesh_input->
filename, curved, read_gf);
3043 MFEM_ABORT(
"NetCDF support requires configuration with"
3044 " MFEM_USE_NETCDF=YES");
3050 MFEM_ABORT(
"Can not determine Cubit mesh filename!"
3051 " Use mfem::named_ifgzstream for input.");
3057 MFEM_ABORT(
"Unknown input mesh format: " << mesh_type);
3085 if (curved && read_gf)
3089 spaceDim = Nodes->VectorDim();
3090 if (ncmesh) { ncmesh->spaceDim = spaceDim; }
3092 for (
int i = 0; i < spaceDim; i++)
3095 Nodes->GetNodalValues(vert_val, i+1);
3096 for (
int j = 0; j < NumOfVertices; j++)
3098 vertices[j](i) = vert_val(j);
3111 MFEM_VERIFY(input.good(),
"Required mesh-end tag not found");
3112 getline(input, line);
3118 if (line ==
"mfem_mesh_end") {
break; }
3120 while (line != parse_tag);
3126 Mesh::Mesh(
Mesh *mesh_array[],
int num_pieces)
3128 int i, j, ie, ib, iv, *v, nv;
3137 if (mesh_array[0]->NURBSext)
3142 NumOfVertices = NURBSext->GetNV();
3143 NumOfElements = NURBSext->GetNE();
3145 NURBSext->GetElementTopo(elements);
3158 NumOfBdrElements = 0;
3159 for (i = 0; i < num_pieces; i++)
3161 NumOfBdrElements += mesh_array[i]->
GetNBE();
3163 boundary.SetSize(NumOfBdrElements);
3164 vertices.SetSize(NumOfVertices);
3166 for (i = 0; i < num_pieces; i++)
3172 for (j = 0; j < m->
GetNE(); j++)
3174 elements[lelem_elem[j]]->SetAttribute(m->
GetAttribute(j));
3177 for (j = 0; j < m->
GetNBE(); j++)
3182 for (
int k = 0; k < nv; k++)
3184 v[k] = lvert_vert[v[k]];
3186 boundary[ib++] = el;
3189 for (j = 0; j < m->
GetNV(); j++)
3199 NumOfBdrElements = 0;
3201 for (i = 0; i < num_pieces; i++)
3204 NumOfElements += m->
GetNE();
3205 NumOfBdrElements += m->
GetNBE();
3206 NumOfVertices += m->
GetNV();
3208 elements.SetSize(NumOfElements);
3209 boundary.SetSize(NumOfBdrElements);
3210 vertices.SetSize(NumOfVertices);
3212 for (i = 0; i < num_pieces; i++)
3216 for (j = 0; j < m->
GetNE(); j++)
3221 for (
int k = 0; k < nv; k++)
3225 elements[ie++] = el;
3228 for (j = 0; j < m->
GetNBE(); j++)
3233 for (
int k = 0; k < nv; k++)
3237 boundary[ib++] = el;
3240 for (j = 0; j < m->
GetNV(); j++)
3254 for (i = 0; i < num_pieces; i++)
3256 gf_array[i] = mesh_array[i]->
GetNodes();
3263 CheckElementOrientation(
false);
3264 CheckBdrElementOrientation(
false);
3268 Mesh::Mesh(
Mesh *orig_mesh,
int ref_factor,
int ref_type)
3271 MFEM_VERIFY(ref_factor >= 1,
"the refinement factor must be >= 1");
3272 MFEM_VERIFY(ref_type == BasisType::ClosedUniform ||
3273 ref_type == BasisType::GaussLobatto,
"invalid refinement type");
3274 MFEM_VERIFY(Dim == 1 || Dim == 2 || Dim == 3,
3275 "only implemented for Segment, Quadrilateral and Hexahedron "
3276 "elements in 1D/2D/3D");
3278 "meshes with mixed elements are not supported");
3285 int r_bndr_factor = pow(ref_factor, Dim - 1);
3286 int r_elem_factor = ref_factor * r_bndr_factor;
3289 int r_num_elem = orig_mesh->
GetNE() * r_elem_factor;
3290 int r_num_bndr = orig_mesh->
GetNBE() * r_bndr_factor;
3292 InitMesh(Dim, orig_mesh->
SpaceDimension(), r_num_vert, r_num_elem,
3296 NumOfVertices = r_num_vert;
3301 for (
int el = 0; el < orig_mesh->
GetNE(); el++)
3305 int nvert = Geometry::NumVerts[
geom];
3308 max_nv = std::max(max_nv, nvert);
3314 const int *c2h_map = rfec.GetDofMap(geom);
3315 for (
int i = 0; i < phys_pts.
Width(); i++)
3317 vertices[rdofs[i]].SetCoords(spaceDim, phys_pts.
GetColumn(i));
3321 Element *elem = NewElement(geom);
3324 for (
int k = 0; k < nvert; k++)
3327 v[k] = rdofs[c2h_map[cid]];
3333 for (
int el = 0; el < orig_mesh->
GetNBE(); el++)
3337 int nvert = Geometry::NumVerts[
geom];
3348 Element *elem = NewElement(geom);
3351 v[0] = rdofs[RG.
RefGeoms[nvert*j]];
3352 AddBdrElement(elem);
3357 const int *c2h_map = rfec.GetDofMap(geom);
3360 Element *elem = NewElement(geom);
3363 for (
int k = 0; k < nvert; k++)
3366 v[k] = rdofs[c2h_map[cid]];
3368 AddBdrElement(elem);
3375 last_operation = Mesh::REFINE;
3378 CoarseFineTr.embeddings.SetSize(GetNE());
3379 if (orig_mesh->
GetNE() > 0)
3383 CoarseFineTr.point_matrices[
geom].SetSize(Dim, max_nv, r_elem_factor);
3384 int nvert = Geometry::NumVerts[
geom];
3386 const int *c2h_map = rfec.GetDofMap(geom);
3391 for (
int k = 0; k < nvert; k++)
3399 for (
int el = 0; el < GetNE(); el++)
3401 Embedding &emb = CoarseFineTr.embeddings[el];
3402 emb.
parent = el / r_elem_factor;
3403 emb.
matrix = el % r_elem_factor;
3406 MFEM_ASSERT(CheckElementOrientation(
false) == 0,
"");
3407 MFEM_ASSERT(CheckBdrElementOrientation(
false) == 0,
"");
3412 if (NURBSext == NULL)
3414 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
3417 if (kv.
Size() != NURBSext->GetNKV())
3419 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
3422 NURBSext->ConvertToPatches(*Nodes);
3424 NURBSext->KnotInsert(kv);
3426 last_operation = Mesh::NONE;
3432 void Mesh::NURBSUniformRefinement()
3435 NURBSext->ConvertToPatches(*Nodes);
3437 NURBSext->UniformRefinement();
3439 last_operation = Mesh::NONE;
3445 void Mesh::DegreeElevate(
int rel_degree,
int degree)
3447 if (NURBSext == NULL)
3449 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
3452 NURBSext->ConvertToPatches(*Nodes);
3454 NURBSext->DegreeElevate(rel_degree, degree);
3456 last_operation = Mesh::NONE;
3462 void Mesh::UpdateNURBS()
3464 NURBSext->SetKnotsFromPatches();
3466 Dim = NURBSext->Dimension();
3469 if (NumOfElements != NURBSext->GetNE())
3471 for (
int i = 0; i < elements.Size(); i++)
3473 FreeElement(elements[i]);
3475 NumOfElements = NURBSext->GetNE();
3476 NURBSext->GetElementTopo(elements);
3479 if (NumOfBdrElements != NURBSext->GetNBE())
3481 for (
int i = 0; i < boundary.Size(); i++)
3483 FreeElement(boundary[i]);
3485 NumOfBdrElements = NURBSext->GetNBE();
3486 NURBSext->GetBdrElementTopo(boundary);
3489 Nodes->FESpace()->Update();
3491 NURBSext->SetCoordsFromPatches(*Nodes);
3493 if (NumOfVertices != NURBSext->GetNV())
3495 NumOfVertices = NURBSext->GetNV();
3496 vertices.SetSize(NumOfVertices);
3497 int vd = Nodes->VectorDim();
3498 for (
int i = 0; i < vd; i++)
3501 Nodes->GetNodalValues(vert_val, i+1);
3502 for (
int j = 0; j < NumOfVertices; j++)
3504 vertices[j](i) = vert_val(j);
3511 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3520 GetElementToFaceTable();
3525 void Mesh::LoadPatchTopo(std::istream &input,
Array<int> &edge_to_knot)
3541 input >> NumOfElements;
3542 elements.SetSize(NumOfElements);
3543 for (
int j = 0; j < NumOfElements; j++)
3545 elements[j] = ReadElement(input);
3551 input >> NumOfBdrElements;
3552 boundary.SetSize(NumOfBdrElements);
3553 for (
int j = 0; j < NumOfBdrElements; j++)
3555 boundary[j] = ReadElement(input);
3561 input >> NumOfEdges;
3562 edge_vertex =
new Table(NumOfEdges, 2);
3563 edge_to_knot.
SetSize(NumOfEdges);
3564 for (
int j = 0; j < NumOfEdges; j++)
3566 int *v = edge_vertex->GetRow(j);
3567 input >> edge_to_knot[j] >> v[0] >> v[1];
3570 edge_to_knot[j] = -1 - edge_to_knot[j];
3577 input >> NumOfVertices;
3578 vertices.SetSize(0);
3581 CheckBdrElementOrientation();
3588 for (
int d = 0; d < v.
Size(); d++)
3596 for (d = 0; d < p.
Size(); d++)
3600 for ( ; d < v.
Size(); d++)
3609 if (Nodes == NULL || Nodes->FESpace() != nodes.
FESpace())
3624 SetNodalGridFunction(nodes,
true);
3627 void Mesh::EnsureNodes()
3629 if (Nodes) {
return; }
3630 SetCurvature(1,
false, -1, Ordering::byVDIM);
3636 NewNodes(*nodes, make_owner);
3641 return ((Nodes) ? Nodes->FESpace() : NULL);
3644 void Mesh::SetCurvature(
int order,
bool discont,
int space_dim,
int ordering)
3646 space_dim = (space_dim == -1) ? spaceDim : space_dim;
3659 SetNodalFESpace(nfes);
3660 Nodes->MakeOwner(nfec);
3663 int Mesh::GetNumFaces()
const
3667 case 1:
return GetNV();
3668 case 2:
return GetNEdges();
3669 case 3:
return GetNFaces();
3674 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3675 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
3678 int Mesh::CheckElementOrientation(
bool fix_it)
3680 int i, j, k, wo = 0, fo = 0, *vi = 0;
3683 if (Dim == 2 && spaceDim == 2)
3687 for (i = 0; i < NumOfElements; i++)
3691 vi = elements[i]->GetVertices();
3692 for (j = 0; j < 3; j++)
3694 v[j] = vertices[vi[j]]();
3696 for (j = 0; j < 2; j++)
3697 for (k = 0; k < 2; k++)
3699 J(j, k) = v[j+1][k] - v[0][k];
3705 GetElementJacobian(i, J);
3711 switch (GetElementType(i))
3713 case Element::TRIANGLE:
3716 case Element::QUADRILATERAL:
3720 MFEM_ABORT(
"Invalid 2D element type \""
3721 << GetElementType(i) <<
"\"");
3735 for (i = 0; i < NumOfElements; i++)
3737 vi = elements[i]->GetVertices();
3738 switch (GetElementType(i))
3740 case Element::TETRAHEDRON:
3743 for (j = 0; j < 4; j++)
3745 v[j] = vertices[vi[j]]();
3747 for (j = 0; j < 3; j++)
3748 for (k = 0; k < 3; k++)
3750 J(j, k) = v[j+1][k] - v[0][k];
3756 GetElementJacobian(i, J);
3769 case Element::WEDGE:
3771 GetElementJacobian(i, J);
3782 case Element::HEXAHEDRON:
3784 GetElementJacobian(i, J);
3796 MFEM_ABORT(
"Invalid 3D element type \""
3797 << GetElementType(i) <<
"\"");
3802 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3804 mfem::out <<
"Elements with wrong orientation: " << wo <<
" / "
3805 << NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
3811 int Mesh::GetTriOrientation(
const int *base,
const int *test)
3819 if (test[0] == base[0])
3820 if (test[1] == base[1])
3828 else if (test[0] == base[1])
3829 if (test[1] == base[0])
3838 if (test[1] == base[0])
3848 const int *aor = tri_t::Orient[orient];
3849 for (
int j = 0; j < 3; j++)
3850 if (test[aor[j]] != base[j])
3859 int Mesh::GetQuadOrientation(
const int *base,
const int *test)
3863 for (i = 0; i < 4; i++)
3864 if (test[i] == base[0])
3871 if (test[(i+1)%4] == base[1])
3879 const int *aor = quad_t::Orient[orient];
3880 for (
int j = 0; j < 4; j++)
3881 if (test[aor[j]] != base[j])
3883 mfem::err <<
"Mesh::GetQuadOrientation(...)" << endl;
3885 for (
int k = 0; k < 4; k++)
3890 for (
int k = 0; k < 4; k++)
3899 if (test[(i+1)%4] == base[1])
3907 int Mesh::CheckBdrElementOrientation(
bool fix_it)
3913 if (el_to_edge == NULL)
3915 el_to_edge =
new Table;
3916 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3919 for (
int i = 0; i < NumOfBdrElements; i++)
3921 if (faces_info[be_to_edge[i]].Elem2No < 0)
3923 int *bv = boundary[i]->GetVertices();
3924 int *fv = faces[be_to_edge[i]]->GetVertices();
3929 mfem::Swap<int>(bv[0], bv[1]);
3939 for (
int i = 0; i < NumOfBdrElements; i++)
3941 const int fi = be_to_face[i];
3943 if (faces_info[fi].Elem2No >= 0) {
continue; }
3946 int *bv = boundary[i]->GetVertices();
3948 MFEM_ASSERT(fi < faces.Size(),
"internal error");
3949 const int *fv = faces[fi]->GetVertices();
3955 case Element::TRIANGLE:
3957 orientation = GetTriOrientation(fv, bv);
3960 case Element::QUADRILATERAL:
3962 orientation = GetQuadOrientation(fv, bv);
3966 MFEM_ABORT(
"Invalid 2D boundary element type \""
3967 << bdr_type <<
"\"");
3972 if (orientation % 2 == 0) {
continue; }
3974 if (!fix_it) {
continue; }
3978 case Element::TRIANGLE:
3982 mfem::Swap<int>(bv[0], bv[1]);
3985 int *be = bel_to_edge->GetRow(i);
3986 mfem::Swap<int>(be[1], be[2]);
3990 case Element::QUADRILATERAL:
3992 mfem::Swap<int>(bv[0], bv[2]);
3995 int *be = bel_to_edge->GetRow(i);
3996 mfem::Swap<int>(be[0], be[1]);
3997 mfem::Swap<int>(be[2], be[3]);
4010 mfem::out <<
"Boundary elements with wrong orientation: " << wo <<
" / "
4011 << NumOfBdrElements <<
" (" << fixed_or_not[fix_it ? 0 : 1]
4018 int Mesh::GetNumGeometries(
int dim)
const
4020 MFEM_ASSERT(0 <= dim && dim <= Dim,
"invalid dim: " << dim);
4022 for (
int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
4031 MFEM_ASSERT(0 <= dim && dim <= Dim,
"invalid dim: " << dim);
4033 for (
int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
4046 el_to_edge->GetRow(i, edges);
4050 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
4051 "is not generated.");
4054 const int *v = elements[i]->GetVertices();
4055 const int ne = elements[i]->GetNEdges();
4057 for (
int j = 0; j < ne; j++)
4059 const int *e = elements[i]->GetEdgeVertices(j);
4060 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4070 edges[0] = be_to_edge[i];
4071 const int *v = boundary[i]->GetVertices();
4072 cor[0] = (v[0] < v[1]) ? (1) : (-1);
4078 bel_to_edge->GetRow(i, edges);
4085 const int *v = boundary[i]->GetVertices();
4086 const int ne = boundary[i]->GetNEdges();
4088 for (
int j = 0; j < ne; j++)
4090 const int *e = boundary[i]->GetEdgeVertices(j);
4091 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4103 const int *v = faces[i]->GetVertices();
4104 o[0] = (v[0] < v[1]) ? (1) : (-1);
4114 face_edge->GetRow(i, edges);
4116 const int *v = faces[i]->GetVertices();
4117 const int ne = faces[i]->GetNEdges();
4119 for (
int j = 0; j < ne; j++)
4121 const int *e = faces[i]->GetEdgeVertices(j);
4122 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
4131 if (!edge_vertex) { GetEdgeVertexTable(); }
4132 edge_vertex->GetRow(i, vert);
4148 if (faces.Size() != NumOfFaces)
4150 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
4154 DSTable v_to_v(NumOfVertices);
4155 GetVertexToVertexTable(v_to_v);
4157 face_edge =
new Table;
4158 GetElementArrayEdgeTable(faces, v_to_v, *face_edge);
4170 DSTable v_to_v(NumOfVertices);
4171 GetVertexToVertexTable(v_to_v);
4174 edge_vertex =
new Table(nedges, 2);
4175 for (
int i = 0; i < NumOfVertices; i++)
4180 edge_vertex->Push(j, i);
4181 edge_vertex->Push(j, it.Column());
4184 edge_vertex->Finalize();
4195 vert_elem->
MakeI(NumOfVertices);
4197 for (i = 0; i < NumOfElements; i++)
4199 nv = elements[i]->GetNVertices();
4200 v = elements[i]->GetVertices();
4201 for (j = 0; j < nv; j++)
4209 for (i = 0; i < NumOfElements; i++)
4211 nv = elements[i]->GetNVertices();
4212 v = elements[i]->GetVertices();
4213 for (j = 0; j < nv; j++)
4224 Table *Mesh::GetFaceToElementTable()
const
4228 face_elem->
MakeI(faces_info.Size());
4230 for (
int i = 0; i < faces_info.Size(); i++)
4232 if (faces_info[i].Elem2No >= 0)
4244 for (
int i = 0; i < faces_info.Size(); i++)
4247 if (faces_info[i].Elem2No >= 0)
4265 el_to_face->
GetRow(i, fcs);
4269 mfem_error(
"Mesh::GetElementFaces(...) : el_to_face not generated.");
4274 for (j = 0; j < n; j++)
4275 if (faces_info[fcs[j]].Elem1No == i)
4277 cor[j] = faces_info[fcs[j]].Elem1Inf % 64;
4280 else if (faces_info[fcs[j]].Elem2No == i)
4282 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
4286 mfem_error(
"Mesh::GetElementFaces(...) : 2");
4291 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
4296 void Mesh::GetBdrElementFace(
int i,
int *f,
int *o)
const
4301 bv = boundary[i]->GetVertices();
4302 fv = faces[be_to_face[i]]->GetVertices();
4306 switch (GetBdrElementType(i))
4308 case Element::TRIANGLE:
4309 *o = GetTriOrientation(fv, bv);
4311 case Element::QUADRILATERAL:
4312 *o = GetQuadOrientation(fv, bv);
4315 mfem_error(
"Mesh::GetBdrElementFace(...) 2");
4319 int Mesh::GetBdrElementEdgeIndex(
int i)
const
4323 case 1:
return boundary[i]->GetVertices()[0];
4324 case 2:
return be_to_edge[i];
4325 case 3:
return be_to_face[i];
4326 default:
mfem_error(
"Mesh::GetBdrElementEdgeIndex: invalid dimension!");
4331 void Mesh::GetBdrElementAdjacentElement(
int bdr_el,
int &el,
int &info)
const
4333 int fid = GetBdrElementEdgeIndex(bdr_el);
4334 const FaceInfo &fi = faces_info[fid];
4335 MFEM_ASSERT(fi.
Elem1Inf%64 == 0,
"internal error");
4336 const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
4337 const int *bv = boundary[bdr_el]->GetVertices();
4339 switch (GetBdrElementBaseGeometry(bdr_el))
4341 case Geometry::POINT: ori = 0;
break;
4342 case Geometry::SEGMENT: ori = (fv[0] == bv[0]) ? 0 : 1;
break;
4343 case Geometry::TRIANGLE: ori = GetTriOrientation(fv, bv);
break;
4344 case Geometry::SQUARE: ori = GetQuadOrientation(fv, bv);
break;
4345 default: MFEM_ABORT(
"boundary element type not implemented"); ori = 0;
4353 return elements[i]->GetType();
4358 return boundary[i]->GetType();
4366 v = elements[i]->GetVertices();
4367 nv = elements[i]->GetNVertices();
4369 pointmat.
SetSize(spaceDim, nv);
4370 for (k = 0; k < spaceDim; k++)
4371 for (j = 0; j < nv; j++)
4373 pointmat(k, j) = vertices[v[j]](k);
4382 v = boundary[i]->GetVertices();
4383 nv = boundary[i]->GetNVertices();
4385 pointmat.
SetSize(spaceDim, nv);
4386 for (k = 0; k < spaceDim; k++)
4387 for (j = 0; j < nv; j++)
4389 pointmat(k, j) = vertices[v[j]](k);
4393 double Mesh::GetLength(
int i,
int j)
const
4395 const double *vi = vertices[i]();
4396 const double *vj = vertices[j]();
4399 for (
int k = 0; k < spaceDim; k++)
4401 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
4404 return sqrt(length);
4412 for (
int i = 0; i < elem_array.
Size(); i++)
4417 for (
int i = 0; i < elem_array.
Size(); i++)
4419 const int *v = elem_array[i]->GetVertices();
4420 const int ne = elem_array[i]->GetNEdges();
4421 for (
int j = 0; j < ne; j++)
4423 const int *e = elem_array[i]->GetEdgeVertices(j);
4430 void Mesh::GetVertexToVertexTable(
DSTable &v_to_v)
const
4434 for (
int i = 0; i < edge_vertex->Size(); i++)
4436 const int *v = edge_vertex->GetRow(i);
4437 v_to_v.
Push(v[0], v[1]);
4442 for (
int i = 0; i < NumOfElements; i++)
4444 const int *v = elements[i]->GetVertices();
4445 const int ne = elements[i]->GetNEdges();
4446 for (
int j = 0; j < ne; j++)
4448 const int *e = elements[i]->GetEdgeVertices(j);
4449 v_to_v.
Push(v[e[0]], v[e[1]]);
4457 int i, NumberOfEdges;
4459 DSTable v_to_v(NumOfVertices);
4460 GetVertexToVertexTable(v_to_v);
4465 GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
4470 be_to_f.
SetSize(NumOfBdrElements);
4471 for (i = 0; i < NumOfBdrElements; i++)
4473 const int *v = boundary[i]->GetVertices();
4474 be_to_f[i] = v_to_v(v[0], v[1]);
4479 if (bel_to_edge == NULL)
4481 bel_to_edge =
new Table;
4483 GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
4487 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
4491 return NumberOfEdges;
4494 const Table & Mesh::ElementToElementTable()
4502 MFEM_ASSERT(faces_info.Size() >= GetNumFaces(),
"faces were not generated!");
4505 conn.
Reserve(2*faces_info.Size());
4507 for (
int i = 0; i < faces_info.Size(); i++)
4509 const FaceInfo &fi = faces_info[i];
4517 int nbr_elem_idx = NumOfElements - 1 - fi.
Elem2No;
4525 el_to_el =
new Table(NumOfElements, conn);
4530 const Table & Mesh::ElementToFaceTable()
const
4532 if (el_to_face == NULL)
4539 const Table & Mesh::ElementToEdgeTable()
const
4541 if (el_to_edge == NULL)
4548 void Mesh::AddPointFaceElement(
int lf,
int gf,
int el)
4550 if (faces_info[gf].Elem1No == -1)
4553 faces_info[gf].Elem1No = el;
4554 faces_info[gf].Elem1Inf = 64 * lf;
4555 faces_info[gf].Elem2No = -1;
4556 faces_info[gf].Elem2Inf = -1;
4560 faces_info[gf].Elem2No = el;
4561 faces_info[gf].Elem2Inf = 64 * lf + 1;
4565 void Mesh::AddSegmentFaceElement(
int lf,
int gf,
int el,
int v0,
int v1)
4567 if (faces[gf] == NULL)
4569 faces[gf] =
new Segment(v0, v1);
4570 faces_info[gf].Elem1No = el;
4571 faces_info[gf].Elem1Inf = 64 * lf;
4572 faces_info[gf].Elem2No = -1;
4573 faces_info[gf].Elem2Inf = -1;
4577 int *v = faces[gf]->GetVertices();
4578 faces_info[gf].Elem2No = el;
4579 if ( v[1] == v0 && v[0] == v1 )
4581 faces_info[gf].Elem2Inf = 64 * lf + 1;
4583 else if ( v[0] == v0 && v[1] == v1 )
4589 faces_info[gf].Elem2Inf = 64 * lf;
4593 MFEM_ABORT(
"internal error");
4598 void Mesh::AddTriangleFaceElement(
int lf,
int gf,
int el,
4599 int v0,
int v1,
int v2)
4601 if (faces[gf] == NULL)
4603 faces[gf] =
new Triangle(v0, v1, v2);
4604 faces_info[gf].Elem1No = el;
4605 faces_info[gf].Elem1Inf = 64 * lf;
4606 faces_info[gf].Elem2No = -1;
4607 faces_info[gf].Elem2Inf = -1;
4611 int orientation, vv[3] = { v0, v1, v2 };
4612 orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
4617 faces_info[gf].Elem2No = el;
4618 faces_info[gf].Elem2Inf = 64 * lf + orientation;
4622 void Mesh::AddQuadFaceElement(
int lf,
int gf,
int el,
4623 int v0,
int v1,
int v2,
int v3)
4625 if (faces_info[gf].Elem1No < 0)
4628 faces_info[gf].Elem1No = el;
4629 faces_info[gf].Elem1Inf = 64 * lf;
4630 faces_info[gf].Elem2No = -1;
4631 faces_info[gf].Elem2Inf = -1;
4635 int vv[4] = { v0, v1, v2, v3 };
4636 int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
4640 faces_info[gf].Elem2No = el;
4641 faces_info[gf].Elem2Inf = 64 * lf + oo;
4645 void Mesh::GenerateFaces()
4647 int i, nfaces = GetNumFaces();
4649 for (i = 0; i < faces.Size(); i++)
4651 FreeElement(faces[i]);
4655 faces.SetSize(nfaces);
4656 faces_info.SetSize(nfaces);
4657 for (i = 0; i < nfaces; i++)
4660 faces_info[i].Elem1No = -1;
4661 faces_info[i].NCFace = -1;
4663 for (i = 0; i < NumOfElements; i++)
4665 const int *v = elements[i]->GetVertices();
4669 AddPointFaceElement(0, v[0], i);
4670 AddPointFaceElement(1, v[1], i);
4674 ef = el_to_edge->GetRow(i);
4675 const int ne = elements[i]->GetNEdges();
4676 for (
int j = 0; j < ne; j++)
4678 const int *e = elements[i]->GetEdgeVertices(j);
4679 AddSegmentFaceElement(j, ef[j], i, v[e[0]], v[e[1]]);
4684 ef = el_to_face->GetRow(i);
4685 switch (GetElementType(i))
4687 case Element::TETRAHEDRON:
4689 for (
int j = 0; j < 4; j++)
4691 const int *fv = tet_t::FaceVert[j];
4692 AddTriangleFaceElement(j, ef[j], i,
4693 v[fv[0]], v[fv[1]], v[fv[2]]);
4697 case Element::WEDGE:
4699 for (
int j = 0; j < 2; j++)
4701 const int *fv = pri_t::FaceVert[j];
4702 AddTriangleFaceElement(j, ef[j], i,
4703 v[fv[0]], v[fv[1]], v[fv[2]]);
4705 for (
int j = 2; j < 5; j++)
4707 const int *fv = pri_t::FaceVert[j];
4708 AddQuadFaceElement(j, ef[j], i,
4709 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4713 case Element::HEXAHEDRON:
4715 for (
int j = 0; j < 6; j++)
4717 const int *fv = hex_t::FaceVert[j];
4718 AddQuadFaceElement(j, ef[j], i,
4719 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4724 MFEM_ABORT(
"Unexpected type of Element.");
4730 void Mesh::GenerateNCFaceInfo()
4732 MFEM_VERIFY(ncmesh,
"missing NCMesh.");
4734 for (
int i = 0; i < faces_info.Size(); i++)
4736 faces_info[i].NCFace = -1;
4740 (Dim == 2) ? ncmesh->GetEdgeList() : ncmesh->GetFaceList();
4742 nc_faces_info.SetSize(0);
4743 nc_faces_info.Reserve(list.
masters.size() + list.
slaves.size());
4745 int nfaces = GetNumFaces();
4748 for (
unsigned i = 0; i < list.
masters.size(); i++)
4751 if (master.
index >= nfaces) {
continue; }
4753 faces_info[master.
index].NCFace = nc_faces_info.Size();
4759 for (
unsigned i = 0; i < list.
slaves.size(); i++)
4762 if (slave.
index >= nfaces || slave.
master >= nfaces) {
continue; }
4768 slave_fi.
NCFace = nc_faces_info.Size();
4780 for (
int i = 0; i < NumOfElements; i++)
4782 const int *v = elements[i]->GetVertices();
4783 switch (GetElementType(i))
4785 case Element::TETRAHEDRON:
4787 for (
int j = 0; j < 4; j++)
4789 const int *fv = tet_t::FaceVert[j];
4790 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
4794 case Element::WEDGE:
4796 for (
int j = 0; j < 2; j++)
4798 const int *fv = pri_t::FaceVert[j];
4799 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
4801 for (
int j = 2; j < 5; j++)
4803 const int *fv = pri_t::FaceVert[j];
4804 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4808 case Element::HEXAHEDRON:
4812 for (
int j = 0; j < 6; j++)
4814 const int *fv = hex_t::FaceVert[j];
4815 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4820 MFEM_ABORT(
"Unexpected type of Element.");
4831 if (el_to_face != NULL)
4835 el_to_face =
new Table(NumOfElements, 6);
4836 faces_tbl =
new STable3D(NumOfVertices);
4837 for (i = 0; i < NumOfElements; i++)
4839 v = elements[i]->GetVertices();
4840 switch (GetElementType(i))
4842 case Element::TETRAHEDRON:
4844 for (
int j = 0; j < 4; j++)
4846 const int *fv = tet_t::FaceVert[j];
4848 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
4852 case Element::WEDGE:
4854 for (
int j = 0; j < 2; j++)
4856 const int *fv = pri_t::FaceVert[j];
4858 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
4860 for (
int j = 2; j < 5; j++)
4862 const int *fv = pri_t::FaceVert[j];
4864 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
4868 case Element::HEXAHEDRON:
4872 for (
int j = 0; j < 6; j++)
4874 const int *fv = hex_t::FaceVert[j];
4876 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
4881 MFEM_ABORT(
"Unexpected type of Element.");
4884 el_to_face->Finalize();
4886 be_to_face.SetSize(NumOfBdrElements);
4887 for (i = 0; i < NumOfBdrElements; i++)
4889 v = boundary[i]->GetVertices();
4890 switch (GetBdrElementType(i))
4892 case Element::TRIANGLE:
4894 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
4897 case Element::QUADRILATERAL:
4899 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
4903 MFEM_ABORT(
"Unexpected type of boundary Element.");
4915 void Mesh::ReorientTetMesh()
4919 if (Dim != 3 || !(meshgen & 1))
4925 Table *old_elem_vert = NULL;
4929 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
4932 for (
int i = 0; i < NumOfElements; i++)
4934 if (GetElementType(i) == Element::TETRAHEDRON)
4936 v = elements[i]->GetVertices();
4938 Rotate3(v[0], v[1], v[2]);
4941 Rotate3(v[1], v[2], v[3]);
4945 ShiftL2R(v[0], v[1], v[3]);
4950 for (
int i = 0; i < NumOfBdrElements; i++)
4952 if (GetBdrElementType(i) == Element::TRIANGLE)
4954 v = boundary[i]->GetVertices();
4956 Rotate3(v[0], v[1], v[2]);
4962 GetElementToFaceTable();
4966 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
4971 DoNodeReorder(old_v_to_v, old_elem_vert);
4972 delete old_elem_vert;
4977 int *Mesh::CartesianPartitioning(
int nxyz[])
4983 for (
int vi = 0; vi < NumOfVertices; vi++)
4985 const double *p = vertices[vi]();
4986 for (
int i = 0; i < spaceDim; i++)
4988 if (p[i] < pmin[i]) { pmin[i] = p[i]; }
4989 if (p[i] > pmax[i]) { pmax[i] = p[i]; }
4993 partitioning =
new int[NumOfElements];
4997 Vector pt(ppt, spaceDim);
4998 for (
int el = 0; el < NumOfElements; el++)
5000 GetElementTransformation(el)->Transform(
5003 for (
int i = spaceDim-1; i >= 0; i--)
5005 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
5006 if (idx < 0) { idx = 0; }
5007 if (idx >= nxyz[i]) { idx = nxyz[i]-1; }
5008 part = part * nxyz[i] + idx;
5010 partitioning[el] = part;
5013 return partitioning;
5016 int *Mesh::GeneratePartitioning(
int nparts,
int part_method)
5018 #ifdef MFEM_USE_METIS
5019 int i, *partitioning;
5021 ElementToElementTable();
5023 partitioning =
new int[NumOfElements];
5027 for (i = 0; i < NumOfElements; i++)
5029 partitioning[i] = 0;
5032 else if (NumOfElements <= nparts)
5034 for (i = 0; i < NumOfElements; i++)
5036 partitioning[i] = i;
5042 #ifndef MFEM_USE_METIS_5
5054 bool freedata =
false;
5056 idx_t *mpartitioning;
5059 if (
sizeof(
idx_t) ==
sizeof(int))
5061 I = (
idx_t*) el_to_el->GetI();
5062 J = (
idx_t*) el_to_el->GetJ();
5063 mpartitioning = (
idx_t*) partitioning;
5067 int *iI = el_to_el->GetI();
5068 int *iJ = el_to_el->GetJ();
5072 for (
int k = 0; k < n+1; k++) { I[k] = iI[k]; }
5073 for (
int k = 0; k < m; k++) { J[k] = iJ[k]; }
5074 mpartitioning =
new idx_t[n];
5077 #ifndef MFEM_USE_METIS_5
5080 METIS_SetDefaultOptions(options);
5081 options[METIS_OPTION_CONTIG] = 1;
5085 if (part_method >= 0 && part_method <= 2)
5087 for (i = 0; i < n; i++)
5093 std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
5099 if (part_method == 0 || part_method == 3)
5101 #ifndef MFEM_USE_METIS_5
5129 " error in METIS_PartGraphRecursive!");
5135 if (part_method == 1 || part_method == 4)
5137 #ifndef MFEM_USE_METIS_5
5165 " error in METIS_PartGraphKway!");
5171 if (part_method == 2 || part_method == 5)
5173 #ifndef MFEM_USE_METIS_5
5186 options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_VOL;
5202 " error in METIS_PartGraphKway!");
5207 mfem::out <<
"Mesh::GeneratePartitioning(...): edgecut = "
5210 nparts = (int) mparts;
5211 if (mpartitioning != (
idx_t*)partitioning)
5213 for (
int k = 0; k<NumOfElements; k++) { partitioning[k] = mpartitioning[k]; }
5219 delete[] mpartitioning;
5232 for (i = 0; i < nparts; i++)
5238 for (i = 0; i < NumOfElements; i++)
5240 psize[partitioning[i]].one++;
5243 int empty_parts = 0;
5244 for (i = 0; i < nparts; i++)
5245 if (psize[i].one == 0)
5254 mfem::err <<
"Mesh::GeneratePartitioning returned " << empty_parts
5255 <<
" empty parts!" << endl;
5257 SortPairs<int,int>(psize, nparts);
5259 for (i = nparts-1; i > nparts-1-empty_parts; i--)
5264 for (
int j = 0; j < NumOfElements; j++)
5265 for (i = nparts-1; i > nparts-1-empty_parts; i--)
5266 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
5272 partitioning[j] = psize[nparts-1-i].two;
5278 return partitioning;
5282 mfem_error(
"Mesh::GeneratePartitioning(...): "
5283 "MFEM was compiled without Metis.");
5297 int num_elem, *i_elem_elem, *j_elem_elem;
5299 num_elem = elem_elem.
Size();
5300 i_elem_elem = elem_elem.
GetI();
5301 j_elem_elem = elem_elem.
GetJ();
5306 int stack_p, stack_top_p, elem;
5310 for (i = 0; i < num_elem; i++)
5312 if (partitioning[i] > num_part)
5314 num_part = partitioning[i];
5321 for (i = 0; i < num_part; i++)
5328 for (elem = 0; elem < num_elem; elem++)
5330 if (component[elem] >= 0)
5335 component[elem] = num_comp[partitioning[elem]]++;
5337 elem_stack[stack_top_p++] = elem;
5339 for ( ; stack_p < stack_top_p; stack_p++)
5341 i = elem_stack[stack_p];
5342 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
5345 if (partitioning[k] == partitioning[i])
5347 if (component[k] < 0)
5349 component[k] = component[i];
5350 elem_stack[stack_top_p++] = k;
5352 else if (component[k] != component[i])
5362 void Mesh::CheckPartitioning(
int *partitioning)
5364 int i, n_empty, n_mcomp;
5366 const Array<int> _partitioning(partitioning, GetNE());
5368 ElementToElementTable();
5372 n_empty = n_mcomp = 0;
5373 for (i = 0; i < num_comp.
Size(); i++)
5374 if (num_comp[i] == 0)
5378 else if (num_comp[i] > 1)
5385 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
5386 <<
"The following subdomains are empty :\n";
5387 for (i = 0; i < num_comp.
Size(); i++)
5388 if (num_comp[i] == 0)
5396 mfem::out <<
"Mesh::CheckPartitioning(...) :\n"
5397 <<
"The following subdomains are NOT connected :\n";
5398 for (i = 0; i < num_comp.
Size(); i++)
5399 if (num_comp[i] > 1)
5405 if (n_empty == 0 && n_mcomp == 0)
5406 mfem::out <<
"Mesh::CheckPartitioning(...) : "
5407 "All subdomains are connected." << endl;
5421 const double *a = A.
Data();
5422 const double *b = B.
Data();
5431 c(0) = a[0]*a[3]-a[1]*a[2];
5432 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
5433 c(2) = b[0]*b[3]-b[1]*b[2];
5454 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
5455 a[1] * (a[5] * a[6] - a[3] * a[8]) +
5456 a[2] * (a[3] * a[7] - a[4] * a[6]));
5458 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
5459 b[1] * (a[5] * a[6] - a[3] * a[8]) +
5460 b[2] * (a[3] * a[7] - a[4] * a[6]) +
5462 a[0] * (b[4] * a[8] - b[5] * a[7]) +
5463 a[1] * (b[5] * a[6] - b[3] * a[8]) +
5464 a[2] * (b[3] * a[7] - b[4] * a[6]) +
5466 a[0] * (a[4] * b[8] - a[5] * b[7]) +
5467 a[1] * (a[5] * b[6] - a[3] * b[8]) +
5468 a[2] * (a[3] * b[7] - a[4] * b[6]));
5470 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
5471 a[1] * (b[5] * b[6] - b[3] * b[8]) +
5472 a[2] * (b[3] * b[7] - b[4] * b[6]) +
5474 b[0] * (a[4] * b[8] - a[5] * b[7]) +
5475 b[1] * (a[5] * b[6] - a[3] * b[8]) +
5476 b[2] * (a[3] * b[7] - a[4] * b[6]) +
5478 b[0] * (b[4] * a[8] - b[5] * a[7]) +
5479 b[1] * (b[5] * a[6] - b[3] * a[8]) +
5480 b[2] * (b[3] * a[7] - b[4] * a[6]));
5482 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
5483 b[1] * (b[5] * b[6] - b[3] * b[8]) +
5484 b[2] * (b[3] * b[7] - b[4] * b[6]));
5530 double a = z(2), b = z(1), c = z(0);
5531 double D = b*b-4*a*c;
5538 x(0) = x(1) = -0.5 * b / a;
5543 x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
5551 t = -0.5 * (b + sqrt(D));
5555 t = -0.5 * (b - sqrt(D));
5561 Swap<double>(x(0), x(1));
5569 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
5572 double Q = (a * a - 3 * b) / 9;
5573 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
5574 double Q3 = Q * Q * Q;
5581 x(0) = x(1) = x(2) = - a / 3;
5585 double sqrtQ = sqrt(Q);
5589 x(0) = -2 * sqrtQ - a / 3;
5590 x(1) = x(2) = sqrtQ - a / 3;
5594 x(0) = x(1) = - sqrtQ - a / 3;
5595 x(2) = 2 * sqrtQ - a / 3;
5602 double theta = acos(R / sqrt(Q3));
5603 double A = -2 * sqrt(Q);
5605 x0 = A * cos(theta / 3) - a / 3;
5606 x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
5607 x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
5612 Swap<double>(x0, x1);
5616 Swap<double>(x1, x2);
5619 Swap<double>(x0, x1);
5632 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
5636 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
5638 x(0) = A + Q / A - a / 3;
5647 const double factor,
const int Dim)
5649 const double c0 = c(0);
5650 c(0) = c0 * (1.0 - pow(factor, -Dim));
5652 for (
int j = 0; j < nr; j++)
5664 c(0) = c0 * (1.0 - pow(factor, Dim));
5666 for (
int j = 0; j < nr; j++)
5680 void Mesh::CheckDisplacements(
const Vector &displacements,
double &tmax)
5682 int nvs = vertices.Size();
5683 DenseMatrix P, V, DS, PDS(spaceDim), VDS(spaceDim);
5684 Vector c(spaceDim+1), x(spaceDim);
5685 const double factor = 2.0;
5692 for (
int i = 0; i < NumOfElements; i++)
5699 for (
int j = 0; j < spaceDim; j++)
5700 for (
int k = 0; k < nv; k++)
5702 P(j, k) = vertices[v[k]](j);
5703 V(j, k) = displacements(v[k]+j*nvs);
5707 GetTransformationFEforElementType(el->
GetType());
5711 case Element::TRIANGLE:
5712 case Element::TETRAHEDRON:
5730 case Element::QUADRILATERAL:
5733 for (
int j = 0; j < nv; j++)
5757 void Mesh::MoveVertices(
const Vector &displacements)
5759 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
5760 for (
int j = 0; j < spaceDim; j++)
5762 vertices[i](j) += displacements(j*nv+i);
5766 void Mesh::GetVertices(
Vector &vert_coord)
const
5768 int nv = vertices.Size();
5769 vert_coord.
SetSize(nv*spaceDim);
5770 for (
int i = 0; i < nv; i++)
5771 for (
int j = 0; j < spaceDim; j++)
5773 vert_coord(j*nv+i) = vertices[i](j);
5777 void Mesh::SetVertices(
const Vector &vert_coord)
5779 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
5780 for (
int j = 0; j < spaceDim; j++)
5782 vertices[i](j) = vert_coord(j*nv+i);
5786 void Mesh::GetNode(
int i,
double *coord)
5791 for (
int j = 0; j < spaceDim; j++)
5793 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
5798 for (
int j = 0; j < spaceDim; j++)
5800 coord[j] = vertices[i](j);
5805 void Mesh::SetNode(
int i,
const double *coord)
5810 for (
int j = 0; j < spaceDim; j++)
5812 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
5817 for (
int j = 0; j < spaceDim; j++)
5819 vertices[i](j) = coord[j];
5825 void Mesh::MoveNodes(
const Vector &displacements)
5829 (*Nodes) += displacements;
5833 MoveVertices(displacements);
5837 void Mesh::GetNodes(
Vector &node_coord)
const
5841 node_coord = (*Nodes);
5845 GetVertices(node_coord);
5849 void Mesh::SetNodes(
const Vector &node_coord)
5853 (*Nodes) = node_coord;
5857 SetVertices(node_coord);
5863 if (own_nodes) {
delete Nodes; }
5866 own_nodes = (int)make_owner;
5877 mfem::Swap<GridFunction*>(Nodes, nodes);
5878 mfem::Swap<int>(own_nodes, own_nodes_);
5885 void Mesh::AverageVertices(
const int *indexes,
int n,
int result)
5889 for (k = 0; k < spaceDim; k++)
5891 vertices[result](k) = vertices[indexes[0]](k);
5894 for (j = 1; j < n; j++)
5895 for (k = 0; k < spaceDim; k++)
5897 vertices[result](k) += vertices[indexes[j]](k);
5900 for (k = 0; k < spaceDim; k++)
5902 vertices[result](k) *= (1.0 / n);
5906 void Mesh::UpdateNodes()
5910 Nodes->FESpace()->Update();
5915 void Mesh::UniformRefinement2D()
5919 if (el_to_edge == NULL)
5921 el_to_edge =
new Table;
5922 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5925 int quad_counter = 0;
5926 for (
int i = 0; i < NumOfElements; i++)
5928 if (elements[i]->GetType() == Element::QUADRILATERAL)
5934 const int oedge = NumOfVertices;
5935 const int oelem = oedge + NumOfEdges;
5937 vertices.SetSize(oelem + quad_counter);
5938 elements.SetSize(4 * NumOfElements);
5940 for (
int i = 0; i < NumOfElements; i++)
5943 const int attr = elements[i]->GetAttribute();
5944 int *v = elements[i]->GetVertices();
5945 const int *e = el_to_edge->GetRow(i);
5946 const int j = NumOfElements + 3 * i;
5949 if (el_type == Element::TRIANGLE)
5951 for (
int ei = 0; ei < 3; ei++)
5953 for (
int k = 0; k < 2; k++)
5955 vv[k] = v[tri_t::Edges[ei][k]];
5957 AverageVertices(vv, 2, oedge+e[ei]);
5960 elements[j+0] =
new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
5961 elements[j+1] =
new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
5962 elements[j+2] =
new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
5967 else if (el_type == Element::QUADRILATERAL)
5969 const int qe = quad_counter;
5971 AverageVertices(v, 4, oelem+qe);
5973 for (
int ei = 0; ei < 4; ei++)
5975 for (
int k = 0; k < 2; k++)
5977 vv[k] = v[quad_t::Edges[ei][k]];
5979 AverageVertices(vv, 2, oedge+e[ei]);
5982 elements[j+0] =
new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
5985 v[2], oedge+e[2], attr);
5987 oedge+e[2], v[3], attr);
5995 MFEM_ABORT(
"unknown element type: " << el_type);
5999 boundary.SetSize(2 * NumOfBdrElements);
6000 for (
int i = 0; i < NumOfBdrElements; i++)
6002 const int attr = boundary[i]->GetAttribute();
6003 int *v = boundary[i]->GetVertices();
6004 const int j = NumOfBdrElements + i;
6006 boundary[j] =
new Segment(oedge+be_to_edge[i], v[1], attr);
6008 v[1] = oedge+be_to_edge[i];
6011 static const double A = 0.0, B = 0.5, C = 1.0;
6012 static double tri_children[2*3*4] =
6019 static double quad_children[2*4*4] =
6027 CoarseFineTr.point_matrices[Geometry::TRIANGLE].
6028 UseExternalData(tri_children, 2, 3, 4);
6029 CoarseFineTr.point_matrices[Geometry::SQUARE].
6030 UseExternalData(quad_children, 2, 4, 4);
6031 CoarseFineTr.embeddings.SetSize(elements.Size());
6033 for (
int i = 0; i < elements.Size(); i++)
6035 Embedding &emb = CoarseFineTr.embeddings[i];
6036 emb.
parent = (i < NumOfElements) ? i : (i - NumOfElements) / 3;
6037 emb.
matrix = (i < NumOfElements) ? 0 : (i - NumOfElements) % 3 + 1;
6040 NumOfVertices = vertices.Size();
6041 NumOfElements = 4 * NumOfElements;
6042 NumOfBdrElements = 2 * NumOfBdrElements;
6045 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6048 last_operation = Mesh::REFINE;
6054 CheckElementOrientation(
false);
6055 CheckBdrElementOrientation(
false);
6059 static inline double sqr(
const double &x)
6068 if (el_to_edge == NULL)
6070 el_to_edge =
new Table;
6071 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6074 if (el_to_face == NULL)
6076 GetElementToFaceTable();
6080 Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
6083 int NumOfQuadFaces = 0;
6084 if (HasGeometry(Geometry::SQUARE))
6086 if (HasGeometry(Geometry::TRIANGLE))
6089 for (
int i = 0; i < faces.Size(); i++)
6091 if (faces[i]->GetType() == Element::QUADRILATERAL)
6093 f2qf[i] = NumOfQuadFaces;
6100 NumOfQuadFaces = faces.
Size();
6104 int hex_counter = 0;
6105 if (HasGeometry(Geometry::CUBE))
6107 for (
int i = 0; i < elements.Size(); i++)
6109 if (elements[i]->GetType() == Element::HEXAHEDRON)
6119 if (HasGeometry(Geometry::TETRAHEDRON))
6123 DSTable *v_to_v_ptr = v_to_v_p;
6126 v_to_v_ptr =
new DSTable(NumOfVertices);
6127 GetVertexToVertexTable(*v_to_v_ptr);
6132 for (
int i = 0; i < NumOfVertices; i++)
6139 std::sort(row_start, J_v2v.
end());
6142 for (
int i = 0; i < J_v2v.
Size(); i++)
6144 e2v[J_v2v[i].two] = i;
6153 for (
int i = 0; i < NumOfVertices; i++)
6157 it.SetIndex(e2v[it.Index()]);
6165 const int oedge = NumOfVertices;
6166 const int oface = oedge + NumOfEdges;
6167 const int oelem = oface + NumOfQuadFaces;
6169 vertices.SetSize(oelem + hex_counter);
6170 elements.SetSize(8 * NumOfElements);
6171 CoarseFineTr.embeddings.SetSize(elements.Size());
6173 for (
int i = 0; i < NumOfElements; i++)
6176 const int attr = elements[i]->GetAttribute();
6177 int *v = elements[i]->GetVertices();
6178 const int *e = el_to_edge->GetRow(i);
6179 const int j = NumOfElements + 7 * i;
6184 const int ne = el_to_edge->RowSize(i);
6185 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
6191 case Element::TETRAHEDRON:
6193 for (
int ei = 0; ei < 6; ei++)
6195 for (
int k = 0; k < 2; k++)
6197 vv[k] = v[tet_t::Edges[ei][k]];
6199 AverageVertices(vv, 2, oedge+e[ei]);
6205 const int rt_algo = 1;
6216 double len_sqr, min_len;
6218 min_len = sqr(J(0,0)-J(0,1)-J(0,2)) +
6219 sqr(J(1,0)-J(1,1)-J(1,2)) +
6220 sqr(J(2,0)-J(2,1)-J(2,2));
6223 len_sqr = sqr(J(0,1)-J(0,0)-J(0,2)) +
6224 sqr(J(1,1)-J(1,0)-J(1,2)) +
6225 sqr(J(2,1)-J(2,0)-J(2,2));
6226 if (len_sqr < min_len) { min_len = len_sqr; rt = 1; }
6228 len_sqr = sqr(J(0,2)-J(0,0)-J(0,1)) +
6229 sqr(J(1,2)-J(1,0)-J(1,1)) +
6230 sqr(J(2,2)-J(2,0)-J(2,1));
6231 if (len_sqr < min_len) { rt = 2; }
6236 double Em_data[18], Js_data[9], Jp_data[9];
6239 double ar1, ar2,
kappa, kappa_min;
6241 for (
int s = 0; s < 3; s++)
6243 for (
int t = 0; t < 3; t++)
6245 Em(t,s) = 0.5*J(t,s);
6248 for (
int t = 0; t < 3; t++)
6250 Em(t,3) = 0.5*(J(t,0)+J(t,1));
6251 Em(t,4) = 0.5*(J(t,0)+J(t,2));
6252 Em(t,5) = 0.5*(J(t,1)+J(t,2));
6256 for (
int t = 0; t < 3; t++)
6258 Js(t,0) = Em(t,5)-Em(t,0);
6259 Js(t,1) = Em(t,1)-Em(t,0);
6260 Js(t,2) = Em(t,2)-Em(t,0);
6264 for (
int t = 0; t < 3; t++)
6266 Js(t,0) = Em(t,5)-Em(t,0);
6267 Js(t,1) = Em(t,2)-Em(t,0);
6268 Js(t,2) = Em(t,4)-Em(t,0);
6272 kappa_min = std::max(ar1, ar2);
6276 for (
int t = 0; t < 3; t++)
6278 Js(t,0) = Em(t,0)-Em(t,1);
6279 Js(t,1) = Em(t,4)-Em(t,1);
6280 Js(t,2) = Em(t,2)-Em(t,1);
6284 for (
int t = 0; t < 3; t++)
6286 Js(t,0) = Em(t,2)-Em(t,1);
6287 Js(t,1) = Em(t,4)-Em(t,1);
6288 Js(t,2) = Em(t,5)-Em(t,1);
6292 kappa = std::max(ar1, ar2);
6293 if (kappa < kappa_min) { kappa_min =
kappa; rt = 1; }
6296 for (
int t = 0; t < 3; t++)
6298 Js(t,0) = Em(t,0)-Em(t,2);
6299 Js(t,1) = Em(t,1)-Em(t,2);
6300 Js(t,2) = Em(t,3)-Em(t,2);
6304 for (
int t = 0; t < 3; t++)
6306 Js(t,0) = Em(t,1)-Em(t,2);
6307 Js(t,1) = Em(t,5)-Em(t,2);
6308 Js(t,2) = Em(t,3)-Em(t,2);
6312 kappa = std::max(ar1, ar2);
6313 if (kappa < kappa_min) { rt = 2; }
6316 static const int mv_all[3][4][4] =
6318 { {0,5,1,2}, {0,5,2,4}, {0,5,4,3}, {0,5,3,1} },
6319 { {1,0,4,2}, {1,2,4,5}, {1,5,4,3}, {1,3,4,0} },
6320 { {2,0,1,3}, {2,1,5,3}, {2,5,4,3}, {2,4,0,3} }
6322 const int (&mv)[4][4] = mv_all[rt];
6324 #ifndef MFEM_USE_MEMALLOC
6326 oedge+e[3], oedge+e[4], attr);
6327 elements[j+1] =
new Tetrahedron(oedge+e[1], oedge+e[3],
6328 v[2], oedge+e[5], attr);
6329 elements[j+2] =
new Tetrahedron(oedge+e[2], oedge+e[4],
6330 oedge+e[5], v[3], attr);
6331 for (
int k = 0; k < 4; k++)
6334 new Tetrahedron(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
6335 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
6339 elements[j+0] = tet = TetMemory.Alloc();
6340 tet->
Init(oedge+e[0], v[1], oedge+e[3], oedge+e[4], attr);
6341 elements[j+1] = tet = TetMemory.Alloc();
6342 tet->
Init(oedge+e[1], oedge+e[3], v[2], oedge+e[5], attr);
6343 elements[j+2] = tet = TetMemory.Alloc();
6344 tet->
Init(oedge+e[2], oedge+e[4], oedge+e[5], v[3], attr);
6345 for (
int k = 0; k < 4; k++)
6347 elements[j+k+3] = tet = TetMemory.Alloc();
6348 tet->
Init(oedge+e[mv[k][0]], oedge+e[mv[k][1]],
6349 oedge+e[mv[k][2]], oedge+e[mv[k][3]], attr);
6356 ((
Tetrahedron*)elements[i])->SetRefinementFlag(0);
6358 CoarseFineTr.embeddings[i].parent = i;
6359 CoarseFineTr.embeddings[i].matrix = 0;
6360 for (
int k = 0; k < 3; k++)
6362 CoarseFineTr.embeddings[j+k].parent = i;
6363 CoarseFineTr.embeddings[j+k].matrix = k+1;
6365 for (
int k = 0; k < 4; k++)
6367 CoarseFineTr.embeddings[j+k+3].parent = i;
6368 CoarseFineTr.embeddings[j+k+3].matrix = 4*(rt+1)+k;
6373 case Element::WEDGE:
6375 const int *f = el_to_face->GetRow(i);
6377 for (
int fi = 2; fi < 5; fi++)
6379 for (
int k = 0; k < 4; k++)
6381 vv[k] = v[pri_t::FaceVert[fi][k]];
6383 AverageVertices(vv, 4, oface + f2qf[f[fi]]);
6386 for (
int ei = 0; ei < 9; ei++)
6388 for (
int k = 0; k < 2; k++)
6390 vv[k] = v[pri_t::Edges[ei][k]];
6392 AverageVertices(vv, 2, oedge+e[ei]);
6395 const int qf2 = f2qf[f[2]];
6396 const int qf3 = f2qf[f[3]];
6397 const int qf4 = f2qf[f[4]];
6399 elements[j+0] =
new Wedge(oedge+e[1], oedge+e[2], oedge+e[0],
6400 oface+qf3, oface+qf4, oface+qf2,
6402 elements[j+1] =
new Wedge(oedge+e[0], v[1], oedge+e[1],
6403 oface+qf2, oedge+e[7], oface+qf3,
6405 elements[j+2] =
new Wedge(oedge+e[2], oedge+e[1], v[2],
6406 oface+qf4, oface+qf3, oedge+e[8],
6408 elements[j+3] =
new Wedge(oedge+e[6], oface+qf2, oface+qf4,
6409 v[3], oedge+e[3], oedge+e[5],
6411 elements[j+4] =
new Wedge(oface+qf3, oface+qf4, oface+qf2,
6412 oedge+e[4], oedge+e[5], oedge+e[3],
6414 elements[j+5] =
new Wedge(oface+qf2, oedge+e[7], oface+qf3,
6415 oedge+e[3], v[4], oedge+e[4],
6417 elements[j+6] =
new Wedge(oface+qf4, oface+qf3, oedge+e[8],
6418 oedge+e[5], oedge+e[4], v[5],
6429 case Element::HEXAHEDRON:
6431 const int *f = el_to_face->GetRow(i);
6432 const int he = hex_counter;
6437 if (f2qf.
Size() == 0)
6443 for (
int k = 0; k < 6; k++) { qf_data[k] = f2qf[f[k]]; }
6447 AverageVertices(v, 8, oelem+he);
6449 for (
int fi = 0; fi < 6; fi++)
6451 for (
int k = 0; k < 4; k++)
6453 vv[k] = v[hex_t::FaceVert[fi][k]];
6455 AverageVertices(vv, 4, oface + qf[fi]);
6458 for (
int ei = 0; ei < 12; ei++)
6460 for (
int k = 0; k < 2; k++)
6462 vv[k] = v[hex_t::Edges[ei][k]];
6464 AverageVertices(vv, 2, oedge+e[ei]);
6467 elements[j+0] =
new Hexahedron(oedge+e[0], v[1], oedge+e[1],
6468 oface+qf[0], oface+qf[1], oedge+e[9],
6469 oface+qf[2], oelem+he, attr);
6470 elements[j+1] =
new Hexahedron(oface+qf[0], oedge+e[1], v[2],
6471 oedge+e[2], oelem+he, oface+qf[2],
6472 oedge+e[10], oface+qf[3], attr);
6473 elements[j+2] =
new Hexahedron(oedge+e[3], oface+qf[0], oedge+e[2],
6474 v[3], oface+qf[4], oelem+he,
6475 oface+qf[3], oedge+e[11], attr);
6476 elements[j+3] =
new Hexahedron(oedge+e[8], oface+qf[1], oelem+he,
6477 oface+qf[4], v[4], oedge+e[4],
6478 oface+qf[5], oedge+e[7], attr);
6479 elements[j+4] =
new Hexahedron(oface+qf[1], oedge+e[9], oface+qf[2],
6480 oelem+he, oedge+e[4], v[5],
6481 oedge+e[5], oface+qf[5], attr);
6482 elements[j+5] =
new Hexahedron(oelem+he, oface+qf[2], oedge+e[10],
6483 oface+qf[3], oface+qf[5], oedge+e[5],
6484 v[6], oedge+e[6], attr);
6485 elements[j+6] =
new Hexahedron(oface+qf[4], oelem+he, oface+qf[3],
6486 oedge+e[11], oedge+e[7], oface+qf[5],
6487 oedge+e[6], v[7], attr);
6500 MFEM_ABORT(
"Unknown 3D element type \"" << el_type <<
"\"");
6505 boundary.SetSize(4 * NumOfBdrElements);
6506 for (
int i = 0; i < NumOfBdrElements; i++)
6509 const int attr = boundary[i]->GetAttribute();
6510 int *v = boundary[i]->GetVertices();
6511 const int *e = bel_to_edge->GetRow(i);
6512 const int j = NumOfBdrElements + 3 * i;
6517 const int ne = bel_to_edge->RowSize(i);
6518 for (
int k = 0; k < ne; k++) { ev[k] = e2v[e[k]]; }
6522 if (bdr_el_type == Element::TRIANGLE)
6524 boundary[j+0] =
new Triangle(oedge+e[1], oedge+e[2], oedge+e[0], attr);
6525 boundary[j+1] =
new Triangle(oedge+e[0], v[1], oedge+e[1], attr);
6526 boundary[j+2] =
new Triangle(oedge+e[2], oedge+e[1], v[2], attr);
6531 else if (bdr_el_type == Element::QUADRILATERAL)
6534 (f2qf.
Size() == 0) ? be_to_face[i] : f2qf[be_to_face[i]];
6536 boundary[j+0] =
new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
6538 boundary[j+1] =
new Quadrilateral(oface+qf, oedge+e[1], v[2],
6541 oedge+e[2], v[3], attr);
6549 MFEM_ABORT(
"boundary Element is not a triangle or a quad!");
6553 static const double A = 0.0, B = 0.5, C = 1.0;
6554 static double tet_children[3*4*16] =
6556 A,A,A, B,A,A, A,B,A, A,A,B,
6557 B,A,A, C,A,A, B,B,A, B,A,B,
6558 A,B,A, B,B,A, A,C,A, A,B,B,
6559 A,A,B, B,A,B, A,B,B, A,A,C,
6564 B,A,A, A,B,B, A,B,A, A,A,B,
6565 B,A,A, A,B,B, A,A,B, B,A,B,
6566 B,A,A, A,B,B, B,A,B, B,B,A,
6567 B,A,A, A,B,B, B,B,A, A,B,A,
6569 A,B,A, B,A,A, B,A,B, A,A,B,
6570 A,B,A, A,A,B, B,A,B, A,B,B,
6571 A,B,A, A,B,B, B,A,B, B,B,A,
6572 A,B,A, B,B,A, B,A,B, B,A,A,
6574 A,A,B, B,A,A, A,B,A, B,B,A,
6575 A,A,B, A,B,A, A,B,B, B,B,A,
6576 A,A,B, A,B,B, B,A,B, B,B,A,
6577 A,A,B, B,A,B, B,A,A, B,B,A
6579 static double pri_children[3*6*8] =
6581 A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
6582 B,B,A, A,B,A, B,A,A, B,B,B, A,B,B, B,A,B,
6583 B,A,A, C,A,A, B,B,A, B,A,B, C,A,B, B,B,B,
6584 A,B,A, B,B,A, A,C,A, A,B,B, B,B,B, A,C,B,
6585 A,A,B, B,A,B, A,B,B, A,A,C, B,A,C, A,B,C,
6586 B,B,B, A,B,B, B,A,B, B,B,C, A,B,C, B,A,C,
6587 B,A,B, C,A,B, B,B,B, B,A,C, C,A,C, B,B,C,
6588 A,B,B, B,B,B, A,C,B, A,B,C, B,B,C, A,C,C
6590 static double hex_children[3*8*8] =
6592 A,A,A, B,A,A, B,B,A, A,B,A, A,A,B, B,A,B, B,B,B, A,B,B,
6593 B,A,A, C,A,A, C,B,A, B,B,A, B,A,B, C,A,B, C,B,B, B,B,B,
6594 B,B,A, C,B,A, C,C,A, B,C,A, B,B,B, C,B,B, C,C,B, B,C,B,
6595 A,B,A, B,B,A, B,C,A, A,C,A, A,B,B, B,B,B, B,C,B, A,C,B,
6596 A,A,B, B,A,B, B,B,B, A,B,B, A,A,C, B,A,C, B,B,C, A,B,C,
6597 B,A,B, C,A,B, C,B,B, B,B,B, B,A,C, C,A,C, C,B,C, B,B,C,
6598 B,B,B, C,B,B, C,C,B, B,C,B, B,B,C, C,B,C, C,C,C, B,C,C,
6599 A,B,B, B,B,B, B,C,B, A,C,B, A,B,C, B,B,C, B,C,C, A,C,C
6602 CoarseFineTr.point_matrices[Geometry::TETRAHEDRON].
6603 UseExternalData(tet_children, 3, 4, 16);
6604 CoarseFineTr.point_matrices[Geometry::PRISM].
6605 UseExternalData(pri_children, 3, 6, 8);
6606 CoarseFineTr.point_matrices[Geometry::CUBE].
6607 UseExternalData(hex_children, 3, 8, 8);
6609 for (
int i = 0; i < elements.Size(); i++)
6612 if (elements[i]->GetType() == Element::TETRAHEDRON) {
continue; }
6613 Embedding &emb = CoarseFineTr.embeddings[i];
6614 emb.
parent = (i < NumOfElements) ? i : (i - NumOfElements) / 7;
6615 emb.
matrix = (i < NumOfElements) ? 0 : (i - NumOfElements) % 7 + 1;
6618 NumOfVertices = vertices.Size();
6619 NumOfElements = 8 * NumOfElements;
6620 NumOfBdrElements = 4 * NumOfBdrElements;
6622 GetElementToFaceTable();
6626 CheckBdrElementOrientation(
false);
6629 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6631 last_operation = Mesh::REFINE;
6637 void Mesh::LocalRefinement(
const Array<int> &marked_el,
int type)
6639 int i, j, ind, nedges;
6646 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
6649 InitRefinementTransforms();
6653 int cne = NumOfElements, cnv = NumOfVertices;
6654 NumOfVertices += marked_el.
Size();
6655 NumOfElements += marked_el.
Size();
6656 vertices.SetSize(NumOfVertices);
6657 elements.SetSize(NumOfElements);
6658 CoarseFineTr.embeddings.SetSize(NumOfElements);
6660 for (j = 0; j < marked_el.
Size(); j++)
6665 int new_v = cnv + j, new_e = cne + j;
6666 AverageVertices(vert, 2, new_v);
6667 elements[new_e] =
new Segment(new_v, vert[1], attr);
6670 CoarseFineTr.embeddings[i] =
Embedding(i, 1);
6671 CoarseFineTr.embeddings[new_e] =
Embedding(i, 2);
6674 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
6675 CoarseFineTr.point_matrices[Geometry::SEGMENT].
6676 UseExternalData(seg_children, 1, 2, 3);
6684 DSTable v_to_v(NumOfVertices);
6685 GetVertexToVertexTable(v_to_v);
6689 int *edge1 =
new int[nedges];
6690 int *edge2 =
new int[nedges];
6691 int *middle =
new int[nedges];
6693 for (i = 0; i < nedges; i++)
6695 edge1[i] = edge2[i] = middle[i] = -1;
6698 for (i = 0; i < NumOfElements; i++)
6700 elements[i]->GetVertices(v);
6701 for (j = 1; j < v.
Size(); j++)
6703 ind = v_to_v(v[j-1], v[j]);
6704 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
6706 ind = v_to_v(v[0], v[v.
Size()-1]);
6707 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
6711 for (i = 0; i < marked_el.
Size(); i++)
6713 RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
6717 int need_refinement;
6720 need_refinement = 0;
6721 for (i = 0; i < nedges; i++)
6723 if (middle[i] != -1 && edge1[i] != -1)
6725 need_refinement = 1;
6726 GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
6730 while (need_refinement == 1);
6733 int v1[2], v2[2], bisect, temp;
6734 temp = NumOfBdrElements;
6735 for (i = 0; i < temp; i++)
6737 boundary[i]->GetVertices(v);
6738 bisect = v_to_v(v[0], v[1]);
6739 if (middle[bisect] != -1)
6741 if (boundary[i]->GetType() == Element::SEGMENT)
6743 v1[0] = v[0]; v1[1] = middle[bisect];
6744 v2[0] = middle[bisect]; v2[1] = v[1];
6746 boundary[i]->SetVertices(v1);
6747 boundary.
Append(
new Segment(v2, boundary[i]->GetAttribute()));
6750 mfem_error(
"Only bisection of segment is implemented"
6754 NumOfBdrElements = boundary.Size();
6761 if (el_to_edge != NULL)
6763 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6774 MFEM_VERIFY(GetNE() == 0 ||
6775 ((
Tetrahedron*)elements[0])->GetRefinementFlag() != 0,
6776 "tetrahedral mesh is not marked for refinement:"
6777 " call Finalize(true)");
6784 for (i = 0; i < marked_el.
Size(); i++)
6786 Bisection(marked_el[i], v_to_v);
6790 for (i = 0; i < marked_el.
Size(); i++)
6792 Bisection(marked_el[i], v_to_v);
6794 Bisection(NumOfElements - 1, v_to_v);
6795 Bisection(marked_el[i], v_to_v);
6799 for (i = 0; i < marked_el.
Size(); i++)
6801 Bisection(marked_el[i], v_to_v);
6803 ii = NumOfElements - 1;
6804 Bisection(ii, v_to_v);
6805 Bisection(NumOfElements - 1, v_to_v);
6806 Bisection(ii, v_to_v);
6808 Bisection(marked_el[i], v_to_v);
6809 Bisection(NumOfElements-1, v_to_v);
6810 Bisection(marked_el[i], v_to_v);
6816 int need_refinement;
6821 need_refinement = 0;
6824 for (i = 0; i < NumOfElements; i++)
6829 if (elements[i]->NeedRefinement(v_to_v))
6831 need_refinement = 1;
6832 Bisection(i, v_to_v);
6836 while (need_refinement == 1);
6843 need_refinement = 0;
6844 for (i = 0; i < NumOfBdrElements; i++)
6845 if (boundary[i]->NeedRefinement(v_to_v))
6847 need_refinement = 1;
6848 BdrBisection(i, v_to_v);
6851 while (need_refinement == 1);
6853 NumOfVertices = vertices.Size();
6854 NumOfBdrElements = boundary.Size();
6857 if (el_to_edge != NULL)
6859 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
6861 if (el_to_face != NULL)
6863 GetElementToFaceTable();
6869 last_operation = Mesh::REFINE;
6875 CheckElementOrientation(
false);
6882 MFEM_VERIFY(!NURBSext,
"Nonconforming refinement of NURBS meshes is "
6883 "not supported. Project the NURBS to Nodes first.");
6890 ncmesh =
new NCMesh(
this);
6893 if (!refinements.
Size())
6895 last_operation = Mesh::NONE;
6900 ncmesh->MarkCoarseLevel();
6901 ncmesh->Refine(refinements);
6905 ncmesh->LimitNCLevel(nc_limit);
6910 ncmesh->OnMeshUpdated(mesh2);
6914 Swap(*mesh2,
false);
6917 GenerateNCFaceInfo();
6919 last_operation = Mesh::REFINE;
6924 Nodes->FESpace()->Update();
6930 const int *fine,
int nfine,
int op)
6933 for (
int i = 0; i < nfine; i++)
6935 MFEM_VERIFY(fine[i] < elem_error.
Size(),
"");
6937 double err_fine = elem_error[fine[i]];
6940 case 0: error = std::min(error, err_fine);
break;
6941 case 1: error += err_fine;
break;
6942 case 2: error = std::max(error, err_fine);
break;
6949 double threshold,
int nc_limit,
int op)
6951 MFEM_VERIFY(ncmesh,
"Only supported for non-conforming meshes.");
6952 MFEM_VERIFY(!NURBSext,
"Derefinement of NURBS meshes is not supported. "
6953 "Project the NURBS to Nodes first.");
6957 const Table &dt = ncmesh->GetDerefinementTable();
6962 ncmesh->CheckDerefinementNCLevel(dt, level_ok, nc_limit);
6966 for (
int i = 0; i < dt.
Size(); i++)
6968 if (nc_limit > 0 && !level_ok[i]) {
continue; }
6971 AggregateError(elem_error, dt.
GetRow(i), dt.
RowSize(i), op);
6973 if (error < threshold) { derefs.
Append(i); }
6976 if (!derefs.
Size()) {
return false; }
6978 ncmesh->Derefine(derefs);
6981 ncmesh->OnMeshUpdated(mesh2);
6983 Swap(*mesh2,
false);
6986 GenerateNCFaceInfo();
6988 last_operation = Mesh::DEREFINE;
6993 Nodes->FESpace()->Update();
7001 int nc_limit,
int op)
7005 if (Nonconforming())
7007 return NonconformingDerefinement(elem_error, threshold, nc_limit, op);
7011 MFEM_ABORT(
"Derefinement is currently supported for non-conforming "
7017 bool Mesh::DerefineByError(
const Vector &elem_error,
double threshold,
7018 int nc_limit,
int op)
7021 for (
int i = 0; i < tmp.Size(); i++)
7023 tmp[i] = elem_error(i);
7025 return DerefineByError(tmp, threshold, nc_limit, op);
7029 void Mesh::InitFromNCMesh(
const NCMesh &ncmesh)
7038 NumOfVertices = vertices.Size();
7039 NumOfElements = elements.Size();
7040 NumOfBdrElements = boundary.Size();
7044 NumOfEdges = NumOfFaces = 0;
7048 el_to_edge =
new Table;
7049 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
7053 GetElementToFaceTable();
7057 CheckBdrElementOrientation(
false);
7068 InitFromNCMesh(ncmesh);
7121 const int nv = Geometry::NumVerts[
geom];
7123 for (
int i = 0; i < elem_array.
Size(); i++)
7125 if (elem_array[i]->GetGeometryType() ==
geom)
7130 elem_vtx.
SetSize(nv*num_elems);
7134 for (
int i = 0; i < elem_array.
Size(); i++)
7140 elem_vtx.
Append(loc_vtx);
7145 void Mesh::UniformRefinement(
int ref_algo)
7149 NURBSUniformRefinement();
7151 else if (ref_algo == 0 && Dim == 3 && meshgen == 1)
7153 UniformRefinement3D();
7155 else if (meshgen == 1 || ncmesh)
7158 for (
int i = 0; i < elem_to_refine.
Size(); i++)
7160 elem_to_refine[i] = i;
7167 LocalRefinement(elem_to_refine);
7171 GeneralRefinement(elem_to_refine, 1);
7178 case 2: UniformRefinement2D();
break;
7179 case 3: UniformRefinement3D();
break;
7180 default: MFEM_ABORT(
"internal error");
7186 int nonconforming,
int nc_limit)
7192 else if (Dim == 1 || (Dim == 3 && (meshgen & 1)))
7196 else if (nonconforming < 0)
7212 NonconformingRefinement(refinements, nc_limit);
7217 for (
int i = 0; i < refinements.
Size(); i++)
7219 el_to_refine[i] = refinements[i].index;
7223 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
7224 if (rt == 1 || rt == 2 || rt == 4)
7228 else if (rt == 3 || rt == 5 || rt == 6)
7238 LocalRefinement(el_to_refine, type);
7242 void Mesh::GeneralRefinement(
const Array<int> &el_to_refine,
int nonconforming,
7246 for (
int i = 0; i < el_to_refine.
Size(); i++)
7248 refinements[i] =
Refinement(el_to_refine[i]);
7250 GeneralRefinement(refinements, nonconforming, nc_limit);
7253 void Mesh::EnsureNCMesh(
bool triangles_nonconforming)
7255 MFEM_VERIFY(!NURBSext,
"Cannot convert a NURBS mesh to an NC mesh. "
7256 "Project the NURBS to Nodes first.");
7260 if ((meshgen & 2) ||
7261 (triangles_nonconforming && Dim == 2 && (meshgen & 1)))
7263 MFEM_VERIFY(GetNumGeometries(Dim) <= 1,
7264 "mixed meshes are not supported");
7265 ncmesh =
new NCMesh(
this);
7266 ncmesh->OnMeshUpdated(
this);
7267 GenerateNCFaceInfo();
7272 void Mesh::RandomRefinement(
double prob,
bool aniso,
int nonconforming,
7276 for (
int i = 0; i < GetNE(); i++)
7278 if ((
double) rand() / RAND_MAX < prob)
7283 type = (Dim == 3) ? (rand() % 7 + 1) : (rand() % 3 + 1);
7288 GeneralRefinement(refs, nonconforming, nc_limit);
7291 void Mesh::RefineAtVertex(
const Vertex& vert,
double eps,
int nonconforming)
7295 for (
int i = 0; i < GetNE(); i++)
7297 GetElementVertices(i, v);
7298 bool refine =
false;
7299 for (
int j = 0; j < v.
Size(); j++)
7302 for (
int l = 0; l < spaceDim; l++)
7304 double d = vert(l) - vertices[v[j]](l);
7307 if (dist <= eps*eps) { refine =
true;
break; }
7314 GeneralRefinement(refs, nonconforming);
7318 int nonconforming,
int nc_limit)
7320 MFEM_VERIFY(elem_error.
Size() == GetNE(),
"");
7322 for (
int i = 0; i < GetNE(); i++)
7324 if (elem_error[i] > threshold)
7329 if (ReduceInt(refs.Size()))
7331 GeneralRefinement(refs, nonconforming, nc_limit);
7337 bool Mesh::RefineByError(
const Vector &elem_error,
double threshold,
7338 int nonconforming,
int nc_limit)
7342 return RefineByError(tmp, threshold, nonconforming, nc_limit);
7346 void Mesh::Bisection(
int i,
const DSTable &v_to_v,
7347 int *edge1,
int *edge2,
int *middle)
7350 int v[2][4], v_new, bisect, t;
7355 if (t == Element::TRIANGLE)
7362 bisect = v_to_v(vert[0], vert[1]);
7363 MFEM_ASSERT(bisect >= 0,
"");
7365 if (middle[bisect] == -1)
7367 v_new = NumOfVertices++;
7368 for (
int d = 0; d < spaceDim; d++)
7370 V(d) = 0.5 * (vertices[vert[0]](d) + vertices[vert[1]](d));
7376 if (edge1[bisect] == i)
7378 edge1[bisect] = edge2[bisect];
7381 middle[bisect] = v_new;
7385 v_new = middle[bisect];
7393 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
7394 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
7399 elements.Append(tri_new);
7408 int coarse = FindCoarseElement(i);
7409 CoarseFineTr.embeddings[i].parent = coarse;
7410 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7415 bisect = v_to_v(v[1][0], v[1][1]);
7416 MFEM_ASSERT(bisect >= 0,
"");
7418 if (edge1[bisect] == i)
7420 edge1[bisect] = NumOfElements;
7422 else if (edge2[bisect] == i)
7424 edge2[bisect] = NumOfElements;
7431 MFEM_ABORT(
"Bisection for now works only for triangles.");
7438 int v[2][4], v_new, bisect, t;
7443 if (t == Element::TETRAHEDRON)
7445 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
7449 "TETRAHEDRON element is not marked for refinement.");
7454 bisect = v_to_v.
FindId(vert[0], vert[1]);
7457 v_new = NumOfVertices + v_to_v.
GetId(vert[0],vert[1]);
7458 for (j = 0; j < 3; j++)
7460 V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
7466 v_new = NumOfVertices + bisect;
7475 new_redges[0][0] = 2;
7476 new_redges[0][1] = 1;
7477 new_redges[1][0] = 2;
7478 new_redges[1][1] = 1;
7479 int tr1 = -1, tr2 = -1;
7480 switch (old_redges[0])
7483 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
7484 if (type == Tetrahedron::TYPE_PF) { new_redges[0][1] = 4; }
7488 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
7492 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
7495 switch (old_redges[1])
7498 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
7499 if (type == Tetrahedron::TYPE_PF) { new_redges[1][0] = 3; }
7503 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
7507 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
7514 #ifdef MFEM_USE_MEMALLOC
7522 elements.Append(tet2);
7528 int coarse = FindCoarseElement(i);
7529 CoarseFineTr.embeddings[i].parent = coarse;
7530 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7535 case Tetrahedron::TYPE_PU:
7536 new_type = Tetrahedron::TYPE_PF;
break;
7537 case Tetrahedron::TYPE_PF:
7538 new_type = Tetrahedron::TYPE_A;
break;
7540 new_type = Tetrahedron::TYPE_PU;
7550 MFEM_ABORT(
"Bisection with HashTable for now works only for tetrahedra.");
7557 int v[2][3], v_new, bisect, t;
7558 Element *bdr_el = boundary[i];
7561 if (t == Element::TRIANGLE)
7568 bisect = v_to_v.
FindId(vert[0], vert[1]);
7569 MFEM_ASSERT(bisect >= 0,
"");
7570 v_new = NumOfVertices + bisect;
7571 MFEM_ASSERT(v_new != -1,
"");
7575 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
7576 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
7586 MFEM_ABORT(
"Bisection of boundary elements with HashTable works only for"
7591 void Mesh::UniformRefinement(
int i,
const DSTable &v_to_v,
7592 int *edge1,
int *edge2,
int *middle)
7595 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
7598 if (elements[i]->GetType() == Element::TRIANGLE)
7604 bisect[0] = v_to_v(v[0],v[1]);
7605 bisect[1] = v_to_v(v[1],v[2]);
7606 bisect[2] = v_to_v(v[0],v[2]);
7607 MFEM_ASSERT(bisect[0] >= 0 && bisect[1] >= 0 && bisect[2] >= 0,
"");
7609 for (j = 0; j < 3; j++)
7611 if (middle[bisect[j]] == -1)
7613 v_new[j] = NumOfVertices++;
7614 for (
int d = 0; d < spaceDim; d++)
7616 V(d) = (vertices[v[j]](d) + vertices[v[(j+1)%3]](d))/2.;
7622 if (edge1[bisect[j]] == i)
7624 edge1[bisect[j]] = edge2[bisect[j]];
7627 middle[bisect[j]] = v_new[j];
7631 v_new[j] = middle[bisect[j]];
7634 edge1[bisect[j]] = -1;
7640 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
7641 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
7642 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
7643 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
7649 elements.Append(tri1);
7650 elements.Append(tri2);
7651 elements.Append(tri3);
7658 tri2->ResetTransform(code);
7659 tri3->ResetTransform(code);
7663 tri2->PushTransform(1);
7664 tri3->PushTransform(2);
7667 int coarse = FindCoarseElement(i);
7668 CoarseFineTr.embeddings[i] =
Embedding(coarse);
7669 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7670 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7671 CoarseFineTr.embeddings.Append(
Embedding(coarse));
7677 MFEM_ABORT(
"Uniform refinement for now works only for triangles.");
7681 void Mesh::InitRefinementTransforms()
7684 map<Geometry::Type,DenseTensor> &pms = CoarseFineTr.point_matrices;
7685 map<Geometry::Type,DenseTensor>::iterator pms_iter;
7686 for (pms_iter = pms.begin(); pms_iter != pms.end(); ++pms_iter)
7688 pms_iter->second.SetSize(0, 0, 0);
7690 CoarseFineTr.embeddings.SetSize(NumOfElements);
7691 for (
int i = 0; i < NumOfElements; i++)
7693 elements[i]->ResetTransform(0);
7694 CoarseFineTr.embeddings[i] =
Embedding(i);
7698 int Mesh::FindCoarseElement(
int i)
7701 while ((coarse = CoarseFineTr.embeddings[i].parent) != i)
7710 MFEM_VERIFY(GetLastOperation() == Mesh::REFINE,
"");
7714 return ncmesh->GetRefinementTransforms();
7718 for (
int i = 0; i < elem_geoms.
Size(); i++)
7721 if (CoarseFineTr.point_matrices[geom].SizeK()) {
continue; }
7723 if (geom == Geometry::TRIANGLE ||
7724 geom == Geometry::TETRAHEDRON)
7726 std::map<unsigned, int> mat_no;
7730 for (
int i = 0; i < elements.Size(); i++)
7733 unsigned code = elements[i]->GetTransform();
7736 int &matrix = mat_no[code];
7737 if (!matrix) { matrix = mat_no.size(); }
7740 CoarseFineTr.embeddings[i].matrix = index;
7744 pmats.
SetSize(Dim, Dim+1, mat_no.size());
7747 std::map<unsigned, int>::iterator it;
7748 for (it = mat_no.begin(); it != mat_no.end(); ++it)
7750 if (geom == Geometry::TRIANGLE)
7752 Triangle::GetPointMatrix(it->first, pmats(it->second-1));
7756 Tetrahedron::GetPointMatrix(it->first, pmats(it->second-1));
7762 MFEM_ABORT(
"Don't know how to construct CoarseFineTransformations for"
7763 " geom = " << geom);
7768 return CoarseFineTr;
7771 void Mesh::PrintXG(std::ostream &
out)
const
7773 MFEM_ASSERT(Dim==spaceDim,
"2D Manifold meshes not supported");
7782 out <<
"areamesh2\n\n";
7786 out <<
"curved_areamesh2\n\n";
7790 out << NumOfBdrElements <<
'\n';
7791 for (i = 0; i < NumOfBdrElements; i++)
7793 boundary[i]->GetVertices(v);
7795 out << boundary[i]->GetAttribute();
7796 for (j = 0; j < v.
Size(); j++)
7798 out <<
' ' << v[j] + 1;
7804 out << NumOfElements <<
'\n';
7805 for (i = 0; i < NumOfElements; i++)
7807 elements[i]->GetVertices(v);
7809 out << elements[i]->GetAttribute() <<
' ' << v.
Size();
7810 for (j = 0; j < v.
Size(); j++)
7812 out <<
' ' << v[j] + 1;
7820 out << NumOfVertices <<
'\n';
7821 for (i = 0; i < NumOfVertices; i++)
7823 out << vertices[i](0);
7824 for (j = 1; j < Dim; j++)
7826 out <<
' ' << vertices[i](j);
7833 out << NumOfVertices <<
'\n';
7841 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
7849 out <<
"NETGEN_Neutral_Format\n";
7851 out << NumOfVertices <<
'\n';
7852 for (i = 0; i < NumOfVertices; i++)
7854 for (j = 0; j < Dim; j++)
7856 out <<
' ' << vertices[i](j);
7862 out << NumOfElements <<
'\n';
7863 for (i = 0; i < NumOfElements; i++)
7865 nv = elements[i]->GetNVertices();
7866 ind = elements[i]->GetVertices();
7867 out << elements[i]->GetAttribute();
7868 for (j = 0; j < nv; j++)
7870 out <<
' ' << ind[j]+1;
7876 out << NumOfBdrElements <<
'\n';
7877 for (i = 0; i < NumOfBdrElements; i++)
7879 nv = boundary[i]->GetNVertices();
7880 ind = boundary[i]->GetVertices();
7881 out << boundary[i]->GetAttribute();
7882 for (j = 0; j < nv; j++)
7884 out <<
' ' << ind[j]+1;
7889 else if (meshgen == 2)
7895 <<
"1 " << NumOfVertices <<
" " << NumOfElements
7896 <<
" 0 0 0 0 0 0 0\n"
7897 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
7898 <<
"0 0 " << NumOfBdrElements <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
7899 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
7900 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
7902 for (i = 0; i < NumOfVertices; i++)
7903 out << i+1 <<
" 0.0 " << vertices[i](0) <<
' ' << vertices[i](1)
7904 <<
' ' << vertices[i](2) <<
" 0.0\n";
7906 for (i = 0; i < NumOfElements; i++)
7908 nv = elements[i]->GetNVertices();
7909 ind = elements[i]->GetVertices();
7910 out << i+1 <<
' ' << elements[i]->GetAttribute();
7911 for (j = 0; j < nv; j++)
7913 out <<
' ' << ind[j]+1;
7918 for (i = 0; i < NumOfBdrElements; i++)
7920 nv = boundary[i]->GetNVertices();
7921 ind = boundary[i]->GetVertices();
7922 out << boundary[i]->GetAttribute();
7923 for (j = 0; j < nv; j++)
7925 out <<
' ' << ind[j]+1;
7927 out <<
" 1.0 1.0 1.0 1.0\n";
7935 void Mesh::Printer(std::ostream &
out, std::string section_delimiter)
const
7942 NURBSext->Print(out);
7953 out << (ncmesh ?
"MFEM mesh v1.1\n" :
7954 section_delimiter.empty() ?
"MFEM mesh v1.0\n" :
7955 "MFEM mesh v1.2\n");
7959 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
7964 "# TETRAHEDRON = 4\n"
7969 out <<
"\ndimension\n" << Dim
7970 <<
"\n\nelements\n" << NumOfElements <<
'\n';
7971 for (i = 0; i < NumOfElements; i++)
7973 PrintElement(elements[i], out);
7976 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
7977 for (i = 0; i < NumOfBdrElements; i++)
7979 PrintElement(boundary[i], out);
7984 out <<
"\nvertex_parents\n";
7985 ncmesh->PrintVertexParents(out);
7987 out <<
"\ncoarse_elements\n";
7988 ncmesh->PrintCoarseElements(out);
7991 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
7994 out << spaceDim <<
'\n';
7995 for (i = 0; i < NumOfVertices; i++)
7997 out << vertices[i](0);
7998 for (j = 1; j < spaceDim; j++)
8000 out <<
' ' << vertices[i](j);
8012 if (!ncmesh && !section_delimiter.empty())
8014 out << section_delimiter << endl;
8023 out <<
"MFEM NURBS mesh v1.0\n";
8027 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8033 out <<
"\ndimension\n" << Dim
8034 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8035 for (i = 0; i < NumOfElements; i++)
8037 PrintElement(elements[i], out);
8040 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
8041 for (i = 0; i < NumOfBdrElements; i++)
8043 PrintElement(boundary[i], out);
8046 out <<
"\nedges\n" << NumOfEdges <<
'\n';
8047 for (i = 0; i < NumOfEdges; i++)
8049 edge_vertex->GetRow(i, vert);
8055 out << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
8057 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8060 void Mesh::PrintVTK(std::ostream &
out)
8063 "# vtk DataFile Version 3.0\n"
8064 "Generated by MFEM\n"
8066 "DATASET UNSTRUCTURED_GRID\n";
8070 out <<
"POINTS " << NumOfVertices <<
" double\n";
8071 for (
int i = 0; i < NumOfVertices; i++)
8073 out << vertices[i](0);
8075 for (j = 1; j < spaceDim; j++)
8077 out <<
' ' << vertices[i](j);
8089 out <<
"POINTS " << Nodes->FESpace()->GetNDofs() <<
" double\n";
8090 for (
int i = 0; i < Nodes->FESpace()->GetNDofs(); i++)
8094 Nodes->FESpace()->DofsToVDofs(vdofs);
8095 out << (*Nodes)(vdofs[0]);
8097 for (j = 1; j < spaceDim; j++)
8099 out <<
' ' << (*Nodes)(vdofs[j]);
8113 for (
int i = 0; i < NumOfElements; i++)
8115 size += elements[i]->GetNVertices() + 1;
8117 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
8118 for (
int i = 0; i < NumOfElements; i++)
8120 const int *v = elements[i]->GetVertices();
8121 const int nv = elements[i]->GetNVertices();
8123 for (
int j = 0; j < nv; j++)
8135 for (
int i = 0; i < NumOfElements; i++)
8137 Nodes->FESpace()->GetElementDofs(i, dofs);
8138 MFEM_ASSERT(Dim != 0 || dofs.
Size() == 1,
8139 "Point meshes should have a single dof per element");
8140 size += dofs.
Size() + 1;
8142 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
8143 const char *fec_name = Nodes->FESpace()->FEColl()->Name();
8145 if (!strcmp(fec_name,
"Linear") ||
8146 !strcmp(fec_name,
"H1_0D_P1") ||
8147 !strcmp(fec_name,
"H1_1D_P1") ||
8148 !strcmp(fec_name,
"H1_2D_P1") ||
8149 !strcmp(fec_name,
"H1_3D_P1"))
8153 else if (!strcmp(fec_name,
"Quadratic") ||
8154 !strcmp(fec_name,
"H1_1D_P2") ||
8155 !strcmp(fec_name,
"H1_2D_P2") ||
8156 !strcmp(fec_name,
"H1_3D_P2"))
8162 mfem::err <<
"Mesh::PrintVTK : can not save '"
8163 << fec_name <<
"' elements!" << endl;
8166 for (
int i = 0; i < NumOfElements; i++)
8168 Nodes->FESpace()->GetElementDofs(i, dofs);
8172 for (
int j = 0; j < dofs.
Size(); j++)
8174 out <<
' ' << dofs[j];
8177 else if (order == 2)
8179 const int *vtk_mfem;
8180 switch (elements[i]->GetGeometryType())
8182 case Geometry::SEGMENT:
8183 case Geometry::TRIANGLE:
8184 case Geometry::SQUARE:
8185 vtk_mfem = vtk_quadratic_hex;
break;
8186 case Geometry::TETRAHEDRON:
8187 vtk_mfem = vtk_quadratic_tet;
break;
8188 case Geometry::PRISM:
8189 vtk_mfem = vtk_quadratic_wedge;
break;
8190 case Geometry::CUBE:
8192 vtk_mfem = vtk_quadratic_hex;
break;
8194 for (
int j = 0; j < dofs.
Size(); j++)
8196 out <<
' ' << dofs[vtk_mfem[j]];
8203 out <<
"CELL_TYPES " << NumOfElements <<
'\n';
8204 for (
int i = 0; i < NumOfElements; i++)
8206 int vtk_cell_type = 5;
8212 case Geometry::POINT: vtk_cell_type = 1;
break;
8213 case Geometry::SEGMENT: vtk_cell_type = 3;
break;
8214 case Geometry::TRIANGLE: vtk_cell_type = 5;
break;
8215 case Geometry::SQUARE: vtk_cell_type = 9;
break;
8216 case Geometry::TETRAHEDRON: vtk_cell_type = 10;
break;
8217 case Geometry::CUBE: vtk_cell_type = 12;
break;
8218 case Geometry::PRISM: vtk_cell_type = 13;
break;
8222 else if (order == 2)
8226 case Geometry::SEGMENT: vtk_cell_type = 21;
break;
8227 case Geometry::TRIANGLE: vtk_cell_type = 22;
break;
8228 case Geometry::SQUARE: vtk_cell_type = 28;
break;
8229 case Geometry::TETRAHEDRON: vtk_cell_type = 24;
break;
8230 case Geometry::CUBE: vtk_cell_type = 29;
break;
8231 case Geometry::PRISM: vtk_cell_type = 32;
break;
8236 out << vtk_cell_type <<
'\n';
8240 out <<
"CELL_DATA " << NumOfElements <<
'\n'
8241 <<
"SCALARS material int\n"
8242 <<
"LOOKUP_TABLE default\n";
8243 for (
int i = 0; i < NumOfElements; i++)
8245 out << elements[i]->GetAttribute() <<
'\n';
8250 void Mesh::PrintVTK(std::ostream &
out,
int ref,
int field_data)
8257 "# vtk DataFile Version 3.0\n"
8258 "Generated by MFEM\n"
8260 "DATASET UNSTRUCTURED_GRID\n";
8265 out <<
"FIELD FieldData 1\n"
8266 <<
"MaterialIds " << 1 <<
" " << attributes.
Size() <<
" int\n";
8267 for (
int i = 0; i < attributes.Size(); i++)
8269 out <<
' ' << attributes[i];
8276 for (
int i = 0; i < GetNE(); i++)
8285 out <<
"POINTS " << np <<
" double\n";
8287 for (
int i = 0; i < GetNE(); i++)
8290 GetElementBaseGeometry(i), ref, 1);
8292 GetElementTransformation(i)->Transform(RefG->
RefPts, pmat);
8294 for (
int j = 0; j < pmat.
Width(); j++)
8296 out << pmat(0, j) <<
' ';
8299 out << pmat(1, j) <<
' ';
8311 out << 0.0 <<
' ' << 0.0;
8318 out <<
"CELLS " << nc <<
' ' << size <<
'\n';
8320 for (
int i = 0; i < GetNE(); i++)
8327 for (
int j = 0; j < RG.
Size(); )
8330 for (
int k = 0; k < nv; k++, j++)
8332 out <<
' ' << np + RG[j];
8338 out <<
"CELL_TYPES " << nc <<
'\n';
8339 for (
int i = 0; i < GetNE(); i++)
8345 int vtk_cell_type = 5;
8349 case Geometry::POINT: vtk_cell_type = 1;
break;
8350 case Geometry::SEGMENT: vtk_cell_type = 3;
break;
8351 case Geometry::TRIANGLE: vtk_cell_type = 5;
break;
8352 case Geometry::SQUARE: vtk_cell_type = 9;
break;
8353 case Geometry::TETRAHEDRON: vtk_cell_type = 10;
break;
8354 case Geometry::CUBE: vtk_cell_type = 12;
break;
8355 case Geometry::PRISM: vtk_cell_type = 13;
break;
8357 MFEM_ABORT(
"Unrecognized VTK element type \"" << geom <<
"\"");
8361 for (
int j = 0; j < RG.
Size(); j += nv)
8363 out << vtk_cell_type <<
'\n';
8367 out <<
"CELL_DATA " << nc <<
'\n'
8368 <<
"SCALARS material int\n"
8369 <<
"LOOKUP_TABLE default\n";
8370 for (
int i = 0; i < GetNE(); i++)
8375 int attr = GetAttribute(i);
8378 out << attr <<
'\n';
8385 srand((
unsigned)time(0));
8386 double a = double(rand()) / (double(RAND_MAX) + 1.);
8387 int el0 = (int)floor(a * GetNE());
8388 GetElementColoring(coloring, el0);
8389 out <<
"SCALARS element_coloring int\n"
8390 <<
"LOOKUP_TABLE default\n";
8391 for (
int i = 0; i < GetNE(); i++)
8398 out << coloring[i] + 1 <<
'\n';
8404 out <<
"POINT_DATA " << np <<
'\n' << flush;
8409 int delete_el_to_el = (el_to_el) ? (0) : (1);
8410 const Table &el_el = ElementToElementTable();
8411 int num_el = GetNE(), stack_p, stack_top_p, max_num_col;
8414 const int *i_el_el = el_el.
GetI();
8415 const int *j_el_el = el_el.
GetJ();
8420 stack_p = stack_top_p = 0;
8421 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
8423 if (colors[el] != -2)
8429 el_stack[stack_top_p++] = el;
8431 for ( ; stack_p < stack_top_p; stack_p++)
8433 int i = el_stack[stack_p];
8434 int num_nb = i_el_el[i+1] - i_el_el[i];
8435 if (max_num_col < num_nb + 1)
8437 max_num_col = num_nb + 1;
8439 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
8442 if (colors[k] == -2)
8445 el_stack[stack_top_p++] = k;
8453 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
8455 int i = el_stack[stack_p], col;
8457 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
8459 col = colors[j_el_el[j]];
8462 col_marker[col] = 1;
8466 for (col = 0; col < max_num_col; col++)
8467 if (col_marker[col] == 0)
8475 if (delete_el_to_el)
8482 void Mesh::PrintWithPartitioning(
int *partitioning, std::ostream &
out,
8483 int elem_attr)
const
8485 if (Dim != 3 && Dim != 2) {
return; }
8487 int i, j, k, l, nv, nbe, *v;
8489 out <<
"MFEM mesh v1.0\n";
8493 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8498 "# TETRAHEDRON = 4\n"
8503 out <<
"\ndimension\n" << Dim
8504 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8505 for (i = 0; i < NumOfElements; i++)
8507 out << int((elem_attr) ? partitioning[i]+1 : elements[i]->GetAttribute())
8508 <<
' ' << elements[i]->GetGeometryType();
8509 nv = elements[i]->GetNVertices();
8510 v = elements[i]->GetVertices();
8511 for (j = 0; j < nv; j++)
8518 for (i = 0; i < faces_info.Size(); i++)
8520 if ((l = faces_info[i].Elem2No) >= 0)
8522 k = partitioning[faces_info[i].Elem1No];
8523 l = partitioning[l];
8527 if (!Nonconforming() || !IsSlaveFace(faces_info[i]))
8538 out <<
"\nboundary\n" << nbe <<
'\n';
8539 for (i = 0; i < faces_info.Size(); i++)
8541 if ((l = faces_info[i].Elem2No) >= 0)
8543 k = partitioning[faces_info[i].Elem1No];
8544 l = partitioning[l];
8547 nv = faces[i]->GetNVertices();
8548 v = faces[i]->GetVertices();
8549 out << k+1 <<
' ' << faces[i]->GetGeometryType();
8550 for (j = 0; j < nv; j++)
8555 if (!Nonconforming() || !IsSlaveFace(faces_info[i]))
8557 out << l+1 <<
' ' << faces[i]->GetGeometryType();
8558 for (j = nv-1; j >= 0; j--)
8568 k = partitioning[faces_info[i].Elem1No];
8569 nv = faces[i]->GetNVertices();
8570 v = faces[i]->GetVertices();
8571 out << k+1 <<
' ' << faces[i]->GetGeometryType();
8572 for (j = 0; j < nv; j++)
8579 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8582 out << spaceDim <<
'\n';
8583 for (i = 0; i < NumOfVertices; i++)
8585 out << vertices[i](0);
8586 for (j = 1; j < spaceDim; j++)
8588 out <<
' ' << vertices[i](j);
8601 void Mesh::PrintElementsWithPartitioning(
int *partitioning,
8605 MFEM_ASSERT(Dim == spaceDim,
"2D Manifolds not supported\n");
8606 if (Dim != 3 && Dim != 2) {
return; }
8613 int *vcount =
new int[NumOfVertices];
8614 for (i = 0; i < NumOfVertices; i++)
8618 for (i = 0; i < NumOfElements; i++)
8620 nv = elements[i]->GetNVertices();
8621 ind = elements[i]->GetVertices();
8622 for (j = 0; j < nv; j++)
8628 int *voff =
new int[NumOfVertices+1];
8630 for (i = 1; i <= NumOfVertices; i++)
8632 voff[i] = vcount[i-1] + voff[i-1];
8635 int **vown =
new int*[NumOfVertices];
8636 for (i = 0; i < NumOfVertices; i++)
8638 vown[i] =
new int[vcount[i]];
8648 Transpose(ElementToEdgeTable(), edge_el);
8651 for (i = 0; i < NumOfElements; i++)
8653 nv = elements[i]->GetNVertices();
8654 ind = elements[i]->GetVertices();
8655 for (j = 0; j < nv; j++)
8658 vown[ind[j]][vcount[ind[j]]] = i;
8662 for (i = 0; i < NumOfVertices; i++)
8664 vcount[i] = voff[i+1] - voff[i];
8668 for (i = 0; i < edge_el.
Size(); i++)
8670 const int *el = edge_el.
GetRow(i);
8673 k = partitioning[el[0]];
8674 l = partitioning[el[1]];
8675 if (interior_faces || k != l)
8687 out <<
"areamesh2\n\n" << nbe <<
'\n';
8689 for (i = 0; i < edge_el.
Size(); i++)
8691 const int *el = edge_el.
GetRow(i);
8694 k = partitioning[el[0]];
8695 l = partitioning[el[1]];
8696 if (interior_faces || k != l)
8699 GetEdgeVertices(i,ev);
8701 for (j = 0; j < 2; j++)
8702 for (s = 0; s < vcount[ev[j]]; s++)
8703 if (vown[ev[j]][s] == el[0])
8705 out <<
' ' << voff[ev[j]]+s+1;
8709 for (j = 1; j >= 0; j--)
8710 for (s = 0; s < vcount[ev[j]]; s++)
8711 if (vown[ev[j]][s] == el[1])
8713 out <<
' ' << voff[ev[j]]+s+1;
8720 k = partitioning[el[0]];
8722 GetEdgeVertices(i,ev);
8724 for (j = 0; j < 2; j++)
8725 for (s = 0; s < vcount[ev[j]]; s++)
8726 if (vown[ev[j]][s] == el[0])
8728 out <<
' ' << voff[ev[j]]+s+1;
8735 out << NumOfElements <<
'\n';
8736 for (i = 0; i < NumOfElements; i++)
8738 nv = elements[i]->GetNVertices();
8739 ind = elements[i]->GetVertices();
8740 out << partitioning[i]+1 <<
' ';
8742 for (j = 0; j < nv; j++)
8744 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
8745 vown[ind[j]][vcount[ind[j]]] = i;
8750 for (i = 0; i < NumOfVertices; i++)
8752 vcount[i] = voff[i+1] - voff[i];
8756 out << voff[NumOfVertices] <<
'\n';
8757 for (i = 0; i < NumOfVertices; i++)
8758 for (k = 0; k < vcount[i]; k++)
8760 for (j = 0; j < Dim; j++)
8762 out << vertices[i](j) <<
' ';
8768 else if (meshgen == 1)
8770 out <<
"NETGEN_Neutral_Format\n";
8772 out << voff[NumOfVertices] <<
'\n';
8773 for (i = 0; i < NumOfVertices; i++)
8774 for (k = 0; k < vcount[i]; k++)
8776 for (j = 0; j < Dim; j++)
8778 out <<
' ' << vertices[i](j);
8784 out << NumOfElements <<
'\n';
8785 for (i = 0; i < NumOfElements; i++)
8787 nv = elements[i]->GetNVertices();
8788 ind = elements[i]->GetVertices();
8789 out << partitioning[i]+1;
8790 for (j = 0; j < nv; j++)
8792 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
8793 vown[ind[j]][vcount[ind[j]]] = i;
8798 for (i = 0; i < NumOfVertices; i++)
8800 vcount[i] = voff[i+1] - voff[i];
8806 for (i = 0; i < NumOfFaces; i++)
8807 if ((l = faces_info[i].Elem2No) >= 0)
8809 k = partitioning[faces_info[i].Elem1No];
8810 l = partitioning[l];
8811 if (interior_faces || k != l)
8822 for (i = 0; i < NumOfFaces; i++)
8823 if ((l = faces_info[i].Elem2No) >= 0)
8825 k = partitioning[faces_info[i].Elem1No];
8826 l = partitioning[l];
8827 if (interior_faces || k != l)
8829 nv = faces[i]->GetNVertices();
8830 ind = faces[i]->GetVertices();
8832 for (j = 0; j < nv; j++)
8833 for (s = 0; s < vcount[ind[j]]; s++)
8834 if (vown[ind[j]][s] == faces_info[i].Elem1No)
8836 out <<
' ' << voff[ind[j]]+s+1;
8840 for (j = nv-1; j >= 0; j--)
8841 for (s = 0; s < vcount[ind[j]]; s++)
8842 if (vown[ind[j]][s] == faces_info[i].Elem2No)
8844 out <<
' ' << voff[ind[j]]+s+1;
8851 k = partitioning[faces_info[i].Elem1No];
8852 nv = faces[i]->GetNVertices();
8853 ind = faces[i]->GetVertices();
8855 for (j = 0; j < nv; j++)
8856 for (s = 0; s < vcount[ind[j]]; s++)
8857 if (vown[ind[j]][s] == faces_info[i].Elem1No)
8859 out <<
' ' << voff[ind[j]]+s+1;
8865 else if (meshgen == 2)
8870 for (i = 0; i < NumOfFaces; i++)
8871 if ((l = faces_info[i].Elem2No) >= 0)
8873 k = partitioning[faces_info[i].Elem1No];
8874 l = partitioning[l];
8875 if (interior_faces || k != l)
8887 <<
"1 " << voff[NumOfVertices] <<
" " << NumOfElements
8888 <<
" 0 0 0 0 0 0 0\n"
8889 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
8890 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
8891 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
8892 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
8894 for (i = 0; i < NumOfVertices; i++)
8895 for (k = 0; k < vcount[i]; k++)
8896 out << voff[i]+k <<
" 0.0 " << vertices[i](0) <<
' '
8897 << vertices[i](1) <<
' ' << vertices[i](2) <<
" 0.0\n";
8899 for (i = 0; i < NumOfElements; i++)
8901 nv = elements[i]->GetNVertices();
8902 ind = elements[i]->GetVertices();
8903 out << i+1 <<
' ' << partitioning[i]+1;
8904 for (j = 0; j < nv; j++)
8906 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
8907 vown[ind[j]][vcount[ind[j]]] = i;
8912 for (i = 0; i < NumOfVertices; i++)
8914 vcount[i] = voff[i+1] - voff[i];
8918 for (i = 0; i < NumOfFaces; i++)
8919 if ((l = faces_info[i].Elem2No) >= 0)
8921 k = partitioning[faces_info[i].Elem1No];
8922 l = partitioning[l];
8923 if (interior_faces || k != l)
8925 nv = faces[i]->GetNVertices();
8926 ind = faces[i]->GetVertices();
8928 for (j = 0; j < nv; j++)
8929 for (s = 0; s < vcount[ind[j]]; s++)
8930 if (vown[ind[j]][s] == faces_info[i].Elem1No)
8932 out <<
' ' << voff[ind[j]]+s+1;
8934 out <<
" 1.0 1.0 1.0 1.0\n";
8936 for (j = nv-1; j >= 0; j--)
8937 for (s = 0; s < vcount[ind[j]]; s++)
8938 if (vown[ind[j]][s] == faces_info[i].Elem2No)
8940 out <<
' ' << voff[ind[j]]+s+1;
8942 out <<
" 1.0 1.0 1.0 1.0\n";
8947 k = partitioning[faces_info[i].Elem1No];
8948 nv = faces[i]->GetNVertices();
8949 ind = faces[i]->GetVertices();
8951 for (j = 0; j < nv; j++)
8952 for (s = 0; s < vcount[ind[j]]; s++)
8953 if (vown[ind[j]][s] == faces_info[i].Elem1No)
8955 out <<
' ' << voff[ind[j]]+s+1;
8957 out <<
" 1.0 1.0 1.0 1.0\n";
8963 for (i = 0; i < NumOfVertices; i++)
8973 void Mesh::PrintSurfaces(
const Table & Aface_face, std::ostream &
out)
const
8980 " NURBS mesh is not supported!");
8984 out <<
"MFEM mesh v1.0\n";
8988 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8993 "# TETRAHEDRON = 4\n"
8998 out <<
"\ndimension\n" << Dim
8999 <<
"\n\nelements\n" << NumOfElements <<
'\n';
9000 for (i = 0; i < NumOfElements; i++)
9002 PrintElement(elements[i], out);
9006 const int *
const i_AF_f = Aface_face.
GetI();
9007 const int *
const j_AF_f = Aface_face.
GetJ();
9009 for (
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
9010 for (
const int * iface = j_AF_f + i_AF_f[iAF];
9011 iface < j_AF_f + i_AF_f[iAF+1];
9014 out << iAF+1 <<
' ';
9015 PrintElementWithoutAttr(faces[*iface],out);
9018 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
9021 out << spaceDim <<
'\n';
9022 for (i = 0; i < NumOfVertices; i++)
9024 out << vertices[i](0);
9025 for (j = 1; j < spaceDim; j++)
9027 out <<
' ' << vertices[i](j);
9040 void Mesh::ScaleSubdomains(
double sf)
9045 int na = attributes.
Size();
9046 double *cg =
new double[na*spaceDim];
9047 int *nbea =
new int[na];
9049 int *vn =
new int[NumOfVertices];
9050 for (i = 0; i < NumOfVertices; i++)
9054 for (i = 0; i < na; i++)
9056 for (j = 0; j < spaceDim; j++)
9058 cg[i*spaceDim+j] = 0.0;
9063 for (i = 0; i < NumOfElements; i++)
9065 GetElementVertices(i, vert);
9066 for (k = 0; k < vert.
Size(); k++)
9072 for (i = 0; i < NumOfElements; i++)
9074 int bea = GetAttribute(i)-1;
9075 GetPointMatrix(i, pointmat);
9076 GetElementVertices(i, vert);
9078 for (k = 0; k < vert.
Size(); k++)
9079 if (vn[vert[k]] == 1)
9082 for (j = 0; j < spaceDim; j++)
9084 cg[bea*spaceDim+j] += pointmat(j,k);
9090 for (i = 0; i < NumOfElements; i++)
9092 int bea = GetAttribute(i)-1;
9093 GetElementVertices (i, vert);
9095 for (k = 0; k < vert.
Size(); k++)
9098 for (j = 0; j < spaceDim; j++)
9099 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
9100 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
9110 void Mesh::ScaleElements(
double sf)
9115 int na = NumOfElements;
9116 double *cg =
new double[na*spaceDim];
9117 int *nbea =
new int[na];
9119 int *vn =
new int[NumOfVertices];
9120 for (i = 0; i < NumOfVertices; i++)
9124 for (i = 0; i < na; i++)
9126 for (j = 0; j < spaceDim; j++)
9128 cg[i*spaceDim+j] = 0.0;
9133 for (i = 0; i < NumOfElements; i++)
9135 GetElementVertices(i, vert);
9136 for (k = 0; k < vert.
Size(); k++)
9142 for (i = 0; i < NumOfElements; i++)
9145 GetPointMatrix(i, pointmat);
9146 GetElementVertices(i, vert);
9148 for (k = 0; k < vert.
Size(); k++)
9149 if (vn[vert[k]] == 1)
9152 for (j = 0; j < spaceDim; j++)
9154 cg[bea*spaceDim+j] += pointmat(j,k);
9160 for (i = 0; i < NumOfElements; i++)
9163 GetElementVertices(i, vert);
9165 for (k = 0; k < vert.
Size(); k++)
9168 for (j = 0; j < spaceDim; j++)
9169 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
9170 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
9185 Vector vold(spaceDim), vnew(NULL, spaceDim);
9186 for (
int i = 0; i < vertices.Size(); i++)
9188 for (
int j = 0; j < spaceDim; j++)
9190 vold(j) = vertices[i](j);
9200 xnew.ProjectCoefficient(f_pert);
9207 MFEM_VERIFY(spaceDim == deformation.
GetVDim(),
9208 "incompatible vector dimensions");
9215 for (
int i = 0; i < NumOfVertices; i++)
9216 for (
int d = 0; d < spaceDim; d++)
9218 vertices[i](d) = xnew(d + spaceDim*i);
9229 void Mesh::RemoveUnusedVertices()
9231 if (NURBSext || ncmesh) {
return; }
9235 for (
int i = 0; i < GetNE(); i++)
9240 for (
int j = 0; j < nv; j++)
9245 for (
int i = 0; i < GetNBE(); i++)
9247 Element *el = GetBdrElement(i);
9250 for (
int j = 0; j < nv; j++)
9256 for (
int i = 0; i < v2v.
Size(); i++)
9260 vertices[num_vert] = vertices[i];
9261 v2v[i] = num_vert++;
9265 if (num_vert == v2v.
Size()) {
return; }
9272 for (
int i = 0; i < GetNE(); i++)
9274 Nodes->FESpace()->GetElementVDofs(i, vdofs);
9279 for (
int i = 0; i < GetNE(); i++)
9281 Nodes->FESpace()->GetElementVDofs(i, vdofs);
9282 Nodes->GetSubVector(vdofs, &nodes_by_element(s));
9286 vertices.SetSize(num_vert);
9287 NumOfVertices = num_vert;
9288 for (
int i = 0; i < GetNE(); i++)
9293 for (
int j = 0; j < nv; j++)
9298 for (
int i = 0; i < GetNBE(); i++)
9300 Element *el = GetBdrElement(i);
9303 for (
int j = 0; j < nv; j++)
9312 el_to_edge =
new Table;
9313 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
9318 GetElementToFaceTable();
9324 Nodes->FESpace()->Update();
9327 for (
int i = 0; i < GetNE(); i++)
9329 Nodes->FESpace()->GetElementVDofs(i, vdofs);
9330 Nodes->SetSubVector(vdofs, &nodes_by_element(s));
9336 void Mesh::RemoveInternalBoundaries()
9338 if (NURBSext || ncmesh) {
return; }
9340 int num_bdr_elem = 0;
9341 int new_bel_to_edge_nnz = 0;
9342 for (
int i = 0; i < GetNBE(); i++)
9344 if (FaceIsInterior(GetBdrElementEdgeIndex(i)))
9346 FreeElement(boundary[i]);
9353 new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
9358 if (num_bdr_elem == GetNBE()) {
return; }
9362 Table *new_bel_to_edge = NULL;
9366 new_be_to_edge.
Reserve(num_bdr_elem);
9370 new_be_to_face.
Reserve(num_bdr_elem);
9371 new_bel_to_edge =
new Table;
9372 new_bel_to_edge->
SetDims(num_bdr_elem, new_bel_to_edge_nnz);
9374 for (
int i = 0; i < GetNBE(); i++)
9376 if (!FaceIsInterior(GetBdrElementEdgeIndex(i)))
9378 new_boundary.
Append(boundary[i]);
9381 new_be_to_edge.
Append(be_to_edge[i]);
9385 int row = new_be_to_face.
Size();
9386 new_be_to_face.
Append(be_to_face[i]);
9387 int *e = bel_to_edge->GetRow(i);
9388 int ne = bel_to_edge->RowSize(i);
9389 int *new_e = new_bel_to_edge->
GetRow(row);
9390 for (
int j = 0; j < ne; j++)
9394 new_bel_to_edge->
GetI()[row+1] = new_bel_to_edge->
GetI()[row] + ne;
9399 NumOfBdrElements = new_boundary.
Size();
9410 bel_to_edge = new_bel_to_edge;
9414 for (
int i = 0; i < attribs.
Size(); i++)
9416 attribs[i] = GetBdrAttribute(i);
9420 bdr_attributes.DeleteAll();
9421 attribs.
Copy(bdr_attributes);
9426 #ifdef MFEM_USE_MEMALLOC
9429 if (E->
GetType() == Element::TETRAHEDRON)
9453 const int npts = point_mat.
Width();
9454 if (!npts) {
return 0; }
9455 MFEM_VERIFY(point_mat.
Height() == spaceDim,
"Invalid points matrix");
9459 if (!GetNE()) {
return 0; }
9461 double *data = point_mat.
GetData();
9468 min_dist = std::numeric_limits<double>::max();
9472 for (
int i = 0; i < GetNE(); i++)
9474 GetElementTransformation(i)->Transform(
9476 for (
int k = 0; k < npts; k++)
9478 double dist = pt.
DistanceTo(data+k*spaceDim);
9479 if (dist < min_dist(k))
9490 for (
int k = 0; k < npts; k++)
9494 int res = inv_tr->
Transform(pt, ips[k]);
9495 if (res == InverseElementTransformation::Inside)
9497 elem_ids[k] = e_idx[k];
9501 if (pts_found != npts)
9504 Table *vtoel = GetVertexToElementTable();
9505 for (
int k = 0; k < npts; k++)
9507 if (elem_ids[k] != -1) {
continue; }
9510 GetElementVertices(e_idx[k], vertices);
9511 for (
int v = 0; v < vertices.
Size(); v++)
9513 int vv = vertices[v];
9515 const int* els = vtoel->
GetRow(vv);
9516 for (
int e = 0; e < ne; e++)
9518 if (els[e] == e_idx[k]) {
continue; }
9520 int res = inv_tr->
Transform(pt, ips[k]);
9521 if (res == InverseElementTransformation::Inside)
9523 elem_ids[k] = els[e];
9533 int le = ncmesh->leaf_elements[e_idx[k]];
9534 ncmesh->FindNeighbors(le,neigh);
9535 for (
int e = 0; e < neigh.
Size(); e++)
9538 if (ncmesh->IsGhost(ncmesh->elements[nn])) {
continue; }
9539 int el = ncmesh->elements[nn].index;
9541 int res = inv_tr->
Transform(pt, ips[k]);
9542 if (res == InverseElementTransformation::Inside)
9554 if (inv_trans == NULL) {
delete inv_tr; }
9556 if (warn && pts_found != npts)
9558 MFEM_WARNING((npts-pts_found) <<
" points were not found");
9569 computed_factors = flags;
9574 const int vdim = fespace->
GetVDim();
9575 const int NE = fespace->
GetNE();
9576 const int ND = fe->
GetDof();
9579 Vector Enodes(vdim*ND*NE);
9582 ElementDofOrdering::NATIVE);
9583 elem_restr->
Mult(*nodes, Enodes);
9585 unsigned eval_flags = 0;
9586 if (flags & GeometricFactors::COORDINATES)
9588 X.SetSize(vdim*NQ*NE);
9589 eval_flags |= QuadratureInterpolator::VALUES;
9591 if (flags & GeometricFactors::JACOBIANS)
9593 J.SetSize(vdim*vdim*NQ*NE);
9594 eval_flags |= QuadratureInterpolator::DERIVATIVES;
9596 if (flags & GeometricFactors::DETERMINANTS)
9598 detJ.SetSize(NQ*NE);
9599 eval_flags |= QuadratureInterpolator::DETERMINANTS;
9605 qi->
Mult(Enodes, eval_flags, X, J, detJ);
9609 NodeExtrudeCoefficient::NodeExtrudeCoefficient(
const int dim,
const int _n,
9623 V(1) = s * ((ip.
y + layer) / n);
9628 V(2) = s * ((ip.
z + layer) / n);
9637 mfem::err <<
"Extrude1D : Not a 1D mesh!" << endl;
9641 int nvy = (closed) ? (ny) : (ny + 1);
9642 int nvt = mesh->
GetNV() * nvy;
9651 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
9656 for (
int i = 0; i < mesh->
GetNV(); i++)
9659 for (
int j = 0; j < nvy; j++)
9661 vc[1] = sy * (double(j) / ny);
9667 for (
int i = 0; i < mesh->
GetNE(); i++)
9672 for (
int j = 0; j < ny; j++)
9675 qv[0] = vert[0] * nvy + j;
9676 qv[1] = vert[1] * nvy + j;
9677 qv[2] = vert[1] * nvy + (j + 1) % nvy;
9678 qv[3] = vert[0] * nvy + (j + 1) % nvy;
9684 for (
int i = 0; i < mesh->
GetNBE(); i++)
9689 for (
int j = 0; j < ny; j++)
9692 sv[0] = vert[0] * nvy + j;
9693 sv[1] = vert[0] * nvy + (j + 1) % nvy;
9697 Swap<int>(sv[0], sv[1]);
9709 for (
int i = 0; i < mesh->
GetNE(); i++)
9715 sv[0] = vert[0] * nvy;
9716 sv[1] = vert[1] * nvy;
9720 sv[0] = vert[1] * nvy + ny;
9721 sv[1] = vert[0] * nvy + ny;
9737 string cname = name;
9738 if (cname ==
"Linear")
9742 else if (cname ==
"Quadratic")
9746 else if (cname ==
"Cubic")
9750 else if (!strncmp(name,
"H1_", 3))
9754 else if (!strncmp(name,
"L2_T", 4))
9758 else if (!strncmp(name,
"L2_", 3))
9765 mfem::err <<
"Extrude1D : The mesh uses unknown FE collection : "
9777 for (
int i = 0; i < mesh->
GetNE(); i++)
9780 for (
int j = ny-1; j >= 0; j--)
9797 mfem::err <<
"Extrude2D : Not a 2D mesh!" << endl;
9802 int nvt = mesh->
GetNV() * nvz;
9807 bool wdgMesh =
false;
9808 bool hexMesh =
false;
9812 for (
int i = 0; i < mesh->
GetNV(); i++)
9816 for (
int j = 0; j < nvz; j++)
9818 vc[2] = sz * (double(j) / nz);
9824 for (
int i = 0; i < mesh->
GetNE(); i++)
9834 for (
int j = 0; j < nz; j++)
9837 pv[0] = vert[0] * nvz + j;
9838 pv[1] = vert[1] * nvz + j;
9839 pv[2] = vert[2] * nvz + j;
9840 pv[3] = vert[0] * nvz + (j + 1) % nvz;
9841 pv[4] = vert[1] * nvz + (j + 1) % nvz;
9842 pv[5] = vert[2] * nvz + (j + 1) % nvz;
9849 for (
int j = 0; j < nz; j++)
9852 hv[0] = vert[0] * nvz + j;
9853 hv[1] = vert[1] * nvz + j;
9854 hv[2] = vert[2] * nvz + j;
9855 hv[3] = vert[3] * nvz + j;
9856 hv[4] = vert[0] * nvz + (j + 1) % nvz;
9857 hv[5] = vert[1] * nvz + (j + 1) % nvz;
9858 hv[6] = vert[2] * nvz + (j + 1) % nvz;
9859 hv[7] = vert[3] * nvz + (j + 1) % nvz;
9861 mesh3d->
AddHex(hv, attr);
9865 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
9866 << geom <<
"\'" << endl;
9872 for (
int i = 0; i < mesh->
GetNBE(); i++)
9877 for (
int j = 0; j < nz; j++)
9880 qv[0] = vert[0] * nvz + j;
9881 qv[1] = vert[1] * nvz + j;
9882 qv[2] = vert[1] * nvz + (j + 1) % nvz;
9883 qv[3] = vert[0] * nvz + (j + 1) % nvz;
9892 for (
int i = 0; i < mesh->
GetNE(); i++)
9903 tv[0] = vert[0] * nvz;
9904 tv[1] = vert[2] * nvz;
9905 tv[2] = vert[1] * nvz;
9909 tv[0] = vert[0] * nvz + nz;
9910 tv[1] = vert[1] * nvz + nz;
9911 tv[2] = vert[2] * nvz + nz;
9919 qv[0] = vert[0] * nvz;
9920 qv[1] = vert[3] * nvz;
9921 qv[2] = vert[2] * nvz;
9922 qv[3] = vert[1] * nvz;
9926 qv[0] = vert[0] * nvz + nz;
9927 qv[1] = vert[1] * nvz + nz;
9928 qv[2] = vert[2] * nvz + nz;
9929 qv[3] = vert[3] * nvz + nz;
9935 mfem::err <<
"Extrude2D : Invalid 2D element type \'"
9936 << geom <<
"\'" << endl;
9942 if ( hexMesh && wdgMesh )
9963 string cname = name;
9964 if (cname ==
"Linear")
9968 else if (cname ==
"Quadratic")
9972 else if (cname ==
"Cubic")
9976 else if (!strncmp(name,
"H1_", 3))
9980 else if (!strncmp(name,
"L2_T", 4))
9984 else if (!strncmp(name,
"L2_", 3))
9991 mfem::err <<
"Extrude3D : The mesh uses unknown FE collection : "
10003 for (
int i = 0; i < mesh->
GetNE(); i++)
10006 for (
int j = nz-1; j >= 0; j--)
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)
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.
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 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 GetMeshComponents(Array< mfem::Vertex > &mvertices, Array< mfem::Element * > &melements, Array< mfem::Element * > &mboundary) const
Return the basic Mesh arrays for the current finest level.
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 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
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)
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.
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
void AddVertex(const double *)
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
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.
Data type hexahedron element.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void GetVertexDofs(int i, Array< int > &dofs) const
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 local
local number within 'element'
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)
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)
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.
const Element * GetElement(int i) const
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
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.
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
FiniteElementSpace * FESpace()
void GetColumn(int c, Vector &col) const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
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 unsigned GetTransform() const
Return current coarse-fine transformation.
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
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
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
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 FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
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...
void Init(int ind1, int ind2, int ind3, int ind4, int attr=1)
Initialize the vertex indices and the attribute of a Tetrahedron.
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.
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
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.
NCMesh * ncmesh
Optional non-conforming mesh extension.
Geometry::Type GetBdrElementBaseGeometry(int i) const
void filter_dos(std::string &line)
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.
int GetId(int p1, int p2)
Get id of item whose parents are p1, p2... Create it if it doesn't exist.
Arbitrary order H1-conforming (continuous) finite elements.
void XYZ_VectorFunction(const Vector &p, Vector &v)
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
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.
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