Skip to content

Commit

Permalink
0 indexed meshes(#114) (#120)
Browse files Browse the repository at this point in the history 
* adding darcy equation

* darcy equation can be built but generates a trivially zero solution

* commenting files to find mesh indexing error

* correcting gN, gE, and IEN after reading meshes

* removing IEN resetting in vtk_xml.cpp since this functionality is moved into load_msh.cpp

* adding general parameters for unsteady darcy problem

* added final version of the hash map for gN, gE, and IEN setting for face meshes

* adding binary mask permutation for non-tet elements (#114)

* adding boundary node counter-clockwise ordering to element matching (#114)

* completed 2D darcy general model formulation for (#71)

* adding scientific notation and 16 precision to global node matching (#114)

* Revert "Merge branch 'Tissue-perfusion-model-#71' into 0-indexed-meshes(#114)"

This reverts commit 7cbf43a0b21cdb5706b5292f8c5f62cd8b0d4070, reversing
changes made to 7a97ef62e0b7b2e1d380230b081f62e508dce471.

* encapsulating hash map generation (#114)

* remove commented lines performing un-encapsulated hash map generation

* adding error checking for duplicate nodes during hash map generation (#114)

---------

Co-authored-by: Dave Parker <[email protected]>
  • Loading branch information
@zasexton @ktbolt
zasexton and ktbolt authored Nov 7, 2023
1 parent b128d28 commit 86f7f0c
Show file tree
Hide file tree
Showing 4 changed files with 224 additions and 89 deletions.
3 changes: 3 additions & 0 deletions Code/Source/svFSI/ComMod.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -863,6 +863,9 @@ class mshType
/// @brief IB: Mesh size parameter
double dx = 0.0;

/// @breif ordering: node ordering for boundaries
std::vector<std::vector<int>> ordering;

/// @brief Element distribution between processors
Vector<int> eDist;

Expand Down
258 changes: 177 additions & 81 deletions Code/Source/svFSI/load_msh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,71 @@
#include <iostream>
#include <fstream>
#include <sstream>
#include <chrono>
#include <unordered_map>
#include <string>
#include <iomanip>
#include <algorithm>
#include <math.h>

namespace load_msh {

/// @brief Generate hash maps for mesh nodes and elements

class MeshHashMaps {
public:
MeshHashMaps() {}

/// @brief Create a hash map for the global mesh nodes
///
/// \code {.cpp}
/// mesh
/// nsd
/// \endcode
///
std::unordered_map<std::string, int> createNodeHashMap(const mshType& mesh, const int& nsd) {
std::unordered_map<std::string, int> mesh_node_map;
for (int i = 0; i < mesh.gnNo; i++) {
std::ostringstream key;
key << std::scientific << std::setprecision(16);
for (int j = 0; j < nsd; j++) {
key << mesh.x(j,i) <<",";
}
mesh_node_map[key.str()] = i;
}
if (mesh_node_map.size() != mesh.gnNo) {
throw std::runtime_error("There may be a duplicate nodes within the mesh.");
}
return mesh_node_map;
}

/// @brief Create a hash map for the global mesh elements
///
/// \code {.cpp}
/// mesh
/// \endcode
///
std::unordered_map<std::string, int> createElementHashMap(const mshType& mesh) {
std::unordered_map<std::string, int> mesh_element_set;
for (int i = 0; i < mesh.gnEl; i++) {
for (unsigned int j = 0; j < mesh.ordering.size(); j++) {
std::vector<int> element_nodes;
for (unsigned int k = 0; k < mesh.ordering[j].size(); k++) {
element_nodes.push_back(mesh.gIEN(mesh.ordering[j][k], i));
}
std::sort(element_nodes.begin(), element_nodes.end());
std::string key = "";
for (int node: element_nodes) {
key += std::to_string(node) + ",";
}
mesh_element_set[key] = i;
}
}
return mesh_element_set;
}
};


#define ndbg_load_msh

/// @brief Read mesh position coordinates and connectivity.
Expand All @@ -54,7 +116,7 @@ void read_ccne(Simulation* simulation, mshType& mesh, const MeshParameters* mesh
{
auto mesh_path = mesh_param->mesh_file_path();
auto mesh_name = mesh_param->name();
throw std::runtime_error("[read_ccne] read_ccne() is not implemented.");
throw std::runtime_error("[read_ccne] read_ccne() is not implemented.");
}

/// @brief Read list of end nodes from a file into the face data structure.
Expand Down Expand Up @@ -113,93 +175,127 @@ void read_ndnlff(const std::string& file_name, faceType& face)
///
/// SUBROUTINE READSV(list, lM)
//
void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_param)
{
auto mesh_path = mesh_param->mesh_file_path();
auto mesh_name = mesh_param->get_name();
#define n_dbg_read_sv
#ifdef dbg_read_sv
DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm());
dmsg.banner();
dmsg << "Mesh name: " << mesh_name;
dmsg << "Mesh path: " << mesh_path;
#endif

// Read in volume mesh.
vtk_xml::read_vtu(mesh_path, mesh);

// Set mesh element properites for the input element type.
nn::select_ele(simulation->com_mod, mesh);

// Check the mesh element node ordering.
//
// Note: This may change element node ordering.
//
auto& com_mod = simulation->get_com_mod();
if (com_mod.ichckIEN) {
read_msh_ns::check_ien(simulation, mesh);
}
void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_param) {
auto mesh_path = mesh_param->mesh_file_path();
auto mesh_name = mesh_param->get_name();
#define n_dbg_read_sv
#ifdef dbg_read_sv
DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm());
dmsg.banner();
dmsg << "Mesh name: " << mesh_name;
dmsg << "Mesh path: " << mesh_path;
#endif

// Read face meshes.
//
// Creates nodal coordinates and element connectivity data.
//
mesh.nFa = mesh_param->face_parameters.size();
mesh.fa.resize(mesh.nFa);
// Read in volume mesh.
vtk_xml::read_vtu(mesh_path, mesh);

if (mesh.lFib && (mesh.nFa > 1)) {
throw std::runtime_error("There are " + std::to_string(mesh.nFa) + " faces defined for the '" +
mesh.name + "' mesh. Only one face is allowed for a 1D fiber-based mesh.");
}
// Set mesh element properites for the input element type.
nn::select_ele(simulation->com_mod, mesh);

for (int i = 0; i < mesh.nFa; i++) {
auto face_param = mesh_param->face_parameters[i];
auto& face = mesh.fa[i];
face.name = face_param->name();
#ifdef dbg_read_sv
dmsg << "face.name: " << face.name;
#endif

face.qmTRI3 = face_param->quadrature_modifier_TRI3();
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].");
}
// Check the mesh element node ordering.
//
// Note: This may change element node ordering.
//
auto &com_mod = simulation->get_com_mod();
if (com_mod.ichckIEN) {
read_msh_ns::check_ien(simulation, mesh);
}
// Read face meshes.
//
// Creates nodal coordinates and element connectivity data.
//
mesh.nFa = mesh_param->face_parameters.size();
mesh.fa.resize(mesh.nFa);

if (mesh.lFib) {
auto face_path = face_param->end_nodes_face_file_path();
#ifdef dbg_read_sv
dmsg << "Read end nodes face file ... " << " ";
dmsg << "face_path: " << face_path;
#endif
if (face_path == "") {
throw std::runtime_error("No end nodes face file path provided.");
}

read_ndnlff(face_path, face);

} else {
auto face_path = face_param->face_file_path();
vtk_xml::read_vtp(face_path, face);

// If node IDs were not read then create them.
if (face.gN.size() == 0) {
read_msh_ns::calc_nbc(mesh, face);

// Reset the connecttivity with the new node IDs?
for (int e = 0; e < face.nEl; e++) {
for (int a = 0; a < face.eNoN; a++) {
int Ac = face.IEN(a,e);
Ac = face.gN(Ac);
face.IEN(a,e) = Ac;
}
if (mesh.lFib && (mesh.nFa > 1)) {
throw std::runtime_error("There are " + std::to_string(mesh.nFa) + " faces defined for the '" +
mesh.name + "' mesh. Only one face is allowed for a 1D fiber-based mesh.");
}
}
}

nn::select_eleb(simulation, mesh, face);
}
for (int i = 0; i < mesh.nFa; i++) {
auto face_param = mesh_param->face_parameters[i];
auto &face = mesh.fa[i];
face.name = face_param->name();
#ifdef dbg_read_sv
dmsg << "face.name: " << face.name;
#endif

}
face.qmTRI3 = face_param->quadrature_modifier_TRI3();
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].");
}

if (mesh.lFib) {
auto face_path = face_param->end_nodes_face_file_path();
#ifdef dbg_read_sv
dmsg << "Read end nodes face file ... " << " ";
dmsg << "face_path: " << face_path;
#endif
if (face_path == "") {
throw std::runtime_error("No end nodes face file path provided.");
}

read_ndnlff(face_path, face);

} else {
auto face_path = face_param->face_file_path();
vtk_xml::read_vtp(face_path, face);

// If node IDs were not read then create them.
if (face.gN.size() == 0) {
read_msh_ns::calc_nbc(mesh, face);
// Reset the connectivity with the new node IDs?
for (int e = 0; e < face.nEl; e++) {
for (int a = 0; a < face.eNoN; a++) {
int Ac = face.IEN(a, e);
Ac = face.gN(Ac);
face.IEN(a, e) = Ac;
}
}
}
}
}
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);

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 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);
}
}
};

17 changes: 9 additions & 8 deletions Code/Source/svFSI/vtk_xml.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -450,15 +450,16 @@ void read_vtp(const std::string& file_name, faceType& face)

if (face.gN.size() == 0) {
std::cout << "[WARNING] No node IDs found in the '" << file_name << "' face file.";
} else {
for (int e = 0; e < face.nEl; e++) {
for (int a = 0; a < face.eNoN; a++) {
int Ac = face.IEN(a,e);
Ac = face.gN(Ac);
face.IEN(a,e) = Ac;
}
}
}
//else {
// for (int e = 0; e < face.nEl; e++) {
// for (int a = 0; a < face.eNoN; a++) {
//int Ac = face.IEN(a,e);
//Ac = face.gN(Ac);
//face.IEN(a,e) = Ac;
// }
// }
//}

// Create essential BC array.
//
Expand Down
35 changes: 35 additions & 0 deletions Code/Source/svFSI/vtk_xml_parser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@

#include <string>
#include <map>
#include <vector>

namespace vtk_xml_parser {

Expand All @@ -68,6 +69,32 @@ std::map<unsigned char,int> vtk_cell_to_elem {
{VTK_WEDGE, 6}
};

std::map<unsigned char, std::vector<std::vector<int>>> vtk_cell_ordering {
{VTK_HEXAHEDRON, {{0,3,2,1},
{4,5,6,7},
{0,1,5,4},
{1,2,6,5},
{2,3,7,6},
{3,0,4,7}}},
{VTK_LINE, {{0},
{1}}},
{VTK_QUAD, {{0,1},
{1,2},
{2,3},
{3,0}}},
{VTK_TETRA, {{0,1,2},
{0,1,3},
{1,2,3},
{2,0,3}}},
{VTK_TRIANGLE, {{0,1},
{1,2},
{2,0}}},
{VTK_WEDGE, {{0,1,2},
{3,4,5},
{0,1,4,3},
{1,2,5,4},
{2,0,3,5}}}
};
/// Names of data arrays store in VTK mesh files.
const std::string NODE_IDS_NAME("GlobalNodeID");
const std::string ELEMENT_IDS_NAME("GlobalElementID");
Expand Down Expand Up @@ -224,18 +251,25 @@ void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType&
#endif

int np_elem = 0;
std::vector<std::vector<int>> ordering;
if (num_line != 0) {
np_elem = vtk_cell_to_elem[VTK_LINE];
ordering = vtk_cell_ordering[VTK_LINE];
} if (num_hex != 0) {
np_elem = vtk_cell_to_elem[VTK_HEXAHEDRON];
ordering = vtk_cell_ordering[VTK_HEXAHEDRON];
} if (num_quad != 0) {
np_elem = vtk_cell_to_elem[VTK_QUAD];
ordering = vtk_cell_ordering[VTK_QUAD];
} if (num_tet != 0) {
np_elem = vtk_cell_to_elem[VTK_TETRA];
ordering = vtk_cell_ordering[VTK_TETRA];
} if (num_tri != 0) {
np_elem = vtk_cell_to_elem[VTK_TRIANGLE];
ordering = vtk_cell_ordering[VTK_TRIANGLE];
} if (num_wedge != 0) {
np_elem = vtk_cell_to_elem[VTK_WEDGE];
ordering = vtk_cell_ordering[VTK_WEDGE];
}

// For higher-order elements with mid-side nodes.
Expand All @@ -249,6 +283,7 @@ void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType&
mesh.gnEl = num_elems;
mesh.eNoN = np_elem;
mesh.gIEN = Array<int>(np_elem, num_elems);
mesh.ordering = ordering;

#ifdef debug_store_element_conn
std::cout << "[store_element_conn] np_elem: " << np_elem << std::endl;
Expand Down

0 comments on commit 86f7f0c

Please sign in to comment.