13 #include "../fem/fem.hpp"
14 #include "../general/text.hpp"
19 #ifdef MFEM_USE_NETCDF
28 bool Mesh::remove_unused_vertices =
true;
30 void Mesh::ReadMFEMMesh(std::istream &input,
bool mfem_v11,
int &curved)
39 MFEM_VERIFY(ident ==
"dimension",
"invalid mesh file");
45 MFEM_VERIFY(ident ==
"elements",
"invalid mesh file");
46 input >> NumOfElements;
47 elements.SetSize(NumOfElements);
48 for (
int j = 0; j < NumOfElements; j++)
50 elements[j] = ReadElement(input);
56 MFEM_VERIFY(ident ==
"boundary",
"invalid mesh file");
57 input >> NumOfBdrElements;
58 boundary.SetSize(NumOfBdrElements);
59 for (
int j = 0; j < NumOfBdrElements; j++)
61 boundary[j] = ReadElement(input);
67 if (mfem_v11 && ident ==
"vertex_parents")
69 ncmesh =
new NCMesh(
this, &input);
75 if (ident ==
"coarse_elements")
77 ncmesh->LoadCoarseElements(input);
84 MFEM_VERIFY(ident ==
"vertices",
"invalid mesh file");
85 input >> NumOfVertices;
86 vertices.SetSize(NumOfVertices);
92 spaceDim = atoi(ident.c_str());
93 for (
int j = 0; j < NumOfVertices; j++)
95 for (
int i = 0; i < spaceDim; i++)
97 input >> vertices[j](i);
102 if (ncmesh) { ncmesh->SetVertexPositions(vertices); }
113 if (remove_unused_vertices) { RemoveUnusedVertices(); }
116 void Mesh::ReadLineMesh(std::istream &input)
122 input >> NumOfVertices;
123 vertices.SetSize(NumOfVertices);
125 for (j = 0; j < NumOfVertices; j++)
127 input >> vertices[j](0);
130 input >> NumOfElements;
131 elements.SetSize(NumOfElements);
133 for (j = 0; j < NumOfElements; j++)
135 input >> a >> p1 >> p2;
136 elements[j] =
new Segment(p1-1, p2-1, a);
140 input >> NumOfBdrElements;
141 boundary.SetSize(NumOfBdrElements);
142 for (j = 0; j < NumOfBdrElements; j++)
144 input >> a >> ind[0];
146 boundary[j] =
new Point(ind,a);
150 void Mesh::ReadNetgen2DMesh(std::istream &input,
int &curved)
152 int ints[32], attr, n;
158 input >> NumOfBdrElements;
159 boundary.SetSize(NumOfBdrElements);
160 for (
int i = 0; i < NumOfBdrElements; i++)
163 >> ints[0] >> ints[1];
164 ints[0]--; ints[1]--;
165 boundary[i] =
new Segment(ints, attr);
169 input >> NumOfElements;
170 elements.SetSize(NumOfElements);
171 for (
int i = 0; i < NumOfElements; i++)
174 for (
int j = 0; j < n; j++)
182 elements[i] =
new Segment(ints, attr);
185 elements[i] =
new Triangle(ints, attr);
196 input >> NumOfVertices;
197 vertices.SetSize(NumOfVertices);
198 for (
int i = 0; i < NumOfVertices; i++)
199 for (
int j = 0; j < Dim; j++)
201 input >> vertices[i](j);
206 input >> NumOfVertices;
207 vertices.SetSize(NumOfVertices);
212 void Mesh::ReadNetgen3DMesh(std::istream &input)
220 input >> NumOfVertices;
222 vertices.SetSize(NumOfVertices);
223 for (
int i = 0; i < NumOfVertices; i++)
224 for (
int j = 0; j < Dim; j++)
226 input >> vertices[i](j);
230 input >> NumOfElements;
231 elements.SetSize(NumOfElements);
232 for (
int i = 0; i < NumOfElements; i++)
235 for (
int j = 0; j < 4; j++)
240 #ifdef MFEM_USE_MEMALLOC
242 tet = TetMemory.Alloc();
252 input >> NumOfBdrElements;
253 boundary.SetSize(NumOfBdrElements);
254 for (
int i = 0; i < NumOfBdrElements; i++)
257 for (
int j = 0; j < 3; j++)
262 boundary[i] =
new Triangle(ints, attr);
266 void Mesh::ReadTrueGridMesh(std::istream &input)
268 int i, j, ints[32], attr;
269 const int buflen = 1024;
280 input >> vari >> NumOfVertices >> vari >> vari >> NumOfElements;
281 input.getline(buf, buflen);
282 input.getline(buf, buflen);
284 input.getline(buf, buflen);
285 input.getline(buf, buflen);
286 input.getline(buf, buflen);
289 vertices.SetSize(NumOfVertices);
290 for (i = 0; i < NumOfVertices; i++)
292 input >> vari >> varf >> vertices[i](0) >> vertices[i](1);
293 input.getline(buf, buflen);
297 elements.SetSize(NumOfElements);
298 for (i = 0; i < NumOfElements; i++)
300 input >> vari >> attr;
301 for (j = 0; j < 4; j++)
306 input.getline(buf, buflen);
307 input.getline(buf, buflen);
315 input >> vari >> NumOfVertices >> NumOfElements;
316 input.getline(buf, buflen);
317 input.getline(buf, buflen);
318 input >> vari >> vari >> NumOfBdrElements;
319 input.getline(buf, buflen);
320 input.getline(buf, buflen);
321 input.getline(buf, buflen);
323 vertices.SetSize(NumOfVertices);
324 for (i = 0; i < NumOfVertices; i++)
326 input >> vari >> varf >> vertices[i](0) >> vertices[i](1)
328 input.getline(buf, buflen);
331 elements.SetSize(NumOfElements);
332 for (i = 0; i < NumOfElements; i++)
334 input >> vari >> attr;
335 for (j = 0; j < 8; j++)
340 input.getline(buf, buflen);
344 boundary.SetSize(NumOfBdrElements);
345 for (i = 0; i < NumOfBdrElements; i++)
348 for (j = 0; j < 4; j++)
353 input.getline(buf, buflen);
360 const int Mesh::vtk_quadratic_tet[10] =
361 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
365 const int Mesh::vtk_quadratic_wedge[18] =
366 { 0, 2, 1, 3, 5, 4, 8, 7, 6, 11, 10, 9, 12, 14, 13, 17, 16, 15};
369 const int Mesh::vtk_quadratic_hex[27] =
371 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
372 24, 22, 21, 23, 20, 25, 26
375 void Mesh::ReadVTKMesh(std::istream &input,
int &curved,
int &read_gf,
387 getline(input, buff);
388 getline(input, buff);
392 MFEM_ABORT(
"VTK mesh is not in ASCII format!");
395 getline(input, buff);
397 if (buff !=
"DATASET UNSTRUCTURED_GRID")
399 MFEM_ABORT(
"VTK mesh is not UNSTRUCTURED_GRID!");
410 MFEM_ABORT(
"VTK mesh does not have POINTS data!");
413 while (buff !=
"POINTS");
419 getline(input, buff);
420 for (i = 0; i < points.
Size(); i++)
427 NumOfElements = n = 0;
432 input >> NumOfElements >> n >> ws;
434 for (i = 0; i < n; i++)
436 input >> cells_data[i];
444 if (buff ==
"CELL_TYPES")
446 input >> NumOfElements;
447 elements.
SetSize(NumOfElements);
448 for (j = i = 0; i < NumOfElements; i++)
450 int ct, elem_dim, elem_order = 1;
456 elements[i] =
new Triangle(&cells_data[j+1]);
464 #ifdef MFEM_USE_MEMALLOC
465 elements[i] = TetMemory.Alloc();
466 elements[i]->SetVertices(&cells_data[j+1]);
473 elements[i] =
new Hexahedron(&cells_data[j+1]);
480 new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
481 cells_data[j+4], cells_data[j+6], cells_data[j+5]);
487 elements[i] =
new Triangle(&cells_data[j+1]);
497 #ifdef MFEM_USE_MEMALLOC
498 elements[i] = TetMemory.Alloc();
499 elements[i]->SetVertices(&cells_data[j+1]);
510 new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
511 cells_data[j+4], cells_data[j+6], cells_data[j+5]);
516 elements[i] =
new Hexahedron(&cells_data[j+1]);
519 MFEM_ABORT(
"VTK mesh : cell type " << ct <<
" is not supported!");
522 MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
523 "elements with different dimensions are not supported");
524 MFEM_VERIFY(order == -1 || order == elem_order,
525 "elements with different orders are not supported");
528 j += cells_data[j] + 1;
533 streampos sp = input.tellg();
535 if (buff ==
"CELL_DATA")
538 getline(input, buff);
541 if (!strncmp(buff.c_str(),
"SCALARS material", 16))
543 getline(input, buff);
544 for (i = 0; i < NumOfElements; i++)
547 elements[i]->SetAttribute(attr);
564 vertices.SetSize(np);
565 for (i = 0; i < np; i++)
567 vertices[i](0) = points(3*i+0);
568 vertices[i](1) = points(3*i+1);
569 vertices[i](2) = points(3*i+2);
574 NumOfBdrElements = 0;
583 for (n = i = 0; i < NumOfElements; i++)
585 int *v = elements[i]->GetVertices();
586 int nv = elements[i]->GetNVertices();
587 for (j = 0; j < nv; j++)
588 if (pts_dof[v[j]] == -1)
594 for (n = i = 0; i < np; i++)
595 if (pts_dof[i] != -1)
600 for (i = 0; i < NumOfElements; i++)
602 int *v = elements[i]->GetVertices();
603 int nv = elements[i]->GetNVertices();
604 for (j = 0; j < nv; j++)
606 v[j] = pts_dof[v[j]];
612 for (i = 0; i < np; i++)
614 if ((j = pts_dof[i]) != -1)
616 vertices[j](0) = points(3*i+0);
617 vertices[j](1) = points(3*i+1);
618 vertices[j](2) = points(3*i+2);
623 NumOfBdrElements = 0;
628 finalize_topo =
false;
634 Nodes->MakeOwner(fec);
639 for (n = i = 0; i < NumOfElements; i++)
643 switch (elements[i]->GetGeometryType())
645 case Geometry::TRIANGLE:
646 case Geometry::SQUARE:
647 vtk_mfem = vtk_quadratic_hex;
break;
648 case Geometry::TETRAHEDRON:
649 vtk_mfem = vtk_quadratic_tet;
break;
651 vtk_mfem = vtk_quadratic_hex;
break;
652 case Geometry::PRISM:
653 vtk_mfem = vtk_quadratic_wedge;
break;
659 for (n++, j = 0; j < dofs.
Size(); j++, n++)
661 if (pts_dof[cells_data[n]] == -1)
663 pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
667 if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
669 MFEM_ABORT(
"VTK mesh : inconsistent quadratic mesh!");
676 for (i = 0; i < np; i++)
679 if ((dofs[0] = pts_dof[i]) != -1)
682 for (j = 0; j < dofs.
Size(); j++)
684 (*Nodes)(dofs[j]) = points(3*i+j);
693 void Mesh::ReadNURBSMesh(std::istream &input,
int &curved,
int &read_gf)
697 Dim = NURBSext->Dimension();
698 NumOfVertices = NURBSext->GetNV();
699 NumOfElements = NURBSext->GetNE();
700 NumOfBdrElements = NURBSext->GetNBE();
702 NURBSext->GetElementTopo(elements);
703 NURBSext->GetBdrElementTopo(boundary);
705 vertices.SetSize(NumOfVertices);
707 if (NURBSext->HavePatches())
713 Nodes->MakeOwner(fec);
714 NURBSext->SetCoordsFromPatches(*Nodes);
717 int vd = Nodes->VectorDim();
718 for (
int i = 0; i < vd; i++)
721 Nodes->GetNodalValues(vert_val, i+1);
722 for (
int j = 0; j < NumOfVertices; j++)
724 vertices[j](i) = vert_val(j);
734 void Mesh::ReadInlineMesh(std::istream &input,
bool generate_edges)
762 MFEM_VERIFY(input.get() ==
'=',
763 "Inline mesh expected '=' after keyword " << name);
770 else if (name ==
"ny")
774 else if (name ==
"nz")
778 else if (name ==
"sx")
782 else if (name ==
"sy")
786 else if (name ==
"sz")
790 else if (name ==
"type")
794 if (eltype ==
"segment")
796 type = Element::SEGMENT;
798 else if (eltype ==
"quad")
800 type = Element::QUADRILATERAL;
802 else if (eltype ==
"tri")
804 type = Element::TRIANGLE;
806 else if (eltype ==
"hex")
808 type = Element::HEXAHEDRON;
810 else if (eltype ==
"wedge")
812 type = Element::WEDGE;
814 else if (eltype ==
"tet")
816 type = Element::TETRAHEDRON;
820 MFEM_ABORT(
"unrecognized element type (read '" << eltype
821 <<
"') in inline mesh format. "
822 "Allowed: segment, tri, quad, tet, hex, wedge");
827 MFEM_ABORT(
"unrecognized keyword (" << name
828 <<
") in inline mesh format. "
829 "Allowed: nx, ny, nz, type, sx, sy, sz");
834 if (input.peek() ==
';')
847 if (type == Element::SEGMENT)
849 MFEM_VERIFY(nx > 0 && sx > 0.0,
850 "invalid 1D inline mesh format, all values must be "
852 <<
" nx = " << nx <<
"\n"
853 <<
" sx = " << sx <<
"\n");
856 else if (type == Element::TRIANGLE || type == Element::QUADRILATERAL)
858 MFEM_VERIFY(nx > 0 && ny > 0 && sx > 0.0 && sy > 0.0,
859 "invalid 2D inline mesh format, all values must be "
861 <<
" nx = " << nx <<
"\n"
862 <<
" ny = " << ny <<
"\n"
863 <<
" sx = " << sx <<
"\n"
864 <<
" sy = " << sy <<
"\n");
865 Make2D(nx, ny, type, sx, sy, generate_edges,
true);
867 else if (type == Element::TETRAHEDRON || type == Element::WEDGE ||
868 type == Element::HEXAHEDRON)
870 MFEM_VERIFY(nx > 0 && ny > 0 && nz > 0 &&
871 sx > 0.0 && sy > 0.0 && sz > 0.0,
872 "invalid 3D inline mesh format, all values must be "
874 <<
" nx = " << nx <<
"\n"
875 <<
" ny = " << ny <<
"\n"
876 <<
" nz = " << nz <<
"\n"
877 <<
" sx = " << sx <<
"\n"
878 <<
" sy = " << sy <<
"\n"
879 <<
" sz = " << sz <<
"\n");
880 Make3D(nx, ny, nz, type, sx, sy, sz,
true);
885 MFEM_ABORT(
"For inline mesh, must specify an element type ="
886 " [segment, tri, quad, tet, hex, wedge]");
890 void Mesh::ReadGmshMesh(std::istream &input)
895 input >> version >> binary >> dsize;
898 MFEM_ABORT(
"Gmsh file version < 2.2");
900 if (dsize !=
sizeof(
double))
902 MFEM_ABORT(
"Gmsh file : dsize != sizeof(double)");
904 getline(input, buff);
909 input.read(reinterpret_cast<char*>(&one),
sizeof(one));
912 MFEM_ABORT(
"Gmsh file : wrong binary format");
919 map<int, int> vertices_map;
922 while (input >> buff)
924 if (buff ==
"$Nodes")
926 input >> NumOfVertices;
927 getline(input, buff);
928 vertices.SetSize(NumOfVertices);
930 const int gmsh_dim = 3;
931 double coord[gmsh_dim];
932 for (
int ver = 0; ver < NumOfVertices; ++ver)
936 input.read(reinterpret_cast<char*>(&serial_number),
sizeof(
int));
937 input.read(reinterpret_cast<char*>(coord), gmsh_dim*
sizeof(
double));
941 input >> serial_number;
942 for (
int ci = 0; ci < gmsh_dim; ++ci)
947 vertices[ver] =
Vertex(coord, gmsh_dim);
948 vertices_map[serial_number] = ver;
950 if (static_cast<int>(vertices_map.size()) != NumOfVertices)
952 MFEM_ABORT(
"Gmsh file : vertices indices are not unique");
955 else if (buff ==
"$Elements")
957 int num_of_all_elements;
958 input >> num_of_all_elements;
960 getline(input, buff);
971 int nodes_of_gmsh_element[] =
1028 vector<Element*> elements_0D, elements_1D, elements_2D, elements_3D;
1029 elements_0D.reserve(num_of_all_elements);
1030 elements_1D.reserve(num_of_all_elements);
1031 elements_2D.reserve(num_of_all_elements);
1032 elements_3D.reserve(num_of_all_elements);
1036 int n_elem_part = 0;
1037 const int header_size = 3;
1040 int header[header_size];
1041 int n_elem_one_type;
1043 while (n_elem_part < num_of_all_elements)
1045 input.read(reinterpret_cast<char*>(header),
1046 header_size*
sizeof(
int));
1047 type_of_element = header[0];
1048 n_elem_one_type = header[1];
1051 n_elem_part += n_elem_one_type;
1053 const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
1054 vector<int> data(1+n_tags+n_elem_nodes);
1055 for (
int el = 0; el < n_elem_one_type; ++el)
1057 input.read(reinterpret_cast<char*>(&data[0]),
1058 data.size()*
sizeof(int));
1060 serial_number = data[dd++];
1063 phys_domain = (n_tags > 0) ? data[dd++] : 1;
1066 elem_domain = (n_tags > 1) ? data[dd++] : 0;
1069 n_partitions = (n_tags > 2) ? data[dd++] : 0;
1072 vector<int> vert_indices(n_elem_nodes);
1073 for (
int vi = 0; vi < n_elem_nodes; ++vi)
1075 map<int, int>::const_iterator it =
1076 vertices_map.find(data[1+n_tags+vi]);
1077 if (it == vertices_map.end())
1079 MFEM_ABORT(
"Gmsh file : vertex index doesn't exist");
1081 vert_indices[vi] = it->second;
1085 if (phys_domain <= 0)
1087 MFEM_ABORT(
"Non-positive element attribute in Gmsh mesh!");
1091 switch (type_of_element)
1095 elements_1D.push_back(
1096 new Segment(&vert_indices[0], phys_domain));
1101 elements_2D.push_back(
1102 new Triangle(&vert_indices[0], phys_domain));
1107 elements_2D.push_back(
1113 elements_3D.push_back(
1119 elements_3D.push_back(
1120 new Hexahedron(&vert_indices[0], phys_domain));
1125 elements_0D.push_back(
1126 new Point(&vert_indices[0], phys_domain));
1130 MFEM_WARNING(
"Unsupported Gmsh element type.");
1139 for (
int el = 0; el < num_of_all_elements; ++el)
1141 input >> serial_number >> type_of_element >> n_tags;
1142 vector<int> data(n_tags);
1143 for (
int i = 0; i < n_tags; ++i) { input >> data[i]; }
1146 phys_domain = (n_tags > 0) ? data[0] : 1;
1149 elem_domain = (n_tags > 1) ? data[1] : 0;
1152 n_partitions = (n_tags > 2) ? data[2] : 0;
1155 const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
1156 vector<int> vert_indices(n_elem_nodes);
1158 for (
int vi = 0; vi < n_elem_nodes; ++vi)
1161 map<int, int>::const_iterator it = vertices_map.find(index);
1162 if (it == vertices_map.end())
1164 MFEM_ABORT(
"Gmsh file : vertex index doesn't exist");
1166 vert_indices[vi] = it->second;
1170 if (phys_domain <= 0)
1172 MFEM_ABORT(
"Non-positive element attribute in Gmsh mesh!");
1176 switch (type_of_element)
1180 elements_1D.push_back(
1181 new Segment(&vert_indices[0], phys_domain));
1186 elements_2D.push_back(
1187 new Triangle(&vert_indices[0], phys_domain));
1192 elements_2D.push_back(
1198 elements_3D.push_back(
1204 elements_3D.push_back(
1205 new Hexahedron(&vert_indices[0], phys_domain));
1210 elements_0D.push_back(
1211 new Point(&vert_indices[0], phys_domain));
1215 MFEM_WARNING(
"Unsupported Gmsh element type.");
1222 if (!elements_3D.empty())
1225 NumOfElements = elements_3D.size();
1226 elements.SetSize(NumOfElements);
1227 for (
int el = 0; el < NumOfElements; ++el)
1229 elements[el] = elements_3D[el];
1231 NumOfBdrElements = elements_2D.size();
1232 boundary.SetSize(NumOfBdrElements);
1233 for (
int el = 0; el < NumOfBdrElements; ++el)
1235 boundary[el] = elements_2D[el];
1238 for (
size_t el = 0; el < elements_1D.size(); ++el)
1240 delete elements_1D[el];
1242 for (
size_t el = 0; el < elements_0D.size(); ++el)
1244 delete elements_0D[el];
1247 else if (!elements_2D.empty())
1250 NumOfElements = elements_2D.size();
1251 elements.SetSize(NumOfElements);
1252 for (
int el = 0; el < NumOfElements; ++el)
1254 elements[el] = elements_2D[el];
1256 NumOfBdrElements = elements_1D.size();
1257 boundary.SetSize(NumOfBdrElements);
1258 for (
int el = 0; el < NumOfBdrElements; ++el)
1260 boundary[el] = elements_1D[el];
1263 for (
size_t el = 0; el < elements_0D.size(); ++el)
1265 delete elements_0D[el];
1268 else if (!elements_1D.empty())
1271 NumOfElements = elements_1D.size();
1272 elements.SetSize(NumOfElements);
1273 for (
int el = 0; el < NumOfElements; ++el)
1275 elements[el] = elements_1D[el];
1277 NumOfBdrElements = elements_0D.size();
1278 boundary.SetSize(NumOfBdrElements);
1279 for (
int el = 0; el < NumOfBdrElements; ++el)
1281 boundary[el] = elements_0D[el];
1286 MFEM_ABORT(
"Gmsh file : no elements found");
1290 MFEM_CONTRACT_VAR(n_partitions);
1291 MFEM_CONTRACT_VAR(elem_domain);
1298 #ifdef MFEM_USE_NETCDF
1299 void Mesh::ReadCubit(
const char *filename,
int &curved,
int &read_gf)
1306 const int sideMapTri3[3][2] =
1313 const int sideMapQuad4[4][2] =
1321 const int sideMapTri6[3][3] =
1328 const int sideMapQuad9[4][3] =
1336 const int sideMapTet4[4][3] =
1344 const int sideMapTet10[4][6] =
1352 const int sideMapHex8[6][4] =
1362 const int sideMapHex27[6][9] =
1364 {1,2,6,5,9,14,17,13,26},
1365 {2,3,7,6,10,15,18,14,25},
1366 {4,3,7,8,11,15,19,16,27},
1367 {1,4,8,5,12,16,20,13,24},
1368 {1,4,3,2,12,11,10,9,22},
1369 {5,8,7,6,20,19,18,17,23}
1374 const int mfemToGenesisTet10[10] = {1,2,3,4,5,7,8,6,9,10};
1377 const int mfemToGenesisHex27[27] = {1,2,3,4,5,6,7,8,9,10,11,
1379 12,17,18,19,20,13,14,15,
1381 16,22,26,25,27,24,23,21
1384 const int mfemToGenesisTri6[6] = {1,2,3,4,5,6};
1385 const int mfemToGenesisQuad9[9] = {1,2,3,4,5,6,7,8,9};
1392 char str_dummy[256];
1399 if ((retval = nc_open(filename, NC_NOWRITE, &ncid)))
1401 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1407 size_t num_dim=0, num_nodes=0, num_elem=0, num_el_blk=0, num_side_sets=0;
1409 if ((retval = nc_inq_dimid(ncid,
"num_dim", &
id)) ||
1410 (retval = nc_inq_dim(ncid,
id, str_dummy, &num_dim)) ||
1412 (retval = nc_inq_dimid(ncid,
"num_nodes", &
id)) ||
1413 (retval = nc_inq_dim(ncid,
id, str_dummy, &num_nodes)) ||
1415 (retval = nc_inq_dimid(ncid,
"num_elem", &
id)) ||
1416 (retval = nc_inq_dim(ncid,
id, str_dummy, &num_elem)) ||
1418 (retval = nc_inq_dimid(ncid,
"num_el_blk", &
id)) ||
1419 (retval = nc_inq_dim(ncid,
id, str_dummy, &num_el_blk)))
1421 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1423 if ((retval = nc_inq_dimid(ncid,
"num_side_sets", &
id)) ||
1424 (retval = nc_inq_dim(ncid,
id, str_dummy, &num_side_sets)))
1432 size_t *num_el_in_blk =
new size_t[num_el_blk];
1433 size_t num_node_per_el;
1435 int previous_num_node_per_el = 0;
1436 for (
int i = 0; i < (int) num_el_blk; i++)
1438 sprintf(temp_str,
"num_el_in_blk%d", i+1);
1439 if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
1440 (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_el_in_blk[i])))
1442 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1445 sprintf(temp_str,
"num_nod_per_el%d", i+1);
1446 if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
1447 (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_node_per_el)))
1449 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1456 if ((
int) num_node_per_el != previous_num_node_per_el)
1458 MFEM_ABORT(
"Element blocks of different element types not supported");
1461 previous_num_node_per_el = num_node_per_el;
1465 enum CubitElementType
1487 CubitElementType cubit_element_type = ELEMENT_TRI3;
1488 CubitFaceType cubit_face_type = FACE_EDGE2;
1489 int num_element_linear_nodes = 0;
1493 switch (num_node_per_el)
1497 cubit_element_type = ELEMENT_TRI3;
1498 cubit_face_type = FACE_EDGE2;
1499 num_element_linear_nodes = 3;
1504 cubit_element_type = ELEMENT_TRI6;
1505 cubit_face_type = FACE_EDGE3;
1506 num_element_linear_nodes = 3;
1511 cubit_element_type = ELEMENT_QUAD4;
1512 cubit_face_type = FACE_EDGE2;
1513 num_element_linear_nodes = 4;
1518 cubit_element_type = ELEMENT_QUAD9;
1519 cubit_face_type = FACE_EDGE3;
1520 num_element_linear_nodes = 4;
1525 MFEM_ABORT(
"Don't know what to do with a " << num_node_per_el <<
1526 " node 2D element\n");
1530 else if (num_dim == 3)
1532 switch (num_node_per_el)
1536 cubit_element_type = ELEMENT_TET4;
1537 cubit_face_type = FACE_TRI3;
1538 num_element_linear_nodes = 4;
1543 cubit_element_type = ELEMENT_TET10;
1544 cubit_face_type = FACE_TRI6;
1545 num_element_linear_nodes = 4;
1550 cubit_element_type = ELEMENT_HEX8;
1551 cubit_face_type = FACE_QUAD4;
1552 num_element_linear_nodes = 8;
1557 cubit_element_type = ELEMENT_HEX27;
1558 cubit_face_type = FACE_QUAD9;
1559 num_element_linear_nodes = 8;
1564 MFEM_ABORT(
"Don't know what to do with a " << num_node_per_el <<
1565 " node 3D element\n");
1571 MFEM_ABORT(
"Invalid dimension: num_dim = " << num_dim);
1576 if (cubit_element_type == ELEMENT_TRI3 || cubit_element_type == ELEMENT_QUAD4 ||
1577 cubit_element_type == ELEMENT_TET4 || cubit_element_type == ELEMENT_HEX8)
1581 else if (cubit_element_type == ELEMENT_TRI6 ||
1582 cubit_element_type == ELEMENT_QUAD9 ||
1583 cubit_element_type == ELEMENT_TET10 || cubit_element_type == ELEMENT_HEX27)
1589 size_t *num_side_in_ss =
new size_t[num_side_sets];
1590 for (
int i = 0; i < (int) num_side_sets; i++)
1592 sprintf(temp_str,
"num_side_ss%d", i+1);
1593 if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
1594 (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_side_in_ss[i])))
1596 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1601 double *coordx =
new double[num_nodes];
1602 double *coordy =
new double[num_nodes];
1603 double *coordz =
new double[num_nodes];
1605 if ((retval = nc_inq_varid(ncid,
"coordx", &
id)) ||
1606 (retval = nc_get_var_double(ncid,
id, coordx)) ||
1607 (retval = nc_inq_varid(ncid,
"coordy", &
id)) ||
1608 (retval = nc_get_var_double(ncid,
id, coordy)))
1610 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1615 if ((retval = nc_inq_varid(ncid,
"coordz", &
id)) ||
1616 (retval = nc_get_var_double(ncid,
id, coordz)))
1618 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1623 int **elem_blk =
new int*[num_el_blk];
1624 for (
int i = 0; i < (int) num_el_blk; i++)
1626 elem_blk[i] =
new int[num_el_in_blk[i] * num_node_per_el];
1627 sprintf(temp_str,
"connect%d", i+1);
1628 if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
1629 (retval = nc_get_var_int(ncid, temp_id, elem_blk[i])))
1631 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1634 int *ebprop =
new int[num_el_blk];
1635 if ((retval = nc_inq_varid(ncid,
"eb_prop1", &
id)) ||
1636 (retval = nc_get_var_int(ncid,
id, ebprop)))
1638 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1643 int **elem_ss =
new int*[num_side_sets];
1644 int **side_ss =
new int*[num_side_sets];
1646 for (
int i = 0; i < (int) num_side_sets; i++)
1648 elem_ss[i] =
new int[num_side_in_ss[i]];
1649 side_ss[i] =
new int[num_side_in_ss[i]];
1651 sprintf(temp_str,
"elem_ss%d", i+1);
1652 if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
1653 (retval = nc_get_var_int(ncid, temp_id, elem_ss[i])))
1655 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1658 sprintf(temp_str,
"side_ss%d",i+1);
1659 if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
1660 (retval = nc_get_var_int(ncid, temp_id, side_ss[i])))
1662 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1666 int *ssprop =
new int[num_side_sets];
1667 if ((num_side_sets > 0) &&
1668 ((retval = nc_inq_varid(ncid,
"ss_prop1", &
id)) ||
1669 (retval = nc_get_var_int(ncid,
id, ssprop))))
1671 MFEM_ABORT(
"Fatal NetCDF error: " << nc_strerror(retval));
1677 int num_face_nodes = 0;
1678 int num_face_linear_nodes = 0;
1680 switch (cubit_face_type)
1685 num_face_linear_nodes = 2;
1691 num_face_linear_nodes = 2;
1697 num_face_linear_nodes = 3;
1703 num_face_linear_nodes = 3;
1709 num_face_linear_nodes = 4;
1715 num_face_linear_nodes = 4;
1722 int *start_of_block =
new int[num_el_blk+1];
1723 start_of_block[0] = 0;
1724 for (
int i = 1; i < (int) num_el_blk+1; i++)
1726 start_of_block[i] = start_of_block[i-1] + num_el_in_blk[i-1];
1729 int **ss_node_id =
new int*[num_side_sets];
1731 for (
int i = 0; i < (int) num_side_sets; i++)
1733 ss_node_id[i] =
new int[num_side_in_ss[i]*num_face_nodes];
1734 for (
int j = 0; j < (int) num_side_in_ss[i]; j++)
1736 int glob_ind = elem_ss[i][j]-1;
1739 while (iblk < (
int) num_el_blk && glob_ind >= start_of_block[iblk+1])
1743 if (iblk >= (
int) num_el_blk)
1745 MFEM_ABORT(
"Sideset element does not exist");
1747 loc_ind = glob_ind - start_of_block[iblk];
1748 int this_side = side_ss[i][j];
1749 int ielem = loc_ind*num_node_per_el;
1751 for (
int k = 0; k < num_face_nodes; k++)
1754 switch (cubit_element_type)
1756 case (ELEMENT_TRI3):
1758 inode = sideMapTri3[this_side-1][k];
1761 case (ELEMENT_TRI6):
1763 inode = sideMapTri6[this_side-1][k];
1766 case (ELEMENT_QUAD4):
1768 inode = sideMapQuad4[this_side-1][k];
1771 case (ELEMENT_QUAD9):
1773 inode = sideMapQuad9[this_side-1][k];
1776 case (ELEMENT_TET4):
1778 inode = sideMapTet4[this_side-1][k];
1781 case (ELEMENT_TET10):
1783 inode = sideMapTet10[this_side-1][k];
1786 case (ELEMENT_HEX8):
1788 inode = sideMapHex8[this_side-1][k];
1791 case (ELEMENT_HEX27):
1793 inode = sideMapHex27[this_side-1][k];
1797 ss_node_id[i][j*num_face_nodes+k] =
1798 elem_blk[iblk][ielem + inode - 1];
1804 std::vector<int> uniqueVertexID;
1806 for (
int iblk = 0; iblk < (int) num_el_blk; iblk++)
1808 for (
int i = 0; i < (int) num_el_in_blk[iblk]; i++)
1810 for (
int j = 0; j < num_element_linear_nodes; j++)
1812 uniqueVertexID.push_back(elem_blk[iblk][i*num_node_per_el + j]);
1816 std::sort(uniqueVertexID.begin(), uniqueVertexID.end());
1817 std::vector<int>::iterator newEnd;
1818 newEnd = std::unique(uniqueVertexID.begin(), uniqueVertexID.end());
1819 uniqueVertexID.resize(std::distance(uniqueVertexID.begin(), newEnd));
1825 std::map<int,int> cubitToMFEMVertMap;
1826 for (
int i = 0; i < (int) uniqueVertexID.size(); i++)
1828 cubitToMFEMVertMap[uniqueVertexID[i]] = i+1;
1830 MFEM_ASSERT(cubitToMFEMVertMap.size() == uniqueVertexID.size(),
1831 "This should never happen\n");
1837 NumOfVertices = uniqueVertexID.size();
1838 vertices.SetSize(NumOfVertices);
1839 for (
int i = 0; i < (int) uniqueVertexID.size(); i++)
1841 vertices[i](0) = coordx[uniqueVertexID[i] - 1];
1842 vertices[i](1) = coordy[uniqueVertexID[i] - 1];
1845 vertices[i](2) = coordz[uniqueVertexID[i] - 1];
1849 NumOfElements = num_elem;
1850 elements.SetSize(num_elem);
1852 int renumberedVertID[8];
1853 for (
int iblk = 0; iblk < (int) num_el_blk; iblk++)
1855 int NumNodePerEl = num_node_per_el;
1856 for (
int i = 0; i < (int) num_el_in_blk[iblk]; i++)
1858 for (
int j = 0; j < num_element_linear_nodes; j++)
1860 renumberedVertID[j] =
1861 cubitToMFEMVertMap[elem_blk[iblk][i*NumNodePerEl+j]]-1;
1864 switch (cubit_element_type)
1866 case (ELEMENT_TRI3):
1867 case (ELEMENT_TRI6):
1869 elements[elcount] =
new Triangle(renumberedVertID,ebprop[iblk]);
1872 case (ELEMENT_QUAD4):
1873 case (ELEMENT_QUAD9):
1875 elements[elcount] =
new Quadrilateral(renumberedVertID,ebprop[iblk]);
1878 case (ELEMENT_TET4):
1879 case (ELEMENT_TET10):
1881 elements[elcount] =
new Tetrahedron(renumberedVertID,ebprop[iblk]);
1884 case (ELEMENT_HEX8):
1885 case (ELEMENT_HEX27):
1887 elements[elcount] =
new Hexahedron(renumberedVertID,ebprop[iblk]);
1897 NumOfBdrElements = 0;
1898 for (
int iss = 0; iss < (int) num_side_sets; iss++)
1900 NumOfBdrElements += num_side_in_ss[iss];
1902 boundary.SetSize(NumOfBdrElements);
1904 for (
int iss = 0; iss < (int) num_side_sets; iss++)
1906 for (
int i = 0; i < (int) num_side_in_ss[iss]; i++)
1908 for (
int j = 0; j < num_face_linear_nodes; j++)
1910 renumberedVertID[j] =
1911 cubitToMFEMVertMap[ss_node_id[iss][i*num_face_nodes+j]] - 1;
1913 switch (cubit_face_type)
1918 boundary[sidecount] =
new Segment(renumberedVertID,ssprop[iss]);
1924 boundary[sidecount] =
new Triangle(renumberedVertID,ssprop[iss]);
1930 boundary[sidecount] =
new Quadrilateral(renumberedVertID,ssprop[iss]);
1943 switch (cubit_element_type)
1945 case (ELEMENT_TRI6):
1947 mymap = (
int *) mfemToGenesisTri6;
1950 case (ELEMENT_QUAD9):
1952 mymap = (
int *) mfemToGenesisQuad9;
1955 case (ELEMENT_TET10):
1957 mymap = (
int *) mfemToGenesisTet10;
1960 case (ELEMENT_HEX27):
1962 mymap = (
int *) mfemToGenesisHex27;
1965 case (ELEMENT_TRI3):
1966 case (ELEMENT_QUAD4):
1967 case (ELEMENT_TET4):
1968 case (ELEMENT_HEX8):
1970 MFEM_ABORT(
"Something went wrong. Linear elements detected when order is 2.");
1982 Nodes->MakeOwner(fec);
1990 for (
int i = 0; i < NumOfElements; i++)
1997 for (
int l = 0; l < dofs.
Size(); l++) { vdofs[l] = dofs[l]; }
2001 while (iblk < (
int) num_el_blk && i >= start_of_block[iblk+1]) { iblk++; }
2002 loc_ind = i - start_of_block[iblk];
2003 for (
int j = 0; j < dofs.
Size(); j++)
2005 int point_id = elem_blk[iblk][loc_ind*num_node_per_el + mymap[j] - 1] - 1;
2006 (*Nodes)(vdofs[j]) = coordx[point_id];
2007 (*Nodes)(vdofs[j]+1) = coordy[point_id];
2010 (*Nodes)(vdofs[j]+2) = coordz[point_id];
2020 for (
int i = 0; i < (int) num_side_sets; i++)
2022 delete [] elem_ss[i];
2023 delete [] side_ss[i];
2028 delete [] num_el_in_blk;
2029 delete [] num_side_in_ss;
2034 for (
int i = 0; i < (int) num_el_blk; i++)
2036 delete [] elem_blk[i];
2040 delete [] start_of_block;
2042 for (
int i = 0; i < (int) num_side_sets; i++)
2044 delete [] ss_node_id[i];
2046 delete [] ss_node_id;
2051 #endif // #ifdef MFEM_USE_NETCDF
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
int Size() const
Logical size of the array.
Class for grid function - Vector with associated FE space.
void SetSize(int s)
Resize the vector to size s.
int Size() const
Returns the size of the vector.
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
Data type quadrilateral element.
void skip_comment_lines(std::istream &is, const char comment_char)
void DeleteAll()
Delete whole array.
Data type hexahedron element.
Type
Constants for the classes derived from Element.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Data type triangle element.
A class for non-conforming AMR on higher-order hexahedral, prismatic, quadrilateral or triangular mes...
Data type tetrahedron element.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void Destroy()
Destroy a vector.
int index(int i, int j, int nx, int ny)
void SetAttribute(const int attr)
Set element's attribute.
Piecewise-(bi)quadratic continuous finite elements.
void filter_dos(std::string &line)
Arbitrary order H1-conforming (continuous) finite elements.
Data type line segment element.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const