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;