Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 16 additions & 50 deletions src/core/io/src/4C_io_exodus.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,8 @@ namespace
/*----------------------------------------------------------------------*
| ctor (public) maf 12/07|
*----------------------------------------------------------------------*/
Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
Core::IO::Exodus::Mesh::Mesh(std::filesystem::path exodus_file, MeshParameters mesh_data)
: mesh_data_(mesh_data)
{
int error;
int CPU_word_size, IO_word_size;
Expand All @@ -118,8 +119,8 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)

// open EXODUS II file
int exo_handle =
ex_open(exofilename.c_str(), EX_READ, &CPU_word_size, &IO_word_size, &exoversion);
if (exo_handle <= 0) FOUR_C_THROW("Error while opening EXODUS II file {}", exofilename);
ex_open(exodus_file.c_str(), EX_READ, &CPU_word_size, &IO_word_size, &exoversion);
if (exo_handle <= 0) FOUR_C_THROW("Error while opening EXODUS II file {}", exodus_file.string());

// read database parameters
int num_elem_blk, num_node_sets, num_side_sets, num_nodes;
Expand All @@ -138,10 +139,12 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
error = ex_get_coord(exo_handle, x.data(), y.data(), z.data());
if (error != 0) FOUR_C_THROW("exo error returned");

FOUR_C_ASSERT_ALWAYS(
mesh_data_.node_start_id >= 0, "Node start id must be greater than or equal to 0");

for (int i = 0; i < num_nodes; ++i)
{
// Node numbering starts at one in 4C
nodes_[i + 1] = {x[i], y[i], z[i]};
nodes_[i + mesh_data_.node_start_id] = {x[i], y[i], z[i]};
}
}

Expand Down Expand Up @@ -175,13 +178,16 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
error = ex_get_conn(exo_handle, EX_ELEM_BLOCK, ebids[i], allconn.data(), nullptr, nullptr);
if (error != 0) FOUR_C_THROW("exo error returned");
std::map<int, std::vector<int>> eleconn;

// Compare the desired start ID to Exodus' one-based indexing to get the offset.
const int node_offset = mesh_data_.node_start_id - 1;
for (int j = 0; j < num_el_in_blk; ++j)
{
std::vector<int>& actconn = eleconn[j];
actconn.reserve(num_nod_per_elem);
for (int k = 0; k < num_nod_per_elem; ++k)
{
actconn.push_back(allconn[k + j * num_nod_per_elem]);
actconn.push_back(allconn[k + j * num_nod_per_elem] + node_offset);
}
reorder_nodes_exodus_to_four_c(actconn, cell_type_from_exodus_string(ele_type));
}
Expand Down Expand Up @@ -217,10 +223,12 @@ Core::IO::Exodus::Mesh::Mesh(const std::string exofilename)
else if (error < 0)
FOUR_C_THROW("error reading node set");
std::set<int> nodes_in_set;
for (int j = 0; j < num_nodes_in_set; ++j) nodes_in_set.insert(node_set_node_list[j]);
// Compare the desired start ID to Exodus' one-based indexing to get the offset.
const int node_offset = mesh_data_.node_start_id - 1;
for (int j = 0; j < num_nodes_in_set; ++j)
nodes_in_set.insert(node_set_node_list[j] + node_offset);
NodeSet actNodeSet(nodes_in_set, nodesetname);

// Add this NodeSet into Mesh map (here prelim due to pro names)
node_sets_.insert(std::pair<int, NodeSet>(npropID[i], actNodeSet));
}
} // end of nodeset section
Expand Down Expand Up @@ -362,48 +370,6 @@ void Core::IO::Exodus::Mesh::set_nsd(const int nsd)

four_c_dim_ = nsd;
}
/*----------------------------------------------------------------------*/
/*----------------------------------------------------------------------*/


/*----------------------------------------------------------------------*/
/*----------------------------------------------------------------------*/
std::vector<double> Core::IO::Exodus::Mesh::normal(
const int head1, const int origin, const int head2) const
{
std::vector<double> normal(3);
std::vector<double> h1 = get_node(head1);
std::vector<double> h2 = get_node(head2);
std::vector<double> o = get_node(origin);

normal[0] = ((h1[1] - o[1]) * (h2[2] - o[2]) - (h1[2] - o[2]) * (h2[1] - o[1]));
normal[1] = -((h1[0] - o[0]) * (h2[2] - o[2]) - (h1[2] - o[2]) * (h2[0] - o[0]));
normal[2] = ((h1[0] - o[0]) * (h2[1] - o[1]) - (h1[1] - o[1]) * (h2[0] - o[0]));

double length = sqrt(normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]);
normal[0] = normal[0] / length;
normal[1] = normal[1] / length;
normal[2] = normal[2] / length;

return normal;
}

/*----------------------------------------------------------------------*/
/*----------------------------------------------------------------------*/
std::vector<double> Core::IO::Exodus::Mesh::node_vec(const int tail, const int head) const
{
std::vector<double> nv(3);
std::vector<double> t = get_node(tail);
std::vector<double> h = get_node(head);
nv[0] = h[0] - t[0];
nv[1] = h[1] - t[1];
nv[2] = h[2] - t[2];
double length = sqrt(nv[0] * nv[0] + nv[1] * nv[1] + nv[2] * nv[2]);
nv[0] = nv[0] / length;
nv[1] = nv[1] / length;
nv[2] = nv[2] / length;
return nv;
}


/*----------------------------------------------------------------------*/
Expand Down
30 changes: 22 additions & 8 deletions src/core/io/src/4C_io_exodus.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@

#include "4C_fem_general_element.hpp"

#include <filesystem>
#include <set>
#include <string>
#include <vector>
Expand All @@ -25,14 +26,31 @@ namespace Core::IO::Exodus
class NodeSet;
class SideSet;

/**
* Additional parameters that are used in the constructor.
*/
struct MeshParameters
{
/**
* The ID of the first node in the mesh. This defaults to 1, since this is the default in
* the Exodus II mesh format.
*/
int node_start_id{1};
};

/**
* A class that stores the mesh information read from an Exodus file.
*/
class Mesh
{
public:
//! constructor
Mesh(std::string exofilename);
/**
* @brief Read a mesh from an Exodus file.
*
* Read the data from the Exodus file and store it in the class. The optional @p
* mesh_data can be used to set options documented in the MeshData struct.
*/
Mesh(std::filesystem::path exodus_file, MeshParameters mesh_data = {});

//! Print mesh info
void print(std::ostream& os, bool verbose = false) const;
Expand Down Expand Up @@ -76,12 +94,6 @@ namespace Core::IO::Exodus
//! Get one SideSet
SideSet get_side_set(const int id) const;

//! Get edge Normal at node
std::vector<double> normal(const int head1, const int origin, const int head2) const;

//! Get normalized Vector between 2 nodes
std::vector<double> node_vec(const int tail, const int head) const;

//! Get Node map
const std::map<int, std::vector<double>>& get_nodes() const { return nodes_; }

Expand All @@ -95,6 +107,8 @@ namespace Core::IO::Exodus
void set_nsd(const int nsd);

private:
MeshParameters mesh_data_;

std::map<int, std::vector<double>> nodes_;

std::map<int, ElementBlock> element_blocks_;
Expand Down
7 changes: 5 additions & 2 deletions src/core/io/src/4C_io_meshreader.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -822,8 +822,11 @@ namespace

const auto& geometry_data = data.group(exodus_reader.section_name);
const auto& exodus_file = geometry_data.get<std::filesystem::path>("FILE");
exodus_reader.mesh_on_rank_zero =
std::make_unique<Core::IO::Exodus::Mesh>(exodus_file.string());
exodus_reader.mesh_on_rank_zero = std::make_unique<Core::IO::Exodus::Mesh>(
exodus_file.string(), Core::IO::Exodus::MeshParameters{
// We internally depend on node numbers starting at 0.
.node_start_id = 0,
});
const auto& mesh = *exodus_reader.mesh_on_rank_zero;

// Initial implementation:
Expand Down
9 changes: 9 additions & 0 deletions src/core/io/tests/4C_io_exodus_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,4 +29,13 @@ namespace

mesh.print(std::cout, true);
}

TEST(Exodus, NodeOffset)
{
Core::IO::Exodus::Mesh mesh(TESTING::get_support_file_path("test_files/exodus/cube.exo"),
Core::IO::Exodus::MeshParameters{.node_start_id = 100});
EXPECT_EQ(mesh.get_node(100), (std::vector<double>{-5.0, 0.0, 0.0}));
EXPECT_EQ(mesh.get_element_block(1).get_ele_nodes(1),
(std::vector<int>{108, 100, 103, 109, 110, 104, 107, 111}));
}
} // namespace
4 changes: 2 additions & 2 deletions tests/input_files/contact3D_exodus_in.4C.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,13 @@ DESIGN SURF DIRICH CONDITIONS:
RESULT DESCRIPTION:
- STRUCTURE:
DIS: structure
NODE: 1270
NODE: 1269
QUANTITY: dispx
VALUE: -8.31974725251016095e-03
TOLERANCE: 1e-8
- STRUCTURE:
DIS: structure
NODE: 1372
NODE: 1371
QUANTITY: dispx
VALUE: -8.36015173005985428e-03
TOLERANCE: 1e-8
24 changes: 12 additions & 12 deletions tests/input_files/fluid_exodus_in.4C.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -65,73 +65,73 @@ DESIGN POINT DIRICH CONDITIONS:
RESULT DESCRIPTION:
- FLUID:
DIS: "fluid"
NODE: 2912
NODE: 2911
QUANTITY: "velx"
VALUE: 0
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 2912
NODE: 2911
QUANTITY: "vely"
VALUE: 0
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 2912
NODE: 2911
QUANTITY: "pressure"
VALUE: 2.316272235592883
TOLERANCE: 1e-10
- FLUID:
DIS: "fluid"
NODE: 4034
NODE: 4033
QUANTITY: "velx"
VALUE: 0
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 4034
NODE: 4033
QUANTITY: "vely"
VALUE: 0
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 4034
NODE: 4033
QUANTITY: "pressure"
VALUE: 1.5110014506671248
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 6667
NODE: 6666
QUANTITY: "velx"
VALUE: 0.8140348125434049
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 6667
NODE: 6666
QUANTITY: "vely"
VALUE: -0.022743581282554975
TOLERANCE: 1e-10
- FLUID:
DIS: "fluid"
NODE: 6667
NODE: 6666
QUANTITY: "pressure"
VALUE: 2.0581385928541125
TOLERANCE: 1e-10
- FLUID:
DIS: "fluid"
NODE: 6998
NODE: 6997
QUANTITY: "velx"
VALUE: 0.9708173796517694
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 6998
NODE: 6997
QUANTITY: "vely"
VALUE: -7.977045433472414e-16
TOLERANCE: 1e-08
- FLUID:
DIS: "fluid"
NODE: 6998
NODE: 6997
QUANTITY: "pressure"
VALUE: 1.8450141345624358
TOLERANCE: 1e-10