15 #include "../fem/fem.hpp"
16 #include "../general/sort_pairs.hpp"
33 int geom = GetElementBaseGeometry(i);
42 GetElementJacobian(i, J);
44 return pow(fabs(J.
Det()), 1./Dim);
55 GetElementJacobian(i, J);
57 return sqrt((d_hat * d_hat) / (dir * dir));
80 double h_min, h_max, kappa_min, kappa_max, h,
kappa;
82 cout <<
"Mesh Characteristics:" << flush;
87 if (Vh) Vh->
SetSize(NumOfElements);
88 if (Vk) Vk->
SetSize(NumOfElements);
90 h_min = kappa_min = numeric_limits<double>::infinity();
91 h_max = kappa_max = -h_min;
92 for (i = 0; i < NumOfElements; i++)
94 GetElementJacobian(i, J);
95 h = pow(fabs(J.
Det()), 1.0/
double(dim));
98 if (Vk) (*Vk)(i) = kappa;
100 if (h < h_min) h_min = h;
101 if (h > h_max) h_max = h;
102 if (kappa < kappa_min) kappa_min =
kappa;
103 if (kappa > kappa_max) kappa_max =
kappa;
108 <<
"Number of vertices : " << GetNV() << endl
109 <<
"Number of elements : " << GetNE() << endl
110 <<
"Number of bdr elem : " << GetNBE() << endl
111 <<
"h_min : " << h_min << endl
112 <<
"h_max : " << h_max << endl
116 <<
"Number of vertices : " << GetNV() << endl
117 <<
"Number of edges : " << GetNEdges() << endl
118 <<
"Number of elements : " << GetNE() << endl
119 <<
"Number of bdr elem : " << GetNBE() << endl
120 <<
"Euler Number : " << EulerNumber2D() << endl
121 <<
"h_min : " << h_min << endl
122 <<
"h_max : " << h_max << endl
123 <<
"kappa_min : " << kappa_min << endl
124 <<
"kappa_max : " << kappa_max << endl
128 <<
"Number of vertices : " << GetNV() << endl
129 <<
"Number of edges : " << GetNEdges() << endl
130 <<
"Number of faces : " << GetNFaces() << endl
131 <<
"Number of elements : " << GetNE() << endl
132 <<
"Number of bdr elem : " << GetNBE() << endl
133 <<
"Euler Number : " << EulerNumber() << endl
134 <<
"h_min : " << h_min << endl
135 <<
"h_max : " << h_max << endl
136 <<
"kappa_min : " << kappa_min << endl
137 <<
"kappa_max : " << kappa_max << endl
152 mfem_error(
"Mesh::GetTransformationFEforElement - unknown ElementType");
164 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
170 Nodes->FESpace()->GetElementVDofs(i, vdofs);
172 int n = vdofs.
Size()/spaceDim;
174 for (
int k = 0; k < spaceDim; k++)
175 for (
int j = 0; j < n; j++)
176 pm(k,j) = (*Nodes)(vdofs[n*k+j]);
177 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
189 int nv = elements[i]->GetNVertices();
190 const int *v = elements[i]->GetVertices();
191 int n = vertices.
Size();
193 for (
int k = 0; k < spaceDim; k++)
194 for (
int j = 0; j < nv; j++)
195 pm(k, j) = nodes(k*n+v[j]);
196 ElTr->
SetFE(GetTransformationFEforElementType(GetElementType(i)));
201 Nodes->FESpace()->GetElementVDofs(i, vdofs);
202 int n = vdofs.
Size()/spaceDim;
204 for (
int k = 0; k < spaceDim; k++)
205 for (
int j = 0; j < n; j++)
206 pm(k,j) = nodes(vdofs[n*k+j]);
207 ElTr->
SetFE(Nodes->FESpace()->GetFE(i));
213 GetElementTransformation(i, &Transformation);
215 return &Transformation;
220 GetBdrElementTransformation(i, &FaceTransformation);
221 return &FaceTransformation;
232 GetTransformationFEforElementType(GetBdrElementType(i)));
238 Nodes->FESpace()->GetBdrElementVDofs(i, vdofs);
239 int n = vdofs.
Size()/spaceDim;
241 for (
int k = 0; k < spaceDim; k++)
242 for (
int j = 0; j < n; j++)
243 pm(k,j) = (*Nodes)(vdofs[n*k+j]);
244 ElTr->
SetFE(Nodes->FESpace()->GetBE(i));
250 FTr->
Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
255 const int *v = (Dim == 1) ? &FaceNo : faces[FaceNo]->GetVertices();
256 const int nv = (Dim == 1) ? 1 : faces[FaceNo]->GetNVertices();
258 for (
int i = 0; i < spaceDim; i++)
259 for (
int j = 0; j < nv; j++)
260 pm(i, j) = vertices[v[j]](i);
261 FTr->
SetFE(GetTransformationFEforElementType(
266 const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
270 Nodes->FESpace()->GetFaceVDofs(FaceNo, vdofs);
271 int n = vdofs.
Size()/spaceDim;
273 for (
int i = 0; i < spaceDim; i++)
274 for (
int j = 0; j < n; j++)
275 pm(i, j) = (*Nodes)(vdofs[n*i+j]);
282 FaceInfo &face_info = faces_info[FaceNo];
284 face_el = Nodes->FESpace()->GetTraceElement(face_info.
Elem1No,
290 GetLocalPtToSegTransformation(FaceElemTr.Loc1.Transf,
295 GetLocalSegToTriTransformation(FaceElemTr.Loc1.Transf,
298 GetLocalSegToQuadTransformation(FaceElemTr.Loc1.Transf,
303 GetLocalTriToTetTransformation(FaceElemTr.Loc1.Transf,
308 GetLocalQuadToHexTransformation(FaceElemTr.Loc1.Transf,
314 FaceElemTr.Loc1.Transform(face_el->
GetNodes(), eir);
316 Nodes->GetVectorValues(Transformation, eir, pm);
325 GetFaceTransformation(FaceNo, &FaceTransformation);
326 return &FaceTransformation;
333 GetFaceTransformation(EdgeNo, EdTr);
337 mfem_error(
"Mesh::GetEdgeTransformation not defined in 1D \n");
345 GetEdgeVertices(EdgeNo, v);
348 for (
int i = 0; i < spaceDim; i++)
349 for (
int j = 0; j < nv; j++)
350 pm(i, j) = vertices[v[j]](i);
357 Nodes->FESpace()->GetEdgeVDofs(EdgeNo, vdofs);
358 int n = vdofs.
Size()/spaceDim;
360 for (
int i = 0; i < spaceDim; i++)
361 for (
int j = 0; j < n; j++)
362 pm(i, j) = (*Nodes)(vdofs[n*i+j]);
369 GetEdgeTransformation(EdgeNo, &EdgeTransformation);
370 return &EdgeTransformation;
392 static const int tri_faces[3][2] = {{0, 1}, {1, 2}, {2, 0}};
393 static const int seg_inv_orient[2][2] = {{0, 1}, {1, 0}};
400 tv = tri_faces[i/64];
401 so = seg_inv_orient[i%64];
404 for (j = 0; j < 2; j++)
406 locpm(0, so[j]) = TriVert->
IntPoint(tv[j]).
x;
407 locpm(1, so[j]) = TriVert->
IntPoint(tv[j]).
y;
415 static const int quad_faces[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
416 static const int seg_inv_orient[2][2] = {{0, 1}, {1, 0}};
423 qv = quad_faces[i/64];
424 so = seg_inv_orient[i%64];
427 for (j = 0; j < 2; j++)
429 locpm(0, so[j]) = QuadVert->
IntPoint(qv[j]).
x;
430 locpm(1, so[j]) = QuadVert->
IntPoint(qv[j]).
y;
435 {{1, 2, 3}, {0, 3, 2},
436 {0, 1, 3}, {0, 2, 1}};
440 {{3, 2, 1, 0}, {0, 1, 5, 4},
441 {1, 2, 6, 5}, {2, 3, 7, 6},
442 {3, 0, 4, 7}, {4, 5, 6, 7}};
445 {{0, 1, 2}, {1, 0, 2},
446 {2, 0, 1}, {2, 1, 0},
447 {1, 2, 0}, {0, 2, 1}};
450 {{0, 1, 2, 3}, {0, 3, 2, 1},
451 {1, 2, 3, 0}, {1, 0, 3, 2},
452 {2, 3, 0, 1}, {2, 1, 0, 3},
453 {3, 0, 1, 2}, {3, 2, 1, 0}};
462 const int *tv = tet_faces[i/64];
465 const int *to = tri_orientations[i%64];
469 for (
int j = 0; j < 3; j++)
472 locpm(0, j) = vert.
x;
473 locpm(1, j) = vert.
y;
474 locpm(2, j) = vert.
z;
485 const int *hv = hex_faces[i/64];
487 const int *qo = quad_orientations[i%64];
490 for (
int j = 0; j < 4; j++)
493 locpm(0, j) = vert.
x;
494 locpm(1, j) = vert.
y;
495 locpm(2, j) = vert.
z;
502 FaceInfo &face_info = faces_info[FaceNo];
508 GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
509 FaceElemTr.Elem1 = &Transformation;
512 FaceElemTr.Elem1 = NULL;
517 if ( (FaceElemTr.Elem2No = face_info.
Elem2No) >= 0 && (mask & 2))
520 if (NURBSext && (mask & 1))
521 mfem_error(
"Mesh::GetFaceElementTransformations :"
522 " NURBS mesh is not supported!");
524 GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
525 FaceElemTr.Elem2 = &Transformation2;
528 FaceElemTr.Elem2 = NULL;
533 FaceElemTr.FaceGeom = faces[FaceNo]->GetGeometryType();
537 FaceElemTr.Face = GetFaceTransformation(FaceNo);
539 FaceElemTr.Face = NULL;
542 int face_type = (Dim == 1) ?
Element::POINT : faces[FaceNo]->GetType();
547 GetLocalPtToSegTransformation(FaceElemTr.Loc1.Transf,
550 if (FaceElemTr.Elem2No >= 0 && (mask & 8))
551 GetLocalPtToSegTransformation(FaceElemTr.Loc2.Transf,
558 GetLocalSegToTriTransformation(FaceElemTr.Loc1.Transf,
561 GetLocalSegToQuadTransformation(FaceElemTr.Loc1.Transf,
565 if (FaceElemTr.Elem2No >= 0 && (mask & 8))
568 GetLocalSegToTriTransformation(FaceElemTr.Loc2.Transf,
571 GetLocalSegToQuadTransformation(FaceElemTr.Loc2.Transf,
578 GetLocalTriToTetTransformation(FaceElemTr.Loc1.Transf,
580 if (FaceElemTr.Elem2No >= 0 && (mask & 8))
581 GetLocalTriToTetTransformation(FaceElemTr.Loc2.Transf,
587 GetLocalQuadToHexTransformation(FaceElemTr.Loc1.Transf,
589 if (FaceElemTr.Elem2No >= 0 && (mask & 8))
590 GetLocalQuadToHexTransformation(FaceElemTr.Loc2.Transf,
603 fn = be_to_face[BdrElemNo];
605 fn = be_to_edge[BdrElemNo];
607 fn = boundary[BdrElemNo]->GetVertices()[0];
608 if (FaceIsTrueInterior(fn))
610 tr = GetFaceElementTransformations(fn);
617 *Elem1 = faces_info[Face].Elem1No;
618 *Elem2 = faces_info[Face].Elem2No;
623 *Inf1 = faces_info[Face].Elem1Inf;
624 *Inf2 = faces_info[Face].Elem2Inf;
629 NumOfVertices = NumOfElements = NumOfBdrElements = NumOfEdges = -1;
630 WantTwoLevelState = 0;
636 nc_coarse_level = NULL;
641 el_to_edge = el_to_face = el_to_el =
642 bel_to_edge = face_edge = edge_vertex = NULL;
647 if (el_to_edge != NULL)
650 if (el_to_face != NULL)
653 if (el_to_el != NULL)
656 if (Dim == 3 && bel_to_edge != NULL)
659 if (face_edge != NULL)
662 if (edge_vertex != NULL)
667 delete nc_coarse_level;
668 nc_coarse_level = NULL;
677 el_to_el = face_edge = edge_vertex = NULL;
679 delete nc_coarse_level;
680 nc_coarse_level = NULL;
689 for (i = 0; i < attribs.
Size(); i++)
690 attribs[i] = GetBdrAttribute(i);
693 if (attribs.
Size() > 0)
697 for (i = 1; i < attribs.
Size(); i++)
698 if (attribs[i] != attribs[i-1])
701 bdr_attributes.SetSize(nattr);
704 bdr_attributes[0] = attribs[0];
705 for (i = j = 1; i < attribs.
Size(); i++)
706 if (attribs[i] != attribs[i-1])
707 bdr_attributes[j++] = attribs[i];
709 cout <<
"Mesh::SetAttributes(): "
710 "Non-positive attributes on the boundary!"
716 for (i = 0; i < attribs.
Size(); i++)
717 attribs[i] = GetAttribute(i);
720 if (attribs.
Size() > 0)
724 for (i = 1; i < attribs.
Size(); i++)
725 if (attribs[i] != attribs[i-1])
728 attributes.SetSize(nattr);
731 attributes[0] = attribs[0];
732 for (i = j = 1; i < attribs.
Size(); i++)
733 if (attribs[i] != attribs[i-1])
734 attributes[j++] = attribs[i];
736 cout <<
"Mesh::SetAttributes(): "
737 "Non-positive attributes in the domain!"
745 spaceDim = _spaceDim;
751 vertices.SetSize(NVert);
754 elements.SetSize(NElem);
756 NumOfBdrElements = 0;
757 boundary.SetSize(NBdrElem);
762 double *y = vertices[NumOfVertices]();
764 for (
int i = 0; i < spaceDim; i++)
771 elements[NumOfElements++] =
new Triangle(vi, attr);
776 elements[NumOfElements++] =
new Triangle(vi, attr);
786 #ifdef MFEM_USE_MEMALLOC
788 tet = TetMemory.Alloc();
791 elements[NumOfElements++] = tet;
793 elements[NumOfElements++] =
new Tetrahedron(vi, attr);
799 elements[NumOfElements++] =
new Hexahedron(vi, attr);
804 static const int hex_to_tet[6][4] =
805 { { 0, 1, 2, 6 }, { 0, 5, 1, 6 }, { 0, 4, 5, 6 },
806 { 0, 2, 3, 6 }, { 0, 3, 7, 6 }, { 0, 7, 4, 6 } };
809 for (
int i = 0; i < 6; i++)
811 for (
int j = 0; j < 4; j++)
812 ti[j] = vi[hex_to_tet[i][j]];
819 boundary[NumOfBdrElements++] =
new Segment(vi, attr);
824 boundary[NumOfBdrElements++] =
new Triangle(vi, attr);
834 static const int quad_to_tri[2][3] = { { 0, 1, 2 }, { 0, 2, 3 } };
837 for (
int i = 0; i < 2; i++)
839 for (
int j = 0; j < 3; j++)
840 ti[j] = vi[quad_to_tri[i][j]];
841 AddBdrTriangle(ti, attr);
848 Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
852 for (i = 0; i < boundary.Size(); i++)
853 FreeElement(boundary[i]);
862 NumOfBdrElements = 0;
863 for (i = 0; i < faces_info.Size(); i++)
864 if (faces_info[i].Elem2No < 0)
867 boundary.
SetSize(NumOfBdrElements);
868 be2face.
SetSize(NumOfBdrElements);
869 for (j = i = 0; i < faces_info.Size(); i++)
870 if (faces_info[i].Elem2No < 0)
872 boundary[j] = faces[i]->Duplicate(
this);
884 static int edge_compare(
const void *ii,
const void *jj)
886 edge_length *i = (edge_length *)ii, *j = (edge_length *)jj;
887 if (i->length > j->length)
return (1);
888 if (i->length < j->length)
return (-1);
894 CheckElementOrientation(fix_orientation);
897 MarkTriMeshForRefinement();
901 el_to_edge =
new Table;
902 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
904 CheckBdrElementOrientation();
917 bool fix_orientation)
920 CheckElementOrientation(fix_orientation);
924 el_to_edge =
new Table;
925 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
927 CheckBdrElementOrientation();
944 MarkTriMeshForRefinement();
946 MarkTetMeshForRefinement();
955 for (
int i = 0; i < NumOfElements; i++)
958 GetPointMatrix(i, pmat);
959 elements[i]->MarkEdge(pmat);
966 edge_length *length =
new edge_length[NumOfEdges];
967 for (
int i = 0; i < NumOfVertices; i++)
972 length[j].length = GetLength(i, it.Column());
978 qsort(length, NumOfEdges,
sizeof(edge_length), edge_compare);
981 for (
int i = 0; i < NumOfEdges; i++)
982 order[length[i].edge] = i;
992 GetVertexToVertexTable(v_to_v);
995 GetEdgeOrdering(v_to_v, order);
997 for (
int i = 0; i < NumOfElements; i++)
999 elements[i]->MarkEdge(v_to_v, order);
1001 for (
int i = 0; i < NumOfBdrElements; i++)
1003 boundary[i]->MarkEdge(v_to_v, order);
1008 if (*old_v_to_v && *old_elem_vert)
1014 if (*old_v_to_v == NULL)
1017 if (num_edge_dofs > 0)
1019 *old_v_to_v =
new DSTable(NumOfVertices);
1020 GetVertexToVertexTable(*(*old_v_to_v));
1023 if (*old_elem_vert == NULL)
1026 int num_elem_dofs = fec->
DofForGeometry(GetElementBaseGeometry(0));
1027 if (num_elem_dofs > 1)
1029 *old_elem_vert =
new Table;
1030 (*old_elem_vert)->
MakeI(GetNE());
1031 for (
int i = 0; i < GetNE(); i++)
1033 (*old_elem_vert)->AddColumnsInRow(i, elements[i]->GetNVertices());
1035 (*old_elem_vert)->MakeJ();
1036 for (
int i = 0; i < GetNE(); i++)
1038 (*old_elem_vert)->AddConnections(i, elements[i]->GetVertices(),
1039 elements[i]->GetNVertices());
1041 (*old_elem_vert)->ShiftUpI();
1055 int num_elem_dofs = fec->
DofForGeometry(GetElementBaseGeometry(0));
1071 if (num_edge_dofs > 0)
1073 DSTable new_v_to_v(NumOfVertices);
1074 GetVertexToVertexTable(new_v_to_v);
1076 for (
int i = 0; i < NumOfVertices; i++)
1080 int old_i = (*old_v_to_v)(i, it.Column());
1081 int new_i = it.Index();
1086 old_dofs.
SetSize(num_edge_dofs);
1087 new_dofs.
SetSize(num_edge_dofs);
1088 for (
int j = 0; j < num_edge_dofs; j++)
1090 old_dofs[j] = offset + old_i * num_edge_dofs + j;
1091 new_dofs[j] = offset + new_i * num_edge_dofs + j;
1095 for (
int j = 0; j < old_dofs.
Size(); j++)
1096 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1099 offset += NumOfEdges * num_edge_dofs;
1102 cout <<
"Mesh::DoNodeReorder : redges = " << redges << endl;
1107 if (num_face_dofs > 0)
1110 Table old_face_vertex;
1111 old_face_vertex.
MakeI(NumOfFaces);
1112 for (
int i = 0; i < NumOfFaces; i++)
1114 old_face_vertex.
MakeJ();
1115 for (
int i = 0; i < NumOfFaces; i++)
1117 faces[i]->GetNVertices());
1121 STable3D *faces_tbl = GetElementToFaceTable(1);
1125 for (
int i = 0; i < NumOfFaces; i++)
1127 int *old_v = old_face_vertex.
GetRow(i), *new_v;
1128 int new_i, new_or, *dof_ord;
1129 switch (old_face_vertex.
RowSize(i))
1132 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2]);
1133 new_v = faces[new_i]->GetVertices();
1134 new_or = GetTriOrientation(old_v, new_v);
1139 new_i = (*faces_tbl)(old_v[0], old_v[1], old_v[2], old_v[3]);
1140 new_v = faces[new_i]->GetVertices();
1141 new_or = GetQuadOrientation(old_v, new_v);
1146 old_dofs.
SetSize(num_face_dofs);
1147 new_dofs.
SetSize(num_face_dofs);
1148 for (
int j = 0; j < num_face_dofs; j++)
1150 old_dofs[j] = offset + i * num_face_dofs + j;
1151 new_dofs[j] = offset + new_i * num_face_dofs + dof_ord[j];
1157 for (
int j = 0; j < old_dofs.
Size(); j++)
1158 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1161 offset += NumOfFaces * num_face_dofs;
1167 if (num_elem_dofs > 1)
1179 for (
int i = 0; i < GetNE(); i++)
1181 int *old_v = old_elem_vert->
GetRow(i);
1182 int *new_v = elements[i]->GetVertices();
1183 int new_or, *dof_ord;
1184 int geom = elements[i]->GetGeometryType();
1188 new_or = (old_v[0] == new_v[0]) ? +1 : -1;
1191 new_or = GetTriOrientation(old_v, new_v);
1194 new_or = GetQuadOrientation(old_v, new_v);
1199 <<
" elements (" << fec->
Name()
1200 <<
" FE collection) are not supported yet!" << endl;
1205 if (dof_ord == NULL)
1207 cerr <<
"Mesh::DoNodeReorder : FE collection '" << fec->
Name()
1209 <<
" elements!" << endl;
1212 old_dofs.
SetSize(num_elem_dofs);
1213 new_dofs.
SetSize(num_elem_dofs);
1214 for (
int j = 0; j < num_elem_dofs; j++)
1217 old_dofs[j] = offset + dof_ord[j];
1218 new_dofs[j] = offset + j;
1222 for (
int j = 0; j < old_dofs.
Size(); j++)
1223 (*Nodes)(new_dofs[j]) = onodes(old_dofs[j]);
1225 offset += num_elem_dofs;
1232 if (num_face_dofs == 0)
1236 GetElementToFaceTable();
1239 CheckBdrElementOrientation();
1244 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1249 CheckBdrElementOrientation();
1256 CheckElementOrientation(fix_orientation);
1258 if (NumOfBdrElements == 0)
1260 GetElementToFaceTable();
1262 GenerateBoundaryElements();
1267 MarkTetMeshForRefinement();
1270 GetElementToFaceTable();
1273 CheckBdrElementOrientation();
1275 if (generate_edges == 1)
1277 el_to_edge =
new Table;
1278 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1294 CheckElementOrientation(fix_orientation);
1296 GetElementToFaceTable();
1299 if (NumOfBdrElements == 0)
1300 GenerateBoundaryElements();
1302 CheckBdrElementOrientation();
1306 el_to_edge =
new Table;
1307 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1318 int generate_edges,
double sx,
double sy,
double sz)
1322 int NVert, NElem, NBdrElem;
1324 NVert = (nx+1) * (ny+1) * (nz+1);
1325 NElem = nx * ny * nz;
1326 NBdrElem = 2*(nx*ny+nx*nz+ny*nz);
1333 InitMesh(3, 3, NVert, NElem, NBdrElem);
1339 for (z = 0; z <= nz; z++)
1341 coord[2] = ((double) z / nz) * sz;
1342 for (y = 0; y <= ny; y++)
1344 coord[1] = ((double) y / ny) * sy;
1345 for (x = 0; x <= nx; x++)
1347 coord[0] = ((double) x / nx) * sx;
1353 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
1356 for (z = 0; z < nz; z++)
1358 for (y = 0; y < ny; y++)
1360 for (x = 0; x < nx; x++)
1362 ind[0] = VTX(x , y , z );
1363 ind[1] = VTX(x+1, y , z );
1364 ind[2] = VTX(x+1, y+1, z );
1365 ind[3] = VTX(x , y+1, z );
1366 ind[4] = VTX(x , y , z+1);
1367 ind[5] = VTX(x+1, y , z+1);
1368 ind[6] = VTX(x+1, y+1, z+1);
1369 ind[7] = VTX(x , y+1, z+1);
1371 AddHexAsTets(ind, 1);
1380 for (y = 0; y < ny; y++)
1381 for (x = 0; x < nx; x++)
1383 ind[0] = VTX(x , y , 0);
1384 ind[1] = VTX(x , y+1, 0);
1385 ind[2] = VTX(x+1, y+1, 0);
1386 ind[3] = VTX(x+1, y , 0);
1388 AddBdrQuadAsTriangles(ind, 1);
1393 for (y = 0; y < ny; y++)
1394 for (x = 0; x < nx; x++)
1396 ind[0] = VTX(x , y , nz);
1397 ind[1] = VTX(x+1, y , nz);
1398 ind[2] = VTX(x+1, y+1, nz);
1399 ind[3] = VTX(x , y+1, nz);
1401 AddBdrQuadAsTriangles(ind, 6);
1406 for (z = 0; z < nz; z++)
1407 for (y = 0; y < ny; y++)
1409 ind[0] = VTX(0 , y , z );
1410 ind[1] = VTX(0 , y , z+1);
1411 ind[2] = VTX(0 , y+1, z+1);
1412 ind[3] = VTX(0 , y+1, z );
1414 AddBdrQuadAsTriangles(ind, 5);
1419 for (z = 0; z < nz; z++)
1420 for (y = 0; y < ny; y++)
1422 ind[0] = VTX(nx, y , z );
1423 ind[1] = VTX(nx, y+1, z );
1424 ind[2] = VTX(nx, y+1, z+1);
1425 ind[3] = VTX(nx, y , z+1);
1427 AddBdrQuadAsTriangles(ind, 3);
1432 for (x = 0; x < nx; x++)
1433 for (z = 0; z < nz; z++)
1435 ind[0] = VTX(x , 0, z );
1436 ind[1] = VTX(x+1, 0, z );
1437 ind[2] = VTX(x+1, 0, z+1);
1438 ind[3] = VTX(x , 0, z+1);
1440 AddBdrQuadAsTriangles(ind, 2);
1445 for (x = 0; x < nx; x++)
1446 for (z = 0; z < nz; z++)
1448 ind[0] = VTX(x , ny, z );
1449 ind[1] = VTX(x , ny, z+1);
1450 ind[2] = VTX(x+1, ny, z+1);
1451 ind[3] = VTX(x+1, ny, z );
1453 AddBdrQuadAsTriangles(ind, 4);
1459 ofstream test_stream(
"debug.mesh");
1461 test_stream.close();
1465 bool fix_orientation =
true;
1468 FinalizeTetMesh(generate_edges, refine, fix_orientation);
1470 FinalizeHexMesh(generate_edges, refine, fix_orientation);
1474 double sx,
double sy)
1487 NumOfVertices = (nx+1) * (ny+1);
1488 NumOfElements = nx * ny;
1489 NumOfBdrElements = 2 * nx + 2 * ny;
1490 vertices.SetSize(NumOfVertices);
1491 elements.SetSize(NumOfElements);
1492 boundary.SetSize(NumOfBdrElements);
1498 for (j = 0; j < ny+1; j++)
1500 cy = ((double) j / ny) * sy;
1501 for (i = 0; i < nx+1; i++)
1503 cx = ((double) i / nx) * sx;
1504 vertices[k](0) = cx;
1505 vertices[k](1) = cy;
1512 for (j = 0; j < ny; j++)
1514 for (i = 0; i < nx; i++)
1516 ind[0] = i + j*(nx+1);
1517 ind[1] = i + 1 +j*(nx+1);
1518 ind[2] = i + 1 + (j+1)*(nx+1);
1519 ind[3] = i + (j+1)*(nx+1);
1527 for (i = 0; i < nx; i++)
1529 boundary[i] =
new Segment(i, i+1, 1);
1530 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
1533 for (j = 0; j < ny; j++)
1535 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
1536 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
1544 NumOfVertices = (nx+1) * (ny+1);
1545 NumOfElements = 2 * nx * ny;
1546 NumOfBdrElements = 2 * nx + 2 * ny;
1547 vertices.SetSize(NumOfVertices);
1548 elements.SetSize(NumOfElements);
1549 boundary.SetSize(NumOfBdrElements);
1555 for (j = 0; j < ny+1; j++)
1557 cy = ((double) j / ny) * sy;
1558 for (i = 0; i < nx+1; i++)
1560 cx = ((double) i / nx) * sx;
1561 vertices[k](0) = cx;
1562 vertices[k](1) = cy;
1569 for (j = 0; j < ny; j++)
1571 for (i = 0; i < nx; i++)
1573 ind[0] = i + j*(nx+1);
1574 ind[1] = i + 1 + (j+1)*(nx+1);
1575 ind[2] = i + (j+1)*(nx+1);
1578 ind[1] = i + 1 + j*(nx+1);
1579 ind[2] = i + 1 + (j+1)*(nx+1);
1587 for (i = 0; i < nx; i++)
1589 boundary[i] =
new Segment(i, i+1, 1);
1590 boundary[nx+i] =
new Segment(m+i+1, m+i, 3);
1593 for (j = 0; j < ny; j++)
1595 boundary[2*nx+j] =
new Segment((j+1)*m, j*m, 4);
1596 boundary[2*nx+ny+j] =
new Segment(j*m+nx, (j+1)*m+nx, 2);
1599 MarkTriMeshForRefinement();
1602 CheckElementOrientation();
1604 if (generate_edges == 1)
1606 el_to_edge =
new Table;
1607 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
1609 CheckBdrElementOrientation();
1616 attributes.Append(1);
1617 bdr_attributes.Append(1); bdr_attributes.Append(2);
1618 bdr_attributes.Append(3); bdr_attributes.Append(4);
1633 NumOfVertices = n + 1;
1635 NumOfBdrElements = 2;
1636 vertices.SetSize(NumOfVertices);
1637 elements.SetSize(NumOfElements);
1638 boundary.SetSize(NumOfBdrElements);
1641 for (j = 0; j < n+1; j++)
1642 vertices[j](0) = ((double) j / n) * sx;
1645 for (j = 0; j < n; j++)
1646 elements[j] =
new Segment(j, j+1, 1);
1650 boundary[0] =
new Point(ind, 1);
1652 boundary[1] =
new Point(ind, 2);
1659 attributes.Append(1);
1660 bdr_attributes.Append(1); bdr_attributes.Append(2);
1663 Mesh::Mesh(std::istream &input,
int generate_edges,
int refine,
bool fix_orientation)
1667 Load(input, generate_edges, refine, fix_orientation);
1680 #ifdef MFEM_USE_MEMALLOC
1681 return TetMemory.Alloc();
1696 el = NewElement(geom);
1699 for (
int i = 0; i < nv; i++)
1710 for (
int j = 0; j < nv; j++)
1721 el = ReadElementWithoutAttr(input);
1730 PrintElementWithoutAttr(el, out);
1736 for (
int i = 0; i < NumOfElements; i++)
1738 switch (elements[i]->GetType())
1743 meshgen |= 1;
break;
1753 static const int vtk_quadratic_tet[10] =
1754 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
1757 static const int vtk_quadratic_hex[27] =
1758 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
1759 24, 22, 21, 23, 20, 25, 26 };
1766 if (is.peek() != comment_char)
1768 is.ignore(numeric_limits<streamsize>::max(),
'\n');
1772 void Mesh::Load(std::istream &input,
int generate_edges,
int refine,
1773 bool fix_orientation)
1775 int i, j, ints[32], n, attr, curved = 0, read_gf = 1;
1776 const int buflen = 1024;
1780 MFEM_ABORT(
"Input stream is not open");
1782 if (NumOfVertices != -1)
1785 for (i = 0; i < NumOfElements; i++)
1787 FreeElement(elements[i]);
1788 elements.DeleteAll();
1791 vertices.DeleteAll();
1794 for (i = 0; i < NumOfBdrElements; i++)
1796 FreeElement(boundary[i]);
1797 boundary.DeleteAll();
1800 for (i = 0; i < faces.Size(); i++)
1801 FreeElement(faces[i]);
1804 faces_info.DeleteAll();
1808 be_to_edge.DeleteAll();
1809 be_to_face.DeleteAll();
1814 if (own_nodes)
delete Nodes;
1820 getline(input, mesh_type);
1822 if (mesh_type ==
"MFEM mesh v1.0")
1836 input >> NumOfElements;
1837 elements.SetSize(NumOfElements);
1838 for (j = 0; j < NumOfElements; j++)
1839 elements[j] = ReadElement(input);
1844 input >> NumOfBdrElements;
1845 boundary.SetSize(NumOfBdrElements);
1846 for (j = 0; j < NumOfBdrElements; j++)
1847 boundary[j] = ReadElement(input);
1852 input >> NumOfVertices;
1853 vertices.SetSize(NumOfVertices);
1855 input >> ws >> ident;
1856 if (ident !=
"nodes")
1859 spaceDim = atoi(ident.c_str());
1860 for (j = 0; j < NumOfVertices; j++)
1861 for (i = 0; i < spaceDim; i++)
1862 input >> vertices[j](i);
1871 else if (mesh_type ==
"linemesh")
1877 input >> NumOfVertices;
1878 vertices.SetSize(NumOfVertices);
1880 for (j = 0; j < NumOfVertices; j++)
1881 input >> vertices[j](0);
1883 input >> NumOfElements;
1884 elements.SetSize(NumOfElements);
1886 for (j = 0; j < NumOfElements; j++)
1888 input >> a >> p1 >> p2;
1889 elements[j] =
new Segment(p1-1, p2-1, a);
1893 input >> NumOfBdrElements;
1894 boundary.SetSize(NumOfBdrElements);
1895 for (j = 0; j < NumOfBdrElements; j++)
1897 input >> a >> ind[0];
1899 boundary[j] =
new Point(ind,a);
1902 else if (mesh_type ==
"areamesh2" || mesh_type ==
"curved_areamesh2")
1907 if (mesh_type ==
"curved_areamesh2")
1911 input >> NumOfBdrElements;
1912 boundary.SetSize(NumOfBdrElements);
1913 for (i = 0; i < NumOfBdrElements; i++)
1916 >> ints[0] >> ints[1];
1917 ints[0]--; ints[1]--;
1918 boundary[i] =
new Segment(ints, attr);
1922 input >> NumOfElements;
1923 elements.SetSize(NumOfElements);
1924 for (i = 0; i < NumOfElements; i++)
1927 for (j = 0; j < n; j++)
1935 elements[i] =
new Segment(ints, attr);
1938 elements[i] =
new Triangle(ints, attr);
1949 input >> NumOfVertices;
1950 vertices.SetSize(NumOfVertices);
1951 for (i = 0; i < NumOfVertices; i++)
1952 for (j = 0; j < Dim; j++)
1953 input >> vertices[i](j);
1957 input >> NumOfVertices;
1958 vertices.SetSize(NumOfVertices);
1962 else if (mesh_type ==
"NETGEN" || mesh_type ==
"NETGEN_Neutral_Format")
1968 input >> NumOfVertices;
1970 vertices.SetSize(NumOfVertices);
1971 for (i = 0; i < NumOfVertices; i++)
1972 for (j = 0; j < Dim; j++)
1973 input >> vertices[i](j);
1976 input >> NumOfElements;
1977 elements.SetSize(NumOfElements);
1978 for (i = 0; i < NumOfElements; i++)
1981 for (j = 0; j < 4; j++)
1986 #ifdef MFEM_USE_MEMALLOC
1988 tet = TetMemory.Alloc();
1998 input >> NumOfBdrElements;
1999 boundary.SetSize(NumOfBdrElements);
2000 for (i = 0; i < NumOfBdrElements; i++)
2003 for (j = 0; j < 3; j++)
2008 boundary[i] =
new Triangle(ints, attr);
2011 else if (mesh_type ==
"TrueGrid")
2023 input >> vari >> NumOfVertices >> vari >> vari >> NumOfElements;
2024 input.getline(buf, buflen);
2025 input.getline(buf, buflen);
2027 input.getline(buf, buflen);
2028 input.getline(buf, buflen);
2029 input.getline(buf, buflen);
2032 vertices.SetSize(NumOfVertices);
2033 for (i = 0; i < NumOfVertices; i++)
2035 input >> vari >> varf >> vertices[i](0) >> vertices[i](1);
2036 input.getline(buf, buflen);
2040 elements.SetSize(NumOfElements);
2041 for (i = 0; i < NumOfElements; i++)
2043 input >> vari >> attr;
2044 for (j = 0; j < 4; j++)
2049 input.getline(buf, buflen);
2050 input.getline(buf, buflen);
2058 input >> vari >> NumOfVertices >> NumOfElements;
2059 input.getline(buf, buflen);
2060 input.getline(buf, buflen);
2061 input >> vari >> vari >> NumOfBdrElements;
2062 input.getline(buf, buflen);
2063 input.getline(buf, buflen);
2064 input.getline(buf, buflen);
2066 vertices.SetSize(NumOfVertices);
2067 for (i = 0; i < NumOfVertices; i++)
2069 input >> vari >> varf >> vertices[i](0) >> vertices[i](1)
2071 input.getline(buf, buflen);
2074 elements.SetSize(NumOfElements);
2075 for (i = 0; i < NumOfElements; i++)
2077 input >> vari >> attr;
2078 for (j = 0; j < 8; j++)
2083 input.getline(buf, buflen);
2087 boundary.SetSize(NumOfBdrElements);
2088 for (i = 0; i < NumOfBdrElements; i++)
2091 for (j = 0; j < 4; j++)
2096 input.getline(buf, buflen);
2101 else if (mesh_type ==
"# vtk DataFile Version 3.0" ||
2102 mesh_type ==
"# vtk DataFile Version 2.0")
2107 getline(input, buff);
2108 getline(input, buff);
2109 if (buff !=
"ASCII")
2111 mfem_error(
"Mesh::Load : VTK mesh is not in ASCII format!");
2114 getline(input, buff);
2115 if (buff !=
"DATASET UNSTRUCTURED_GRID")
2117 mfem_error(
"Mesh::Load : VTK mesh is not UNSTRUCTURED_GRID!");
2127 mfem_error(
"Mesh::Load : VTK mesh does not have POINTS data!");
2129 while (buff !=
"POINTS");
2135 getline(input, buff);
2136 for (i = 0; i < points.
Size(); i++)
2141 NumOfElements = n = 0;
2143 input >> ws >> buff;
2144 if (buff ==
"CELLS")
2146 input >> NumOfElements >> n >> ws;
2148 for (i = 0; i < n; i++)
2149 input >> cells_data[i];
2155 input >> ws >> buff;
2156 if (buff ==
"CELL_TYPES")
2158 input >> NumOfElements;
2159 elements.
SetSize(NumOfElements);
2160 for (j = i = 0; i < NumOfElements; i++)
2168 elements[i] =
new Triangle(&cells_data[j+1]);
2176 #ifdef MFEM_USE_MEMALLOC
2177 elements[i] = TetMemory.Alloc();
2178 elements[i]->SetVertices(&cells_data[j+1]);
2185 elements[i] =
new Hexahedron(&cells_data[j+1]);
2191 elements[i] =
new Triangle(&cells_data[j+1]);
2201 #ifdef MFEM_USE_MEMALLOC
2202 elements[i] = TetMemory.Alloc();
2203 elements[i]->SetVertices(&cells_data[j+1]);
2211 elements[i] =
new Hexahedron(&cells_data[j+1]);
2214 cerr <<
"Mesh::Load : VTK mesh : cell type " << ct
2215 <<
" is not supported!" << endl;
2219 j += cells_data[j] + 1;
2224 streampos sp = input.tellg();
2225 input >> ws >> buff;
2226 if (buff ==
"CELL_DATA")
2229 getline(input, buff);
2230 if (buff ==
"SCALARS material int" || buff ==
"SCALARS material float")
2232 getline(input, buff);
2233 for (i = 0; i < NumOfElements; i++)
2236 elements[i]->SetAttribute(attr);
2249 vertices.SetSize(np);
2250 for (i = 0; i < np; i++)
2252 vertices[i](0) = points(3*i+0);
2253 vertices[i](1) = points(3*i+1);
2254 vertices[i](2) = points(3*i+2);
2259 NumOfBdrElements = 0;
2261 else if (order == 2)
2268 for (n = i = 0; i < NumOfElements; i++)
2270 int *v = elements[i]->GetVertices();
2271 int nv = elements[i]->GetNVertices();
2272 for (j = 0; j < nv; j++)
2273 if (pts_dof[v[j]] == -1)
2274 pts_dof[v[j]] = n++;
2277 for (n = i = 0; i < np; i++)
2278 if (pts_dof[i] != -1)
2281 for (i = 0; i < NumOfElements; i++)
2283 int *v = elements[i]->GetVertices();
2284 int nv = elements[i]->GetNVertices();
2285 for (j = 0; j < nv; j++)
2286 v[j] = pts_dof[v[j]];
2291 for (i = 0; i < np; i++)
2293 if ((j = pts_dof[i]) != -1)
2295 vertices[j](0) = points(3*i+0);
2296 vertices[j](1) = points(3*i+1);
2297 vertices[j](2) = points(3*i+2);
2302 NumOfBdrElements = 0;
2310 GetElementToFaceTable();
2317 el_to_edge =
new Table;
2318 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2326 Nodes->MakeOwner(fec);
2331 for (n = i = 0; i < NumOfElements; i++)
2334 const int *vtk_mfem;
2335 switch (elements[i]->GetGeometryType())
2339 vtk_mfem = vtk_quadratic_hex;
break;
2341 vtk_mfem = vtk_quadratic_tet;
break;
2344 vtk_mfem = vtk_quadratic_hex;
break;
2347 for (n++, j = 0; j < dofs.
Size(); j++, n++)
2349 if (pts_dof[cells_data[n]] == -1)
2351 pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
2355 if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
2357 "inconsistent quadratic mesh!");
2363 for (i = 0; i < np; i++)
2366 if ((dofs[0] = pts_dof[i]) != -1)
2369 for (j = 0; j < dofs.
Size(); j++)
2370 (*Nodes)(dofs[j]) = points(3*i+j);
2377 else if (mesh_type ==
"MFEM NURBS mesh v1.0")
2381 Dim = NURBSext->Dimension();
2382 NumOfVertices = NURBSext->GetNV();
2383 NumOfElements = NURBSext->GetNE();
2384 NumOfBdrElements = NURBSext->GetNBE();
2386 NURBSext->GetElementTopo(elements);
2387 NURBSext->GetBdrElementTopo(boundary);
2389 vertices.SetSize(NumOfVertices);
2391 if (NURBSext->HavePatches())
2397 Nodes->MakeOwner(fec);
2398 NURBSext->SetCoordsFromPatches(*Nodes);
2401 int vd = Nodes->VectorDim();
2402 for (i = 0; i < vd; i++)
2405 Nodes->GetNodalValues(vert_val, i+1);
2406 for (j = 0; j < NumOfVertices; j++)
2407 vertices[j](i) = vert_val(j);
2413 else if (mesh_type ==
"MFEM INLINE mesh v1.0")
2441 if (input.get() !=
'=')
2443 std::ostringstream os;
2444 os <<
"Mesh::Load : Inline mesh expected '=' after keyword " << name;
2453 else if (name ==
"ny")
2457 else if (name ==
"nz")
2461 else if (name ==
"sx")
2465 else if (name ==
"sy")
2469 else if (name ==
"sz")
2473 else if (name ==
"type")
2477 if (eltype ==
"segment")
2481 else if (eltype ==
"quad")
2485 else if (eltype ==
"tri")
2489 else if (eltype ==
"hex")
2493 else if (eltype ==
"tet")
2499 std::ostringstream os;
2500 os <<
"Mesh::Load : unrecognized element type (read '" << eltype
2501 <<
"') in inline mesh format. Allowed: segment, tri, tet, quad, hex";
2507 std::ostringstream os;
2508 os <<
"Mesh::Load : unrecognized keyword (" << name
2509 <<
") in inline mesh format. Allowed: nx, ny, nz, type, sx, sy, sz";
2515 if (input.peek() ==
';')
2530 if (nx < 0 || sx < 0.0)
2532 std::ostringstream os;
2533 os <<
"Mesh::Load : invalid 1D inline mesh format, all values must be "
2535 <<
" nx = " << nx <<
"\n"
2536 <<
" sx = " << sx <<
"\n";
2543 if (nx < 0 || ny < 0 || sx < 0.0 || sy < 0.0)
2545 std::ostringstream os;
2546 os <<
"Mesh::Load : invalid 2D inline mesh format, all values must be "
2548 <<
" nx = " << nx <<
"\n"
2549 <<
" ny = " << ny <<
"\n"
2550 <<
" sx = " << sx <<
"\n"
2551 <<
" sy = " << sy <<
"\n";
2554 Make2D(nx, ny, type, generate_edges, sx, sy);
2558 if (nx < 0 || ny < 0 || nz < 0 || sx < 0.0 || sy < 0.0 || sz < 0.0)
2560 std::ostringstream os;
2561 os <<
"Mesh::Load : invalid 3D inline mesh format, all values must be "
2563 <<
" nx = " << nx <<
"\n"
2564 <<
" ny = " << ny <<
"\n"
2565 <<
" nz = " << nz <<
"\n"
2566 <<
" sx = " << sx <<
"\n"
2567 <<
" sy = " << sy <<
"\n"
2568 <<
" sz = " << sz <<
"\n";
2571 Make3D(nx, ny, nz, type, generate_edges, sx, sy, sz);
2575 mfem_error(
"Mesh::Load : For inline mesh, must specify an "
2576 "element type = [segment, tri, quad, tet, hex]");
2582 MFEM_ABORT(
"Unknown input mesh format");
2604 if (NumOfBdrElements == 0 && Dim > 2)
2607 GetElementToFaceTable();
2609 GenerateBoundaryElements();
2615 CheckElementOrientation(fix_orientation);
2618 MarkForRefinement();
2627 GetElementToFaceTable();
2630 if ( !(curved && (meshgen & 1)) )
2631 CheckBdrElementOrientation();
2637 if (Dim > 1 && generate_edges == 1)
2639 el_to_edge =
new Table;
2640 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2644 if (NumOfBdrElements == 0)
2645 GenerateBoundaryElements();
2647 if ( !(curved && (meshgen & 1)) )
2648 CheckBdrElementOrientation();
2650 c_el_to_edge = NULL;
2664 spaceDim = Nodes->VectorDim();
2666 for (i = 0; i < spaceDim; i++)
2669 Nodes->GetNodalValues(vert_val, i+1);
2670 for (j = 0; j < NumOfVertices; j++)
2671 vertices[j](i) = vert_val(j);
2679 Table *old_elem_vert = NULL;
2680 if (fix_orientation || refine)
2682 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
2687 CheckElementOrientation(fix_orientation);
2689 MarkForRefinement();
2691 if (fix_orientation || refine)
2693 DoNodeReorder(old_v_to_v, old_elem_vert);
2694 delete old_elem_vert;
2703 int i, j, ie, ib, iv, *v, nv;
2713 if (mesh_array[0]->NURBSext)
2718 NumOfVertices = NURBSext->GetNV();
2719 NumOfElements = NURBSext->GetNE();
2721 NURBSext->GetElementTopo(elements);
2734 NumOfBdrElements = 0;
2735 for (i = 0; i < num_pieces; i++)
2736 NumOfBdrElements += mesh_array[i]->GetNBE();
2737 boundary.SetSize(NumOfBdrElements);
2738 vertices.SetSize(NumOfVertices);
2740 for (i = 0; i < num_pieces; i++)
2746 for (j = 0; j < m->
GetNE(); j++)
2747 elements[lelem_elem[j]]->SetAttribute(m->
GetAttribute(j));
2749 for (j = 0; j < m->
GetNBE(); j++)
2754 for (
int k = 0; k < nv; k++)
2755 v[k] = lvert_vert[v[k]];
2756 boundary[ib++] = el;
2759 for (j = 0; j < m->
GetNV(); j++)
2760 vertices[lvert_vert[j]].SetCoords(m->
GetVertex(j));
2766 NumOfBdrElements = 0;
2768 for (i = 0; i < num_pieces; i++)
2771 NumOfElements += m->
GetNE();
2772 NumOfBdrElements += m->
GetNBE();
2773 NumOfVertices += m->
GetNV();
2775 elements.SetSize(NumOfElements);
2776 boundary.SetSize(NumOfBdrElements);
2777 vertices.SetSize(NumOfVertices);
2779 for (i = 0; i < num_pieces; i++)
2783 for (j = 0; j < m->
GetNE(); j++)
2788 for (
int k = 0; k < nv; k++)
2790 elements[ie++] = el;
2793 for (j = 0; j < m->
GetNBE(); j++)
2798 for (
int k = 0; k < nv; k++)
2800 boundary[ib++] = el;
2803 for (j = 0; j < m->
GetNV(); j++)
2804 vertices[iv++].SetCoords(m->
GetVertex(j));
2810 for (i = 0; i < num_pieces; i++)
2811 meshgen |= mesh_array[i]->MeshGenerator();
2816 GetElementToFaceTable();
2825 el_to_edge =
new Table;
2826 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2841 for (i = 0; i < num_pieces; i++)
2842 gf_array[i] = mesh_array[i]->GetNodes();
2848 CheckElementOrientation(
false);
2849 CheckBdrElementOrientation(
false);
2855 if (NURBSext == NULL)
2856 mfem_error(
"Mesh::KnotInsert : Not a NURBS mesh!");
2858 if (kv.
Size() != NURBSext->GetNKV())
2859 mfem_error(
"Mesh::KnotInsert : KnotVector array size mismatch!");
2861 NURBSext->ConvertToPatches(*Nodes);
2863 NURBSext->KnotInsert(kv);
2871 NURBSext->ConvertToPatches(*Nodes);
2873 NURBSext->UniformRefinement();
2880 if (NURBSext == NULL)
2881 mfem_error(
"Mesh::DegreeElevate : Not a NURBS mesh!");
2883 NURBSext->ConvertToPatches(*Nodes);
2885 NURBSext->DegreeElevate(t);
2898 NURBSext->SetKnotsFromPatches();
2900 Dim = NURBSext->Dimension();
2903 if (NumOfElements != NURBSext->GetNE())
2905 for (
int i = 0; i < elements.Size(); i++)
2906 FreeElement(elements[i]);
2907 NumOfElements = NURBSext->GetNE();
2908 NURBSext->GetElementTopo(elements);
2911 if (NumOfBdrElements != NURBSext->GetNBE())
2913 for (
int i = 0; i < boundary.Size(); i++)
2914 FreeElement(boundary[i]);
2915 NumOfBdrElements = NURBSext->GetNBE();
2916 NURBSext->GetBdrElementTopo(boundary);
2919 Nodes->FESpace()->Update();
2921 NURBSext->SetCoordsFromPatches(*Nodes);
2923 if (NumOfVertices != NURBSext->GetNV())
2925 NumOfVertices = NURBSext->GetNV();
2926 vertices.SetSize(NumOfVertices);
2927 int vd = Nodes->VectorDim();
2928 for (
int i = 0; i < vd; i++)
2931 Nodes->GetNodalValues(vert_val, i+1);
2932 for (
int j = 0; j < NumOfVertices; j++)
2933 vertices[j](i) = vert_val(j);
2939 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
2946 GetElementToFaceTable();
2970 input >> NumOfElements;
2971 elements.SetSize(NumOfElements);
2972 for (j = 0; j < NumOfElements; j++)
2973 elements[j] = ReadElement(input);
2978 input >> NumOfBdrElements;
2979 boundary.SetSize(NumOfBdrElements);
2980 for (j = 0; j < NumOfBdrElements; j++)
2981 boundary[j] = ReadElement(input);
2986 input >> NumOfEdges;
2987 edge_vertex =
new Table(NumOfEdges, 2);
2988 edge_to_knot.
SetSize(NumOfEdges);
2989 for (j = 0; j < NumOfEdges; j++)
2991 int *v = edge_vertex->GetRow(j);
2992 input >> edge_to_knot[j] >> v[0] >> v[1];
2994 edge_to_knot[j] = -1 - edge_to_knot[j];
3000 input >> NumOfVertices;
3001 vertices.SetSize(0);
3008 GetElementToFaceTable();
3010 if (NumOfBdrElements == 0)
3011 GenerateBoundaryElements();
3012 CheckBdrElementOrientation();
3020 el_to_edge =
new Table;
3021 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
3025 if (NumOfBdrElements == 0)
3026 GenerateBoundaryElements();
3027 CheckBdrElementOrientation();
3041 for (
int d = 0; d < v.
Size(); d++)
3047 for (d = 0; d < p.
Size(); d++)
3049 for ( ; d < v.
Size(); d++)
3056 if (Nodes == NULL || Nodes->FESpace() != nodes.
FESpace())
3069 SetNodalGridFunction(nodes,
true);
3075 NewNodes(*nodes, make_owner);
3080 return ((Nodes) ? Nodes->FESpace() : NULL);
3087 case 1:
return GetNV();
3088 case 2:
return GetNEdges();
3089 case 3:
return GetNFaces();
3094 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3095 static const char *fixed_or_not[] = {
"fixed",
"NOT FIXED" };
3100 int i, j, k, wo = 0, fo = 0, *vi = 0;
3103 if (Dim == 2 && spaceDim == 2)
3107 for (i = 0; i < NumOfElements; i++)
3111 vi = elements[i]->GetVertices();
3112 for (j = 0; j < 3; j++)
3113 v[j] = vertices[vi[j]]();
3114 for (j = 0; j < 2; j++)
3115 for (k = 0; k < 2; k++)
3116 J(j, k) = v[j+1][k] - v[0][k];
3121 GetElementJacobian(i, J);
3127 switch (GetElementType(i))
3147 for (i = 0; i < NumOfElements; i++)
3149 vi = elements[i]->GetVertices();
3150 switch (GetElementType(i))
3155 for (j = 0; j < 4; j++)
3156 v[j] = vertices[vi[j]]();
3157 for (j = 0; j < 3; j++)
3158 for (k = 0; k < 3; k++)
3159 J(j, k) = v[j+1][k] - v[0][k];
3164 GetElementJacobian(i, J);
3179 GetElementJacobian(i, J);
3192 #if (!defined(MFEM_USE_MPI) || defined(MFEM_DEBUG))
3194 cout <<
"Elements with wrong orientation: " << wo <<
" / "
3195 << NumOfElements <<
" (" << fixed_or_not[(wo == fo) ? 0 : 1]
3204 if (test[0] == base[0])
3205 if (test[1] == base[1])
3209 else if (test[0] == base[1])
3210 if (test[1] == base[0])
3215 if (test[1] == base[0])
3221 const int *aor = tri_orientations[orient];
3222 for (
int j = 0; j < 3; j++)
3223 if (test[aor[j]] != base[j])
3234 for (i = 0; i < 4; i++)
3235 if (test[i] == base[0])
3240 if (test[(i+1)%4] == base[1])
3244 const int *aor = quad_orientations[orient];
3245 for (
int j = 0; j < 4; j++)
3246 if (test[aor[j]] != base[j])
3248 cerr <<
"Mesh::GetQuadOrientation(...)" << endl;
3249 cerr <<
" base = [";
3250 for (
int k = 0; k < 4; k++)
3251 cerr <<
" " << base[k];
3252 cerr <<
" ]\n test = [";
3253 for (
int k = 0; k < 4; k++)
3254 cerr <<
" " << test[k];
3255 cerr <<
" ]" << endl;
3260 if (test[(i+1)%4] == base[1])
3272 for (i = 0; i < NumOfBdrElements; i++)
3274 if (faces_info[be_to_edge[i]].Elem2No < 0)
3276 int *bv = boundary[i]->GetVertices();
3277 int *fv = faces[be_to_edge[i]]->GetVertices();
3281 mfem::Swap<int>(bv[0], bv[1]);
3293 for (i = 0; i < NumOfBdrElements; i++)
3295 if (faces_info[be_to_face[i]].Elem2No < 0)
3297 bv = boundary[i]->GetVertices();
3298 el = faces_info[be_to_face[i]].Elem1No;
3299 ev = elements[el]->GetVertices();
3300 switch (GetElementType(el))
3304 int *fv = faces[be_to_face[i]]->GetVertices();
3307 orientation = GetTriOrientation(fv, bv);
3308 if (orientation % 2)
3313 mfem::Swap<int>(bv[0], bv[1]);
3321 int lf = faces_info[be_to_face[i]].Elem1Inf/64;
3322 for (
int j = 0; j < 4; j++)
3323 v[j] = ev[hex_faces[lf][j]];
3324 if (GetQuadOrientation(v, bv) % 2)
3327 mfem::Swap<int>(bv[0], bv[2]);
3339 cout <<
"Boundary elements with wrong orientation: " << wo <<
" / "
3340 << NumOfBdrElements <<
" (" << fixed_or_not[fix_it ? 0 : 1]
3349 el_to_edge->GetRow(i, edges);
3351 mfem_error(
"Mesh::GetElementEdges(...) element to edge table "
3352 "is not generated.");
3354 const int *v = elements[i]->GetVertices();
3355 const int ne = elements[i]->GetNEdges();
3357 for (
int j = 0; j < ne; j++)
3359 const int *e = elements[i]->GetEdgeVertices(j);
3360 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
3371 edges[0] = be_to_edge[i];
3372 const int *v = boundary[i]->GetVertices();
3373 cor[0] = (v[0] < v[1]) ? (1) : (-1);
3378 bel_to_edge->GetRow(i, edges);
3382 const int *v = boundary[i]->GetVertices();
3383 const int ne = boundary[i]->GetNEdges();
3385 for (
int j = 0; j < ne; j++)
3387 const int *e = boundary[i]->GetEdgeVertices(j);
3388 cor[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
3400 const int *v = faces[i]->GetVertices();
3401 o[0] = (v[0] < v[1]) ? (1) : (-1);
3409 face_edge->GetRow(i, edges);
3411 const int *v = faces[i]->GetVertices();
3412 const int ne = faces[i]->GetNEdges();
3414 for (
int j = 0; j < ne; j++)
3416 const int *e = faces[i]->GetEdgeVertices(j);
3417 o[j] = (v[e[0]] < v[e[1]]) ? (1) : (-1);
3425 GetEdgeVertexTable();
3426 edge_vertex->GetRow(i, vert);
3438 if (faces.Size() != NumOfFaces)
3439 mfem_error(
"Mesh::GetFaceEdgeTable : faces were not generated!");
3442 DSTable v_to_v(NumOfVertices);
3443 GetVertexToVertexTable(v_to_v);
3445 face_edge =
new Table;
3446 GetElementArrayEdgeTable(faces, v_to_v, *face_edge);
3456 DSTable v_to_v(NumOfVertices);
3457 GetVertexToVertexTable(v_to_v);
3460 edge_vertex =
new Table(nedges, 2);
3461 for (
int i = 0; i < NumOfVertices; i++)
3466 edge_vertex->Push(j, i);
3467 edge_vertex->Push(j, it.Column());
3470 edge_vertex->Finalize();
3481 vert_elem->
MakeI(NumOfVertices);
3483 for (i = 0; i < NumOfElements; i++)
3485 nv = elements[i]->GetNVertices();
3486 v = elements[i]->GetVertices();
3487 for (j = 0; j < nv; j++)
3493 for (i = 0; i < NumOfElements; i++)
3495 nv = elements[i]->GetNVertices();
3496 v = elements[i]->GetVertices();
3497 for (j = 0; j < nv; j++)
3510 face_elem->
MakeI(faces_info.Size());
3512 for (
int i = 0; i < faces_info.Size(); i++)
3514 if (faces_info[i].Elem2No >= 0)
3522 for (
int i = 0; i < faces_info.Size(); i++)
3525 if (faces_info[i].Elem2No >= 0)
3540 el_to_face->GetRow(i, fcs);
3542 mfem_error(
"Mesh::GetElementFaces(...) : el_to_face not generated.");
3546 for (j = 0; j < n; j++)
3547 if (faces_info[fcs[j]].Elem1No == i)
3548 cor[j] = faces_info[fcs[j]].Elem1Inf % 64;
3550 else if (faces_info[fcs[j]].Elem2No == i)
3551 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
3553 mfem_error(
"Mesh::GetElementFaces(...) : 2");
3556 cor[j] = faces_info[fcs[j]].Elem2Inf % 64;
3571 bv = boundary[i]->GetVertices();
3572 fv = faces[be_to_face[i]]->GetVertices();
3576 switch (GetBdrElementType(i))
3579 *o = GetTriOrientation(fv, bv);
3582 *o = GetQuadOrientation(fv, bv);
3585 mfem_error(
"Mesh::GetBdrElementFace(...) 2");
3592 switch (GetElementType(0))
3606 mfem_error(
"Mesh::GetFaceBaseGeometry(...) #1");
3610 if (faces[i] == NULL)
3611 switch (GetElementType(faces_info[i].Elem1No))
3618 mfem_error(
"Mesh::GetFaceBaseGeometry(...) #2");
3621 return faces[i]->GetGeometryType();
3628 return be_to_edge[i];
3629 return be_to_face[i];
3665 v = elements[i]->GetVertices();
3666 nv = elements[i]->GetNVertices();
3668 pointmat.
SetSize(spaceDim, nv);
3669 for (k = 0; k < spaceDim; k++)
3670 for (j = 0; j < nv; j++)
3671 pointmat(k, j) = vertices[v[j]](k);
3679 v = boundary[i]->GetVertices();
3680 nv = boundary[i]->GetNVertices();
3682 pointmat.
SetSize(spaceDim, nv);
3683 for (k = 0; k < spaceDim; k++)
3684 for (j = 0; j < nv; j++)
3685 pointmat(k, j) = vertices[v[j]](k);
3690 const double *vi = vertices[i]();
3691 const double *vj = vertices[j]();
3694 for (
int k = 0; k < spaceDim; k++)
3695 length += (vi[k]-vj[k])*(vi[k]-vj[k]);
3697 return sqrt(length);
3705 for (
int i = 0; i < elem_array.
Size(); i++)
3710 for (
int i = 0; i < elem_array.
Size(); i++)
3712 const int *v = elem_array[i]->GetVertices();
3713 const int ne = elem_array[i]->GetNEdges();
3714 for (
int j = 0; j < ne; j++)
3716 const int *e = elem_array[i]->GetEdgeVertices(j);
3727 for (
int i = 0; i < edge_vertex->Size(); i++)
3729 const int *v = edge_vertex->GetRow(i);
3730 v_to_v.
Push(v[0], v[1]);
3735 for (
int i = 0; i < NumOfElements; i++)
3737 const int *v = elements[i]->GetVertices();
3738 const int ne = elements[i]->GetNEdges();
3739 for (
int j = 0; j < ne; j++)
3741 const int *e = elements[i]->GetEdgeVertices(j);
3742 v_to_v.
Push(v[e[0]], v[e[1]]);
3750 int i, NumberOfEdges;
3752 DSTable v_to_v(NumOfVertices);
3753 GetVertexToVertexTable(v_to_v);
3758 GetElementArrayEdgeTable(elements, v_to_v, e_to_f);
3763 be_to_f.
SetSize(NumOfBdrElements);
3764 for (i = 0; i < NumOfBdrElements; i++)
3766 const int *v = boundary[i]->GetVertices();
3767 be_to_f[i] = v_to_v(v[0], v[1]);
3772 if (bel_to_edge == NULL)
3773 bel_to_edge =
new Table;
3774 GetElementArrayEdgeTable(boundary, v_to_v, *bel_to_edge);
3777 mfem_error(
"1D GetElementToEdgeTable is not yet implemented.");
3780 return NumberOfEdges;
3793 f_el = GetVertexToElementTable();
3801 el_to_el =
new Table;
3802 el_to_el->
MakeI(NumOfElements);
3803 for (
int i = 0; i < f_el->
Size(); i++)
3807 const int *el = f_el->
GetRow(i);
3808 el_to_el->AddAColumnInRow(el[0]);
3809 el_to_el->AddAColumnInRow(el[1]);
3813 for (
int i = 0; i < f_el->
Size(); i++)
3817 const int *el = f_el->
GetRow(i);
3818 el_to_el->AddConnection(el[0], el[1]);
3819 el_to_el->AddConnection(el[1], el[0]);
3822 el_to_el->ShiftUpI();
3828 el_to_el =
new Table(NumOfElements, 6);
3831 if (faces_info.Size() != NumOfFaces)
3832 mfem_error(
"Mesh::ElementToElementTable : faces were not generated!");
3835 for (
int i = 0; i < faces_info.Size(); i++)
3836 if (faces_info[i].Elem2No >= 0)
3838 el_to_el->Push(faces_info[i].Elem1No, faces_info[i].Elem2No);
3839 el_to_el->Push(faces_info[i].Elem2No, faces_info[i].Elem1No);
3842 el_to_el->Finalize();
3850 if (el_to_face == NULL)
3857 if (el_to_edge == NULL)
3864 if (faces_info[gf].Elem1No == -1)
3867 faces_info[gf].Elem1No = el;
3868 faces_info[gf].Elem1Inf = 64 * lf;
3869 faces_info[gf].Elem2No = -1;
3870 faces_info[gf].Elem2Inf = -1;
3879 faces_info[gf].Elem2No = el;
3880 faces_info[gf].Elem2Inf = 64 * lf + 1;
3886 if (faces[gf] == NULL)
3888 faces[gf] =
new Segment(v0, v1);
3889 faces_info[gf].Elem1No = el;
3890 faces_info[gf].Elem1Inf = 64 * lf;
3891 faces_info[gf].Elem2No = -1;
3892 faces_info[gf].Elem2Inf = -1;
3897 int *v = faces[gf]->GetVertices();
3898 if (v[1] != v0 || v[0] != v1)
3899 mfem_error(
"Mesh::AddSegmentFaceElement(...)");
3901 faces_info[gf].Elem2No = el;
3902 faces_info[gf].Elem2Inf = 64 * lf + 1;
3907 int v0,
int v1,
int v2)
3909 if (faces[gf] == NULL)
3911 faces[gf] =
new Triangle(v0, v1, v2);
3912 faces_info[gf].Elem1No = el;
3913 faces_info[gf].Elem1Inf = 64 * lf;
3914 faces_info[gf].Elem2No = -1;
3915 faces_info[gf].Elem2Inf = -1;
3919 int orientation, vv[3] = { v0, v1, v2 };
3920 orientation = GetTriOrientation(faces[gf]->GetVertices(), vv);
3922 if (orientation % 2 == 0)
3923 mfem_error(
"Mesh::AddTriangleFaceElement(...)");
3925 faces_info[gf].Elem2No = el;
3926 faces_info[gf].Elem2Inf = 64 * lf + orientation;
3931 int v0,
int v1,
int v2,
int v3)
3933 if (faces_info[gf].Elem1No < 0)
3936 faces_info[gf].Elem1No = el;
3937 faces_info[gf].Elem1Inf = 64 * lf;
3938 faces_info[gf].Elem2No = -1;
3939 faces_info[gf].Elem2Inf = -1;
3943 int vv[4] = { v0, v1, v2, v3 };
3944 int oo = GetQuadOrientation(faces[gf]->GetVertices(), vv);
3949 faces_info[gf].Elem2No = el;
3950 faces_info[gf].Elem2Inf = 64 * lf + oo;
3958 nfaces = (Dim == 1) ? NumOfVertices : ((Dim == 2) ? NumOfEdges : NumOfFaces);
3960 for (i = 0; i < faces.Size(); i++)
3961 FreeElement(faces[i]);
3964 faces.SetSize(nfaces);
3965 faces_info.SetSize(nfaces);
3966 for (i = 0; i < nfaces; i++)
3969 faces_info[i].Elem1No = -1;
3971 for (i = 0; i < NumOfElements; i++)
3973 const int *v = elements[i]->GetVertices();
3977 AddPointFaceElement(0, v[0], i);
3978 AddPointFaceElement(1, v[1], i);
3982 ef = el_to_edge->GetRow(i);
3983 const int ne = elements[i]->GetNEdges();
3984 for (
int j = 0; j < ne; j++)
3986 const int *e = elements[i]->GetEdgeVertices(j);
3987 AddSegmentFaceElement(j, ef[j], i, v[e[0]], v[e[1]]);
3992 ef = el_to_face->GetRow(i);
3993 switch (GetElementType(i))
3996 for (
int j = 0; j < 4; j++)
3998 const int *fv = tet_faces[j];
3999 AddTriangleFaceElement(j, ef[j], i,
4000 v[fv[0]], v[fv[1]], v[fv[2]]);
4004 for (
int j = 0; j < 6; j++)
4006 const int *fv = hex_faces[j];
4007 AddQuadFaceElement(j, ef[j], i,
4008 v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4013 MFEM_ABORT(
"Unexpected type of Element.");
4023 for (
int i = 0; i < NumOfElements; i++)
4025 const int *v = elements[i]->GetVertices();
4026 switch (GetElementType(i))
4029 for (
int j = 0; j < 4; j++)
4031 const int *fv = tet_faces[j];
4032 faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
4038 for (
int j = 0; j < 6; j++)
4040 const int *fv = hex_faces[j];
4041 faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
4046 MFEM_ABORT(
"Unexpected type of Element.");
4058 if (el_to_face != NULL)
4060 el_to_face =
new Table(NumOfElements, 6);
4061 faces_tbl =
new STable3D(NumOfVertices);
4062 for (i = 0; i < NumOfElements; i++)
4064 v = elements[i]->GetVertices();
4065 switch (GetElementType(i))
4068 for (
int j = 0; j < 4; j++)
4070 const int *fv = tet_faces[j];
4072 i, faces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]));
4078 for (
int j = 0; j < 6; j++)
4080 const int *fv = hex_faces[j];
4082 i, faces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
4087 MFEM_ABORT(
"Unexpected type of Element.");
4091 el_to_face->Finalize();
4093 be_to_face.SetSize(NumOfBdrElements);
4094 for (i = 0; i < NumOfBdrElements; i++)
4096 v = boundary[i]->GetVertices();
4097 switch (GetBdrElementType(i))
4100 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2]);
4103 be_to_face[i] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
4107 MFEM_ABORT(
"Unexpected type of boundary Element.");
4122 if (Dim != 3 || !(meshgen & 1))
4126 Table *old_elem_vert = NULL;
4129 PrepareNodeReorder(&old_v_to_v, &old_elem_vert);
4131 DeleteCoarseTables();
4133 for (
int i = 0; i < NumOfElements; i++)
4136 v = elements[i]->GetVertices();
4138 Rotate3(v[0], v[1], v[2]);
4140 Rotate3(v[1], v[2], v[3]);
4142 ShiftL2R(v[0], v[1], v[3]);
4145 for (
int i = 0; i < NumOfBdrElements; i++)
4148 v = boundary[i]->GetVertices();
4150 Rotate3(v[0], v[1], v[2]);
4155 GetElementToFaceTable();
4158 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
4162 DoNodeReorder(old_v_to_v, old_elem_vert);
4163 delete old_elem_vert;
4170 static int mfem_less(
const void *x,
const void *y)
4172 if (*(
int*)x < *(
int*)y)
4174 if (*(
int*)x > *(
int*)y)
4178 #ifndef MFEM_USE_METIS_5
4183 int*,
int*,
int*,
int*,
int*,
idxtype*);
4185 int*,
int*,
int*,
int*,
int*,
idxtype*);
4187 int*,
int*,
int*,
int*,
int*,
idxtype*);
4212 double pmin[3] = { numeric_limits<double>::infinity(),
4213 numeric_limits<double>::infinity(),
4214 numeric_limits<double>::infinity()};
4215 double pmax[3] = { -numeric_limits<double>::infinity(),
4216 -numeric_limits<double>::infinity(),
4217 -numeric_limits<double>::infinity()};
4219 for (
int vi = 0; vi < NumOfVertices; vi++)
4221 const double *p = vertices[vi]();
4222 for (
int i = 0; i < spaceDim; i++)
4224 if (p[i] < pmin[i]) pmin[i] = p[i];
4225 if (p[i] > pmax[i]) pmax[i] = p[i];
4229 partitioning =
new int[NumOfElements];
4233 Vector pt(ppt, spaceDim);
4234 for (
int el = 0; el < NumOfElements; el++)
4236 GetElementTransformation(el)->Transform(
4239 for (
int i = spaceDim-1; i >= 0; i--)
4241 int idx = (int)floor(nxyz[i]*((pt(i) - pmin[i])/(pmax[i] - pmin[i])));
4242 if (idx < 0) idx = 0;
4243 if (idx >= nxyz[i]) idx = nxyz[i]-1;
4244 part = part * nxyz[i] + idx;
4246 partitioning[el] = part;
4249 return partitioning;
4255 int i, *partitioning;
4257 ElementToElementTable();
4259 partitioning =
new int[NumOfElements];
4263 for (i = 0; i < NumOfElements; i++)
4264 partitioning[i] = 0;
4269 #ifndef MFEM_USE_METIS_5
4281 I = el_to_el->GetI();
4282 J = el_to_el->GetJ();
4283 #ifndef MFEM_USE_METIS_5
4291 if (part_method >= 0 && part_method <= 2)
4292 for (i = 0; i < n; i++)
4293 qsort(&J[I[i]], I[i+1]-I[i],
sizeof(
int), &mfem_less);
4297 if (part_method == 0 || part_method == 3)
4299 #ifndef MFEM_USE_METIS_5
4327 " error in METIS_PartGraphRecursive!");
4333 if (part_method == 1 || part_method == 4)
4335 #ifndef MFEM_USE_METIS_5
4363 " error in METIS_PartGraphKway!");
4369 if (part_method == 2 || part_method == 5)
4371 #ifndef MFEM_USE_METIS_5
4400 " error in METIS_PartGraphKway!");
4405 cout <<
"Mesh::GeneratePartitioning(...): edgecut = "
4417 for (i = 0; i < nparts; i++)
4423 for (i = 0; i < NumOfElements; i++)
4424 psize[partitioning[i]].one++;
4426 int empty_parts = 0;
4427 for (i = 0; i < nparts; i++)
4428 if (psize[i].one == 0)
4435 cerr <<
"Mesh::GeneratePartitioning returned " << empty_parts
4436 <<
" empty parts!" << endl;
4440 for (i = nparts-1; i > nparts-1-empty_parts; i--)
4443 for (
int j = 0; j < NumOfElements; j++)
4444 for (i = nparts-1; i > nparts-1-empty_parts; i--)
4445 if (psize[i].one == 0 || partitioning[j] != psize[i].two)
4449 partitioning[j] = psize[nparts-1-i].two;
4455 return partitioning;
4459 mfem_error(
"Mesh::GeneratePartitioning(...): "
4460 "MFEM was compiled without Metis.");
4474 int num_elem, *i_elem_elem, *j_elem_elem;
4476 num_elem = elem_elem.
Size();
4477 i_elem_elem = elem_elem.
GetI();
4478 j_elem_elem = elem_elem.
GetJ();
4483 int stack_p, stack_top_p, elem;
4487 for (i = 0; i < num_elem; i++)
4489 if (partitioning[i] > num_part)
4490 num_part = partitioning[i];
4496 for (i = 0; i < num_part; i++)
4501 for (elem = 0; elem < num_elem; elem++)
4503 if (component[elem] >= 0)
4506 component[elem] = num_comp[partitioning[elem]]++;
4508 elem_stack[stack_top_p++] = elem;
4510 for ( ; stack_p < stack_top_p; stack_p++)
4512 i = elem_stack[stack_p];
4513 for (j = i_elem_elem[i]; j < i_elem_elem[i+1]; j++)
4516 if (partitioning[k] == partitioning[i])
4518 if (component[k] < 0)
4520 component[k] = component[i];
4521 elem_stack[stack_top_p++] = k;
4523 else if (component[k] != component[i])
4535 int i, n_empty, n_mcomp;
4537 const Array<int> _partitioning(partitioning, GetNE());
4539 ElementToElementTable();
4543 n_empty = n_mcomp = 0;
4544 for (i = 0; i < num_comp.
Size(); i++)
4545 if (num_comp[i] == 0)
4547 else if (num_comp[i] > 1)
4552 cout <<
"Mesh::CheckPartitioning(...) :\n"
4553 <<
"The following subdomains are empty :\n";
4554 for (i = 0; i < num_comp.
Size(); i++)
4555 if (num_comp[i] == 0)
4561 cout <<
"Mesh::CheckPartitioning(...) :\n"
4562 <<
"The following subdomains are NOT connected :\n";
4563 for (i = 0; i < num_comp.
Size(); i++)
4564 if (num_comp[i] > 1)
4568 if (n_empty == 0 && n_mcomp == 0)
4569 cout <<
"Mesh::CheckPartitioning(...) : "
4570 "All subdomains are connected." << endl;
4582 const double *a = A.
Data();
4583 const double *b = B.
Data();
4592 c(0) = a[0]*a[3]-a[1]*a[2];
4593 c(1) = a[0]*b[3]-a[1]*b[2]+b[0]*a[3]-b[1]*a[2];
4594 c(2) = b[0]*b[3]-b[1]*b[2];
4615 c(0) = (a[0] * (a[4] * a[8] - a[5] * a[7]) +
4616 a[1] * (a[5] * a[6] - a[3] * a[8]) +
4617 a[2] * (a[3] * a[7] - a[4] * a[6]));
4619 c(1) = (b[0] * (a[4] * a[8] - a[5] * a[7]) +
4620 b[1] * (a[5] * a[6] - a[3] * a[8]) +
4621 b[2] * (a[3] * a[7] - a[4] * a[6]) +
4623 a[0] * (b[4] * a[8] - b[5] * a[7]) +
4624 a[1] * (b[5] * a[6] - b[3] * a[8]) +
4625 a[2] * (b[3] * a[7] - b[4] * a[6]) +
4627 a[0] * (a[4] * b[8] - a[5] * b[7]) +
4628 a[1] * (a[5] * b[6] - a[3] * b[8]) +
4629 a[2] * (a[3] * b[7] - a[4] * b[6]));
4631 c(2) = (a[0] * (b[4] * b[8] - b[5] * b[7]) +
4632 a[1] * (b[5] * b[6] - b[3] * b[8]) +
4633 a[2] * (b[3] * b[7] - b[4] * b[6]) +
4635 b[0] * (a[4] * b[8] - a[5] * b[7]) +
4636 b[1] * (a[5] * b[6] - a[3] * b[8]) +
4637 b[2] * (a[3] * b[7] - a[4] * b[6]) +
4639 b[0] * (b[4] * a[8] - b[5] * a[7]) +
4640 b[1] * (b[5] * a[6] - b[3] * a[8]) +
4641 b[2] * (b[3] * a[7] - b[4] * a[6]));
4643 c(3) = (b[0] * (b[4] * b[8] - b[5] * b[7]) +
4644 b[1] * (b[5] * b[6] - b[3] * b[8]) +
4645 b[2] * (b[3] * b[7] - b[4] * b[6]));
4687 double a = z(2), b = z(1), c = z(0);
4688 double D = b*b-4*a*c;
4695 x(0) = x(1) = -0.5 * b / a;
4700 x(0) = -(x(1) = fabs(0.5 * sqrt(D) / a));
4707 t = -0.5 * (b + sqrt(D));
4709 t = -0.5 * (b - sqrt(D));
4713 Swap<double>(x(0), x(1));
4720 double a = z(2)/z(3), b = z(1)/z(3), c = z(0)/z(3);
4723 double Q = (a * a - 3 * b) / 9;
4724 double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
4725 double Q3 = Q * Q * Q;
4732 x(0) = x(1) = x(2) = - a / 3;
4736 double sqrtQ = sqrt(Q);
4740 x(0) = -2 * sqrtQ - a / 3;
4741 x(1) = x(2) = sqrtQ - a / 3;
4745 x(0) = x(1) = - sqrtQ - a / 3;
4746 x(2) = 2 * sqrtQ - a / 3;
4753 double theta = acos(R / sqrt(Q3));
4754 double A = -2 * sqrt(Q);
4756 x0 = A * cos(theta / 3) - a / 3;
4757 x1 = A * cos((theta + 2.0 * M_PI) / 3) - a / 3;
4758 x2 = A * cos((theta - 2.0 * M_PI) / 3) - a / 3;
4762 Swap<double>(x0, x1);
4765 Swap<double>(x1, x2);
4767 Swap<double>(x0, x1);
4778 A = -pow(sqrt(R2 - Q3) + R, 1.0/3.0);
4780 A = pow(sqrt(R2 - Q3) - R, 1.0/3.0);
4781 x(0) = A + Q / A - a / 3;
4790 const double factor,
const int Dim)
4792 const double c0 = c(0);
4793 c(0) = c0 * (1.0 - pow(factor, -Dim));
4795 for (
int j = 0; j < nr; j++)
4805 c(0) = c0 * (1.0 - pow(factor, Dim));
4807 for (
int j = 0; j < nr; j++)
4821 int nvs = vertices.Size();
4822 DenseMatrix P, V, DS, PDS(spaceDim), VDS(spaceDim);
4823 Vector c(spaceDim+1), x(spaceDim);
4824 const double factor = 2.0;
4829 for (
int i = 0; i < NumOfElements; i++)
4836 for (
int j = 0; j < spaceDim; j++)
4837 for (
int k = 0; k < nv; k++)
4839 P(j, k) = vertices[v[k]](j);
4840 V(j, k) = displacements(v[k]+j*nvs);
4844 GetTransformationFEforElementType(el->
GetType());
4866 for (
int j = 0; j < nv; j++)
4888 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
4889 for (
int j = 0; j < spaceDim; j++)
4890 vertices[i](j) += displacements(j*nv+i);
4895 int nv = vertices.Size();
4896 vert_coord.
SetSize(nv*spaceDim);
4897 for (
int i = 0; i < nv; i++)
4898 for (
int j = 0; j < spaceDim; j++)
4899 vert_coord(j*nv+i) = vertices[i](j);
4904 for (
int i = 0, nv = vertices.Size(); i < nv; i++)
4905 for (
int j = 0; j < spaceDim; j++)
4906 vertices[i](j) = vert_coord(j*nv+i);
4914 for (
int j = 0; j < spaceDim; j++)
4915 coord[j] = (*Nodes)(fes->
DofToVDof(i, j));
4919 for (
int j = 0; j < spaceDim; j++)
4920 coord[j] = vertices[i](j);
4929 for (
int j = 0; j < spaceDim; j++)
4930 (*Nodes)(fes->
DofToVDof(i, j)) = coord[j];
4934 for (
int j = 0; j < spaceDim; j++)
4935 vertices[i](j) = coord[j];
4943 (*Nodes) += displacements;
4945 MoveVertices(displacements);
4951 node_coord = (*Nodes);
4953 GetVertices(node_coord);
4959 (*Nodes) = node_coord;
4961 SetVertices(node_coord);
4966 if (own_nodes)
delete Nodes;
4969 own_nodes = (int)make_owner;
4980 mfem::Swap<GridFunction*>(Nodes, nodes);
4981 mfem::Swap<int>(own_nodes, own_nodes_);
4992 for (k = 0; k < spaceDim; k++)
4993 vertices[result](k) = vertices[indexes[0]](k);
4995 for (j = 1; j < n; j++)
4996 for (k = 0; k < spaceDim; k++)
4997 vertices[result](k) += vertices[indexes[j]](k);
4999 for (k = 0; k < spaceDim; k++)
5000 vertices[result](k) *= (1.0 / n);
5005 Nodes->FESpace()->UpdateAndInterpolate(Nodes);
5011 int i, j, *v, vv[2], attr, wtls = WantTwoLevelState;
5016 UseTwoLevelState(1);
5021 if (el_to_edge == NULL)
5023 el_to_edge =
new Table;
5024 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5027 int oedge = NumOfVertices;
5028 int oelem = oedge + NumOfEdges;
5030 DeleteCoarseTables();
5032 vertices.SetSize(oelem + NumOfElements);
5034 for (i = 0; i < NumOfElements; i++)
5036 v = elements[i]->GetVertices();
5038 AverageVertices(v, 4, oelem+i);
5040 e = el_to_edge->GetRow(i);
5042 vv[0] = v[0], vv[1] = v[1]; AverageVertices(vv, 2, oedge+e[0]);
5043 vv[0] = v[1], vv[1] = v[2]; AverageVertices(vv, 2, oedge+e[1]);
5044 vv[0] = v[2], vv[1] = v[3]; AverageVertices(vv, 2, oedge+e[2]);
5045 vv[0] = v[3], vv[1] = v[0]; AverageVertices(vv, 2, oedge+e[3]);
5048 elements.SetSize(4 * NumOfElements);
5049 for (i = 0; i < NumOfElements; i++)
5051 attr = elements[i]->GetAttribute();
5052 v = elements[i]->GetVertices();
5053 e = el_to_edge->GetRow(i);
5054 j = NumOfElements + 3 * i;
5056 elements[j+0] =
new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
5058 elements[j+1] =
new Quadrilateral(oelem+i, oedge+e[1], v[2],
5060 elements[j+2] =
new Quadrilateral(oedge+e[3], oelem+i, oedge+e[2],
5063 if (WantTwoLevelState)
5080 boundary.SetSize(2 * NumOfBdrElements);
5081 for (i = 0; i < NumOfBdrElements; i++)
5084 v = boundary[i]->GetVertices();
5085 j = NumOfBdrElements + i;
5087 boundary[j] =
new Segment(oedge+be_to_edge[i], v[1], attr);
5089 if (WantTwoLevelState)
5091 #ifdef MFEM_USE_MEMALLOC
5102 v[1] = oedge+be_to_edge[i];
5105 if (WantTwoLevelState)
5107 c_NumOfVertices = NumOfVertices;
5108 c_NumOfEdges = NumOfEdges;
5109 c_NumOfElements = NumOfElements;
5110 c_NumOfBdrElements = NumOfBdrElements;
5116 NumOfVertices = oelem + NumOfElements;
5117 NumOfElements = 4 * NumOfElements;
5118 NumOfBdrElements = 2 * NumOfBdrElements;
5121 if (WantTwoLevelState)
5123 f_NumOfVertices = NumOfVertices;
5124 f_NumOfElements = NumOfElements;
5125 f_NumOfBdrElements = NumOfBdrElements;
5128 if (el_to_edge != NULL)
5130 if (WantTwoLevelState)
5132 c_el_to_edge = el_to_edge;
5134 f_el_to_edge =
new Table;
5135 NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5136 el_to_edge = f_el_to_edge;
5137 f_NumOfEdges = NumOfEdges;
5140 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5147 UseTwoLevelState(wtls);
5151 CheckElementOrientation(
false);
5152 CheckBdrElementOrientation(
false);
5158 int i, wtls = WantTwoLevelState;
5165 UseTwoLevelState(1);
5170 if (el_to_edge == NULL)
5172 el_to_edge =
new Table;
5173 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5175 if (el_to_face == NULL)
5176 GetElementToFaceTable();
5178 int oedge = NumOfVertices;
5179 int oface = oedge + NumOfEdges;
5180 int oelem = oface + NumOfFaces;
5182 DeleteCoarseTables();
5184 vertices.SetSize(oelem + NumOfElements);
5185 for (i = 0; i < NumOfElements; i++)
5188 "Element is not a hex!");
5189 v = elements[i]->GetVertices();
5191 AverageVertices(v, 8, oelem+i);
5193 f = el_to_face->GetRow(i);
5195 for (
int j = 0; j < 6; j++)
5197 for (
int k = 0; k < 4; k++)
5198 vv[k] = v[hex_faces[j][k]];
5199 AverageVertices(vv, 4, oface+f[j]);
5202 e = el_to_edge->GetRow(i);
5204 for (
int j = 0; j < 12; j++)
5206 for (
int k = 0; k < 2; k++)
5208 AverageVertices(vv, 2, oedge+e[j]);
5213 elements.SetSize(8 * NumOfElements);
5214 for (i = 0; i < NumOfElements; i++)
5216 attr = elements[i]->GetAttribute();
5217 v = elements[i]->GetVertices();
5218 e = el_to_edge->GetRow(i);
5219 f = el_to_face->GetRow(i);
5220 j = NumOfElements + 7 * i;
5222 elements[j+0] =
new Hexahedron(oedge+e[0], v[1], oedge+e[1], oface+f[0],
5223 oface+f[1], oedge+e[9], oface+f[2],
5225 elements[j+1] =
new Hexahedron(oface+f[0], oedge+e[1], v[2], oedge+e[2],
5226 oelem+i, oface+f[2], oedge+e[10],
5228 elements[j+2] =
new Hexahedron(oedge+e[3], oface+f[0], oedge+e[2], v[3],
5229 oface+f[4], oelem+i, oface+f[3],
5231 elements[j+3] =
new Hexahedron(oedge+e[8], oface+f[1], oelem+i,
5232 oface+f[4], v[4], oedge+e[4], oface+f[5],
5234 elements[j+4] =
new Hexahedron(oface+f[1], oedge+e[9], oface+f[2],
5235 oelem+i, oedge+e[4], v[5], oedge+e[5],
5237 elements[j+5] =
new Hexahedron(oelem+i, oface+f[2], oedge+e[10],
5238 oface+f[3], oface+f[5], oedge+e[5], v[6],
5240 elements[j+6] =
new Hexahedron(oface+f[4], oelem+i, oface+f[3],
5241 oedge+e[11], oedge+e[7], oface+f[5],
5242 oedge+e[6], v[7], attr);
5244 if (WantTwoLevelState)
5250 for (k = 0; k < 7; k++)
5251 oe->
Child[k] = j + k;
5264 boundary.SetSize(4 * NumOfBdrElements);
5265 for (i = 0; i < NumOfBdrElements; i++)
5268 "boundary Element is not a quad!");
5269 attr = boundary[i]->GetAttribute();
5270 v = boundary[i]->GetVertices();
5271 e = bel_to_edge->GetRow(i);
5272 f = & be_to_face[i];
5273 j = NumOfBdrElements + 3 * i;
5275 boundary[j+0] =
new Quadrilateral(oedge+e[0], v[1], oedge+e[1],
5277 boundary[j+1] =
new Quadrilateral(oface+f[0], oedge+e[1], v[2],
5279 boundary[j+2] =
new Quadrilateral(oedge+e[3], oface+f[0], oedge+e[2],
5282 if (WantTwoLevelState)
5299 if (WantTwoLevelState)
5301 c_NumOfVertices = NumOfVertices;
5302 c_NumOfEdges = NumOfEdges;
5303 c_NumOfFaces = NumOfFaces;
5304 c_NumOfElements = NumOfElements;
5305 c_NumOfBdrElements = NumOfBdrElements;
5311 NumOfVertices = oelem + NumOfElements;
5312 NumOfElements = 8 * NumOfElements;
5313 NumOfBdrElements = 4 * NumOfBdrElements;
5315 if (WantTwoLevelState)
5317 f_NumOfVertices = NumOfVertices;
5318 f_NumOfElements = NumOfElements;
5319 f_NumOfBdrElements = NumOfBdrElements;
5322 if (el_to_face != NULL)
5324 if (WantTwoLevelState)
5326 c_el_to_face = el_to_face;
5330 GetElementToFaceTable();
5332 if (WantTwoLevelState)
5334 f_el_to_face = el_to_face;
5335 f_NumOfFaces = NumOfFaces;
5340 CheckBdrElementOrientation(
false);
5343 if (el_to_edge != NULL)
5345 if (WantTwoLevelState)
5347 c_el_to_edge = el_to_edge;
5348 f_el_to_edge =
new Table;
5349 c_bel_to_edge = bel_to_edge;
5351 NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5352 el_to_edge = f_el_to_edge;
5353 f_bel_to_edge = bel_to_edge;
5354 f_NumOfEdges = NumOfEdges;
5357 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5363 UseTwoLevelState(wtls);
5373 int i, j, ind, nedges, wtls = WantTwoLevelState;
5378 mfem_error(
"Local and nonconforming refinements cannot be mixed.");
5383 UseTwoLevelState(1);
5387 DeleteCoarseTables();
5391 if (WantTwoLevelState)
5393 c_NumOfVertices = NumOfVertices;
5394 c_NumOfElements = NumOfElements;
5395 c_NumOfBdrElements = NumOfBdrElements;
5398 int cne = NumOfElements, cnv = NumOfVertices;
5399 NumOfVertices += marked_el.
Size();
5400 NumOfElements += marked_el.
Size();
5401 vertices.SetSize(NumOfVertices);
5402 elements.SetSize(NumOfElements);
5403 for (j = 0; j < marked_el.
Size(); j++)
5408 int new_v = cnv + j, new_e = cne + j;
5409 AverageVertices(vert, 2, new_v);
5410 elements[new_e] =
new Segment(new_v, vert[1], attr);
5411 if (WantTwoLevelState)
5413 #ifdef MFEM_USE_MEMALLOC
5428 if (WantTwoLevelState)
5430 f_NumOfVertices = NumOfVertices;
5431 f_NumOfElements = NumOfElements;
5432 f_NumOfBdrElements = NumOfBdrElements;
5442 if (WantTwoLevelState)
5444 c_NumOfVertices = NumOfVertices;
5445 c_NumOfEdges = NumOfEdges;
5446 c_NumOfElements = NumOfElements;
5447 c_NumOfBdrElements = NumOfBdrElements;
5451 DSTable v_to_v(NumOfVertices);
5452 GetVertexToVertexTable(v_to_v);
5456 int *edge1 =
new int[nedges];
5457 int *edge2 =
new int[nedges];
5458 int *middle =
new int[nedges];
5460 for (i = 0; i < nedges; i++)
5461 edge1[i] = edge2[i] = middle[i] = -1;
5463 for (i = 0; i < NumOfElements; i++)
5465 elements[i]->GetVertices(v);
5466 for (j = 1; j < v.
Size(); j++)
5468 ind = v_to_v(v[j-1], v[j]);
5469 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
5471 ind = v_to_v(v[0], v[v.
Size()-1]);
5472 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
5476 for (i = 0; i < marked_el.
Size(); i++)
5477 RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
5480 int need_refinement;
5483 need_refinement = 0;
5484 for (i = 0; i < nedges; i++)
5485 if (middle[i] != -1 && edge1[i] != -1)
5487 need_refinement = 1;
5488 GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
5491 while (need_refinement == 1);
5494 int v1[2], v2[2], bisect, temp;
5495 temp = NumOfBdrElements;
5496 for (i = 0; i < temp; i++)
5498 boundary[i]->GetVertices(v);
5499 bisect = v_to_v(v[0], v[1]);
5500 if (middle[bisect] != -1)
5504 v1[0] = v[0]; v1[1] = middle[bisect];
5505 v2[0] = middle[bisect]; v2[1] = v[1];
5507 if (WantTwoLevelState)
5509 boundary.
Append(
new Segment(v2, boundary[i]->GetAttribute()));
5510 #ifdef MFEM_USE_MEMALLOC
5517 new Segment(v1, boundary[i]->GetAttribute());
5525 boundary.Append(
new Segment(v2, boundary[i]->GetAttribute()));
5529 mfem_error(
"Only bisection of segment is implemented"
5533 NumOfBdrElements = boundary.Size();
5540 if (WantTwoLevelState)
5542 f_NumOfVertices = NumOfVertices;
5543 f_NumOfElements = NumOfElements;
5544 f_NumOfBdrElements = NumOfBdrElements;
5549 if (el_to_edge != NULL)
5551 if (WantTwoLevelState)
5553 c_el_to_edge = el_to_edge;
5555 f_el_to_edge =
new Table;
5556 NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5557 el_to_edge = f_el_to_edge;
5558 f_NumOfEdges = NumOfEdges;
5561 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5568 if (WantTwoLevelState)
5570 c_NumOfVertices = NumOfVertices;
5571 c_NumOfEdges = NumOfEdges;
5572 c_NumOfFaces = NumOfFaces;
5573 c_NumOfElements = NumOfElements;
5574 c_NumOfBdrElements = NumOfBdrElements;
5578 DSTable v_to_v(NumOfVertices);
5579 GetVertexToVertexTable(v_to_v);
5583 int *middle =
new int[nedges];
5585 for (i = 0; i < nedges; i++)
5593 for (i = 0; i < marked_el.
Size(); i++)
5594 Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5597 for (i = 0; i < marked_el.
Size(); i++)
5599 Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5601 Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
5602 Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5606 for (i = 0; i < marked_el.
Size(); i++)
5608 Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5610 ii = NumOfElements - 1;
5611 Bisection(ii, v_to_v, NULL, NULL, middle);
5612 Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
5613 Bisection(ii, v_to_v, NULL, NULL, middle);
5615 Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5616 Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
5617 Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
5622 if (WantTwoLevelState)
5629 int need_refinement;
5634 need_refinement = 0;
5637 for (i = 0; i < NumOfElements; i++)
5641 if (elements[i]->NeedRefinement(v_to_v, middle))
5643 need_refinement = 1;
5644 Bisection(i, v_to_v, NULL, NULL, middle);
5648 while (need_refinement == 1);
5655 need_refinement = 0;
5656 for (i = 0; i < NumOfBdrElements; i++)
5657 if (boundary[i]->NeedRefinement(v_to_v, middle))
5659 need_refinement = 1;
5660 Bisection(i, v_to_v, middle);
5663 while (need_refinement == 1);
5666 int refinement_edges[2], type, flag;
5667 for (i = 0; i < NumOfElements; i++)
5672 ((
Tetrahedron *)El)->ParseRefinementFlag(refinement_edges, type,
5675 ((
Tetrahedron *)El)->CreateRefinementFlag(refinement_edges,
5680 NumOfBdrElements = boundary.Size();
5685 if (el_to_edge != NULL)
5687 if (WantTwoLevelState)
5689 c_el_to_edge = el_to_edge;
5690 f_el_to_edge =
new Table;
5691 c_bel_to_edge = bel_to_edge;
5693 NumOfEdges = GetElementToEdgeTable(*f_el_to_edge, be_to_edge);
5694 el_to_edge = f_el_to_edge;
5695 f_bel_to_edge = bel_to_edge;
5698 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5700 if (el_to_face != NULL)
5702 if (WantTwoLevelState)
5704 c_el_to_face = el_to_face;
5708 GetElementToFaceTable();
5710 if (WantTwoLevelState)
5712 f_el_to_face = el_to_face;
5716 if (WantTwoLevelState)
5718 f_NumOfVertices = NumOfVertices;
5719 f_NumOfEdges = NumOfEdges;
5720 f_NumOfFaces = NumOfFaces;
5721 f_NumOfElements = NumOfElements;
5722 f_NumOfBdrElements = NumOfBdrElements;
5730 UseTwoLevelState(wtls);
5734 CheckElementOrientation(
false);
5743 mfem_error(
"Mesh::NonconformingRefinement: NURBS meshes are not supported."
5744 " Project the NURBS to Nodes first.");
5747 int wtls = WantTwoLevelState;
5751 UseTwoLevelState(1);
5755 DeleteCoarseTables();
5760 ncmesh =
new NCMesh(
this);
5763 if (WantTwoLevelState)
5764 ncmesh->MarkCoarseLevel();
5767 ncmesh->Refine(refinements);
5770 ncmesh->LimitNCLevel(nc_limit);
5775 ncmesh->SetEdgeIndicesFromMesh(mesh2);
5776 ncmesh->SetFaceIndicesFromMesh(mesh2);
5780 Swap(*mesh2,
false);
5783 if (WantTwoLevelState)
5785 nc_coarse_level = mesh2;
5786 State = TWO_LEVEL_FINE;
5794 UseTwoLevelState(wtls);
5807 NumOfVertices = vertices.Size();
5808 NumOfElements = elements.Size();
5809 NumOfBdrElements = boundary.Size();
5814 NumOfEdges = NumOfFaces = 0;
5818 GetElementToFaceTable();
5821 CheckBdrElementOrientation(
false);
5825 el_to_edge =
new Table;
5826 NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
5827 c_el_to_edge = NULL;
5877 NURBSUniformRefinement();
5878 else if (meshgen == 1)
5882 for (
int i = 0; i < elem_to_refine.
Size(); i++)
5883 elem_to_refine[i] = i;
5884 LocalRefinement(elem_to_refine);
5887 QuadUniformRefinement();
5889 HexUniformRefinement();
5897 if (nonconforming < 0)
5900 int type = elements[0]->GetType();
5910 NonconformingRefinement(refinements, nc_limit);
5915 for (
int i = 0; i < refinements.
Size(); i++)
5916 el_to_refine.
Append(refinements[i].index);
5919 int type, rt = (refinements.
Size() ? refinements[0].ref_type : 7);
5920 if (rt == 1 || rt == 2 || rt == 4)
5922 else if (rt == 3 || rt == 5 || rt == 6)
5928 LocalRefinement(el_to_refine, type);
5936 for (
int i = 0; i < el_to_refine.
Size(); i++)
5939 GeneralRefinement(refinements, nonconforming, nc_limit);
5943 int *edge1,
int *edge2,
int *middle)
5946 int v[2][4], v_new, bisect, t;
5951 if (WantTwoLevelState)
5972 bisect = v_to_v(vert[0], vert[1]);
5975 mfem_error(
"Mesh::Bisection(...) of triangle! #1");
5977 if (middle[bisect] == -1)
5979 v_new = NumOfVertices++;
5980 for (
int d = 0; d < spaceDim; d++)
5981 V(d) = 0.5 * (vertices[vert[0]](d) + vertices[vert[1]](d));
5986 if (edge1[bisect] == i)
5987 edge1[bisect] = edge2[bisect];
5989 middle[bisect] = v_new;
5993 v_new = middle[bisect];
6001 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
6002 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
6004 if (WantTwoLevelState)
6006 #ifdef MFEM_USE_MEMALLOC
6023 bisect = v_to_v(v[1][0], v[1][1]);
6026 mfem_error(
"Mesh::Bisection(...) of triangle! #2");
6028 if (edge1[bisect] == i)
6029 edge1[bisect] = NumOfElements;
6030 else if (edge2[bisect] == i)
6031 edge2[bisect] = NumOfElements;
6037 int j, type, new_type, old_redges[2], new_redges[2][2], flag;
6041 mfem_error(
"Mesh::Bisection : TETRAHEDRON element is not marked for "
6047 bisect = v_to_v(vert[0], vert[1]);
6051 cerr <<
"Error in Bisection(...) of tetrahedron!" << endl
6052 <<
" redge[0] = " << old_redges[0]
6053 <<
" redge[1] = " << old_redges[1]
6054 <<
" type = " << type
6055 <<
" flag = " << flag << endl;
6059 if (middle[bisect] == -1)
6061 v_new = NumOfVertices++;
6062 for (j = 0; j < 3; j++)
6063 V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
6066 middle[bisect] = v_new;
6069 v_new = middle[bisect];
6077 new_redges[0][0] = 2;
6078 new_redges[0][1] = 1;
6079 new_redges[1][0] = 2;
6080 new_redges[1][1] = 1;
6081 switch (old_redges[0])
6084 v[0][0] = vert[0]; v[0][1] = vert[2]; v[0][2] = vert[3];
6088 v[0][0] = vert[3]; v[0][1] = vert[0]; v[0][2] = vert[2];
6091 v[0][0] = vert[2]; v[0][1] = vert[3]; v[0][2] = vert[0];
6093 switch (old_redges[1])
6096 v[1][0] = vert[2]; v[1][1] = vert[1]; v[1][2] = vert[3];
6100 v[1][0] = vert[1]; v[1][1] = vert[3]; v[1][2] = vert[2];
6103 v[1][0] = vert[3]; v[1][1] = vert[2]; v[1][2] = vert[1];
6107 if (WantTwoLevelState)
6109 #ifdef MFEM_USE_MEMALLOC
6112 tet = TetMemory.Alloc();
6128 #ifdef MFEM_USE_MEMALLOC
6132 elements.Append(tet2);
6151 CreateRefinementFlag(new_redges[1], new_type, flag+1);
6156 mfem_error(
"Bisection for now works only for triangles & tetrahedra.");
6162 int v[2][3], v_new, bisect, t;
6166 if (WantTwoLevelState)
6187 bisect = v_to_v(vert[0], vert[1]);
6190 mfem_error(
"Mesh::Bisection(...) of boundary triangle! #1");
6192 v_new = middle[bisect];
6195 mfem_error(
"Mesh::Bisection(...) of boundary triangle! #2");
6200 v[0][0] = vert[2]; v[0][1] = vert[0]; v[0][2] = v_new;
6201 v[1][0] = vert[1]; v[1][1] = vert[2]; v[1][2] = v_new;
6202 if (WantTwoLevelState)
6204 #ifdef MFEM_USE_MEMALLOC
6222 mfem_error(
"Bisection of boundary elements works only for triangles!");
6226 int *edge1,
int *edge2,
int *middle)
6229 int j, v1[3], v2[3], v3[3], v4[3], v_new[3], bisect[3];
6234 elements[i]->GetVertices(v);
6237 bisect[0] = v_to_v(v[0],v[1]);
6238 bisect[1] = v_to_v(v[1],v[2]);
6239 bisect[2] = v_to_v(v[0],v[2]);
6241 if (bisect[0] < 0 || bisect[1] < 0 || bisect[2] < 0)
6242 mfem_error(
"Mesh::UniformRefinement(...): ERROR");
6245 for (j = 0; j < 3; j++)
6246 if (middle[bisect[j]] == -1)
6248 v_new[j] = NumOfVertices++;
6249 for (
int d = 0; d < spaceDim; d++)
6250 V(d) = (vertices[v[j]](d) + vertices[v[(j+1)%3]](d))/2.;
6255 if (edge1[bisect[j]] == i)
6256 edge1[bisect[j]] = edge2[bisect[j]];
6258 middle[bisect[j]] = v_new[j];
6262 v_new[j] = middle[bisect[j]];
6265 edge1[bisect[j]] = -1;
6270 v1[0] = v[0]; v1[1] = v_new[0]; v1[2] = v_new[2];
6271 v2[0] = v_new[0]; v2[1] = v[1]; v2[2] = v_new[1];
6272 v3[0] = v_new[2]; v3[1] = v_new[1]; v3[2] = v[2];
6273 v4[0] = v_new[1]; v4[1] = v_new[2]; v4[2] = v_new[0];
6276 elements.Append(
new Triangle(v2, elements[i]->GetAttribute()));
6277 elements.Append(
new Triangle(v3, elements[i]->GetAttribute()));
6278 if (WantTwoLevelState)
6282 aux->
Child2 = NumOfElements;
6283 aux->
Child3 = NumOfElements+1;
6284 aux->
Child4 = NumOfElements+2;
6295 mfem_error(
"Uniform refinement for now works only for triangles.");
6304 delete nc_coarse_level;
6305 nc_coarse_level = NULL;
6306 ncmesh->ClearCoarseLevel();
6311 this->
Swap(*nc_coarse_level,
false);
6314 else if (State != s)
6326 for (i = 0; i < f_NumOfElements; )
6328 t = elements[i]->GetType();
6342 for (i = 0; i < f_NumOfBdrElements; )
6344 t = boundary[i]->GetType();
6358 if (el_to_edge != NULL)
6360 delete c_el_to_edge;
6361 el_to_edge = f_el_to_edge;
6366 fc_be_to_edge.DeleteAll();
6370 delete c_bel_to_edge;
6371 bel_to_edge = f_bel_to_edge;
6374 if (el_to_face != NULL)
6376 delete c_el_to_face;
6377 el_to_face = f_el_to_face;
6380 fc_faces_info.DeleteAll();
6383 NumOfVertices = f_NumOfVertices;
6384 NumOfEdges = f_NumOfEdges;
6386 NumOfFaces = f_NumOfFaces;
6387 NumOfElements = f_NumOfElements;
6388 NumOfBdrElements = f_NumOfBdrElements;
6394 if (el_to_edge != NULL)
6396 el_to_edge = f_el_to_edge;
6400 bel_to_edge = f_bel_to_edge;
6402 if (el_to_face != NULL)
6404 el_to_face = f_el_to_face;
6407 NumOfVertices = f_NumOfVertices;
6408 NumOfEdges = f_NumOfEdges;
6410 NumOfFaces = f_NumOfFaces;
6411 NumOfElements = f_NumOfElements;
6412 NumOfBdrElements = f_NumOfBdrElements;
6418 if (el_to_edge != NULL)
6420 el_to_edge = c_el_to_edge;
6424 bel_to_edge = c_bel_to_edge;
6426 if (el_to_face != NULL)
6428 el_to_face = c_el_to_face;
6431 NumOfVertices = c_NumOfVertices;
6432 NumOfEdges = c_NumOfEdges;
6434 NumOfFaces = c_NumOfFaces;
6435 NumOfElements = c_NumOfElements;
6436 NumOfBdrElements = c_NumOfBdrElements;
6440 else if (State != s)
6447 int t = elements[i]->GetType();
6506 int L, R, n, s, lb, rb;
6510 R = GetBisectionHierarchy(elements[n]);
6519 n |= ((L & 1) << s);
6527 n |= ((R & 1) << s);
6533 lb = 2 * nlb; rb = 2 * nrb;
6535 while (lb > 0 || rb > 0);
6552 t = elements[i]->GetType();
6564 type = GetBisectionHierarchy(elements[i]);
6572 int redges[2], type, flag;
6592 type |= ( GetBisectionHierarchy(E) << 3 );
6593 if (type < 8) type = 0;
6603 int t = elements[i]->GetType();
6625 case 1:
return aux->
Child2;
6626 case 2:
return aux->
Child3;
6627 case 3:
return aux->
Child4;
6633 case 0:
return aux->
Child2;
6634 case 1:
return aux->
Child3;
6635 case 2:
return aux->
Child4;
6687 if (j == 0)
return i;
6705 np[0] = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6706 np[1] = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6707 pointmat(0,1) = pointmat(0,0); pointmat(1,1) = pointmat(1,0);
6708 pointmat(0,0) = pointmat(0,2); pointmat(1,0) = pointmat(1,2);
6709 pointmat(0,2) = np[0]; pointmat(1,2) = np[1];
6714 np[0] = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6715 np[1] = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6716 pointmat(0,0) = pointmat(0,1); pointmat(1,0) = pointmat(1,1);
6717 pointmat(0,1) = pointmat(0,2); pointmat(1,1) = pointmat(1,2);
6718 pointmat(0,2) = np[0]; pointmat(1,2) = np[1];
6724 int i, j, redges[2], type, flag, ind[4];
6732 pointmat(0,1) = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6733 pointmat(1,1) = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6734 pointmat(2,1) = 0.5 * ( pointmat(2,0) + pointmat(2,1) );
6738 case 2: ind[0] = 0; ind[1] = 3; ind[2] = 1; ind[3] = 2;
break;
6739 case 3: ind[0] = 1; ind[1] = 3; ind[2] = 2; ind[3] = 0;
break;
6741 default: ind[0] = 2; ind[1] = 3; ind[2] = 0; ind[3] = 1;
6747 pointmat(0,0) = 0.5 * ( pointmat(0,0) + pointmat(0,1) );
6748 pointmat(1,0) = 0.5 * ( pointmat(1,0) + pointmat(1,1) );
6749 pointmat(2,0) = 0.5 * ( pointmat(2,0) + pointmat(2,1) );
6753 case 1: ind[0] = 3; ind[1] = 1; ind[2] = 0; ind[3] = 2;
break;
6754 case 4: ind[0] = 3; ind[1] = 0; ind[2] = 2; ind[3] = 1;
break;
6756 default: ind[0] = 3; ind[1] = 2; ind[2] = 1; ind[3] = 0;
6760 for (i = 0; i < 3; i++)
6762 for (j = 0; j < 4; j++)
6763 t[j] = pointmat(i,j);
6764 for (j = 0; j < 4; j++)
6765 pointmat(i,ind[j]) = t[j];
6781 return ((k << 1)+1) << l;
6810 Transformation.Attribute = 0;
6811 Transformation.ElementNo = 0;
6817 case 0: pm(0,0) = 0.0; pm(0,1) = 0.5;
break;
6818 default: pm(0,0) = 0.5; pm(0,1) = 1.0;
break;
6823 pm(0,0) = 0.0; pm(0,1) = 1.0;
6829 Transformation.Attribute = 0;
6830 Transformation.ElementNo = 0;
6831 t = elements[i]->GetType();
6843 pm(0,0) = 0.0; pm(1,0) = 0.0;
6844 pm(0,1) = 0.5; pm(1,1) = 0.0;
6845 pm(0,2) = 0.5; pm(1,2) = 0.5;
6846 pm(0,3) = 0.0; pm(1,3) = 0.5;
6849 pm(0,0) = 0.5; pm(1,0) = 0.0;
6850 pm(0,1) = 1.0; pm(1,1) = 0.0;
6851 pm(0,2) = 1.0; pm(1,2) = 0.5;
6852 pm(0,3) = 0.5; pm(1,3) = 0.5;
6855 pm(0,0) = 0.5; pm(1,0) = 0.5;
6856 pm(0,1) = 1.0; pm(1,1) = 0.5;
6857 pm(0,2) = 1.0; pm(1,2) = 1.0;
6858 pm(0,3) = 0.5; pm(1,3) = 1.0;
6861 pm(0,0) = 0.0; pm(1,0) = 0.5;
6862 pm(0,1) = 0.5; pm(1,1) = 0.5;
6863 pm(0,2) = 0.5; pm(1,2) = 1.0;
6864 pm(0,3) = 0.0; pm(1,3) = 1.0;
6878 pm(0,0) = 0.0; pm(0,1) = 0.5; pm(0,2) = 0.0;
6879 pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 0.5;
6882 pm(0,0) = 0.5; pm(0,1) = 1.0; pm(0,2) = 0.5;
6883 pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 0.5;
6886 pm(0,0) = 0.0; pm(0,1) = 0.5; pm(0,2) = 0.0;
6887 pm(1,0) = 0.5; pm(1,1) = 0.5; pm(1,2) = 1.0;
6890 pm(0,0) = 0.5; pm(0,1) = 0.0; pm(0,2) = 0.5;
6891 pm(1,0) = 0.5; pm(1,1) = 0.5; pm(1,2) = 0.0;
6909 pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0;
6910 pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0;
6912 path = GetFineElemPath(i, j);
6929 pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0;
6930 pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0;
6932 return &Transformation;
6942 Transformation.Attribute = 0;
6943 Transformation.ElementNo = 0;
6945 if (j < 4) dz = 0.0;
6948 if (jj < 2) dy = 0.0;
6950 if (jj == 0 || jj == 3) dx = 0.0;
6952 pm(0,0) = dx; pm(1,0) = dy; pm(2,0) = dz;
6953 pm(0,1) = 0.5 + dx; pm(1,1) = dy; pm(2,1) = dz;
6954 pm(0,2) = 0.5 + dx; pm(1,2) = 0.5 + dy; pm(2,2) = dz;
6955 pm(0,3) = dx; pm(1,3) = 0.5 + dy; pm(2,3) = dz;
6956 pm(0,4) = dx; pm(1,4) = dy; pm(2,4) = 0.5 + dz;
6957 pm(0,5) = 0.5 + dx; pm(1,5) = dy; pm(2,5) = 0.5 + dz;
6958 pm(0,6) = 0.5 + dx; pm(1,6) = 0.5 + dy; pm(2,6) = 0.5 + dz;
6959 pm(0,7) = dx; pm(1,7) = 0.5 + dy; pm(2,7) = 0.5 + dz;
6960 return &Transformation;
6966 Transformation.Attribute = 0;
6967 Transformation.ElementNo = 0;
6972 pm(0,0) = 0.0; pm(0,1) = 1.0; pm(0,2) = 0.0; pm(0,3) = 0.0;
6973 pm(1,0) = 0.0; pm(1,1) = 0.0; pm(1,2) = 1.0; pm(1,3) = 0.0;
6974 pm(2,0) = 0.0; pm(2,1) = 0.0; pm(2,2) = 0.0; pm(2,3) = 1.0;
6976 path = GetFineElemPath(i, j);
6989 return &Transformation;
6994 MFEM_ASSERT(Dim==spaceDim,
"2D Manifold meshes not supported");
7002 out <<
"areamesh2\n\n";
7004 out <<
"curved_areamesh2\n\n";
7007 out << NumOfBdrElements <<
'\n';
7008 for (i = 0; i < NumOfBdrElements; i++)
7010 boundary[i]->GetVertices(v);
7012 out << boundary[i]->GetAttribute();
7013 for (j = 0; j < v.
Size(); j++)
7014 out <<
' ' << v[j] + 1;
7019 out << NumOfElements <<
'\n';
7020 for (i = 0; i < NumOfElements; i++)
7022 elements[i]->GetVertices(v);
7024 out << elements[i]->GetAttribute() <<
' ' << v.
Size();
7025 for (j = 0; j < v.
Size(); j++)
7026 out <<
' ' << v[j] + 1;
7033 out << NumOfVertices <<
'\n';
7034 for (i = 0; i < NumOfVertices; i++)
7036 out << vertices[i](0);
7037 for (j = 1; j < Dim; j++)
7038 out <<
' ' << vertices[i](j);
7044 out << NumOfVertices <<
'\n';
7052 mfem_error(
"Mesh::PrintXG(...) : Curved mesh in 3D");
7060 out <<
"NETGEN_Neutral_Format\n";
7062 out << NumOfVertices <<
'\n';
7063 for (i = 0; i < NumOfVertices; i++)
7065 for (j = 0; j < Dim; j++)
7066 out <<
' ' << vertices[i](j);
7071 out << NumOfElements <<
'\n';
7072 for (i = 0; i < NumOfElements; i++)
7074 nv = elements[i]->GetNVertices();
7075 ind = elements[i]->GetVertices();
7076 out << elements[i]->GetAttribute();
7077 for (j = 0; j < nv; j++)
7078 out <<
' ' << ind[j]+1;
7083 out << NumOfBdrElements <<
'\n';
7084 for (i = 0; i < NumOfBdrElements; i++)
7086 nv = boundary[i]->GetNVertices();
7087 ind = boundary[i]->GetVertices();
7088 out << boundary[i]->GetAttribute();
7089 for (j = 0; j < nv; j++)
7090 out <<
' ' << ind[j]+1;
7094 else if (meshgen == 2)
7100 <<
"1 " << NumOfVertices <<
" " << NumOfElements
7101 <<
" 0 0 0 0 0 0 0\n"
7102 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
7103 <<
"0 0 " << NumOfBdrElements <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
7104 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
7105 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
7107 for (i = 0; i < NumOfVertices; i++)
7108 out << i+1 <<
" 0.0 " << vertices[i](0) <<
' ' << vertices[i](1)
7109 <<
' ' << vertices[i](2) <<
" 0.0\n";
7111 for (i = 0; i < NumOfElements; i++)
7113 nv = elements[i]->GetNVertices();
7114 ind = elements[i]->GetVertices();
7115 out << i+1 <<
' ' << elements[i]->GetAttribute();
7116 for (j = 0; j < nv; j++)
7117 out <<
' ' << ind[j]+1;
7121 for (i = 0; i < NumOfBdrElements; i++)
7123 nv = boundary[i]->GetNVertices();
7124 ind = boundary[i]->GetVertices();
7125 out << boundary[i]->GetAttribute();
7126 for (j = 0; j < nv; j++)
7127 out <<
' ' << ind[j]+1;
7128 out <<
" 1.0 1.0 1.0 1.0\n";
7142 NURBSext->Print(out);
7148 out <<
"MFEM mesh v1.0\n";
7152 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
7157 "# TETRAHEDRON = 4\n"
7161 out <<
"\ndimension\n" << Dim
7162 <<
"\n\nelements\n" << NumOfElements <<
'\n';
7163 for (i = 0; i < NumOfElements; i++)
7164 PrintElement(elements[i], out);
7166 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
7167 for (i = 0; i < NumOfBdrElements; i++)
7168 PrintElement(boundary[i], out);
7170 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
7173 out << spaceDim <<
'\n';
7174 for (i = 0; i < NumOfVertices; i++)
7176 out << vertices[i](0);
7177 for (j = 1; j < spaceDim; j++)
7178 out <<
' ' << vertices[i](j);
7194 out <<
"MFEM NURBS mesh v1.0\n";
7198 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
7204 out <<
"\ndimension\n" << Dim
7205 <<
"\n\nelements\n" << NumOfElements <<
'\n';
7206 for (i = 0; i < NumOfElements; i++)
7207 PrintElement(elements[i], out);
7209 out <<
"\nboundary\n" << NumOfBdrElements <<
'\n';
7210 for (i = 0; i < NumOfBdrElements; i++)
7211 PrintElement(boundary[i], out);
7213 out <<
"\nedges\n" << NumOfEdges <<
'\n';
7214 for (i = 0; i < NumOfEdges; i++)
7216 edge_vertex->GetRow(i, vert);
7220 out << ki <<
' ' << vert[0] <<
' ' << vert[1] <<
'\n';
7222 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
7228 "# vtk DataFile Version 3.0\n"
7229 "Generated by MFEM\n"
7231 "DATASET UNSTRUCTURED_GRID\n";
7235 out <<
"POINTS " << NumOfVertices <<
" double\n";
7236 for (
int i = 0; i < NumOfVertices; i++)
7238 out << vertices[i](0);
7240 for (j = 1; j < spaceDim; j++)
7241 out <<
' ' << vertices[i](j);
7250 out <<
"POINTS " << Nodes->FESpace()->GetNDofs() <<
" double\n";
7251 for (
int i = 0; i < Nodes->FESpace()->GetNDofs(); i++)
7255 Nodes->FESpace()->DofsToVDofs(vdofs);
7256 out << (*Nodes)(vdofs[0]);
7258 for (j = 1; j < spaceDim; j++)
7259 out <<
' ' << (*Nodes)(vdofs[j]);
7270 for (
int i = 0; i < NumOfElements; i++)
7271 size += elements[i]->GetNVertices() + 1;
7272 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
7273 for (
int i = 0; i < NumOfElements; i++)
7275 const int *v = elements[i]->GetVertices();
7276 const int nv = elements[i]->GetNVertices();
7278 for (
int j = 0; j < nv; j++)
7288 for (
int i = 0; i < NumOfElements; i++)
7290 Nodes->FESpace()->GetElementDofs(i, dofs);
7291 size += dofs.
Size() + 1;
7293 out <<
"CELLS " << NumOfElements <<
' ' << size <<
'\n';
7294 const char *fec_name = Nodes->FESpace()->FEColl()->Name();
7295 if (!strcmp(fec_name,
"Linear") ||
7296 !strcmp(fec_name,
"H1_2D_P1") ||
7297 !strcmp(fec_name,
"H1_3D_P1"))
7299 else if (!strcmp(fec_name,
"Quadratic") ||
7300 !strcmp(fec_name,
"H1_2D_P2") ||
7301 !strcmp(fec_name,
"H1_3D_P2"))
7305 cerr <<
"Mesh::PrintVTK : can not save '"
7306 << fec_name <<
"' elements!" << endl;
7309 for (
int i = 0; i < NumOfElements; i++)
7311 Nodes->FESpace()->GetElementDofs(i, dofs);
7315 for (
int j = 0; j < dofs.
Size(); j++)
7316 out <<
' ' << dofs[j];
7318 else if (order == 2)
7320 const int *vtk_mfem;
7321 switch (elements[i]->GetGeometryType())
7325 vtk_mfem = vtk_quadratic_hex;
break;
7327 vtk_mfem = vtk_quadratic_tet;
break;
7330 vtk_mfem = vtk_quadratic_hex;
break;
7332 for (
int j = 0; j < dofs.
Size(); j++)
7333 out <<
' ' << dofs[vtk_mfem[j]];
7339 out <<
"CELL_TYPES " << NumOfElements <<
'\n';
7340 for (
int i = 0; i < NumOfElements; i++)
7342 int vtk_cell_type = 5;
7345 switch (elements[i]->GetGeometryType())
7353 else if (order == 2)
7355 switch (elements[i]->GetGeometryType())
7364 out << vtk_cell_type <<
'\n';
7368 out <<
"CELL_DATA " << NumOfElements <<
'\n'
7369 <<
"SCALARS material int\n"
7370 <<
"LOOKUP_TABLE default\n";
7371 for (
int i = 0; i < NumOfElements; i++)
7373 out << elements[i]->GetAttribute() <<
'\n';
7384 "# vtk DataFile Version 3.0\n"
7385 "Generated by MFEM\n"
7387 "DATASET UNSTRUCTURED_GRID\n";
7392 out <<
"FIELD FieldData 1" << endl
7393 <<
"MaterialIds " << 1 <<
" " << attributes.
Size() <<
" int" << endl;
7394 for (
int i = 0; i < attributes.Size(); i++)
7395 out << attributes[i] <<
" ";
7401 for (
int i = 0; i < GetNE(); i++)
7403 int geom = GetElementBaseGeometry(i);
7410 out <<
"POINTS " << np <<
" double\n";
7412 for (
int i = 0; i < GetNE(); i++)
7415 GetElementBaseGeometry(i), ref, 1);
7417 GetElementTransformation(i)->Transform(RefG->
RefPts, pmat);
7419 for (
int j = 0; j < pmat.
Width(); j++)
7421 out << pmat(0, j) <<
' ';
7424 out << pmat(1, j) <<
' ';
7431 out << 0.0 <<
' ' << 0.0;
7437 out <<
"CELLS " << nc <<
' ' << size <<
'\n';
7439 for (
int i = 0; i < GetNE(); i++)
7441 int geom = GetElementBaseGeometry(i);
7446 for (
int j = 0; j < RG.
Size(); )
7449 for (
int k = 0; k < nv; k++, j++)
7450 out <<
' ' << np + RG[j];
7455 out <<
"CELL_TYPES " << nc <<
'\n';
7456 for (
int i = 0; i < GetNE(); i++)
7458 int geom = GetElementBaseGeometry(i);
7462 int vtk_cell_type = 5;
7473 for (
int j = 0; j < RG.
Size(); j += nv)
7475 out << vtk_cell_type <<
'\n';
7479 out <<
"CELL_DATA " << nc <<
'\n'
7480 <<
"SCALARS material int\n"
7481 <<
"LOOKUP_TABLE default\n";
7482 for (
int i = 0; i < GetNE(); i++)
7484 int geom = GetElementBaseGeometry(i);
7487 int attr = GetAttribute(i);
7490 out << attr <<
'\n';
7495 srand((
unsigned)time(0));
7496 double a = double(rand()) / (double(RAND_MAX) + 1.);
7497 int el0 = (int)floor(a * GetNE());
7498 GetElementColoring(coloring, el0);
7499 out <<
"SCALARS element_coloring int\n"
7500 <<
"LOOKUP_TABLE default\n";
7501 for (
int i = 0; i < GetNE(); i++)
7503 int geom = GetElementBaseGeometry(i);
7508 out << coloring[i] + 1 <<
'\n';
7512 out <<
"POINT_DATA " << np <<
'\n';
7517 int delete_el_to_el = (el_to_el) ? (0) : (1);
7518 const Table &el_el = ElementToElementTable();
7519 int num_el = GetNE(), stack_p, stack_top_p, max_num_col;
7522 const int *i_el_el = el_el.
GetI();
7523 const int *j_el_el = el_el.
GetJ();
7528 stack_p = stack_top_p = 0;
7529 for (
int el = el0; stack_top_p < num_el; el=(el+1)%num_el)
7531 if (colors[el] != -2)
7535 el_stack[stack_top_p++] = el;
7537 for ( ; stack_p < stack_top_p; stack_p++)
7539 int i = el_stack[stack_p];
7540 int num_nb = i_el_el[i+1] - i_el_el[i];
7541 if (max_num_col < num_nb + 1)
7542 max_num_col = num_nb + 1;
7543 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
7546 if (colors[k] == -2)
7549 el_stack[stack_top_p++] = k;
7557 for (stack_p = 0; stack_p < stack_top_p; stack_p++)
7559 int i = el_stack[stack_p], col;
7561 for (
int j = i_el_el[i]; j < i_el_el[i+1]; j++)
7563 col = colors[j_el_el[j]];
7565 col_marker[col] = 1;
7568 for (col = 0; col < max_num_col; col++)
7569 if (col_marker[col] == 0)
7575 if (delete_el_to_el)
7583 int elem_attr)
const
7585 if (Dim != 3 && Dim != 2)
return;
7587 int i, j, k, l, nv, nbe, *v;
7589 out <<
"MFEM mesh v1.0\n";
7593 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
7598 "# TETRAHEDRON = 4\n"
7602 out <<
"\ndimension\n" << Dim
7603 <<
"\n\nelements\n" << NumOfElements <<
'\n';
7604 for (i = 0; i < NumOfElements; i++)
7606 out << int((elem_attr) ? partitioning[i]+1 : elements[i]->GetAttribute())
7607 <<
' ' << elements[i]->GetGeometryType();
7608 nv = elements[i]->GetNVertices();
7609 v = elements[i]->GetVertices();
7610 for (j = 0; j < nv; j++)
7615 for (i = 0; i < faces_info.Size(); i++)
7617 if ((l = faces_info[i].Elem2No) >= 0)
7619 k = partitioning[faces_info[i].Elem1No];
7620 l = partitioning[l];
7627 out <<
"\nboundary\n" << nbe <<
'\n';
7628 for (i = 0; i < faces_info.Size(); i++)
7630 if ((l = faces_info[i].Elem2No) >= 0)
7632 k = partitioning[faces_info[i].Elem1No];
7633 l = partitioning[l];
7636 nv = faces[i]->GetNVertices();
7637 v = faces[i]->GetVertices();
7638 out << k+1 <<
' ' << faces[i]->GetGeometryType();
7639 for (j = 0; j < nv; j++)
7642 out << l+1 <<
' ' << faces[i]->GetGeometryType();
7643 for (j = nv-1; j >= 0; j--)
7650 k = partitioning[faces_info[i].Elem1No];
7651 nv = faces[i]->GetNVertices();
7652 v = faces[i]->GetVertices();
7653 out << k+1 <<
' ' << faces[i]->GetGeometryType();
7654 for (j = 0; j < nv; j++)
7659 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
7662 out << spaceDim <<
'\n';
7663 for (i = 0; i < NumOfVertices; i++)
7665 out << vertices[i](0);
7666 for (j = 1; j < spaceDim; j++)
7667 out <<
' ' << vertices[i](j);
7682 MFEM_ASSERT(Dim == spaceDim,
"2D Manifolds not supported\n");
7683 if (Dim != 3 && Dim != 2)
return;
7690 int *vcount =
new int[NumOfVertices];
7691 for (i = 0; i < NumOfVertices; i++)
7693 for (i = 0; i < NumOfElements; i++)
7695 nv = elements[i]->GetNVertices();
7696 ind = elements[i]->GetVertices();
7697 for (j = 0; j < nv; j++)
7701 int *voff =
new int[NumOfVertices+1];
7703 for (i = 1; i <= NumOfVertices; i++)
7704 voff[i] = vcount[i-1] + voff[i-1];
7706 int **vown =
new int*[NumOfVertices];
7707 for (i = 0; i < NumOfVertices; i++)
7708 vown[i] =
new int[vcount[i]];
7717 Transpose(ElementToEdgeTable(), edge_el);
7720 for (i = 0; i < NumOfElements; i++)
7722 nv = elements[i]->GetNVertices();
7723 ind = elements[i]->GetVertices();
7724 for (j = 0; j < nv; j++)
7727 vown[ind[j]][vcount[ind[j]]] = i;
7731 for (i = 0; i < NumOfVertices; i++)
7732 vcount[i] = voff[i+1] - voff[i];
7735 for (i = 0; i < edge_el.
Size(); i++)
7737 const int *el = edge_el.
GetRow(i);
7740 k = partitioning[el[0]];
7741 l = partitioning[el[1]];
7742 if (interior_faces || k != l)
7750 out <<
"areamesh2\n\n" << nbe <<
'\n';
7752 for (i = 0; i < edge_el.
Size(); i++)
7754 const int *el = edge_el.
GetRow(i);
7757 k = partitioning[el[0]];
7758 l = partitioning[el[1]];
7759 if (interior_faces || k != l)
7762 GetEdgeVertices(i,ev);
7764 for (j = 0; j < 2; j++)
7765 for (s = 0; s < vcount[ev[j]]; s++)
7766 if (vown[ev[j]][s] == el[0])
7767 out <<
' ' << voff[ev[j]]+s+1;
7770 for (j = 1; j >= 0; j--)
7771 for (s = 0; s < vcount[ev[j]]; s++)
7772 if (vown[ev[j]][s] == el[1])
7773 out <<
' ' << voff[ev[j]]+s+1;
7779 k = partitioning[el[0]];
7781 GetEdgeVertices(i,ev);
7783 for (j = 0; j < 2; j++)
7784 for (s = 0; s < vcount[ev[j]]; s++)
7785 if (vown[ev[j]][s] == el[0])
7786 out <<
' ' << voff[ev[j]]+s+1;
7792 out << NumOfElements <<
'\n';
7793 for (i = 0; i < NumOfElements; i++)
7795 nv = elements[i]->GetNVertices();
7796 ind = elements[i]->GetVertices();
7797 out << partitioning[i]+1 <<
' ';
7799 for (j = 0; j < nv; j++)
7801 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
7802 vown[ind[j]][vcount[ind[j]]] = i;
7807 for (i = 0; i < NumOfVertices; i++)
7808 vcount[i] = voff[i+1] - voff[i];
7811 out << voff[NumOfVertices] <<
'\n';
7812 for (i = 0; i < NumOfVertices; i++)
7813 for (k = 0; k < vcount[i]; k++)
7815 for (j = 0; j < Dim; j++)
7816 out << vertices[i](j) <<
' ';
7826 out <<
"NETGEN_Neutral_Format\n";
7828 out << voff[NumOfVertices] <<
'\n';
7829 for (i = 0; i < NumOfVertices; i++)
7830 for (k = 0; k < vcount[i]; k++)
7832 for (j = 0; j < Dim; j++)
7833 out <<
' ' << vertices[i](j);
7838 out << NumOfElements <<
'\n';
7839 for (i = 0; i < NumOfElements; i++)
7841 nv = elements[i]->GetNVertices();
7842 ind = elements[i]->GetVertices();
7843 out << partitioning[i]+1;
7844 for (j = 0; j < nv; j++)
7846 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
7847 vown[ind[j]][vcount[ind[j]]] = i;
7852 for (i = 0; i < NumOfVertices; i++)
7853 vcount[i] = voff[i+1] - voff[i];
7858 for (i = 0; i < NumOfFaces; i++)
7859 if ((l = faces_info[i].Elem2No) >= 0)
7861 k = partitioning[faces_info[i].Elem1No];
7862 l = partitioning[l];
7863 if (interior_faces || k != l)
7870 for (i = 0; i < NumOfFaces; i++)
7871 if ((l = faces_info[i].Elem2No) >= 0)
7873 k = partitioning[faces_info[i].Elem1No];
7874 l = partitioning[l];
7875 if (interior_faces || k != l)
7877 nv = faces[i]->GetNVertices();
7878 ind = faces[i]->GetVertices();
7880 for (j = 0; j < nv; j++)
7881 for (s = 0; s < vcount[ind[j]]; s++)
7882 if (vown[ind[j]][s] == faces_info[i].Elem1No)
7883 out <<
' ' << voff[ind[j]]+s+1;
7886 for (j = nv-1; j >= 0; j--)
7887 for (s = 0; s < vcount[ind[j]]; s++)
7888 if (vown[ind[j]][s] == faces_info[i].Elem2No)
7889 out <<
' ' << voff[ind[j]]+s+1;
7895 k = partitioning[faces_info[i].Elem1No];
7896 nv = faces[i]->GetNVertices();
7897 ind = faces[i]->GetVertices();
7899 for (j = 0; j < nv; j++)
7900 for (s = 0; s < vcount[ind[j]]; s++)
7901 if (vown[ind[j]][s] == faces_info[i].Elem1No)
7902 out <<
' ' << voff[ind[j]]+s+1;
7906 for (i = 0; i < NumOfVertices; i++)
7909 else if (meshgen == 2)
7914 for (i = 0; i < NumOfFaces; i++)
7915 if ((l = faces_info[i].Elem2No) >= 0)
7917 k = partitioning[faces_info[i].Elem1No];
7918 l = partitioning[l];
7919 if (interior_faces || k != l)
7927 <<
"1 " << voff[NumOfVertices] <<
" " << NumOfElements
7928 <<
" 0 0 0 0 0 0 0\n"
7929 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
7930 <<
"0 0 " << nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
7931 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
7932 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
7934 for (i = 0; i < NumOfVertices; i++)
7935 for (k = 0; k < vcount[i]; k++)
7936 out << voff[i]+k <<
" 0.0 " << vertices[i](0) <<
' '
7937 << vertices[i](1) <<
' ' << vertices[i](2) <<
" 0.0\n";
7939 for (i = 0; i < NumOfElements; i++)
7941 nv = elements[i]->GetNVertices();
7942 ind = elements[i]->GetVertices();
7943 out << i+1 <<
' ' << partitioning[i]+1;
7944 for (j = 0; j < nv; j++)
7946 out <<
' ' << voff[ind[j]]+vcount[ind[j]]--;
7947 vown[ind[j]][vcount[ind[j]]] = i;
7952 for (i = 0; i < NumOfVertices; i++)
7953 vcount[i] = voff[i+1] - voff[i];
7956 for (i = 0; i < NumOfFaces; i++)
7957 if ((l = faces_info[i].Elem2No) >= 0)
7959 k = partitioning[faces_info[i].Elem1No];
7960 l = partitioning[l];
7961 if (interior_faces || k != l)
7963 nv = faces[i]->GetNVertices();
7964 ind = faces[i]->GetVertices();
7966 for (j = 0; j < nv; j++)
7967 for (s = 0; s < vcount[ind[j]]; s++)
7968 if (vown[ind[j]][s] == faces_info[i].Elem1No)
7969 out <<
' ' << voff[ind[j]]+s+1;
7970 out <<
" 1.0 1.0 1.0 1.0\n";
7972 for (j = nv-1; j >= 0; j--)
7973 for (s = 0; s < vcount[ind[j]]; s++)
7974 if (vown[ind[j]][s] == faces_info[i].Elem2No)
7975 out <<
' ' << voff[ind[j]]+s+1;
7976 out <<
" 1.0 1.0 1.0 1.0\n";
7981 k = partitioning[faces_info[i].Elem1No];
7982 nv = faces[i]->GetNVertices();
7983 ind = faces[i]->GetVertices();
7985 for (j = 0; j < nv; j++)
7986 for (s = 0; s < vcount[ind[j]]; s++)
7987 if (vown[ind[j]][s] == faces_info[i].Elem1No)
7988 out <<
' ' << voff[ind[j]]+s+1;
7989 out <<
" 1.0 1.0 1.0 1.0\n";
8007 " NURBS mesh is not supported!");
8011 out <<
"MFEM mesh v1.0\n";
8015 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
8020 "# TETRAHEDRON = 4\n"
8024 out <<
"\ndimension\n" << Dim
8025 <<
"\n\nelements\n" << NumOfElements <<
'\n';
8026 for (i = 0; i < NumOfElements; i++)
8027 PrintElement(elements[i], out);
8030 const int *
const i_AF_f = Aface_face.
GetI();
8031 const int *
const j_AF_f = Aface_face.
GetJ();
8033 for(
int iAF=0; iAF < Aface_face.
Size(); ++iAF)
8034 for(
const int * iface = j_AF_f + i_AF_f[iAF]; iface < j_AF_f + i_AF_f[iAF+1]; ++iface)
8036 out << iAF+1 <<
' ';
8037 PrintElementWithoutAttr(faces[*iface],out);
8040 out <<
"\nvertices\n" << NumOfVertices <<
'\n';
8043 out << spaceDim <<
'\n';
8044 for (i = 0; i < NumOfVertices; i++)
8046 out << vertices[i](0);
8047 for (j = 1; j < spaceDim; j++)
8048 out <<
' ' << vertices[i](j);
8064 int na = attributes.
Size();
8065 double *cg =
new double[na*spaceDim];
8066 int *nbea =
new int[na];
8068 int *vn =
new int[NumOfVertices];
8069 for (i = 0; i < NumOfVertices; i++)
8071 for (i = 0; i < na; i++)
8073 for (j = 0; j < spaceDim; j++)
8074 cg[i*spaceDim+j] = 0.0;
8078 for (i = 0; i < NumOfElements; i++)
8080 GetElementVertices(i, vert);
8081 for (k = 0; k < vert.
Size(); k++)
8085 for (i = 0; i < NumOfElements; i++)
8087 int bea = GetAttribute(i)-1;
8088 GetPointMatrix(i, pointmat);
8089 GetElementVertices(i, vert);
8091 for (k = 0; k < vert.
Size(); k++)
8092 if (vn[vert[k]] == 1)
8095 for (j = 0; j < spaceDim; j++)
8096 cg[bea*spaceDim+j] += pointmat(j,k);
8101 for (i = 0; i < NumOfElements; i++)
8103 int bea = GetAttribute(i)-1;
8104 GetElementVertices (i, vert);
8106 for (k = 0; k < vert.
Size(); k++)
8109 for (j = 0; j < spaceDim; j++)
8110 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
8111 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
8126 int na = NumOfElements;
8127 double *cg =
new double[na*spaceDim];
8128 int *nbea =
new int[na];
8130 int *vn =
new int[NumOfVertices];
8131 for (i = 0; i < NumOfVertices; i++)
8133 for (i = 0; i < na; i++)
8135 for (j = 0; j < spaceDim; j++)
8136 cg[i*spaceDim+j] = 0.0;
8140 for (i = 0; i < NumOfElements; i++)
8142 GetElementVertices(i, vert);
8143 for (k = 0; k < vert.
Size(); k++)
8147 for (i = 0; i < NumOfElements; i++)
8150 GetPointMatrix(i, pointmat);
8151 GetElementVertices(i, vert);
8153 for (k = 0; k < vert.
Size(); k++)
8154 if (vn[vert[k]] == 1)
8157 for (j = 0; j < spaceDim; j++)
8158 cg[bea*spaceDim+j] += pointmat(j,k);
8163 for (i = 0; i < NumOfElements; i++)
8166 GetElementVertices(i, vert);
8168 for (k = 0; k < vert.
Size(); k++)
8171 for (j = 0; j < spaceDim; j++)
8172 vertices[vert[k]](j) = sf*vertices[vert[k]](j) +
8173 (1-sf)*cg[bea*spaceDim+j]/nbea[bea];
8188 Vector vold(spaceDim), vnew(NULL, spaceDim);
8189 for (
int i = 0; i < vertices.Size(); i++)
8191 for (
int j = 0; j < spaceDim; j++)
8192 vold(j) = vertices[i](j);
8201 xnew.ProjectCoefficient(f_pert);
8208 #ifdef MFEM_USE_MEMALLOC
8214 default:
delete E;
break;
8225 if (own_nodes)
delete Nodes;
8231 for (i = 0; i < NumOfElements; i++)
8232 FreeElement(elements[i]);
8234 for (i = 0; i < NumOfBdrElements; i++)
8235 FreeElement(boundary[i]);
8237 for (i = 0; i < faces.Size(); i++)
8238 FreeElement(faces[i]);
8264 V(1) = s * ((ip.
y + layer) / n);
8269 V(2) = s * ((ip.
z + layer) / n);
8278 cerr <<
"Extrude1D : Not a 1D mesh!" << endl;
8282 int nvy = (closed) ? (ny) : (ny + 1);
8283 int nvt = mesh->
GetNV() * nvy;
8290 mesh2d =
new Mesh(2, nvt, mesh->
GetNE()*ny,
8295 for (
int i = 0; i < mesh->
GetNV(); i++)
8298 for (
int j = 0; j < nvy; j++)
8300 vc[1] = sy * (double(j) / ny);
8306 for (
int i = 0; i < mesh->
GetNE(); i++)
8311 for (
int j = 0; j < ny; j++)
8314 qv[0] = vert[0] * nvy + j;
8315 qv[1] = vert[1] * nvy + j;
8316 qv[2] = vert[1] * nvy + (j + 1) % nvy;
8317 qv[3] = vert[0] * nvy + (j + 1) % nvy;
8323 for (
int i = 0; i < mesh->
GetNBE(); i++)
8328 for (
int j = 0; j < ny; j++)
8331 sv[0] = vert[0] * nvy + j;
8332 sv[1] = vert[0] * nvy + (j + 1) % nvy;
8335 Swap<int>(sv[0], sv[1]);
8346 for (
int i = 0; i < mesh->
GetNE(); i++)
8352 sv[0] = vert[0] * nvy;
8353 sv[1] = vert[1] * nvy;
8357 sv[0] = vert[1] * nvy + ny;
8358 sv[1] = vert[0] * nvy + ny;
8374 string cname = name;
8375 if (cname ==
"Linear")
8377 else if (cname ==
"Quadratic")
8379 else if (cname ==
"Cubic")
8381 else if (!strncmp(name,
"H1_", 3))
8383 else if (!strncmp(name,
"L2_T", 4))
8385 else if (!strncmp(name,
"L2_", 3))
8390 cerr <<
"Extrude1D : The mesh uses unknown FE collection : "
8402 for (
int i = 0; i < mesh->
GetNE(); i++)
8405 for (
int j = ny-1; j >= 0; j--)
int GetNPoints() const
Returns the number of the points in the integration rule.
Abstract class for Finite Elements.
void AddHex(const int *vi, int attr=1)
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
int Size() const
For backward compatibility define Size to be synonym of Width()
void SetState(int s)
Change the mesh state to NORMAL, TWO_LEVEL_COARSE, TWO_LEVEL_FINE.
std::ostream & operator<<(std::ostream &out, const Mesh &mesh)
int Size() const
Logical size of the array.
void GetFaceEdges(int i, Array< int > &, Array< int > &) const
void PrintSurfaces(const Table &Aface_face, std::ostream &out) const
Print set of disjoint surfaces:
void GetPointMatrix(int i, DenseMatrix &pointmat) const
int * CartesianPartitioning(int nxyz[])
void GetVerticesElementsBoundary(Array< mfem::Vertex > &vertices, Array< mfem::Element * > &elements, Array< mfem::Element * > &boundary)
Class for integration rule.
void UpdateOrder(int Order)
Change the order of the collection.
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
void ScaleElements(double sf)
Class for grid function - Vector with associated FE space.
Table * GetEdgeVertexTable() const
Returns the edge-to-vertex Table (3D)
void FreeElement(Element *E)
void GetFaceInfos(int Face, int *Inf1, int *Inf2)
void SetVertices(const Vector &vert_coord)
virtual Element * Duplicate(Mesh *m) const =0
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
void METIS_PartGraphRecursive(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
void AddColumnsInRow(int r, int ncol)
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
void AddHexAsTets(const int *vi, int attr=1)
void GetFaceElements(int Face, int *Elem1, int *Elem2)
Array< Element * > boundary
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
int * GeneratePartitioning(int nparts, int part_method=1)
void BisectTetTrans(DenseMatrix &pointmat, Tetrahedron *tet, int child)
void MoveVertices(const Vector &displacements)
void SetSize(int s)
Resizes the vector if the new size is different.
int GetNBE() const
Returns number of boundary elements.
void GetLocalPtToSegTransformation(IsoparametricTransformation &, int)
Used in GetFaceElementTransformations (...)
void PrintWithPartitioning(int *partitioning, std::ostream &out, int elem_attr=0) const
void MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of 'fec' and 'fes'.
double Det() const
Calculates the determinant of the matrix (for 2x2 or 3x3 matrices)
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
void DegreeElevate(int t)
RefinedGeometry * Refine(int Geom, int Times, int ETimes=1)
void AddTri(const int *vi, int attr=1)
int Push(int r, int c, int f)
Piecewise-(bi)linear continuous finite elements.
int Size() const
Returns the size of the vector.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Mesh * Extrude1D(Mesh *mesh, const int ny, const double sy, const bool closed)
Extrude a 1D mesh.
void Transform(void(*f)(const Vector &, Vector &))
int GetBdrElementEdgeIndex(int i) const
int GetNE() const
Returns number of elements.
void BisectTriTrans(DenseMatrix &pointmat, Triangle *tri, int child)
void GenerateBoundaryElements()
void GetElementEdges(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all edges of element i.
void DeleteCoarseTables()
void GetElementLocalToGlobal(Array< int > &lelem_elem)
ElementTransformation * GetFineElemTrans(int i, int j)
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
static int GetQuadOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
Data type quadrilateral element.
void GetVertices(Vector &vert_coord) const
virtual ~Mesh()
Destroys mesh.
void AddBdrSegment(const int *vi, int attr=1)
const IntegrationPoint & GetCenter(int GeomType)
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
void CheckBdrElementOrientation(bool fix_it=true)
Check the orientation of the boundary elements.
NURBSExtension * GetNURBSext()
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
int Size_of_connections() const
const IntegrationRule * GetVertices(int GeomType)
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
void GetVertexToVertexTable(DSTable &) const
void MarkTriMeshForRefinement()
void skip_comment_lines(std::istream &is, const char comment_char)
int GetGeometryType() const
void DeleteAll()
Delete whole array.
void AddConnections(int r, const int *c, int nc)
Element * ReadElement(std::istream &)
void KnotInsert(Array< KnotVector * > &kv)
int GetFineElemPath(int i, int j)
Element * NewElement(int geom)
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr)
static const int tet_faces[4][3]
ElementTransformation * GetBdrElementTransformation(int i)
Returns the transformation defining the i-th boundary element.
template void SortPairs< int, int >(Pair< int, int > *, int)
void AddVertex(const double *)
void PrintTopo(std::ostream &out, const Array< int > &e_to_k) const
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void Load(std::istream &input, int generate_edges=0, int refine=1, bool fix_orientation=true)
Piecewise-(bi)cubic continuous finite elements.
PointFiniteElement PointFE
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void GetElementJacobian(int i, DenseMatrix &J)
Data type hexahedron element.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
void MoveNodes(const Vector &displacements)
void GetEdgeOrdering(DSTable &v_to_v, Array< int > &order)
IntegrationRules IntRules(0)
A global object with all integration rules (defined in intrules.cpp)
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
static const int quad_orientations[8][4]
void FinalizeTetMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
void DetOfLinComb(const DenseMatrix &A, const DenseMatrix &B, Vector &c)
void AddBdrQuadAsTriangles(const int *vi, int attr=1)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
void AddSegmentFaceElement(int lf, int gf, int el, int v0, int v1)
int Append(const T &el)
Append element to array, resize if necessary.
void MarkTetMeshForRefinement()
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
void GetLocalSegToTriTransformation(IsoparametricTransformation &loc, int i)
void AddQuadFaceElement(int lf, int gf, int el, int v0, int v1, int v2, int v3)
TriLinear3DFiniteElement HexahedronFE
int GetRefinementType(int i)
For a given coarse element i returns its refinement type.
void AddConnection(int r, int c)
static FiniteElement * GetTransformationFEforElementType(int)
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
STable3D * GetElementToFaceTable(int ret_ftbl=0)
void Make2D(int nx, int ny, Element::Type type, int generate_edges, double sx, double sy)
Type
Constants for the classes derived from Element.
void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void GetBdrElementEdges(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all edges of bdr element i.
void CheckDisplacements(const Vector &displacements, double &tmax)
void CheckElementOrientation(bool fix_it=true)
Check the orientation of the elements.
NURBSExtension * StealNURBSext()
const IntegrationRule & GetNodes() const
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
GeometryRefiner GlobGeometryRefiner
void SetLayer(const int l)
int GetAttribute() const
Return element's attribute.
int GetElementToEdgeTable(Table &, Array< int > &)
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
int GetElementType(int i) const
Returns the type of element i.
Data type triangle element.
const Element * GetElement(int i) const
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
void GetLocalQuadToHexTransformation(IsoparametricTransformation &loc, int i)
Used in GetFaceElementTransformations (...)
void Sort()
Sorts the array.
void SetNodalFESpace(FiniteElementSpace *nfes)
int Size() const
Returns the number of TYPE I elements.
int GetVDim() const
Returns vector dimension.
A class for non-conforming AMR on higher-order hexahedral, quadrilateral or triangular meshes...
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
static const char * Name[NumGeom]
FiniteElementSpace * FESpace()
virtual void ReorientTetMesh()
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
virtual void LocalRefinement(const Array< int > &marked_el, int type=3)
This function is not public anymore. Use GeneralRefinement instead.
int METIS_SetDefaultOptions(idx_t *options)
int GetBdrElementType(int i) const
Returns the type of boundary element i.
void FindPartitioningComponents(Table &elem_elem, const Array< int > &partitioning, Array< int > &component, Array< int > &num_comp)
Data type tetrahedron element.
int GetGeomType() const
Returns the geometry type:
double GetElementSize(int i, int type=0)
int SpaceDimension() const
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
void CheckPartitioning(int *partitioning)
virtual void Print(std::ostream &out=std::cout) const
Print the mesh to the given stream using the default MFEM mesh format.
void PrintElementsWithPartitioning(int *partitioning, std::ostream &out, int interior_faces=0)
void AddPointFaceElement(int lf, int gf, int el)
Used in GenerateFaces()
void Make1D(int n, double sx=1.0)
Creates a 1D mesh for the interval [0,sx] divided into n equal intervals.
double * Data() const
Returns vector of the elements.
int GetNumFineElems(int i)
void Swap(Array< T > &, Array< T > &)
Array< int > bdr_attributes
void NonconformingRefinement(const Array< Refinement > &refinements, int nc_limit=0)
This function is not public anymore. Use GeneralRefinement instead.
void PrintVTK(std::ostream &out)
Print the mesh in VTK format (linear and quadratic meshes only).
void SetCoarseElem(Element *ce)
MemAlloc< BisectedElement, 1024 > BEMemory
int GetBisectionHierarchy(Element *E)
Abstract finite element space.
int GetDof() const
Returns the degrees of freedom in the FE space.
static void PrintElement(const Element *, std::ostream &)
virtual void NURBSUniformRefinement()
Refine NURBS mesh.
void METIS_PartGraphKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
static void GetElementArrayEdgeTable(const Array< Element * > &elem_array, const DSTable &v_to_v, Table &el_to_edge)
void AddAColumnInRow(int r)
void mfem_error(const char *msg)
GridFunction * GetNodes()
Return a pointer to the internal node GridFunction (may be NULL).
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert)
static const int edges[12][2]
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
int FindRoots(const Vector &z, Vector &x)
void AddQuad(const int *vi, int attr=1)
int Push4(int r, int c, int f, int t)
void GetElementVDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void Save(std::ostream &out)
Prints array to stream out.
virtual void QuadUniformRefinement()
Refine quadrilateral mesh.
virtual void HexUniformRefinement()
Refine hexahedral mesh.
Linear3DFiniteElement TetrahedronFE
void Swap(Mesh &other, bool non_geometry=false)
Array< Element * > elements
const Table & ElementToElementTable()
double GetLength(int i, int j) const
Return the length of the segment from node i to node j.
Linear2DFiniteElement TriangleFE
void METIS_PartGraphVKway(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *, int *, idxtype *)
virtual int * DofOrderForOrientation(int GeomType, int Or) const =0
void AddBdrTriangle(const int *vi, int attr=1)
Table * GetFaceToElementTable() const
virtual void GetVertices(Array< int > &v) const
Returns the indices of the element's vertices.
void GetLocalTriToTetTransformation(IsoparametricTransformation &loc, int i)
Used in GetFaceElementTransformations (...)
NURBSExtension * NURBSext
void AddBdrQuad(const int *vi, int attr=1)
virtual const char * Name() const
const Table & ElementToFaceTable() const
Class for integration point with weight.
Element * ReadElementWithoutAttr(std::istream &)
void AddTriangle(const int *vi, int attr=1)
static const int hex_faces[6][4]
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
void GetLocalSegToQuadTransformation(IsoparametricTransformation &loc, int i)
const FiniteElementSpace * GetNodalFESpace()
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
NodeExtrudeCoefficient(const int dim, const int _n, const double _s)
void Destroy()
Destroy a vector.
void FindTMax(Vector &c, Vector &x, double &tmax, const double factor, const int Dim)
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i'th element.
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Array< FaceInfo > faces_info
virtual int GetNVertices() const =0
void SetAttribute(const int attr)
Set element's attribute.
void ProjectCoefficient(Coefficient &coeff)
STable3D * GetFacesTable()
Piecewise-(bi)quadratic continuous finite elements.
void Make3D(int nx, int ny, int nz, Element::Type type, int generate_edges, double sx, double sy, double sz)
static void PrintElementWithoutAttr(const Element *, std::ostream &)
int GetFineElem(int i, int j)
int DofToVDof(int dof, int vd) const
void DofsToVDofs(Array< int > &dofs) const
const FiniteElementCollection * FEColl() const
void GetNodes(Vector &node_coord) const
int GetFaceBaseGeometry(int i) const
void AverageVertices(int *indexes, int n, int result)
int NumberOfEntries() const
static int GetTriOrientation(const int *base, const int *test)
Returns the orientation of "test" relative to "base".
void SetNode(int i, const double *coord)
void GetNode(int i, double *coord)
Arbitrary order H1-conforming (continuous) finite elements.
void XYZ_VectorFunction(const Vector &p, Vector &v)
void GetElementColoring(Array< int > &colors, int el0=0)
void AddTet(const int *vi, int attr=1)
Table * GetFaceEdgeTable() const
Returns the face-to-edge Table (3D)
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
void SetNodes(const Vector &node_coord)
virtual int GetRefinementFlag()
virtual int DofForGeometry(int GeomType) const =0
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
void SetSize(int s)
If the matrix is not a square matrix of size s then recreate it.
void ScaleSubdomains(double sf)
virtual void PrintXG(std::ostream &out=std::cout) const
Print the mesh to the given stream using Netgen/Truegrid format.
void Bisection(int i, const DSTable &, int *, int *, int *)
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
BiLinear2DFiniteElement QuadrilateralFE
Abstract data type element.
int GetAttribute(int i) const
Return the attribute of element i.
void AddTriangleFaceElement(int lf, int gf, int el, int v0, int v1, int v2)
const Table & ElementToEdgeTable() const
Data type line segment element.
Linear1DFiniteElement SegmentFE
double GetElementVolume(int i)
void InitMesh(int _Dim, int _spaceDim, int NVert, int NElem, int NBdrElem)
Begin construction of a mesh.
void PrintCharacteristics(Vector *Vh=NULL, Vector *Vk=NULL)
Table * GetVertexToElementTable()
The returned Table must be destroyed by the caller.
const Element * GetBdrElement(int i) const
void GeneralRefinement(Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
virtual int GetType() const =0
Returns element's type.
void UpdateNodes()
Update the nodes of a curved mesh after refinement.
Arbitrary order "L2-conforming" discontinuous finite elements.
static const int tri_orientations[6][3]
Class used to exrude the nodes of a mesh.