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);