Skip to content

Commit

Permalink
add repr for node and element sections
Browse files Browse the repository at this point in the history
  • Loading branch information
akaszynski committed Sep 6, 2024
1 parent b016df0 commit b25f24a
Show file tree
Hide file tree
Showing 4 changed files with 179 additions and 31 deletions.
28 changes: 23 additions & 5 deletions README.md
Original file line number Diff line number Diff line change
@@ -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:
Expand All @@ -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

```
98 changes: 93 additions & 5 deletions src/deck.cpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <math.h>
#include <sstream>
Expand Down Expand Up @@ -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 {
Expand All @@ -347,11 +383,14 @@ struct ElementSection {
NDArray<int, 1> node_ids;
NDArray<int, 1> node_id_offsets;

ElementSection() {}

ElementSection(std::vector<int> eid_vec, std::vector<int> pid_vec,
std::vector<int> node_ids_vec,
std::vector<int> node_id_offsets_vec) {
std::array<int, 1> nel_shape = {static_cast<int>(eid_vec.size())};
std::array<int, 1> node_ids_shape = {static_cast<int>(node_ids_vec.size())};
std::array<int, 1> node_ids_shape = {
static_cast<int>(node_id_offsets_vec.back())};
std::array<int, 1> node_ids_offsets_shape = {
static_cast<int>(eid_vec.size() + 1)};

Expand All @@ -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 {
Expand All @@ -371,10 +452,6 @@ class Deck {
MemoryMappedFile memmap;

public:
// int n_nodes = 0;
// int n_elem = 0;
// NDArray<double, 2> nodes_arr;

std::vector<NodeSection> node_sections;
std::vector<ElementSection> element_solid_sections;

Expand Down Expand Up @@ -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 =
Expand All @@ -526,13 +612,15 @@ class Deck {
NB_MODULE(_deck, m) {
nb::class_<NodeSection>(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)
.def_ro("rc", &NodeSection::rc, nb::rv_policy::automatic);

nb::class_<ElementSection>(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)
Expand Down
42 changes: 41 additions & 1 deletion src/lsdyna_mesh_reader/_deck.pyi
Original file line number Diff line number Diff line change
@@ -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: ...
42 changes: 22 additions & 20 deletions tests/test_deck.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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")
Expand All @@ -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])
Expand All @@ -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)

0 comments on commit b25f24a

Please sign in to comment.