diff --git a/README.md b/README.md index ac4cf9a..0b1b168 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,20 @@ ## LS-DYNA Mesh Reader -This library can be used to read in LS-DYNA meshes stored within `*.k` files, -commonly known as keyword format "input decks". +This library can be used to read in LS-DYNA meshes stored within keyword +(`*.k`, `*.key`, `*.dyn`) files, also known as keyword format "input +decks". Many of these example files were obtained from the excellent documentation at [LS-DYNA Examples](https://www.dynaexamples.com/). +### Motivation + +Despite its popularity, there doesn't appear to be a reader for LS-DYNA keyword +files. I need a reader for a closed source project and hope that this helps +someone else who also wishes to read in these files. It borrows from +[mapdl-archive](https://github.com/akaszynski/mapdl-archive) as MAPDL also +follows many of the same FORTRAN conventions when writing out FEMs. + ## Installation Install the fully featured reader with visualization with: @@ -25,13 +34,22 @@ pip install lsdyna-mesh-reader ### Examples -Let's load the [Contact Eroding +Before going through a basic example, let's talk about how these "decks" are organized. Each keyword file contains "keywords" that describe the start of sections of "cards". This terminology dates back to when DYNA3D was developed in 1976 where programs were written on [punch cards](https://en.wikipedia.org/wiki/Punched_card). + +To read in nodes and elements, we have to read in one or more node and element sections, each starting with `*NODE` or `*ELEMENT_SOLID`. This library loads in those raw sections as well as parsed sections as a higher level abstraction. + +Let's start by loading the [Contact Eroding I](https://www.dynaexamples.com/introduction/intro-by-a.-tabiei/contact/contact-eroding-i) -example mesh and output the node coordinates as a numpy array +example mesh and getting the first node section. ```py + +Load the example file and get the first node section. + >>> import lsdyna_mesh_reader >>> from lsdyna_mesh_reader import examples >>> deck = lsdyna_mesh_reader.Deck(examples.birdball) ->>> deck.nodes +>>> node_section = deck.node_sections[0] +>>> node_section + ``` diff --git a/src/deck.cpp b/src/deck.cpp index 5b764c4..bc9c6a0 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -339,6 +340,41 @@ struct NodeSection { tc = WrapVectorAsNDArray(std::move(tc_vec), nnum_shape); rc = WrapVectorAsNDArray(std::move(rc_vec), nnum_shape); } + + std::string ToString() const { + std::ostringstream oss; + int num_nodes = nnum.size(); + if (num_nodes > 1) { + oss << "NodeSection containing " << num_nodes << " nodes\n\n"; + } else { + oss << "NodeSection containing " << num_nodes << " node\n\n"; + } + + // Header + oss << "| NID | X | Y | Z | tc | " + " rc |\n"; + oss << "|-------|---------------|---------------|---------------|--------|-" + "-------|\n"; + + // Output first 3 nodes in Fortran-like format + for (int i = 0; i < std::min(10, num_nodes); ++i) { + oss << std::setw(8) << nnum(i) << " " // Node ID (nnum) + << std::setw(15) << std::scientific << std::setprecision(8) + << nodes(i, 0) << " " // X + << std::setw(15) << std::scientific << std::setprecision(8) + << nodes(i, 1) << " " // Y + << std::setw(15) << std::scientific << std::setprecision(8) + << nodes(i, 2) << " " // Z + << std::setw(8) << tc(i) << " " // tc + << std::setw(8) << rc(i) << "\n"; // rc + } + + if (num_nodes > 10) { + oss << "..." << std::endl; + } + + return oss.str(); + } }; struct ElementSection { @@ -347,11 +383,14 @@ struct ElementSection { NDArray node_ids; NDArray node_id_offsets; + ElementSection() {} + ElementSection(std::vector eid_vec, std::vector pid_vec, std::vector node_ids_vec, std::vector node_id_offsets_vec) { std::array nel_shape = {static_cast(eid_vec.size())}; - std::array node_ids_shape = {static_cast(node_ids_vec.size())}; + std::array node_ids_shape = { + static_cast(node_id_offsets_vec.back())}; std::array node_ids_offsets_shape = { static_cast(eid_vec.size() + 1)}; @@ -361,6 +400,48 @@ struct ElementSection { node_id_offsets = WrapVectorAsNDArray(std::move(node_id_offsets_vec), node_ids_offsets_shape); } + + std::string ToString() const { + std::ostringstream oss; + int num_elements = eid.size(); + + if (num_elements > 1) { + oss << "ElementSection containing " << num_elements << " elements\n\n"; + } else { + oss << "ElementSection containing " << num_elements << " element\n\n"; + } + + // Header + oss << "| EID | PID | N1 | N2 | N3 | N4 | N5 | N6 | " + "N7 | N8 |\n"; + oss << "|-------|-------|-------|-------|-------|-------|-------|-------|--" + "-----|-------|\n"; + + // Output first 10 elements (or less) in Fortran-like format + for (int i = 0; i < std::min(10, num_elements); ++i) { + // Print element ID and part ID + oss << std::setw(8) << eid(i) << "" // EID + << std::setw(8) << pid(i) << ""; // PID + + // Retrieve node IDs associated with this element + int start = node_id_offsets(i); + int end = node_id_offsets(i + 1); + + // Print node IDs for the element + for (int j = start; j < end; ++j) { + oss << std::setw(8) << node_ids(j) << ""; // NODE_ID + } + + oss << "\n"; + } + + // Add "..." if there are more than 10 elements + if (num_elements > 10) { + oss << "...\n"; + } + + return oss.str(); + } }; class Deck { @@ -371,10 +452,6 @@ class Deck { MemoryMappedFile memmap; public: - // int n_nodes = 0; - // int n_elem = 0; - // NDArray nodes_arr; - std::vector node_sections; std::vector element_solid_sections; @@ -507,12 +584,21 @@ class Deck { // at the moment, we're making the assumption that the entire // element can be defined on this line eid.push_back(fast_atoi(memmap.current, 8)); + memmap += 8; pid.push_back(fast_atoi(memmap.current, 8)); + memmap += 8; // Is this always 8? for (int i = 0; i < 8; i++) { node_ids.push_back(fast_atoi(memmap.current, 8)); + memmap += 8; } + + // append end of node section (and start of next one) in our locator + node_id_offsets.push_back(node_ids.size()); + + // skip remainder of the line + memmap.seek_eol(); } ElementSection *element_section = @@ -526,6 +612,7 @@ class Deck { NB_MODULE(_deck, m) { nb::class_(m, "NodeSection") .def(nb::init()) + .def("__repr__", &NodeSection::ToString) .def_ro("nodes", &NodeSection::nodes, nb::rv_policy::automatic) .def_ro("nnum", &NodeSection::nnum, nb::rv_policy::automatic) .def_ro("tc", &NodeSection::tc, nb::rv_policy::automatic) @@ -533,6 +620,7 @@ NB_MODULE(_deck, m) { nb::class_(m, "ElementSection") .def(nb::init()) + .def("__repr__", &ElementSection::ToString) .def_ro("eid", &ElementSection::eid, nb::rv_policy::automatic) .def_ro("pid", &ElementSection::pid, nb::rv_policy::automatic) .def_ro("node_ids", &ElementSection::node_ids, nb::rv_policy::automatic) diff --git a/src/lsdyna_mesh_reader/_deck.pyi b/src/lsdyna_mesh_reader/_deck.pyi index 70d7a47..8dec18f 100644 --- a/src/lsdyna_mesh_reader/_deck.pyi +++ b/src/lsdyna_mesh_reader/_deck.pyi @@ -1 +1,41 @@ -# class Node +from typing import List + +import numpy as np +from numpy.typing import NDArray + +# shape typing added in 2.1, not adding it here yet +IntArray = NDArray[np.int32] +FloatArray1D = NDArray[np.float64] +FloatArray2D = NDArray[np.float64] + +class NodeSection: + def __init__(self) -> None: ... + @property + def nnum(self) -> IntArray: ... + @property + def nodes(self) -> FloatArray2D: ... + @property + def tc(self) -> IntArray: ... + @property + def rc(self) -> IntArray: ... + +class ElementSection: + def __init__(self) -> None: ... + @property + def eid(self) -> IntArray: ... + @property + def pid(self) -> IntArray: ... + @property + def node_ids(self) -> IntArray: ... + @property + def node_id_offsets(self) -> IntArray: ... + +class Deck: + def __init__(self, fname: str) -> None: ... + @property + def node_sections(self) -> List[NodeSection]: ... + @property + def element_solid_sections(self) -> List[ElementSection]: ... + def read_line(self) -> int: ... + def read_element_section(self) -> None: ... + def read_node_section(self) -> None: ... diff --git a/tests/test_deck.py b/tests/test_deck.py index 9906599..fd0da5f 100644 --- a/tests/test_deck.py +++ b/tests/test_deck.py @@ -23,22 +23,6 @@ ] ) - -def test_node_card(tmp_path: Path) -> None: - filename = str(tmp_path / "tmp.k") - with open(filename, "w") as fid: - fid.write(NODE_CARD_BLOCK) - - deck = Deck(filename) - deck.read_line() - deck.read_node_card() - - assert len(deck.node_cards) == 1 - node_card = deck.node_cards[0] - assert np.allclose(node_card.nodes, NODE_CARD_NODES_EXPECTED) - assert np.allclose(node_card.nnum, range(1, 6)) - - ELEMENT_SECTION = """*ELEMENT_SOLID 1 1 1 2 6 5 17 18 22 21 2 1 2 3 7 6 18 19 23 22 @@ -48,6 +32,14 @@ def test_node_card(tmp_path: Path) -> None: *END """ +ELEMENT_SECTION_ELEMS = [ + [1, 2, 6, 5, 17, 18, 22, 21], + [2, 3, 7, 6, 18, 19, 23, 22], + [3, 4, 8, 7, 19, 20, 24, 23], + [5, 6, 10, 9, 21, 22, 26, 25], + [6, 7, 11, 10, 22, 23, 27, 26], +] + def test_node_section(tmp_path: Path) -> None: filename = str(tmp_path / "tmp.k") @@ -60,6 +52,9 @@ def test_node_section(tmp_path: Path) -> None: assert len(deck.node_sections) == 1 node_section = deck.node_sections[0] + + assert "NodeSection containing 5 nodes" in str(node_section) + assert np.allclose(node_section.nodes, NODE_SECTION_NODES_EXPECTED) assert np.allclose(node_section.nnum, range(1, 6)) assert np.allclose(node_section.tc, [0, 0, 1, 2, 0]) @@ -69,13 +64,20 @@ def test_node_section(tmp_path: Path) -> None: def test_element_section(tmp_path: Path) -> None: filename = str(tmp_path / "tmp.k") with open(filename, "w") as fid: - fid.write(NODE_SECTION) + fid.write(ELEMENT_SECTION) deck = Deck(filename) deck.read_line() deck.read_element_section() - assert len(deck.element_sections) == 1 - element_section = deck.element_sections[0] + assert len(deck.element_solid_sections) == 1 + element_section = deck.element_solid_sections[0] + + assert "ElementSection containing 5 elements" in str(element_section) + + assert np.allclose(element_section.eid, range(1, 6)) + assert np.allclose(element_section.pid, [1] * 5) + assert np.allclose(element_section.node_ids, np.array(ELEMENT_SECTION_ELEMS).ravel()) - # assert np.allclose(element_section.element_id, + offsets = np.cumsum([0] + [len(element) for element in ELEMENT_SECTION_ELEMS]) + assert np.allclose(element_section.node_id_offsets, offsets)