Skip to content

Commit

Permalink
implement wedges and tets
Browse files Browse the repository at this point in the history
  • Loading branch information
@akaszynski
akaszynski committed Sep 7, 2024
1 parent e133ee9 commit 80a5a23
Show file tree
Hide file tree
Showing 2 changed files with 55 additions and 33 deletions.
57 changes: 38 additions & 19 deletions src/deck.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -505,31 +505,55 @@ struct ElementSolidSection : public ElementSection {
// convert cells, offset, and celltypes to vtk style arrays
nb::tuple ToVTK() {
NDArray<uint8_t, 1> celltypes_arr = MakeNDArray<uint8_t, 1>({(int)n_elem});
// NDArray<int64_t, 1> offsets_arr = MakeNDArray<int64_t, 1>({(int)(n_elem +
// 1)}); NDArray<int64_t, 1> cells_arr = MakeNDArray<int64_t,
// 1>({(int)node_ids.size()});
NDArray<int64_t, 1> offsets_arr =
MakeNDArray<int64_t, 1>({(int)(n_elem + 1)});

uint8_t *celltypes = celltypes_arr.data();
// int64_t *offsets = offsets_arr.data();
// int64_t *cells = cells_arr.data();
int64_t *offsets = offsets_arr.data();
int64_t *cells = AllocateArray<int64_t>(node_ids.size());

const int *node_id_offsets_data = node_id_offsets.data();
const int *node_ids_data = node_ids.data();

int c = 0;
// offsets[0] = 0;
offsets[0] = 0;
uint8_t celltype = VTK_EMPTY_CELL;
for (int i; i < n_elem; i++) {
// int offset = node_id_offsets_data[i];
// int el_sz = 8; // assuming hex
celltypes[i] = VTK_HEXAHEDRON;
// offsets[i + 1] = offsets[i] + el_sz; // start of next cell
// for (int j=0; j< el_sz; j++){
// cells[c++] = node_ids_data[offset + j];
// }

int el_sz = 0;
int offset = node_id_offsets_data[i];
if (node_ids_data[offset + 3] == node_ids_data[offset + 4]) {
celltypes[i] = VTK_TETRA;
cells[c++] = node_ids_data[offset + 0];
cells[c++] = node_ids_data[offset + 1];
cells[c++] = node_ids_data[offset + 2];
cells[c++] = node_ids_data[offset + 3];
el_sz = 4;
} else if (node_ids_data[offset + 5] == node_ids_data[offset + 6]) {
celltypes[i] = VTK_WEDGE;
// map to vtk style
cells[c++] = node_ids_data[offset + 0];
cells[c++] = node_ids_data[offset + 1];
cells[c++] = node_ids_data[offset + 4];
cells[c++] = node_ids_data[offset + 3];
cells[c++] = node_ids_data[offset + 2];
cells[c++] = node_ids_data[offset + 5];
el_sz = 6;
} else {
celltypes[i] = VTK_HEXAHEDRON;
// assuming compiler is smart enough to unroll this
for (int j = 0; j < 8; j++) {
cells[c++] = node_ids_data[offset + j];
}
el_sz = 8;
} // cell type branching

offsets[i + 1] = offsets[i] + el_sz; // start of next cell

} // for each cell

return nb::make_tuple(node_ids, node_id_offsets, celltypes_arr);
NDArray<int64_t, 1> cells_arr = WrapNDarray<int64_t, 1>(cells, {c});
return nb::make_tuple(cells_arr, offsets_arr, celltypes_arr);
} // to vtk
};

Expand Down Expand Up @@ -563,18 +587,13 @@ struct ElementShellSection : public ElementSection {
for (int i; i < n_elem; i++) {
// determine if the cell is a quad or triangle
int offset = node_id_offsets_data[i];
// int el_sz = node_id_offsets_data[i+1] - offset;
// if (el_sz == 4){
if (node_ids_data[offset + 2] == node_ids_data[offset + 3]) {
celltypes[i] = VTK_TRIANGLE;
el_sz = 3;
} else {
celltypes[i] = VTK_QUAD;
el_sz = 4;
}
// } else {
// throw std::runtime_error("Unsupported shell cell size");
// }

offsets[i + 1] = offsets[i] + el_sz; // start of next cell
for (int j = 0; j < el_sz; j++) {
Expand Down
31 changes: 17 additions & 14 deletions tests/test_deck.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,19 +28,15 @@

ELEMENT_SOLID_SECTION = """*ELEMENT_SOLID
1 1 1 2 6 5 17 18 22 21
2 1 2 3 7 6 18 19 23 22
3 1 3 4 8 7 19 20 24 23
4 1 5 6 10 9 21 22 26 25
5 1 6 7 11 10 22 23 27 26
2 1 5 6 10 9 9 9 9 9
3 1 6 7 11 10 22 23 23 23
*END
"""

ELEMENT_SOLID_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],
[5, 6, 10, 9, 9, 9, 9, 9],
[6, 7, 11, 10, 22, 23, 23, 23],
]


Expand Down Expand Up @@ -93,16 +89,23 @@ def test_element_solid_section(tmp_path: Path) -> None:

assert len(deck.element_solid_sections) == 1
section = deck.element_solid_sections[0]
assert "ElementSolidSection containing 5 elements" in str(section)
assert "ElementSolidSection containing 3 elements" in str(section)

assert np.allclose(section.eid, range(1, 6))
assert np.allclose(section.pid, [1] * 5)
assert np.allclose(section.eid, range(1, 4))
assert np.allclose(section.pid, [1] * 3)
assert np.allclose(section.node_ids, np.array(ELEMENT_SOLID_SECTION_ELEMS).ravel())

cells, offset, celltypes = section.to_vtk()
assert np.allclose(cells, section.node_ids) # expect no change
assert np.allclose(offset, section.node_id_offsets) # expect no change
assert np.allclose(celltypes, [pv.CellType.HEXAHEDRON] * 5)
expected_cells = np.hstack(
[
ELEMENT_SOLID_SECTION_ELEMS[0], # hex
ELEMENT_SOLID_SECTION_ELEMS[1][:4], # tet
[6, 7, 22, 10, 11, 23], # remapped to vtk style wedge
]
)
assert np.allclose(cells, expected_cells)
assert np.allclose(offset, [0, 8, 12, 18])
assert np.allclose(celltypes, [pv.CellType.HEXAHEDRON, pv.CellType.TETRA, pv.CellType.WEDGE])

offsets = np.cumsum([0] + [len(element) for element in ELEMENT_SOLID_SECTION_ELEMS])
assert np.allclose(section.node_id_offsets, offsets)
Expand Down

0 comments on commit 80a5a23

Please sign in to comment.