From 470872c399bfed6706282a0c6cb0ecd61398facf Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 14:29:13 -0600 Subject: [PATCH 1/8] add node section overwrite --- src/deck.cpp | 116 ++++++++++++++++++++++++++++++- src/lsdyna_mesh_reader/_deck.pyi | 4 ++ src/lsdyna_mesh_reader/deck.py | 63 ++++++++++++++++- tests/test_deck.py | 16 +++++ 4 files changed, 193 insertions(+), 6 deletions(-) diff --git a/src/deck.cpp b/src/deck.cpp index 70ef831..4d03851 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -353,22 +353,55 @@ static inline void ans_strtod(char *raw, int fltsz, #endif } +// FORTRAN-like scientific notation string formatting +void FormatWithExp(char *buffer, size_t buffer_size, double value, int width, + int precision, int num_exp) { + // Format the value using snprintf with given width and precision + snprintf(buffer, buffer_size, "% *.*E", width, precision, value); + + // Find the 'E' in the output, which marks the beginning of the exponent. + char *exponent_pos = strchr(buffer, 'E'); + + // Check the current length of the exponent (after 'E+') + int exp_length = strlen(exponent_pos + 2); + + // If the exponent is shorter than the desired number of digits + if (exp_length < num_exp) { + // Shift the exponent one place to the right and prepend a '0' + for (int i = exp_length + 1; i > 0; i--) { + exponent_pos[2 + i] = exponent_pos[1 + i]; + } + exponent_pos[2] = '0'; + + // Shift the entire buffer to the left to remove the space added by snprintf + memmove( + buffer, buffer + 1, + strlen( + buffer)); // length of buffer not recalculated, using previous value + } +} + struct NodeSection { NDArray nid; NDArray coord; NDArray tc; NDArray rc; int n_nodes = 0; + int fpos = 0; // Default constructor NodeSection() {} NodeSection(std::vector nid_vec, std::vector nodes_vec, - std::vector tc_vec, std::vector rc_vec) { + std::vector tc_vec, std::vector rc_vec, + int file_position) { n_nodes = nid_vec.size(); std::array nid_shape = {static_cast(n_nodes)}; std::array coord_shape = {static_cast(n_nodes), 3}; + // store file position where node block began + fpos = file_position; + // Wrap vectors as NDArrays nid = WrapVectorAsNDArray(std::move(nid_vec), nid_shape); coord = WrapVectorAsNDArray(std::move(nodes_vec), coord_shape); @@ -681,6 +714,8 @@ class Deck { std::vector rc; rc.reserve(NNUM_RESERVE); + int start_pos = memmap.tellg(); + int index = 0; while (memmap[0] != '*') { @@ -743,7 +778,7 @@ class Deck { memmap.seek_eol(); } - NodeSection *node_section = new NodeSection(nid, coord, tc, rc); + NodeSection *node_section = new NodeSection(nid, coord, tc, rc, start_pos); node_sections.push_back(*node_section); return; @@ -862,6 +897,78 @@ class Deck { int ReadLine() { return memmap.read_line(); } }; +void OverwriteNodeSection(const char *filename, int fpos, + const NDArray coord_arr) { + + // Open the file with read and write permissions + FILE *fp = fopen(filename, "r+b"); + if (!fp) { + throw std::runtime_error("Cannot open file for reading and writing."); + } + + // Seek to the start position of the node section + if (fseeko(fp, fpos, SEEK_SET) != 0) { + fclose(fp); + throw std::runtime_error( + "Cannot seek to the start position of the node section."); + } + + int node_idx = 0; + off_t line_start_pos; + char line[256]; + + int n_coord = coord_arr.shape(0); + const double *coord = coord_arr.data(); + + // Prepare buffers for formatted coordinates + char coord_str[17] = {0}; + + while (fgets(line, sizeof(line), fp)) { + // Record the position of the beginning of the line + line_start_pos = ftello(fp) - strlen(line); + + // Skip comment lines + if (line[0] == '$') { + continue; + } + + // Check for end of node section + if (line[0] == '*') { + break; + } + + // Check if we have more nodes to process + if (node_idx >= n_coord) { + break; + } + + // Ensure the line is long enough + size_t line_length = strlen(line); + if (line_length < 56) { + throw std::runtime_error("Insufficient line length."); + } + + // Format the coordinates + for (int ii = 0; ii < 3; ii++) { + // Generate the formatted string + FormatWithExp(coord_str, 17, coord[node_idx * 3 + ii], 16, 9, 2); + memcpy(&line[8 + ii * 16], coord_str, 16); // coordinate field + } + + // Seek back to the beginning of the line and write the updated line + fseeko(fp, line_start_pos, SEEK_SET); + if (fwrite(line, 1, line_length, fp) != line_length) { + fclose(fp); + throw std::runtime_error("Failed to write modified line to file."); + } + + // Move to the next node + node_idx++; + } + + fclose(fp); +} + NB_MODULE(_deck, m) { nb::class_(m, "NodeSection") .def(nb::init()) @@ -870,7 +977,8 @@ NB_MODULE(_deck, m) { .def_ro("coordinates", &NodeSection::coord, nb::rv_policy::automatic) .def_ro("nid", &NodeSection::nid, nb::rv_policy::automatic) .def_ro("tc", &NodeSection::tc, nb::rv_policy::automatic) - .def_ro("rc", &NodeSection::rc, nb::rv_policy::automatic); + .def_ro("rc", &NodeSection::rc, nb::rv_policy::automatic) + .def_ro("fpos", &NodeSection::fpos); nb::class_(m, "ElementSolidSection") .def(nb::init()) @@ -906,4 +1014,6 @@ NB_MODULE(_deck, m) { .def("read_element_solid_section", &Deck::ReadElementSolidSection) .def("read_element_shell_section", &Deck::ReadElementShellSection) .def("read_node_section", &Deck::ReadNodeSection); + + m.def("overwrite_node_section", &OverwriteNodeSection); } diff --git a/src/lsdyna_mesh_reader/_deck.pyi b/src/lsdyna_mesh_reader/_deck.pyi index 4a18df5..2b27e92 100644 --- a/src/lsdyna_mesh_reader/_deck.pyi +++ b/src/lsdyna_mesh_reader/_deck.pyi @@ -21,6 +21,8 @@ class NodeSection: @property def rc(self) -> IntArray: ... def __len__(self) -> int: ... + @property + def fpos(self) -> int: ... class ElementSection: def __init__(self) -> None: ... @@ -51,3 +53,5 @@ class _Deck: def read_element_shell_section(self) -> None: ... def read_node_section(self) -> None: ... def read(self) -> None: ... + +def overwrite_node_section(filename: str, fpos: int, nodes: FloatArray2D) -> None: ... diff --git a/src/lsdyna_mesh_reader/deck.py b/src/lsdyna_mesh_reader/deck.py index 2741b20..17e93b4 100644 --- a/src/lsdyna_mesh_reader/deck.py +++ b/src/lsdyna_mesh_reader/deck.py @@ -1,10 +1,18 @@ +import os +import shutil from pathlib import Path from typing import TYPE_CHECKING, List, Union import numpy as np from numpy.typing import NDArray -from lsdyna_mesh_reader._deck import ElementShellSection, ElementSolidSection, NodeSection, _Deck +from lsdyna_mesh_reader._deck import ( + ElementShellSection, + ElementSolidSection, + NodeSection, + _Deck, + overwrite_node_section, +) if TYPE_CHECKING: try: @@ -19,7 +27,7 @@ class Deck: Parameters ---------- filename : str | pathlib.Path - Path to the keyword file (`*.k`, `*.key`, `*.dyn`). + Path to the keyword file (``*.k``, ``*.key``, ``*.dyn``). Examples -------- @@ -35,8 +43,12 @@ class Deck: def __init__(self, filename: Union[str, Path]) -> None: """Initialize the deck object.""" - self._deck = _Deck(str(filename)) + filename = str(filename) + if not os.path.isfile(filename): + raise FileNotFoundError(f"Invalid file or unable to locate {filename}") + self._deck = _Deck(filename) self._deck.read() + self._filename = filename @property def element_solid_sections(self) -> List[ElementSolidSection]: @@ -235,6 +247,51 @@ def to_grid(self) -> "UnstructuredGrid": return grid + def overwrite_node_section(self, filename: Union[str, Path], nodes: NDArray[float]) -> None: + """ + Create a new deck file with the same content but overwritten node section. + + Parameters + ---------- + filename : str | pathlib.Path + Path to the new file. + nodes : np.ndarray[float] + New node coordinates to write. Number of nodes must match the + number of nodes in the node section. + + Notes + ----- + Overwrites the first node section. + + Examples + -------- + Load an existing LS-DYNA deck and write new nodes to it in a new file. + + >>> import numpy as np + >>> import lsdyna_mesh_reader + >>> from lsdyna_mesh_reader import examples + >>> deck = lsdyna_mesh_reader.Deck(examples.birdball) + >>> nodes = deck.node_sections[0].coordinates + >>> new_nodes = nodes + np.random.random(nodes.shape) + >>> deck.overwrite_node_section("new_deck.k", new_nodes) + + """ + if not self.node_sections: + raise RuntimeError("This deck is missing a node section") + nsec = self.node_sections[0] + + if nodes.ndim != 2 or nodes.shape[1] != 3: + raise ValueError(f"Expected `nodes` to contain 3D coordinates, not {nodes.shape}") + elif nodes.shape[0] != nsec.nid.size: + raise RuntimeError( + f"Number of coordinates in `nodes` ({nodes.shape[0]}) must match number of " + f"nodes in the node section ({nsec.nid.size})" + ) + + filename = str(filename) + shutil.copy(self._filename, filename) + overwrite_node_section(filename, nsec.fpos, nodes) + def __repr__(self) -> str: lines = ["LSDYNA Deck with:"] lines.append(f" Node sections: {len(self.node_sections)}") diff --git a/tests/test_deck.py b/tests/test_deck.py index 575927c..38a7aba 100644 --- a/tests/test_deck.py +++ b/tests/test_deck.py @@ -203,3 +203,19 @@ def test_element_tshell() -> None: assert len(node_section) == 324 grid = deck.to_grid() + + +def test_overwrite_node_section(tmp_path: Path) -> None: + filename = str(tmp_path / "tmp.k") + with open(filename, "w") as fid: + fid.write(NODE_SECTION) + + deck = lsdyna_mesh_reader.Deck(filename) + + new_filename = str(tmp_path / "tmp_overwrite.k") + node_section = deck.node_sections[0] + new_nodes = node_section.coordinates + np.random.random(node_section.coordinates.shape) + deck.overwrite_node_section(new_filename, new_nodes) + + deck_new = lsdyna_mesh_reader.Deck(new_filename) + assert np.allclose(deck_new.node_sections[0].coordinates, new_nodes) From f1931f7a060abcd9299e378ddff8009cf8af22a3 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 14:50:48 -0600 Subject: [PATCH 2/8] win compat --- src/deck.cpp | 4 ++++ src/lsdyna_mesh_reader/deck.py | 4 +++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/deck.cpp b/src/deck.cpp index 4d03851..ab69b53 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -22,7 +22,11 @@ using namespace nb::literals; #if defined(_WIN32) || defined(_WIN64) #define NOMINMAX #define strtok_r strtok_s +#include #include +#define fseeko _fseeki64 +#define ftello _ftelli64 +typedef __int64 off_t; #else #include #include diff --git a/src/lsdyna_mesh_reader/deck.py b/src/lsdyna_mesh_reader/deck.py index 17e93b4..ca07c27 100644 --- a/src/lsdyna_mesh_reader/deck.py +++ b/src/lsdyna_mesh_reader/deck.py @@ -247,7 +247,9 @@ def to_grid(self) -> "UnstructuredGrid": return grid - def overwrite_node_section(self, filename: Union[str, Path], nodes: NDArray[float]) -> None: + def overwrite_node_section( + self, filename: Union[str, Path], nodes: NDArray[np.float64] + ) -> None: """ Create a new deck file with the same content but overwritten node section. From 6b9890effa8912de98a2dc9739a2f3e454688a5b Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 14:58:42 -0600 Subject: [PATCH 3/8] just use seek --- src/deck.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/deck.cpp b/src/deck.cpp index ab69b53..a2a27b5 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -22,11 +22,7 @@ using namespace nb::literals; #if defined(_WIN32) || defined(_WIN64) #define NOMINMAX #define strtok_r strtok_s -#include #include -#define fseeko _fseeki64 -#define ftello _ftelli64 -typedef __int64 off_t; #else #include #include @@ -911,7 +907,7 @@ void OverwriteNodeSection(const char *filename, int fpos, } // Seek to the start position of the node section - if (fseeko(fp, fpos, SEEK_SET) != 0) { + if (fseek(fp, fpos, SEEK_SET) != 0) { fclose(fp); throw std::runtime_error( "Cannot seek to the start position of the node section."); @@ -960,7 +956,7 @@ void OverwriteNodeSection(const char *filename, int fpos, } // Seek back to the beginning of the line and write the updated line - fseeko(fp, line_start_pos, SEEK_SET); + fseek(fp, line_start_pos, SEEK_SET); if (fwrite(line, 1, line_length, fp) != line_length) { fclose(fp); throw std::runtime_error("Failed to write modified line to file."); From 64c82fa292b31d0bc63b30fe25098cee32a08dbc Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 15:02:20 -0600 Subject: [PATCH 4/8] forgot ftell --- src/deck.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/deck.cpp b/src/deck.cpp index a2a27b5..d02205a 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -925,7 +925,7 @@ void OverwriteNodeSection(const char *filename, int fpos, while (fgets(line, sizeof(line), fp)) { // Record the position of the beginning of the line - line_start_pos = ftello(fp) - strlen(line); + line_start_pos = ftell(fp) - strlen(line); // Skip comment lines if (line[0] == '$') { From 91a898103bf5f9ce126db68c4c3552ec02ad7390 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 15:12:21 -0600 Subject: [PATCH 5/8] read and write in binary --- src/deck.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/deck.cpp b/src/deck.cpp index d02205a..63bc133 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -901,7 +901,7 @@ void OverwriteNodeSection(const char *filename, int fpos, const NDArray coord_arr) { // Open the file with read and write permissions - FILE *fp = fopen(filename, "r+b"); + FILE *fp = fopen(filename, "rb+"); if (!fp) { throw std::runtime_error("Cannot open file for reading and writing."); } @@ -957,10 +957,7 @@ void OverwriteNodeSection(const char *filename, int fpos, // Seek back to the beginning of the line and write the updated line fseek(fp, line_start_pos, SEEK_SET); - if (fwrite(line, 1, line_length, fp) != line_length) { - fclose(fp); - throw std::runtime_error("Failed to write modified line to file."); - } + fwrite(line, 1, line_length, fp); // Move to the next node node_idx++; From 608b0a4eba91e1fec63a85101977dc591c89d3db Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 14:59:37 -0700 Subject: [PATCH 6/8] seek to eol after writing line --- src/deck.cpp | 4 +++- tests/test_deck.py | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/deck.cpp b/src/deck.cpp index 63bc133..a322b9f 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -925,7 +925,8 @@ void OverwriteNodeSection(const char *filename, int fpos, while (fgets(line, sizeof(line), fp)) { // Record the position of the beginning of the line - line_start_pos = ftell(fp) - strlen(line); + int line_end_pos = ftell(fp); + line_start_pos = line_end_pos - strlen(line); // Skip comment lines if (line[0] == '$') { @@ -958,6 +959,7 @@ void OverwriteNodeSection(const char *filename, int fpos, // Seek back to the beginning of the line and write the updated line fseek(fp, line_start_pos, SEEK_SET); fwrite(line, 1, line_length, fp); + fseek(fp, line_start_pos, SEEK_SET); // Move to the next node node_idx++; diff --git a/tests/test_deck.py b/tests/test_deck.py index 38a7aba..acf2a6a 100644 --- a/tests/test_deck.py +++ b/tests/test_deck.py @@ -219,3 +219,4 @@ def test_overwrite_node_section(tmp_path: Path) -> None: deck_new = lsdyna_mesh_reader.Deck(new_filename) assert np.allclose(deck_new.node_sections[0].coordinates, new_nodes) + assert np.allclose(deck_new.node_sections[0].nid, deck.node_sections[0].nid) From 1cc9a175bb971bf898a67561c265e2096e98111e Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 16:02:50 -0600 Subject: [PATCH 7/8] turns out seek is only necessary for windows --- src/deck.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/deck.cpp b/src/deck.cpp index a322b9f..5393b83 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -959,7 +959,9 @@ void OverwriteNodeSection(const char *filename, int fpos, // Seek back to the beginning of the line and write the updated line fseek(fp, line_start_pos, SEEK_SET); fwrite(line, 1, line_length, fp); +#if defined(_WIN32) || defined(_WIN64) fseek(fp, line_start_pos, SEEK_SET); +#endif // Move to the next node node_idx++; From 0ad04a75fc3e2807d6b20f1547d912b77f4348b9 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 25 Sep 2024 15:06:48 -0700 Subject: [PATCH 8/8] should be line_end_pos --- src/deck.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/deck.cpp b/src/deck.cpp index 5393b83..822fe32 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -960,7 +960,7 @@ void OverwriteNodeSection(const char *filename, int fpos, fseek(fp, line_start_pos, SEEK_SET); fwrite(line, 1, line_length, fp); #if defined(_WIN32) || defined(_WIN64) - fseek(fp, line_start_pos, SEEK_SET); + fseek(fp, line_end_pos, SEEK_SET); #endif // Move to the next node