From c3fac1e81c774a7b1a7ca70657f43492a368b29e Mon Sep 17 00:00:00 2001 From: Dave Parker Date: Mon, 10 Feb 2025 15:51:00 -0800 Subject: [PATCH] Check spatial dimension (#339) * Change file name. * Add setting the mesh element type from vtk cell type and a check for spatial dimension consistent with mesh element type. * Change to check using element dimension. * Add check for fiber mesh. * Add setting mesh element type for vtp files, reorganize checks. --- Code/Source/solver/VtkData.cpp | 14 +++++ Code/Source/solver/VtkData.h | 3 + Code/Source/solver/consts.cpp | 42 ++++++++++++- Code/Source/solver/consts.h | 2 + Code/Source/solver/load_msh.cpp | 87 ++++++++++++++++++--------- Code/Source/solver/nn.cpp | 5 +- Code/Source/solver/nn_elem_props.h | 9 --- Code/Source/solver/vtk_xml.cpp | 7 ++- Code/Source/solver/vtk_xml_parser.cpp | 59 +++++++++++++++++- 9 files changed, 180 insertions(+), 48 deletions(-) diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index 0c7e6fb6b..207306f0a 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -66,6 +66,7 @@ class VtkVtpData::VtkVtpDataImpl { void write(const std::string& file_name); vtkSmartPointer vtk_polydata; + int elem_type; int num_elems; int np_elem; int num_points; @@ -92,6 +93,7 @@ void VtkVtpData::VtkVtpDataImpl::read_file(const std::string& file_name) auto cell = vtkGenericCell::New(); vtk_polydata->GetCell(0, cell); np_elem = cell->GetNumberOfPoints(); + elem_type = cell->GetCellType(); } void VtkVtpData::VtkVtpDataImpl::set_connectivity(const int nsd, const Array& conn, const int pid) @@ -272,6 +274,7 @@ class VtkVtuData::VtkVtuDataImpl { }; vtkSmartPointer vtk_ugrid; + int elem_type; int num_elems; int np_elem; int num_points; @@ -298,6 +301,7 @@ void VtkVtuData::VtkVtuDataImpl::read_file(const std::string& file_name) auto cell = vtkGenericCell::New(); vtk_ugrid->GetCell(0, cell); np_elem = cell->GetNumberOfPoints(); + elem_type = cell->GetCellType(); } void VtkVtuData::VtkVtuDataImpl::set_connectivity(const int nsd, const Array& conn, const int pid) @@ -745,6 +749,11 @@ bool VtkVtpData::has_point_data(const std::string& data_name) return false; } +int VtkVtpData::elem_type() +{ + return impl->elem_type; +} + int VtkVtpData::num_elems() { return impl->num_elems; @@ -1007,6 +1016,11 @@ Array VtkVtuData::get_points() return points_array; } +int VtkVtuData::elem_type() +{ + return impl->elem_type; +} + int VtkVtuData::num_elems() { return impl->num_elems; diff --git a/Code/Source/solver/VtkData.h b/Code/Source/solver/VtkData.h index b60c38fcf..41a240948 100644 --- a/Code/Source/solver/VtkData.h +++ b/Code/Source/solver/VtkData.h @@ -44,6 +44,7 @@ class VtkData { virtual Array get_connectivity() = 0; virtual Array get_points() = 0; virtual int num_elems() = 0; + virtual int elem_type() = 0; virtual int np_elem() = 0; virtual int num_points() = 0; virtual void read_file(const std::string& file_name) = 0; @@ -79,6 +80,7 @@ class VtkVtpData : public VtkData { virtual Array get_connectivity(); virtual Array get_points(); + virtual int elem_type(); virtual int num_elems(); virtual int np_elem(); virtual int num_points(); @@ -115,6 +117,7 @@ class VtkVtuData : public VtkData { ~VtkVtuData(); virtual Array get_connectivity(); + virtual int elem_type(); virtual int num_elems(); virtual int np_elem(); virtual int num_points(); diff --git a/Code/Source/solver/consts.cpp b/Code/Source/solver/consts.cpp index 4f50991fe..06575b11a 100644 --- a/Code/Source/solver/consts.cpp +++ b/Code/Source/solver/consts.cpp @@ -32,7 +32,6 @@ namespace consts { - /// @brief Reproduces the 'SELECT CASE (lM%eType)' statements in the /// Fortran 'PARTMSH' subroutine. // @@ -53,6 +52,47 @@ const std::map element_type_to_elem_nonb = {ElementType::WDG, 3 } }; +/// @brief Map containing the element type as a string. +/// +const std::map element_type_to_string = +{ + {ElementType::NA, "unknown" }, + {ElementType::PNT, "point" }, + {ElementType::LIN1, "line1" }, + {ElementType::LIN2, "line2" }, + {ElementType::TRI3, "tri3" }, + {ElementType::TRI6, "tri6" }, + {ElementType::QUD4, "quad4" }, + {ElementType::QUD8, "quad8" }, + {ElementType::QUD9, "qud9" }, + {ElementType::TET4, "tet4" }, + {ElementType::TET10, "tet10" }, + {ElementType::HEX8, "hex8" }, + {ElementType::HEX20, "hex20" }, + {ElementType::HEX27, "hex27" }, + {ElementType::WDG, "wedge" } +}; + +/// @brief Map containing the dimension for each element type. +/// +const std::map element_dimension = +{ + {ElementType::PNT, 0 }, + {ElementType::LIN1, 1 }, + {ElementType::LIN2, 1 }, + {ElementType::TRI3, 2 }, + {ElementType::TRI6, 2 }, + {ElementType::QUD4, 2 }, + {ElementType::QUD8, 2 }, + {ElementType::QUD9, 2 }, + {ElementType::TET4, 3 }, + {ElementType::TET10, 3 }, + {ElementType::HEX8, 3 }, + {ElementType::HEX20, 3 }, + {ElementType::HEX27, 3 }, + {ElementType::WDG, 3 } +}; + /// @brief Map for constitutive_model string to ConstitutiveModelType. /// /// Type of constitutive model (isochoric) for structure equation: diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index b4d2ff173..326b471ff 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -271,7 +271,9 @@ enum class ElementType NRB = 115 }; +extern const std::map element_type_to_string; extern const std::map element_type_to_elem_nonb; +extern const std::map element_dimension; // Template for printing ElementType. /* diff --git a/Code/Source/solver/load_msh.cpp b/Code/Source/solver/load_msh.cpp index f70b6b0f4..9d87f7099 100644 --- a/Code/Source/solver/load_msh.cpp +++ b/Code/Source/solver/load_msh.cpp @@ -183,11 +183,31 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p dmsg.banner(); dmsg << "Mesh name: " << mesh_name; dmsg << "Mesh path: " << mesh_path; + dmsg << "mesh.lShl: " << mesh.lShl; #endif - // Read in volume mesh. vtk_xml::read_vtu(mesh_path, mesh); + // Check that the input number of spatial dimensions is consistent + // with the types of elements defined for the simulation mesh. + // + int nsd = simulation->com_mod.nsd; + int elem_dim = consts::element_dimension.at(mesh.eType); + auto elem_type = consts::element_type_to_string.at(mesh.eType); + + if (mesh.lShl) { + if (nsd == 1) { + throw std::runtime_error("The number of spatial dimensions (" + std::to_string(nsd) + + ") is not consistent with the mesh '" + mesh.name + "' which contains shell elements."); + } + + } else if (!mesh.lFib) { + if (elem_dim != nsd) { + throw std::runtime_error("The number of spatial dimensions (" + std::to_string(nsd) + + ") is not consistent with the mesh '" + mesh.name + "' which contains " + elem_type + " elements."); + } + } + // Set mesh element properites for the input element type. nn::select_ele(simulation->com_mod, mesh); @@ -226,6 +246,9 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p if (face.qmTRI3 < (1.0 / 3.0) || face.qmTRI3 > 1.0) { throw std::runtime_error("Quadrature_modifier_TRI3 must be in the range [1/3, 1]."); } +#ifdef dbg_read_sv + dmsg << "face.qmTRI3: " << face.qmTRI3; +#endif if (mesh.lFib) { auto face_path = face_param->end_nodes_face_file_path(); @@ -258,42 +281,46 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p } } if (!mesh.lFib) { - // Create a hash map for nodes and elements. - MeshHashMaps mesh_hash_maps; - auto mesh_node_map = mesh_hash_maps.createNodeHashMap(mesh, com_mod.nsd); - auto mesh_element_set = mesh_hash_maps.createElementHashMap(mesh); + // Create a hash map for nodes and elements. + MeshHashMaps mesh_hash_maps; + auto mesh_node_map = mesh_hash_maps.createNodeHashMap(mesh, com_mod.nsd); + auto mesh_element_set = mesh_hash_maps.createElementHashMap(mesh); - for (int i = 0; i < mesh.nFa; i++) { + for (int i = 0; i < mesh.nFa; i++) { auto &face = mesh.fa[i]; if (!mesh.lFib) { - // Set face global node IDs - for (int a = 0; a < face.nNo; a++) { - std::ostringstream key; - key << std::scientific << std::setprecision(16); - for (int j = 0; j < com_mod.nsd; j++) { - key << face.x(j, a) << ","; - } - face.gN(a) = mesh_node_map[key.str()]; + // Set face global node IDs + for (int a = 0; a < face.nNo; a++) { + std::ostringstream key; + key << std::scientific << std::setprecision(16); + for (int j = 0; j < com_mod.nsd; j++) { + key << face.x(j, a) << ","; } - // Set face element IDs - for (int e = 0; e < face.nEl; e++) { - std::vector element_nodes; - for (int a = 0; a < face.eNoN; a++) { - int Ac = face.IEN(a, e); - Ac = face.gN(Ac); - face.IEN(a, e) = Ac; - element_nodes.push_back(Ac); - } - std::string key = ""; - std::sort(element_nodes.begin(), element_nodes.end()); - for (int a = 0; a < face.eNoN; a++) { - key += std::to_string(element_nodes[a]) + ","; - } - face.gE(e) = mesh_element_set[key]; + face.gN(a) = mesh_node_map[key.str()]; + } + + // Set face element IDs + for (int e = 0; e < face.nEl; e++) { + std::vector element_nodes; + for (int a = 0; a < face.eNoN; a++) { + int Ac = face.IEN(a, e); + Ac = face.gN(Ac); + face.IEN(a, e) = Ac; + element_nodes.push_back(Ac); + } + + std::string key = ""; + std::sort(element_nodes.begin(), element_nodes.end()); + for (int a = 0; a < face.eNoN; a++) { + key += std::to_string(element_nodes[a]) + ","; } + + face.gE(e) = mesh_element_set[key]; + } } + } } - } + for (int i = 0; i> // SetElementPropsMapType set_3d_element_props = { - {2, [](int insd, mshType& mesh) -> void { - mesh.eType = ElementType::LIN1; - mesh.nG = 2; - mesh.vtkType = 3; - mesh.nEf = 2; - mesh.lShpF = true; - } - }, - {4, [](int insd, mshType& mesh) -> void { mesh.eType = ElementType::TET4; mesh.nG = 4; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 536127df6..256e31318 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -440,7 +440,7 @@ void read_vtp(const std::string& file_name, faceType& face) using namespace vtk_xml_parser; if (FILE *file = fopen(file_name.c_str(), "r")) { - fclose(file); + fclose(file); } else { throw std::runtime_error("The VTK face file '" + file_name + "' can't be read."); } @@ -560,6 +560,7 @@ void read_vtp_pdata(const std::string& fName, const std::string& kwrd, const int // mesh.eNoN - number of noders per element // mesh.x - node coordinates // mesh.gIEN - element connectivity +// mesh.eType - element type // // Replicates Fortran READVTU subroutine defined in VTKXML.f. // @@ -570,7 +571,7 @@ void read_vtu(const std::string& file_name, mshType& mesh) using namespace vtk_xml_parser; if (FILE *file = fopen(file_name.c_str(), "r")) { - fclose(file); + fclose(file); } else { throw std::runtime_error("The VTU mesh file '" + file_name + "' can't be read."); } @@ -582,6 +583,7 @@ void read_vtu(const std::string& file_name, mshType& mesh) auto vtk_data = VtkData::create_reader(file_name); int num_elems = vtk_data->num_elems(); int np_elem = vtk_data->np_elem(); + int elem_type = vtk_data->elem_type(); // Set mesh data. mesh.nEl = num_elems; @@ -602,7 +604,6 @@ void read_vtu(const std::string& file_name, mshType& mesh) #endif } - //---------- // read_precomputed_solution_vtu //---------- diff --git a/Code/Source/solver/vtk_xml_parser.cpp b/Code/Source/solver/vtk_xml_parser.cpp index 74a77a1f2..08d562218 100644 --- a/Code/Source/solver/vtk_xml_parser.cpp +++ b/Code/Source/solver/vtk_xml_parser.cpp @@ -80,6 +80,22 @@ std::map vtk_cell_to_elem { {VTK_TRIQUADRATIC_HEXAHEDRON, 27} }; +/// Map used to convert VTK cell types to consts::ElementType. +std::map vtk_cell_to_elem_type { + {VTK_HEXAHEDRON, consts::ElementType::HEX8}, + {VTK_LINE, consts::ElementType::LIN1}, + {VTK_QUAD, consts::ElementType::QUD4}, + {VTK_TETRA, consts::ElementType::TET4}, + {VTK_TRIANGLE, consts::ElementType::TRI3}, + {VTK_WEDGE, consts::ElementType::WDG}, + {VTK_QUADRATIC_TRIANGLE, consts::ElementType::TRI6}, + {VTK_QUADRATIC_QUAD, consts::ElementType::QUD8}, + {VTK_BIQUADRATIC_QUAD, consts::ElementType::QUD9}, + {VTK_QUADRATIC_TETRA, consts::ElementType::TET10}, + {VTK_QUADRATIC_HEXAHEDRON, consts::ElementType::HEX20}, + {VTK_TRIQUADRATIC_HEXAHEDRON, consts::ElementType::HEX27} +}; + std::map>> vtk_cell_ordering { {VTK_HEXAHEDRON, {{0,3,2,1}, {4,5,6,7}, @@ -714,6 +730,17 @@ void load_vtp(const std::string& file_name, faceType& face) } vtkIdType num_elems = vtk_polydata->GetNumberOfCells(); + auto cell = vtkGenericCell::New(); + vtk_polydata->GetCell(0, cell); + int np_elem = cell->GetNumberOfPoints(); + int elem_type = cell->GetCellType(); + + // Set the mesh element type from the vtk cell type. + try { + face.eType = vtk_cell_to_elem_type[elem_type]; + } catch (const std::bad_function_call& exception) { + throw std::runtime_error("[load_vtp] The VTK cell type " + std::to_string(elem_type) + " is not supported."); + } #ifdef debug_load_vtp std::cout << "[load_vtp] Number of nodes: " << num_nodes << std::endl; std::cout << "[load_vtp] Number of elements: " << num_elems << std::endl; @@ -741,7 +768,7 @@ void load_vtp(const std::string& file_name, faceType& face) // void load_vtp(const std::string& file_name, mshType& mesh) { - #ifdef debug_load_vtp + #ifdef n_debug_load_vtp std::cout << "[load_vtp] " << std::endl; std::cout << "[load_vtp] ===== vtk_xml_parser.cpp::load_vtp ===== " << std::endl; #endif @@ -756,6 +783,18 @@ void load_vtp(const std::string& file_name, mshType& mesh) } vtkIdType num_elems = vtk_polydata->GetNumberOfCells(); + auto cell = vtkGenericCell::New(); + vtk_polydata->GetCell(0, cell); + int np_elem = cell->GetNumberOfPoints(); + int elem_type = cell->GetCellType(); + + // Set the mesh element type from the vtk cell type. + try { + mesh.eType = vtk_cell_to_elem_type[elem_type]; + } catch (const std::bad_function_call& exception) { + throw std::runtime_error("[load_vtp] The VTK cell type " + std::to_string(elem_type) + " is not supported."); + } + #ifdef debug_load_vtp std::cout << "[load_vtp] Number of nodes: " << num_nodes << std::endl; std::cout << "[load_vtp] Number of elements: " << num_elems << std::endl; @@ -788,6 +827,7 @@ void load_vtp(const std::string& file_name, mshType& mesh) /// mesh.gnEl - number of elements /// mesh.eNoN - number of noders per element /// mesh.gIEN - element connectivity (num_nodes_per_elem, num_elems) +/// mesh.eType - element type /// /// Replicates 'subroutine loadVTK(vtk,fName,istat)' defined in vtkXMLParser.f90. // @@ -814,11 +854,20 @@ void load_vtu(const std::string& file_name, mshType& mesh) auto cell = vtkGenericCell::New(); vtk_ugrid->GetCell(0, cell); int np_elem = cell->GetNumberOfPoints(); + int elem_type = cell->GetCellType(); + + // Set the mesh element type from the vtk cell type. + try { + mesh.eType = vtk_cell_to_elem_type[elem_type]; + } catch (const std::bad_function_call& exception) { + throw std::runtime_error("[load_vtu] The VTK cell type " + std::to_string(elem_type) + " is not supported."); + } #ifdef debug_load_vtu std::cout << "[load_vtu] Number of nodes: " << num_nodes << std::endl; std::cout << "[load_vtu] Number of elements: " << num_elems << std::endl; std::cout << "[load_vtu] Number of nodes per element: " << np_elem << std::endl; + std::cout << "[load_vtu] Element type: " << elem_type << std::endl; #endif // Store nodal coordinates. @@ -857,6 +906,14 @@ void load_vtu(const std::string& file_name, faceType& face) auto cell = vtkGenericCell::New(); vtk_ugrid->GetCell(0, cell); int np_elem = cell->GetNumberOfPoints(); + int elem_type = cell->GetCellType(); + + // Set the mesh element type from the vtk cell type. + try { + face.eType = vtk_cell_to_elem_type[elem_type]; + } catch (const std::bad_function_call& exception) { + throw std::runtime_error("[load_vtu] The VTK cell type " + std::to_string(elem_type) + " is not supported."); + } #ifdef debug_load_vtu std::cout << "[load_vtu] Number of nodes: " << num_nodes << std::endl;