Skip to content

Commit

Permalink
Check spatial dimension (SimVascular#339)
Browse files Browse the repository at this point in the history 
* 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.
  • Loading branch information
@ktbolt
ktbolt authored Feb 10, 2025
1 parent 6d38248 commit c3fac1e
Show file tree
Hide file tree
Showing 9 changed files with 180 additions and 48 deletions.
14 changes: 14 additions & 0 deletions Code/Source/solver/VtkData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@ class VtkVtpData::VtkVtpDataImpl {
void write(const std::string& file_name);

vtkSmartPointer<vtkPolyData> vtk_polydata;
int elem_type;
int num_elems;
int np_elem;
int num_points;
Expand All @@ -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<int>& conn, const int pid)
Expand Down Expand Up @@ -272,6 +274,7 @@ class VtkVtuData::VtkVtuDataImpl {
};

vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid;
int elem_type;
int num_elems;
int np_elem;
int num_points;
Expand All @@ -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<int>& conn, const int pid)
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -1007,6 +1016,11 @@ Array<double> VtkVtuData::get_points()
return points_array;
}

int VtkVtuData::elem_type()
{
return impl->elem_type;
}

int VtkVtuData::num_elems()
{
return impl->num_elems;
Expand Down
3 changes: 3 additions & 0 deletions Code/Source/solver/VtkData.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ class VtkData {
virtual Array<int> get_connectivity() = 0;
virtual Array<double> 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;
Expand Down Expand Up @@ -79,6 +80,7 @@ class VtkVtpData : public VtkData {

virtual Array<int> get_connectivity();
virtual Array<double> get_points();
virtual int elem_type();
virtual int num_elems();
virtual int np_elem();
virtual int num_points();
Expand Down Expand Up @@ -115,6 +117,7 @@ class VtkVtuData : public VtkData {
~VtkVtuData();

virtual Array<int> get_connectivity();
virtual int elem_type();
virtual int num_elems();
virtual int np_elem();
virtual int num_points();
Expand Down
42 changes: 41 additions & 1 deletion Code/Source/solver/consts.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@

namespace consts {


/// @brief Reproduces the 'SELECT CASE (lM%eType)' statements in the
/// Fortran 'PARTMSH' subroutine.
//
Expand All @@ -53,6 +52,47 @@ const std::map<ElementType,int> element_type_to_elem_nonb =
{ElementType::WDG, 3 }
};

/// @brief Map containing the element type as a string.
///
const std::map<ElementType,std::string> 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<ElementType,int> 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:
Expand Down
2 changes: 2 additions & 0 deletions Code/Source/solver/consts.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,9 @@ enum class ElementType
NRB = 115
};

extern const std::map<ElementType,std::string> element_type_to_string;
extern const std::map<ElementType,int> element_type_to_elem_nonb;
extern const std::map<ElementType,int> element_dimension;

// Template for printing ElementType.
/*
Expand Down
87 changes: 57 additions & 30 deletions Code/Source/solver/load_msh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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<int> 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<int> 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<mesh.nFa; i++){
auto &face = mesh.fa[i];
nn::select_eleb(simulation, mesh, face);
Expand Down
5 changes: 1 addition & 4 deletions Code/Source/solver/nn.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -942,7 +942,7 @@ void gn_nxx(const int l, const int eNoN, const int nsd, const int insd, Array<do
//
void select_ele(const ComMod& com_mod, mshType& mesh)
{
// Set integration dimension.
// Set integration dimension: shell, fiber or solid (2d or 3d)
int insd;
if (mesh.lShl) {
insd = com_mod.nsd - 1;
Expand All @@ -952,9 +952,6 @@ void select_ele(const ComMod& com_mod, mshType& mesh)
insd = com_mod.nsd;
}

if (mesh.eType == ElementType::NRB) {
}

// Set element properties based on integration dimension
// and number of element nodes.
//
Expand Down
9 changes: 0 additions & 9 deletions Code/Source/solver/nn_elem_props.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,15 +38,6 @@ using SetElementPropsMapType = std::map<int, std::function<void(int, mshType&)>>
//
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;
Expand Down
7 changes: 4 additions & 3 deletions Code/Source/solver/vtk_xml.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.");
}
Expand Down Expand Up @@ -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.
//
Expand All @@ -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.");
}
Expand All @@ -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;
Expand All @@ -602,7 +604,6 @@ void read_vtu(const std::string& file_name, mshType& mesh)
#endif
}


//----------
// read_precomputed_solution_vtu
//----------
Expand Down
Loading

0 comments on commit c3fac1e

Please sign in to comment.