From 80a5a2376efefff175ce529dc70e86e0824ded68 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Sat, 7 Sep 2024 11:48:24 -0600 Subject: [PATCH] implement wedges and tets --- src/deck.cpp | 57 ++++++++++++++++++++++++++++++---------------- tests/test_deck.py | 31 +++++++++++++------------ 2 files changed, 55 insertions(+), 33 deletions(-) diff --git a/src/deck.cpp b/src/deck.cpp index 83faa65..590c034 100644 --- a/src/deck.cpp +++ b/src/deck.cpp @@ -505,31 +505,55 @@ struct ElementSolidSection : public ElementSection { // convert cells, offset, and celltypes to vtk style arrays nb::tuple ToVTK() { NDArray celltypes_arr = MakeNDArray({(int)n_elem}); - // NDArray offsets_arr = MakeNDArray({(int)(n_elem + - // 1)}); NDArray cells_arr = MakeNDArray({(int)node_ids.size()}); + NDArray offsets_arr = + MakeNDArray({(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(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 cells_arr = WrapNDarray(cells, {c}); + return nb::make_tuple(cells_arr, offsets_arr, celltypes_arr); } // to vtk }; @@ -563,8 +587,6 @@ 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; @@ -572,9 +594,6 @@ struct ElementShellSection : public ElementSection { 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++) { diff --git a/tests/test_deck.py b/tests/test_deck.py index 81a6e15..4887c73 100644 --- a/tests/test_deck.py +++ b/tests/test_deck.py @@ -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], ] @@ -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)