From 50e7e94ce45632b4f64d834b9631727d3e6b1400 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Thu, 29 Aug 2024 18:05:16 -0600 Subject: [PATCH] add pre-commit --- .github/workflows/testing-and-deployment.yml | 2 +- .pre-commit-config.yaml | 70 + src/archive.cpp | 927 +++++----- src/array_support.h | 82 +- src/lsdyna_mesh_reader/archive.py | 5 +- src/lsdyna_mesh_reader/mesh.py | 5 +- src/reader.cpp | 1684 +++++++++--------- src/vtk_support.cpp | 681 ++++--- src/vtk_support.h | 12 +- 9 files changed, 1733 insertions(+), 1735 deletions(-) create mode 100644 .pre-commit-config.yaml diff --git a/.github/workflows/testing-and-deployment.yml b/.github/workflows/testing-and-deployment.yml index cec0ad1..fe7cfd5 100644 --- a/.github/workflows/testing-and-deployment.yml +++ b/.github/workflows/testing-and-deployment.yml @@ -68,7 +68,7 @@ jobs: url: https://pypi.org/p/lsdyna-mesh-reader permissions: id-token: write # this permission is mandatory for trusted publishing - contents: write # required to create a release + contents: write # required to create a release steps: - uses: actions/download-artifact@v4 - name: Flatten directory structure diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..a8cfc37 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,70 @@ +# See https://pre-commit.ci/ +ci: + autofix_prs: true + autoupdate_schedule: quarterly + +repos: +- repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.6.3 + hooks: + - id: ruff + args: [--fix, --exit-non-zero-on-fix] + exclude: ^(docs/|tests) + - id: ruff-format + +- repo: https://github.com/keewis/blackdoc + rev: v0.3.9 + hooks: + - id: blackdoc + exclude: README.rst + +- repo: https://github.com/codespell-project/codespell + rev: v2.3.0 + hooks: + - id: codespell + args: [--skip=*.vt*] + +- repo: https://github.com/pycqa/pydocstyle + rev: 6.3.0 + hooks: + - id: pydocstyle + additional_dependencies: [tomli==2.0.1] + files: ^src/pyminiply/.*\.py + +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.6.0 + hooks: + - id: check-merge-conflict + - id: debug-statements + - id: trailing-whitespace + exclude: .*\.(cdb|k|dat)$ + +- repo: https://github.com/pre-commit/mirrors-mypy + rev: v1.11.2 + hooks: + - id: mypy + additional_dependencies: [numpy>=2.0, npt-promote==0.1] + exclude: ^tests/ + +- repo: https://github.com/macisamuele/language-formatters-pre-commit-hooks + rev: v2.14.0 + hooks: + - id: pretty-format-toml + args: [--autofix] + - id: pretty-format-yaml + args: [--autofix, --indent, '2'] + +- repo: https://github.com/pre-commit/mirrors-clang-format + rev: v17.0.6 + hooks: + - id: clang-format + +- repo: https://github.com/python-jsonschema/check-jsonschema + rev: 0.29.2 + hooks: + - id: check-github-workflows + +- repo: https://github.com/dzhu/rstfmt + rev: v0.0.14 + hooks: + - id: rstfmt diff --git a/src/archive.cpp b/src/archive.cpp index 9e94221..772e9fa 100644 --- a/src/archive.cpp +++ b/src/archive.cpp @@ -30,137 +30,121 @@ #define VTK_QUADRATIC_PYRAMID 27 char *create_format_str(int fields, int sig_digits) { - static char format_str[64]; // Buffer for our format string - char field_format[16]; // Buffer for a single field format string + static char format_str[64]; // Buffer for our format string + char field_format[16]; // Buffer for a single field format string - // Leaves one whitespace, sign, one before the decimal, decimal, and four - // for the scientific notation - // For example: " -1.233333333333E+09" - int total_char = sig_digits + 7; + // Leaves one whitespace, sign, one before the decimal, decimal, and four + // for the scientific notation + // For example: " -1.233333333333E+09" + int total_char = sig_digits + 7; - // Create a single field format string - sprintf(field_format, "%%%d.%dE", total_char, sig_digits); + // Create a single field format string + sprintf(field_format, "%%%d.%dE", total_char, sig_digits); - // Start the format string with the fixed part - strcpy(format_str, "%8d 0 0"); + // Start the format string with the fixed part + strcpy(format_str, "%8d 0 0"); - // Add the field format string for each field - for (int i = 0; i < fields; i++) { - strcat(format_str, field_format); - } + // Add the field format string for each field + for (int i = 0; i < fields; i++) { + strcat(format_str, field_format); + } - // Add a newline to the end of the format string - strcat(format_str, "\n"); + // Add a newline to the end of the format string + strcat(format_str, "\n"); - return format_str; + return format_str; } // Write node IDs, node coordinates, and angles to file as a NBLOCK template -void WriteNblock( - std::string &filename, - const int max_node_id, - const NDArray node_id_arr, - const NDArray nodes_arr, - const NDArray angles_arr, - int sig_digits, - std::string &mode) { - - FILE *file = fopen(filename.c_str(), mode.c_str()); - if (file == nullptr) { - throw std::runtime_error("Failed to open file: " + filename); +void WriteNblock(std::string &filename, const int max_node_id, + const NDArray node_id_arr, + const NDArray nodes_arr, + const NDArray angles_arr, int sig_digits, + std::string &mode) { + + FILE *file = fopen(filename.c_str(), mode.c_str()); + if (file == nullptr) { + throw std::runtime_error("Failed to open file: " + filename); + } + + const int n_nodes = node_id_arr.size(); + bool has_angles = angles_arr.size() > 0; + + // Header + // Tell ANSYS to start reading the node block with 6 fields, + // associated with a solid, the maximum node number and the number + // of lines in the node block + fprintf(file, "/PREP7\n"); + fprintf(file, "NBLOCK,6,SOLID,%10d,%10d\n", max_node_id, n_nodes); + fprintf(file, "(3i8,6e%d.%d)\n", sig_digits + 7, sig_digits); + + char *format_str; + const int *node_id = node_id_arr.data(); + const double *nodes = nodes_arr.data(); + const double *angles = angles_arr.data(); + + int i; + if (has_angles) { + format_str = create_format_str(6, sig_digits); + for (i = 0; i < n_nodes; i++) { + fprintf(file, format_str, node_id[i], nodes[i * 3 + 0], nodes[i * 3 + 1], + nodes[i * 3 + 2], angles[i * 3 + 0], angles[i * 3 + 1], + angles[i * 3 + 2]); } - - const int n_nodes = node_id_arr.size(); - bool has_angles = angles_arr.size() > 0; - - // Header - // Tell ANSYS to start reading the node block with 6 fields, - // associated with a solid, the maximum node number and the number - // of lines in the node block - fprintf(file, "/PREP7\n"); - fprintf(file, "NBLOCK,6,SOLID,%10d,%10d\n", max_node_id, n_nodes); - fprintf(file, "(3i8,6e%d.%d)\n", sig_digits + 7, sig_digits); - - char *format_str; - const int *node_id = node_id_arr.data(); - const double *nodes = nodes_arr.data(); - const double *angles = angles_arr.data(); - - int i; - if (has_angles) { - format_str = create_format_str(6, sig_digits); - for (i = 0; i < n_nodes; i++) { - fprintf( - file, - format_str, - node_id[i], - nodes[i * 3 + 0], - nodes[i * 3 + 1], - nodes[i * 3 + 2], - angles[i * 3 + 0], - angles[i * 3 + 1], - angles[i * 3 + 2]); - } - } else { - format_str = create_format_str(3, sig_digits); - for (i = 0; i < n_nodes; i++) { - fprintf( - file, - format_str, - node_id[i], - nodes[i * 3 + 0], - nodes[i * 3 + 1], - nodes[i * 3 + 2]); - } + } else { + format_str = create_format_str(3, sig_digits); + for (i = 0; i < n_nodes; i++) { + fprintf(file, format_str, node_id[i], nodes[i * 3 + 0], nodes[i * 3 + 1], + nodes[i * 3 + 2]); } + } - fprintf(file, "N,R5.3,LOC, -1,\n"); + fprintf(file, "N,R5.3,LOC, -1,\n"); - fclose(file); + fclose(file); } // Write an ANSYS EBLOCK to file void WriteEblock( const std::string &filename, - const int n_elem, // number of elements - const NDArray elem_id_arr, // element ID array - const NDArray etype_arr, // element type ID array - const NDArray mtype_arr, // material ID array - const NDArray rcon_arr, // real constant ID array - const NDArray elem_nnodes_arr, // number of nodes per element + const int n_elem, // number of elements + const NDArray elem_id_arr, // element ID array + const NDArray etype_arr, // element type ID array + const NDArray mtype_arr, // material ID array + const NDArray rcon_arr, // real constant ID array + const NDArray elem_nnodes_arr, // number of nodes per element const NDArray celltypes_arr, // VTK celltypes array const NDArray offset_arr, // VTK offset array - const NDArray cells_arr, // VTK cell connectivity array - const NDArray typenum_arr, // ANSYS type number (e.g. 187 for SOLID187) + const NDArray cells_arr, // VTK cell connectivity array + const NDArray + typenum_arr, // ANSYS type number (e.g. 187 for SOLID187) const NDArray nodenum_arr, // ANSYS node numbering std::string &mode) { - FILE *file = fopen(filename.c_str(), mode.c_str()); - - const int *elem_id = elem_id_arr.data(); - const int *etype = etype_arr.data(); - const int *mtype = mtype_arr.data(); - const int *rcon = rcon_arr.data(); - const int *elem_nnodes = elem_nnodes_arr.data(); - const uint8_t *celltypes = celltypes_arr.data(); - const int64_t *offset = offset_arr.data(); - const int64_t *cells = cells_arr.data(); - const int *typenum = typenum_arr.data(); - const int *nodenum = nodenum_arr.data(); - - // Write header - fprintf(file, "EBLOCK,19,SOLID,%10d,%10d\n", elem_id[n_elem - 1], n_elem); - fprintf(file, "(19i8)\n"); - - int c; // position within offset array - for (int i = 0; i < n_elem; i++) { - c = offset[i]; - - // Write cell info - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d", + FILE *file = fopen(filename.c_str(), mode.c_str()); + + const int *elem_id = elem_id_arr.data(); + const int *etype = etype_arr.data(); + const int *mtype = mtype_arr.data(); + const int *rcon = rcon_arr.data(); + const int *elem_nnodes = elem_nnodes_arr.data(); + const uint8_t *celltypes = celltypes_arr.data(); + const int64_t *offset = offset_arr.data(); + const int64_t *cells = cells_arr.data(); + const int *typenum = typenum_arr.data(); + const int *nodenum = nodenum_arr.data(); + + // Write header + fprintf(file, "EBLOCK,19,SOLID,%10d,%10d\n", elem_id[n_elem - 1], n_elem); + fprintf(file, "(19i8)\n"); + + int c; // position within offset array + for (int i = 0; i < n_elem; i++) { + c = offset[i]; + + // Write cell info + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d", mtype[i], // Field 1: material reference number etype[i], // Field 2: element type number rcon[i], // Field 3: real constant reference number @@ -173,277 +157,245 @@ void WriteEblock( 0, // Field 10: Not Used elem_id[i]); // Field 11: Element number - // Write element nodes - switch (celltypes[i]) { - case VTK_QUADRATIC_TETRA: - if (typenum[i] == 187) { - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 3]], // 3, L - nodenum[cells[c + 4]], // 4, M - nodenum[cells[c + 5]], // 5, N - nodenum[cells[c + 6]], // 6, O - nodenum[cells[c + 7]], // 7, P - nodenum[cells[c + 8]], // 8, Q - nodenum[cells[c + 9]]); // 9, R - } else { - // Using SOLID186-like format - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 2]], // 3, L (duplicate of K) - nodenum[cells[c + 3]], // 4, M - nodenum[cells[c + 3]], // 5, N (duplicate of M) - nodenum[cells[c + 3]], // 6, O (duplicate of M) - nodenum[cells[c + 3]], // 7, P (duplicate of M) - nodenum[cells[c + 4]], // 8, Q - nodenum[cells[c + 5]], // 9, R - nodenum[cells[c + 3]], // 10, S (duplicate of K) - nodenum[cells[c + 6]], // 11, T - nodenum[cells[c + 3]], // 12, U (duplicate of M) - nodenum[cells[c + 3]], // 13, V (duplicate of M) - nodenum[cells[c + 3]], // 14, W (duplicate of M) - nodenum[cells[c + 3]], // 15, X (duplicate of M) - nodenum[cells[c + 7]], // 16, Y - nodenum[cells[c + 8]], // 17, Z - nodenum[cells[c + 9]], // 18, A - nodenum[cells[c + 9]]); // 19, B (duplicate of A) - } - break; - case VTK_TETRA: // point - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 2]], // 3, L (duplicate of K) - nodenum[cells[c + 3]], // 4, M - nodenum[cells[c + 3]], // 5, N (duplicate of M) - nodenum[cells[c + 3]], // 6, O (duplicate of M) - nodenum[cells[c + 3]]); // 7, P (duplicate of M) - break; - case VTK_WEDGE: - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 2]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 0]], // 2, K - nodenum[cells[c + 0]], // 3, L (duplicate of K) - nodenum[cells[c + 5]], // 4, M - nodenum[cells[c + 4]], // 5, N - nodenum[cells[c + 3]], // 6, O - nodenum[cells[c + 3]]); // 7, P (duplicate of O) - break; - case VTK_QUADRATIC_WEDGE: - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 2]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 0]], // 2, K - nodenum[cells[c + 0]], // 3, L (duplicate of K) - nodenum[cells[c + 5]], // 4, M - nodenum[cells[c + 4]], // 5, N - nodenum[cells[c + 3]], // 6, O - nodenum[cells[c + 3]], // 7, P (duplicate of O) - nodenum[cells[c + 7]], // 8, Q - nodenum[cells[c + 6]], // 9, R - nodenum[cells[c + 0]], // 10, S (duplicate of K) - nodenum[cells[c + 8]], // 11, T - nodenum[cells[c + 10]], // 12, U - nodenum[cells[c + 9]], // 13, V - nodenum[cells[c + 3]], // 14, W (duplicate of O) - nodenum[cells[c + 11]], // 15, X - nodenum[cells[c + 14]], // 16, Y - nodenum[cells[c + 13]], // 17, Z - nodenum[cells[c + 12]], // 18, A - nodenum[cells[c + 12]]); // 19, B (duplicate of A) - break; - case VTK_QUADRATIC_PYRAMID: - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 3]], // 3, L - nodenum[cells[c + 4]], // 4, M - nodenum[cells[c + 4]], // 5, N (duplicate of M) - nodenum[cells[c + 4]], // 6, O (duplicate of M) - nodenum[cells[c + 4]], // 7, P (duplicate of M) - nodenum[cells[c + 5]], // 8, Q - nodenum[cells[c + 6]], // 9, R - nodenum[cells[c + 7]], // 10, S - nodenum[cells[c + 8]], // 11, T - nodenum[cells[c + 4]], // 12, U (duplicate of M) - nodenum[cells[c + 4]], // 13, V (duplicate of M) - nodenum[cells[c + 4]], // 14, W (duplicate of M) - nodenum[cells[c + 4]], // 15, X (duplicate of M) - nodenum[cells[c + 9]], // 16, Y - nodenum[cells[c + 10]], // 17, Z - nodenum[cells[c + 11]], // 18, A - nodenum[cells[c + 12]]); // 19, B (duplicate of A) - break; - case VTK_PYRAMID: - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 3]], // 3, L - nodenum[cells[c + 4]], // 4, M - nodenum[cells[c + 4]], // 5, N (duplicate of M) - nodenum[cells[c + 4]], // 6, O (duplicate of M) - nodenum[cells[c + 4]]); // 7, P (duplicate of M) - break; - case VTK_VOXEL: - // note the flipped order for nodes (K, L) and (O, P) - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n", + // Write element nodes + switch (celltypes[i]) { + case VTK_QUADRATIC_TETRA: + if (typenum[i] == 187) { + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d\n", nodenum[cells[c + 0]], // 0, I nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 3]], // 2, K - nodenum[cells[c + 2]], // 3, L + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 3]], // 3, L nodenum[cells[c + 4]], // 4, M nodenum[cells[c + 5]], // 5, N - nodenum[cells[c + 7]], // 6, O - nodenum[cells[c + 6]]); // 7, P - break; - case VTK_HEXAHEDRON: - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], - nodenum[cells[c + 1]], - nodenum[cells[c + 2]], - nodenum[cells[c + 3]], - nodenum[cells[c + 4]], - nodenum[cells[c + 5]], - nodenum[cells[c + 6]], - nodenum[cells[c + 7]]); - break; - case VTK_QUADRATIC_HEXAHEDRON: - fprintf( - file, - "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], - nodenum[cells[c + 1]], - nodenum[cells[c + 2]], - nodenum[cells[c + 3]], - nodenum[cells[c + 4]], - nodenum[cells[c + 5]], - nodenum[cells[c + 6]], - nodenum[cells[c + 7]], - nodenum[cells[c + 8]], - nodenum[cells[c + 9]], - nodenum[cells[c + 10]], - nodenum[cells[c + 11]], - nodenum[cells[c + 12]], - nodenum[cells[c + 13]], - nodenum[cells[c + 14]], - nodenum[cells[c + 15]], - nodenum[cells[c + 16]], - nodenum[cells[c + 17]], - nodenum[cells[c + 18]], - nodenum[cells[c + 19]]); - break; - case VTK_TRIANGLE: - fprintf( - file, - "%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 2]]); // 3, L (duplicate of K) - break; - case VTK_QUAD: - fprintf( - file, - "%8d%8d%8d%8d\n", - nodenum[cells[c + 0]], // 0, I - nodenum[cells[c + 1]], // 1, J - nodenum[cells[c + 2]], // 2, K - nodenum[cells[c + 3]]); // 3, L - break; - } + nodenum[cells[c + 6]], // 6, O + nodenum[cells[c + 7]], // 7, P + nodenum[cells[c + 8]], // 8, Q + nodenum[cells[c + 9]]); // 9, R + } else { + // Using SOLID186-like format + fprintf( + file, + "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 2]], // 3, L (duplicate of K) + nodenum[cells[c + 3]], // 4, M + nodenum[cells[c + 3]], // 5, N (duplicate of M) + nodenum[cells[c + 3]], // 6, O (duplicate of M) + nodenum[cells[c + 3]], // 7, P (duplicate of M) + nodenum[cells[c + 4]], // 8, Q + nodenum[cells[c + 5]], // 9, R + nodenum[cells[c + 3]], // 10, S (duplicate of K) + nodenum[cells[c + 6]], // 11, T + nodenum[cells[c + 3]], // 12, U (duplicate of M) + nodenum[cells[c + 3]], // 13, V (duplicate of M) + nodenum[cells[c + 3]], // 14, W (duplicate of M) + nodenum[cells[c + 3]], // 15, X (duplicate of M) + nodenum[cells[c + 7]], // 16, Y + nodenum[cells[c + 8]], // 17, Z + nodenum[cells[c + 9]], // 18, A + nodenum[cells[c + 9]]); // 19, B (duplicate of A) + } + break; + case VTK_TETRA: // point + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 2]], // 3, L (duplicate of K) + nodenum[cells[c + 3]], // 4, M + nodenum[cells[c + 3]], // 5, N (duplicate of M) + nodenum[cells[c + 3]], // 6, O (duplicate of M) + nodenum[cells[c + 3]]); // 7, P (duplicate of M) + break; + case VTK_WEDGE: + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 2]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 0]], // 2, K + nodenum[cells[c + 0]], // 3, L (duplicate of K) + nodenum[cells[c + 5]], // 4, M + nodenum[cells[c + 4]], // 5, N + nodenum[cells[c + 3]], // 6, O + nodenum[cells[c + 3]]); // 7, P (duplicate of O) + break; + case VTK_QUADRATIC_WEDGE: + fprintf( + file, + "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 2]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 0]], // 2, K + nodenum[cells[c + 0]], // 3, L (duplicate of K) + nodenum[cells[c + 5]], // 4, M + nodenum[cells[c + 4]], // 5, N + nodenum[cells[c + 3]], // 6, O + nodenum[cells[c + 3]], // 7, P (duplicate of O) + nodenum[cells[c + 7]], // 8, Q + nodenum[cells[c + 6]], // 9, R + nodenum[cells[c + 0]], // 10, S (duplicate of K) + nodenum[cells[c + 8]], // 11, T + nodenum[cells[c + 10]], // 12, U + nodenum[cells[c + 9]], // 13, V + nodenum[cells[c + 3]], // 14, W (duplicate of O) + nodenum[cells[c + 11]], // 15, X + nodenum[cells[c + 14]], // 16, Y + nodenum[cells[c + 13]], // 17, Z + nodenum[cells[c + 12]], // 18, A + nodenum[cells[c + 12]]); // 19, B (duplicate of A) + break; + case VTK_QUADRATIC_PYRAMID: + fprintf( + file, + "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 3]], // 3, L + nodenum[cells[c + 4]], // 4, M + nodenum[cells[c + 4]], // 5, N (duplicate of M) + nodenum[cells[c + 4]], // 6, O (duplicate of M) + nodenum[cells[c + 4]], // 7, P (duplicate of M) + nodenum[cells[c + 5]], // 8, Q + nodenum[cells[c + 6]], // 9, R + nodenum[cells[c + 7]], // 10, S + nodenum[cells[c + 8]], // 11, T + nodenum[cells[c + 4]], // 12, U (duplicate of M) + nodenum[cells[c + 4]], // 13, V (duplicate of M) + nodenum[cells[c + 4]], // 14, W (duplicate of M) + nodenum[cells[c + 4]], // 15, X (duplicate of M) + nodenum[cells[c + 9]], // 16, Y + nodenum[cells[c + 10]], // 17, Z + nodenum[cells[c + 11]], // 18, A + nodenum[cells[c + 12]]); // 19, B (duplicate of A) + break; + case VTK_PYRAMID: + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 3]], // 3, L + nodenum[cells[c + 4]], // 4, M + nodenum[cells[c + 4]], // 5, N (duplicate of M) + nodenum[cells[c + 4]], // 6, O (duplicate of M) + nodenum[cells[c + 4]]); // 7, P (duplicate of M) + break; + case VTK_VOXEL: + // note the flipped order for nodes (K, L) and (O, P) + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 3]], // 2, K + nodenum[cells[c + 2]], // 3, L + nodenum[cells[c + 4]], // 4, M + nodenum[cells[c + 5]], // 5, N + nodenum[cells[c + 7]], // 6, O + nodenum[cells[c + 6]]); // 7, P + break; + case VTK_HEXAHEDRON: + fprintf(file, "%8d%8d%8d%8d%8d%8d%8d%8d\n", nodenum[cells[c + 0]], + nodenum[cells[c + 1]], nodenum[cells[c + 2]], + nodenum[cells[c + 3]], nodenum[cells[c + 4]], + nodenum[cells[c + 5]], nodenum[cells[c + 6]], + nodenum[cells[c + 7]]); + break; + case VTK_QUADRATIC_HEXAHEDRON: + fprintf( + file, + "%8d%8d%8d%8d%8d%8d%8d%8d\n%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], nodenum[cells[c + 1]], nodenum[cells[c + 2]], + nodenum[cells[c + 3]], nodenum[cells[c + 4]], nodenum[cells[c + 5]], + nodenum[cells[c + 6]], nodenum[cells[c + 7]], nodenum[cells[c + 8]], + nodenum[cells[c + 9]], nodenum[cells[c + 10]], nodenum[cells[c + 11]], + nodenum[cells[c + 12]], nodenum[cells[c + 13]], + nodenum[cells[c + 14]], nodenum[cells[c + 15]], + nodenum[cells[c + 16]], nodenum[cells[c + 17]], + nodenum[cells[c + 18]], nodenum[cells[c + 19]]); + break; + case VTK_TRIANGLE: + fprintf(file, "%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 2]]); // 3, L (duplicate of K) + break; + case VTK_QUAD: + fprintf(file, "%8d%8d%8d%8d\n", + nodenum[cells[c + 0]], // 0, I + nodenum[cells[c + 1]], // 1, J + nodenum[cells[c + 2]], // 2, K + nodenum[cells[c + 3]]); // 3, L + break; } - fprintf(file, " -1\n"); + } + fprintf(file, " -1\n"); - // TODO: consider using C++ ifstream - fclose(file); + // TODO: consider using C++ ifstream + fclose(file); } NDArray CmblockItems(const NDArray array) { - const int *data = array.data(); - - // Use vector for initial operations - std::vector temp(data, data + array.size()); - std::sort(temp.begin(), temp.end()); - temp.erase(std::unique(temp.begin(), temp.end()), temp.end()); - - std::vector items; - items.push_back(temp[0]); - - int last_index = 0; // Tracks the index of the last unique range start - for (size_t i = 0; i < temp.size() - 1; ++i) { - // Check if the next item forms a continuous range - if (temp[i + 1] - temp[i] != 1) { - if (temp[i] != items[last_index]) { - items.push_back(-temp[i]); // End the current range - } - items.push_back(temp[i + 1]); // Start a new range - last_index = items.size() - 1; - } + const int *data = array.data(); + + // Use vector for initial operations + std::vector temp(data, data + array.size()); + std::sort(temp.begin(), temp.end()); + temp.erase(std::unique(temp.begin(), temp.end()), temp.end()); + + std::vector items; + items.push_back(temp[0]); + + int last_index = 0; // Tracks the index of the last unique range start + for (size_t i = 0; i < temp.size() - 1; ++i) { + // Check if the next item forms a continuous range + if (temp[i + 1] - temp[i] != 1) { + if (temp[i] != items[last_index]) { + items.push_back(-temp[i]); // End the current range + } + items.push_back(temp[i + 1]); // Start a new range + last_index = items.size() - 1; } + } - // End the last range if it was continuous - if (temp.back() != items.back()) { - items.push_back(-temp.back()); - } + // End the last range if it was continuous + if (temp.back() != items.back()) { + items.push_back(-temp.back()); + } - // Copy to ndarray - NDArray items_arr = MakeNDArray({(int)items.size()}); - int *items_data = items_arr.data(); - std::copy(items.begin(), items.end(), items_data); + // Copy to ndarray + NDArray items_arr = MakeNDArray({(int)items.size()}); + int *items_data = items_arr.data(); + std::copy(items.begin(), items.end(), items_data); - return items_arr; + return items_arr; } // Resets the midside nodes of the tetrahedral starting at index c. // midside nodes between // (0,1), (1,2), (2,0), (0,3), (1,3), and (2,3) template inline void ResetMidTet(const int64_t *cells, T *points) { - int64_t ind0 = cells[0] * 3; - int64_t ind1 = cells[1] * 3; - int64_t ind2 = cells[2] * 3; - int64_t ind3 = cells[3] * 3; - int64_t ind4 = cells[4] * 3; - int64_t ind5 = cells[5] * 3; - int64_t ind6 = cells[6] * 3; - int64_t ind7 = cells[7] * 3; - int64_t ind8 = cells[8] * 3; - int64_t ind9 = cells[9] * 3; - - for (size_t j = 0; j < 3; j++) { - points[ind4 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; - points[ind5 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; - points[ind6 + j] = (points[ind2 + j] + points[ind0 + j]) * 0.5; - points[ind7 + j] = (points[ind0 + j] + points[ind3 + j]) * 0.5; - points[ind8 + j] = (points[ind1 + j] + points[ind3 + j]) * 0.5; - points[ind9 + j] = (points[ind2 + j] + points[ind3 + j]) * 0.5; - } + int64_t ind0 = cells[0] * 3; + int64_t ind1 = cells[1] * 3; + int64_t ind2 = cells[2] * 3; + int64_t ind3 = cells[3] * 3; + int64_t ind4 = cells[4] * 3; + int64_t ind5 = cells[5] * 3; + int64_t ind6 = cells[6] * 3; + int64_t ind7 = cells[7] * 3; + int64_t ind8 = cells[8] * 3; + int64_t ind9 = cells[9] * 3; + + for (size_t j = 0; j < 3; j++) { + points[ind4 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; + points[ind5 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; + points[ind6 + j] = (points[ind2 + j] + points[ind0 + j]) * 0.5; + points[ind7 + j] = (points[ind0 + j] + points[ind3 + j]) * 0.5; + points[ind8 + j] = (points[ind1 + j] + points[ind3 + j]) * 0.5; + points[ind9 + j] = (points[ind2 + j] + points[ind3 + j]) * 0.5; + } } // Reset pyr midside nodes between: @@ -456,70 +408,70 @@ template inline void ResetMidTet(const int64_t *cells, T *points) { // 11(2, 4) // 12(3, 4) template inline void ResetMidPyr(const int64_t *cells, T *points) { - int64_t ind0 = cells[0] * 3; - int64_t ind1 = cells[1] * 3; - int64_t ind2 = cells[2] * 3; - int64_t ind3 = cells[3] * 3; - int64_t ind4 = cells[4] * 3; - int64_t ind5 = cells[5] * 3; - int64_t ind6 = cells[6] * 3; - int64_t ind7 = cells[7] * 3; - int64_t ind8 = cells[8] * 3; - int64_t ind9 = cells[9] * 3; - int64_t ind10 = cells[10] * 3; - int64_t ind11 = cells[11] * 3; - int64_t ind12 = cells[12] * 3; - - for (size_t j = 0; j < 3; j++) { - points[ind5 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; - points[ind6 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; - points[ind7 + j] = (points[ind2 + j] + points[ind3 + j]) * 0.5; - points[ind8 + j] = (points[ind3 + j] + points[ind0 + j]) * 0.5; - points[ind9 + j] = (points[ind0 + j] + points[ind4 + j]) * 0.5; - points[ind10 + j] = (points[ind1 + j] + points[ind4 + j]) * 0.5; - points[ind11 + j] = (points[ind2 + j] + points[ind4 + j]) * 0.5; - points[ind12 + j] = (points[ind3 + j] + points[ind4 + j]) * 0.5; - } + int64_t ind0 = cells[0] * 3; + int64_t ind1 = cells[1] * 3; + int64_t ind2 = cells[2] * 3; + int64_t ind3 = cells[3] * 3; + int64_t ind4 = cells[4] * 3; + int64_t ind5 = cells[5] * 3; + int64_t ind6 = cells[6] * 3; + int64_t ind7 = cells[7] * 3; + int64_t ind8 = cells[8] * 3; + int64_t ind9 = cells[9] * 3; + int64_t ind10 = cells[10] * 3; + int64_t ind11 = cells[11] * 3; + int64_t ind12 = cells[12] * 3; + + for (size_t j = 0; j < 3; j++) { + points[ind5 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; + points[ind6 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; + points[ind7 + j] = (points[ind2 + j] + points[ind3 + j]) * 0.5; + points[ind8 + j] = (points[ind3 + j] + points[ind0 + j]) * 0.5; + points[ind9 + j] = (points[ind0 + j] + points[ind4 + j]) * 0.5; + points[ind10 + j] = (points[ind1 + j] + points[ind4 + j]) * 0.5; + points[ind11 + j] = (points[ind2 + j] + points[ind4 + j]) * 0.5; + points[ind12 + j] = (points[ind3 + j] + points[ind4 + j]) * 0.5; + } } template inline void ResetMidWeg(const int64_t *cells, T *points) { - // Reset midside nodes of a wedge cell: - // 6 (0,1) - // 7 (1,2) - // 8 (2,0) - // 9 (3,4) - // 10 (4,5) - // 11 (5,3) - // 12 (0,3) - // 13 (1,4) - // 14 (2,5) - int64_t ind0 = cells[0] * 3; - int64_t ind1 = cells[1] * 3; - int64_t ind2 = cells[2] * 3; - int64_t ind3 = cells[3] * 3; - int64_t ind4 = cells[4] * 3; - int64_t ind5 = cells[5] * 3; - int64_t ind6 = cells[6] * 3; - int64_t ind7 = cells[7] * 3; - int64_t ind8 = cells[8] * 3; - int64_t ind9 = cells[9] * 3; - int64_t ind10 = cells[10] * 3; - int64_t ind11 = cells[11] * 3; - int64_t ind12 = cells[12] * 3; - int64_t ind13 = cells[13] * 3; - int64_t ind14 = cells[14] * 3; - - for (size_t j = 0; j < 3; j++) { - points[ind6 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; - points[ind7 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; - points[ind8 + j] = (points[ind2 + j] + points[ind0 + j]) * 0.5; - points[ind9 + j] = (points[ind3 + j] + points[ind4 + j]) * 0.5; - points[ind10 + j] = (points[ind4 + j] + points[ind5 + j]) * 0.5; - points[ind11 + j] = (points[ind5 + j] + points[ind3 + j]) * 0.5; - points[ind12 + j] = (points[ind0 + j] + points[ind3 + j]) * 0.5; - points[ind13 + j] = (points[ind1 + j] + points[ind4 + j]) * 0.5; - points[ind14 + j] = (points[ind2 + j] + points[ind5 + j]) * 0.5; - } + // Reset midside nodes of a wedge cell: + // 6 (0,1) + // 7 (1,2) + // 8 (2,0) + // 9 (3,4) + // 10 (4,5) + // 11 (5,3) + // 12 (0,3) + // 13 (1,4) + // 14 (2,5) + int64_t ind0 = cells[0] * 3; + int64_t ind1 = cells[1] * 3; + int64_t ind2 = cells[2] * 3; + int64_t ind3 = cells[3] * 3; + int64_t ind4 = cells[4] * 3; + int64_t ind5 = cells[5] * 3; + int64_t ind6 = cells[6] * 3; + int64_t ind7 = cells[7] * 3; + int64_t ind8 = cells[8] * 3; + int64_t ind9 = cells[9] * 3; + int64_t ind10 = cells[10] * 3; + int64_t ind11 = cells[11] * 3; + int64_t ind12 = cells[12] * 3; + int64_t ind13 = cells[13] * 3; + int64_t ind14 = cells[14] * 3; + + for (size_t j = 0; j < 3; j++) { + points[ind6 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; + points[ind7 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; + points[ind8 + j] = (points[ind2 + j] + points[ind0 + j]) * 0.5; + points[ind9 + j] = (points[ind3 + j] + points[ind4 + j]) * 0.5; + points[ind10 + j] = (points[ind4 + j] + points[ind5 + j]) * 0.5; + points[ind11 + j] = (points[ind5 + j] + points[ind3 + j]) * 0.5; + points[ind12 + j] = (points[ind0 + j] + points[ind3 + j]) * 0.5; + points[ind13 + j] = (points[ind1 + j] + points[ind4 + j]) * 0.5; + points[ind14 + j] = (points[ind2 + j] + points[ind5 + j]) * 0.5; + } } // Reset midside nodes of a hexahedral cell @@ -536,80 +488,79 @@ template inline void ResetMidWeg(const int64_t *cells, T *points) { // 18 (2,6) // 19 (3,7) template inline void ResetMidHex(const int64_t *cells, T *points) { - int64_t ind0 = cells[0] * 3; - int64_t ind1 = cells[1] * 3; - int64_t ind2 = cells[2] * 3; - int64_t ind3 = cells[3] * 3; - int64_t ind4 = cells[4] * 3; - int64_t ind5 = cells[5] * 3; - int64_t ind6 = cells[6] * 3; - int64_t ind7 = cells[7] * 3; - int64_t ind8 = cells[8] * 3; - int64_t ind9 = cells[9] * 3; - int64_t ind10 = cells[10] * 3; - int64_t ind11 = cells[11] * 3; - int64_t ind12 = cells[12] * 3; - int64_t ind13 = cells[13] * 3; - int64_t ind14 = cells[14] * 3; - int64_t ind15 = cells[15] * 3; - int64_t ind16 = cells[16] * 3; - int64_t ind17 = cells[17] * 3; - int64_t ind18 = cells[18] * 3; - int64_t ind19 = cells[19] * 3; - - for (size_t j = 0; j < 3; j++) { - points[ind8 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; - points[ind9 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; - points[ind10 + j] = (points[ind2 + j] + points[ind3 + j]) * 0.5; - points[ind11 + j] = (points[ind3 + j] + points[ind0 + j]) * 0.5; - points[ind12 + j] = (points[ind4 + j] + points[ind5 + j]) * 0.5; - points[ind13 + j] = (points[ind5 + j] + points[ind6 + j]) * 0.5; - points[ind14 + j] = (points[ind6 + j] + points[ind7 + j]) * 0.5; - points[ind15 + j] = (points[ind7 + j] + points[ind4 + j]) * 0.5; - points[ind16 + j] = (points[ind0 + j] + points[ind4 + j]) * 0.5; - points[ind17 + j] = (points[ind1 + j] + points[ind5 + j]) * 0.5; - points[ind18 + j] = (points[ind2 + j] + points[ind6 + j]) * 0.5; - points[ind19 + j] = (points[ind3 + j] + points[ind7 + j]) * 0.5; - } + int64_t ind0 = cells[0] * 3; + int64_t ind1 = cells[1] * 3; + int64_t ind2 = cells[2] * 3; + int64_t ind3 = cells[3] * 3; + int64_t ind4 = cells[4] * 3; + int64_t ind5 = cells[5] * 3; + int64_t ind6 = cells[6] * 3; + int64_t ind7 = cells[7] * 3; + int64_t ind8 = cells[8] * 3; + int64_t ind9 = cells[9] * 3; + int64_t ind10 = cells[10] * 3; + int64_t ind11 = cells[11] * 3; + int64_t ind12 = cells[12] * 3; + int64_t ind13 = cells[13] * 3; + int64_t ind14 = cells[14] * 3; + int64_t ind15 = cells[15] * 3; + int64_t ind16 = cells[16] * 3; + int64_t ind17 = cells[17] * 3; + int64_t ind18 = cells[18] * 3; + int64_t ind19 = cells[19] * 3; + + for (size_t j = 0; j < 3; j++) { + points[ind8 + j] = (points[ind0 + j] + points[ind1 + j]) * 0.5; + points[ind9 + j] = (points[ind1 + j] + points[ind2 + j]) * 0.5; + points[ind10 + j] = (points[ind2 + j] + points[ind3 + j]) * 0.5; + points[ind11 + j] = (points[ind3 + j] + points[ind0 + j]) * 0.5; + points[ind12 + j] = (points[ind4 + j] + points[ind5 + j]) * 0.5; + points[ind13 + j] = (points[ind5 + j] + points[ind6 + j]) * 0.5; + points[ind14 + j] = (points[ind6 + j] + points[ind7 + j]) * 0.5; + points[ind15 + j] = (points[ind7 + j] + points[ind4 + j]) * 0.5; + points[ind16 + j] = (points[ind0 + j] + points[ind4 + j]) * 0.5; + points[ind17 + j] = (points[ind1 + j] + points[ind5 + j]) * 0.5; + points[ind18 + j] = (points[ind2 + j] + points[ind6 + j]) * 0.5; + points[ind19 + j] = (points[ind3 + j] + points[ind7 + j]) * 0.5; + } } template -void ResetMidside( - NDArray celltypes_arr, - NDArray cells_arr, - NDArray offset_arr, - NDArray points_arr) { - - int n_cells = celltypes_arr.size(); - const int64_t *cells = cells_arr.data(); - const int64_t *offset = offset_arr.data(); - const uint8_t *celltypes = celltypes_arr.data(); - T *points = points_arr.data(); - - for (int i = 0; i < n_cells; i++) { - int64_t c = offset[i]; - - switch (celltypes[i]) { - case VTK_QUADRATIC_TETRA: - ResetMidTet(&cells[c], points); - continue; - case VTK_QUADRATIC_PYRAMID: - ResetMidPyr(&cells[c], points); - continue; - case VTK_QUADRATIC_WEDGE: - ResetMidWeg(&cells[c], points); - continue; - case VTK_QUADRATIC_HEXAHEDRON: - ResetMidHex(&cells[c], points); - } +void ResetMidside(NDArray celltypes_arr, + NDArray cells_arr, + NDArray offset_arr, + NDArray points_arr) { + + int n_cells = celltypes_arr.size(); + const int64_t *cells = cells_arr.data(); + const int64_t *offset = offset_arr.data(); + const uint8_t *celltypes = celltypes_arr.data(); + T *points = points_arr.data(); + + for (int i = 0; i < n_cells; i++) { + int64_t c = offset[i]; + + switch (celltypes[i]) { + case VTK_QUADRATIC_TETRA: + ResetMidTet(&cells[c], points); + continue; + case VTK_QUADRATIC_PYRAMID: + ResetMidPyr(&cells[c], points); + continue; + case VTK_QUADRATIC_WEDGE: + ResetMidWeg(&cells[c], points); + continue; + case VTK_QUADRATIC_HEXAHEDRON: + ResetMidHex(&cells[c], points); } + } } NB_MODULE(_archive, m) { - m.def("write_nblock", &WriteNblock); - m.def("write_nblock", &WriteNblock); - m.def("write_eblock", &WriteEblock); - m.def("cmblock_items_from_array", &CmblockItems); - m.def("reset_midside", &ResetMidside); - m.def("reset_midside", &ResetMidside); + m.def("write_nblock", &WriteNblock); + m.def("write_nblock", &WriteNblock); + m.def("write_eblock", &WriteEblock); + m.def("cmblock_items_from_array", &CmblockItems); + m.def("reset_midside", &ResetMidside); + m.def("reset_midside", &ResetMidside); } diff --git a/src/array_support.h b/src/array_support.h index cdfda54..a2fb312 100644 --- a/src/array_support.h +++ b/src/array_support.h @@ -15,61 +15,63 @@ namespace nb = nanobind; template using NDArray = nb::ndarray, nb::c_contig>; -template T *AllocateArray(size_t total, bool zero_initialize = false) { - T *data = zero_initialize ? new T[total]() : new T[total]; +template +T *AllocateArray(size_t total, bool zero_initialize = false) { + T *data = zero_initialize ? new T[total]() : new T[total]; #ifdef __linux__ - const size_t hugepage_threshold = 1u << 22u; // 4MB threshold - const size_t page_size = 4096u; + const size_t hugepage_threshold = 1u << 22u; // 4MB threshold + const size_t page_size = 4096u; - if (total * sizeof(T) >= hugepage_threshold) { - uintptr_t data_addr = reinterpret_cast(data); - size_t offset = page_size - (data_addr % page_size); - size_t length = total * sizeof(T) - offset; + if (total * sizeof(T) >= hugepage_threshold) { + uintptr_t data_addr = reinterpret_cast(data); + size_t offset = page_size - (data_addr % page_size); + size_t length = total * sizeof(T) - offset; - madvise(reinterpret_cast(data_addr + offset), length, MADV_HUGEPAGE); - } + madvise(reinterpret_cast(data_addr + offset), length, + MADV_HUGEPAGE); + } #endif - return data; + return data; } // wrap an existing array as a numpy ndarray template -NDArray -WrapNDarray(T *data, const std::array shape, bool zero_initialize = false) { - size_t shape_[N]; - for (size_t i = 0; i < N; ++i) { - shape_[i] = shape[i]; - } +NDArray WrapNDarray(T *data, const std::array shape, + bool zero_initialize = false) { + size_t shape_[N]; + for (size_t i = 0; i < N; ++i) { + shape_[i] = shape[i]; + } - nb::capsule owner(data, [](void *p) noexcept { delete[] (T *)p; }); + nb::capsule owner(data, [](void *p) noexcept { delete[] (T *)p; }); - return NDArray(data, N, shape_, owner); + return NDArray(data, N, shape_, owner); } template -NDArray MakeNDArray( - const std::array shape, bool zero_initialize = false, bool use_owner = true) { - - // Calculate the total number of elements in the ndarray - size_t total = 1; - for (size_t i = 0; i < N; i++) { - total *= shape[i]; - } - - size_t shape_[N]; - for (size_t i = 0; i < N; i++) { - shape_[i] = shape[i]; - } - - T *data = AllocateArray(total, zero_initialize); - if (use_owner) { - nb::capsule owner(data, [](void *p) noexcept { delete[] (T *)p; }); - return NDArray(data, N, shape_, owner); - } else { - return NDArray(data, N, shape_, nullptr); - } +NDArray MakeNDArray(const std::array shape, + bool zero_initialize = false, bool use_owner = true) { + + // Calculate the total number of elements in the ndarray + size_t total = 1; + for (size_t i = 0; i < N; i++) { + total *= shape[i]; + } + + size_t shape_[N]; + for (size_t i = 0; i < N; i++) { + shape_[i] = shape[i]; + } + + T *data = AllocateArray(total, zero_initialize); + if (use_owner) { + nb::capsule owner(data, [](void *p) noexcept { delete[] (T *)p; }); + return NDArray(data, N, shape_, owner); + } else { + return NDArray(data, N, shape_, nullptr); + } } #endif // ARRAY_SUPPORT_HEADER_H diff --git a/src/lsdyna_mesh_reader/archive.py b/src/lsdyna_mesh_reader/archive.py index aff32db..0bc56b4 100644 --- a/src/lsdyna_mesh_reader/archive.py +++ b/src/lsdyna_mesh_reader/archive.py @@ -8,11 +8,10 @@ from typing import Any, Dict, List, Optional, Sequence, TypeVar, Union import numpy as np -from numpy.typing import NDArray -from pyvista import ID_TYPE, CellType, UnstructuredGrid - from mapdl_archive import _archive, _reader from mapdl_archive.mesh import Mesh +from numpy.typing import NDArray +from pyvista import ID_TYPE, CellType, UnstructuredGrid # types NPArray_FLOAT32 = NDArray[np.float32] diff --git a/src/lsdyna_mesh_reader/mesh.py b/src/lsdyna_mesh_reader/mesh.py index 342f3d0..17c94f5 100644 --- a/src/lsdyna_mesh_reader/mesh.py +++ b/src/lsdyna_mesh_reader/mesh.py @@ -5,14 +5,13 @@ from typing import Dict, List, Optional, Tuple, TypeVar, Union, cast import numpy as np +from mapdl_archive import _archive, _reader +from mapdl_archive.elements import ETYPE_MAP from numpy.typing import NDArray from pyvista import ID_TYPE, CellArray from pyvista.core.pointset import PolyData, UnstructuredGrid from vtkmodules.util.numpy_support import numpy_to_vtk -from mapdl_archive import _archive, _reader -from mapdl_archive.elements import ETYPE_MAP - COMP_DICT = Dict[str, NDArray[np.int32]] INVALID_ALLOWABLE_TYPES = TypeError( diff --git a/src/reader.cpp b/src/reader.cpp index 35309b8..36f9bc2 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -34,49 +34,50 @@ using namespace nb::literals; // #define DEBUG static const double DIV_OF_TEN[] = { - 1.0e-0, 1.0e-1, 1.0e-2, 1.0e-3, 1.0e-4, 1.0e-5, 1.0e-6, 1.0e-7, 1.0e-8, 1.0e-9, - 1.0e-10, 1.0e-11, 1.0e-12, 1.0e-13, 1.0e-14, 1.0e-15, 1.0e-16, 1.0e-17, 1.0e-18, 1.0e-19, - 1.0e-20, 1.0e-21, 1.0e-22, 1.0e-23, 1.0e-24, 1.0e-25, 1.0e-26, 1.0e-27, + 1.0e-0, 1.0e-1, 1.0e-2, 1.0e-3, 1.0e-4, 1.0e-5, 1.0e-6, + 1.0e-7, 1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12, 1.0e-13, + 1.0e-14, 1.0e-15, 1.0e-16, 1.0e-17, 1.0e-18, 1.0e-19, 1.0e-20, + 1.0e-21, 1.0e-22, 1.0e-23, 1.0e-24, 1.0e-25, 1.0e-26, 1.0e-27, }; static inline double power_of_ten(int exponent) { - double result = 1.0; - double base = (exponent < 0) ? 0.1 : 10.0; - int abs_exponent = abs(exponent); - for (int i = 0; i < abs_exponent; ++i) { - result *= base; - } - return result; + double result = 1.0; + double base = (exponent < 0) ? 0.1 : 10.0; + int abs_exponent = abs(exponent); + for (int i = 0; i < abs_exponent; ++i) { + result *= base; + } + return result; } // Fast string to integer converter static inline int fast_atoi(const char *raw, const int intsz) { - int val = 0; - int c; - char current_char; + int val = 0; + int c; + char current_char; - for (c = 0; c < intsz; ++c) { - current_char = *raw++; + for (c = 0; c < intsz; ++c) { + current_char = *raw++; - // Multiply by 10 only if current_char is a digit - if (current_char >= '0' && current_char <= '9') { - val = val * 10 + (current_char - '0'); - } + // Multiply by 10 only if current_char is a digit + if (current_char >= '0' && current_char <= '9') { + val = val * 10 + (current_char - '0'); } + } - return val; + return val; } // Checks for negative character static inline int checkneg(char *raw, int intsz) { - int c; - for (c = 0; c < intsz; ++c) { - if (raw[0] == '-') { - return 1; - } - ++raw; + int c; + for (c = 0; c < intsz; ++c) { + if (raw[0] == '-') { + return 1; } - return 0; + ++raw; + } + return 0; } // Reads fortran formatted float formats in the form of @@ -85,846 +86,850 @@ static inline int checkneg(char *raw, int intsz) { // // fltsz : Number of characters to read in a floating point number static inline int strtod(char *raw, int fltsz, double *arr) { - char *end = raw + fltsz; - double sign = 1; - - while (raw < end) { - if (*raw == '\0') { - // section empty, value is zero - *arr = 0; - return 1; - } else if (*raw != ' ') { // always skip whitespace - break; - } - raw++; + char *end = raw + fltsz; + double sign = 1; + + while (raw < end) { + if (*raw == '\0') { + // section empty, value is zero + *arr = 0; + return 1; + } else if (*raw != ' ') { // always skip whitespace + break; } - - // either a number or a sign + raw++; + } + + // either a number or a sign + if (*raw == '-') { + sign = -1; + ++raw; + } + + // next value is always a number + // Use integer arithmetric and then convert to a float + uint64_t val_int = *raw++ - '0'; + raw++; // next value is always a "." + + // Read through the rest of the number + int decimal_digits = 0; + while (raw < end) { + if (*raw == 'e' || *raw == 'E') { // incredibly, can be lowercase + break; + } else if (*raw >= '0' && *raw <= '9') { + val_int = val_int * 10 + (*raw++ - '0'); + decimal_digits++; + } + } + + // Compute the floating-point value + double val; + if (decimal_digits < 27) { + val = (double)val_int * DIV_OF_TEN[decimal_digits]; + } else { + val = (double)val_int * power_of_ten(10); + } + + // Might have scientific notation remaining, for example: + // 1.0000000000000E-001 + int evalue = 0; + int esign = 1; + if (*raw == 'e' || *raw == 'E') { + raw++; // skip "E" + // always a sign of some sort if (*raw == '-') { - sign = -1; - ++raw; + esign = -1; } + raw++; - // next value is always a number - // Use integer arithmetric and then convert to a float - uint64_t val_int = *raw++ - '0'; - raw++; // next value is always a "." - - // Read through the rest of the number - int decimal_digits = 0; while (raw < end) { - if (*raw == 'e' || *raw == 'E') { // incredibly, can be lowercase - break; - } else if (*raw >= '0' && *raw <= '9') { - val_int = val_int * 10 + (*raw++ - '0'); - decimal_digits++; - } + // read to whitespace or end of the line + if (*raw == ' ' || *raw == '\0') { + break; + } + evalue = evalue * 10 + (*raw++ - '0'); } - - // Compute the floating-point value - double val; - if (decimal_digits < 27) { - val = (double)val_int * DIV_OF_TEN[decimal_digits]; + if (esign == 1) { + val *= power_of_ten(evalue); } else { - val = (double)val_int * power_of_ten(10); - } - - // Might have scientific notation remaining, for example: - // 1.0000000000000E-001 - int evalue = 0; - int esign = 1; - if (*raw == 'e' || *raw == 'E') { - raw++; // skip "E" - // always a sign of some sort - if (*raw == '-') { - esign = -1; - } - raw++; - - while (raw < end) { - // read to whitespace or end of the line - if (*raw == ' ' || *raw == '\0') { - break; - } - evalue = evalue * 10 + (*raw++ - '0'); - } - if (esign == 1) { - val *= power_of_ten(evalue); - } else { - val /= power_of_ten(evalue); - } + val /= power_of_ten(evalue); } + } - // seek through end of float value - if (sign == -1) { - *arr = -val; - } else { - *arr = val; - } + // seek through end of float value + if (sign == -1) { + *arr = -val; + } else { + *arr = val; + } - // check if at end of line - return 0; + // check if at end of line + return 0; } class MemoryMappedFile { - private: - size_t size; - char *start; +private: + size_t size; + char *start; #ifdef _WIN32 - HANDLE fileHandle; - HANDLE mapHandle; + HANDLE fileHandle; + HANDLE mapHandle; #else - int fd; + int fd; #endif - public: - std::string line; - char *current; +public: + std::string line; + char *current; - MemoryMappedFile(const char *filename) - : start(nullptr), current(nullptr), size(0) + MemoryMappedFile(const char *filename) + : start(nullptr), current(nullptr), size(0) #ifdef _WIN32 - , - fileHandle(INVALID_HANDLE_VALUE), mapHandle(nullptr) + , + fileHandle(INVALID_HANDLE_VALUE), mapHandle(nullptr) #else - , - fd(-1) + , + fd(-1) #endif - { + { #ifdef _WIN32 - fileHandle = CreateFile( - filename, - GENERIC_READ, - FILE_SHARE_READ, - nullptr, - OPEN_EXISTING, - FILE_ATTRIBUTE_NORMAL, - nullptr); - if (fileHandle == INVALID_HANDLE_VALUE) { - throw std::runtime_error("Error opening file"); - } + fileHandle = CreateFile(filename, GENERIC_READ, FILE_SHARE_READ, nullptr, + OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, nullptr); + if (fileHandle == INVALID_HANDLE_VALUE) { + throw std::runtime_error("Error opening file"); + } - LARGE_INTEGER fileSize; - if (!GetFileSizeEx(fileHandle, &fileSize)) { - CloseHandle(fileHandle); - throw std::runtime_error("Error getting file size"); - } + LARGE_INTEGER fileSize; + if (!GetFileSizeEx(fileHandle, &fileSize)) { + CloseHandle(fileHandle); + throw std::runtime_error("Error getting file size"); + } - size = static_cast(fileSize.QuadPart); - mapHandle = CreateFileMapping(fileHandle, nullptr, PAGE_READONLY, 0, 0, nullptr); - if (mapHandle == nullptr) { - CloseHandle(fileHandle); - throw std::runtime_error("Error creating file mapping"); - } + size = static_cast(fileSize.QuadPart); + mapHandle = + CreateFileMapping(fileHandle, nullptr, PAGE_READONLY, 0, 0, nullptr); + if (mapHandle == nullptr) { + CloseHandle(fileHandle); + throw std::runtime_error("Error creating file mapping"); + } - start = static_cast(MapViewOfFile(mapHandle, FILE_MAP_READ, 0, 0, size)); - if (start == nullptr) { - CloseHandle(mapHandle); - CloseHandle(fileHandle); - throw std::runtime_error("Error mapping file"); - } + start = static_cast( + MapViewOfFile(mapHandle, FILE_MAP_READ, 0, 0, size)); + if (start == nullptr) { + CloseHandle(mapHandle); + CloseHandle(fileHandle); + throw std::runtime_error("Error mapping file"); + } #else - fd = open(filename, O_RDONLY); - if (fd == -1) { - throw std::runtime_error("Error opening file"); - } + fd = open(filename, O_RDONLY); + if (fd == -1) { + throw std::runtime_error("Error opening file"); + } - struct stat st; - if (fstat(fd, &st) == -1) { - close(fd); - throw std::runtime_error("Error getting file size"); - } + struct stat st; + if (fstat(fd, &st) == -1) { + close(fd); + throw std::runtime_error("Error getting file size"); + } - size = st.st_size; - start = static_cast(mmap(nullptr, size, PROT_READ, MAP_PRIVATE, fd, 0)); - if (start == MAP_FAILED) { - close(fd); - throw std::runtime_error("Error mapping file"); - } -#endif - current = start; + size = st.st_size; + start = + static_cast(mmap(nullptr, size, PROT_READ, MAP_PRIVATE, fd, 0)); + if (start == MAP_FAILED) { + close(fd); + throw std::runtime_error("Error mapping file"); } +#endif + current = start; + } - ~MemoryMappedFile() { - close_file(); + ~MemoryMappedFile() { + close_file(); #ifdef _WIN32 - if (fileHandle != INVALID_HANDLE_VALUE) { - CloseHandle(fileHandle); - } - if (mapHandle != nullptr) { - CloseHandle(mapHandle); - } + if (fileHandle != INVALID_HANDLE_VALUE) { + CloseHandle(fileHandle); + } + if (mapHandle != nullptr) { + CloseHandle(mapHandle); + } #else - if (fd != -1) { - close(fd); - } -#endif + if (fd != -1) { + close(fd); } +#endif + } - void close_file() { - if (start) { + void close_file() { + if (start) { #ifdef _WIN32 - UnmapViewOfFile(start); + UnmapViewOfFile(start); #else - munmap(start, size); + munmap(start, size); #endif - start = nullptr; - current = nullptr; - } - } - - char &operator[](size_t index) { - // implement bounds checking? - // if (index >= size) { - // throw std::out_of_range("Index out of bounds"); - // } - return current[index]; + start = nullptr; + current = nullptr; } + } - void operator+=(size_t offset) { current += offset; } + char &operator[](size_t index) { + // implement bounds checking? + // if (index >= size) { + // throw std::out_of_range("Index out of bounds"); + // } + return current[index]; + } - // Seek to the end of the line - void seek_eol() { - // check if at end of file - if (current >= start + size) { - // std::cout << "end" << std::endl; - return; - } + void operator+=(size_t offset) { current += offset; } - while (current < start + size && *current != '\n') { - current++; - } + // Seek to the end of the line + void seek_eol() { + // check if at end of file + if (current >= start + size) { + // std::cout << "end" << std::endl; + return; + } - if (current < start + size && *current == '\n') { - current++; - } + while (current < start + size && *current != '\n') { + current++; } - bool eof() { return current >= start + size; } + if (current < start + size && *current == '\n') { + current++; + } + } - bool read_line() { - line.clear(); - if (current >= start + size) { - return false; - } + bool eof() { return current >= start + size; } - char *line_start = current; - while (current < start + size && *current != '\n') { - line += *current++; - } + bool read_line() { + line.clear(); + if (current >= start + size) { + return false; + } - if (current < start + size && *current == '\n') { - current++; - } + char *line_start = current; + while (current < start + size && *current != '\n') { + line += *current++; + } - return line_start != current; + if (current < start + size && *current == '\n') { + current++; } - size_t current_line_length() const { - char *temp = current; - size_t length = 0; - while (temp < start + size && *temp != '\n') { - length++; - temp++; - } - return length; + return line_start != current; + } + + size_t current_line_length() const { + char *temp = current; + size_t length = 0; + while (temp < start + size && *temp != '\n') { + length++; + temp++; } + return length; + } - off_t tellg() const { return current - start; } + off_t tellg() const { return current - start; } }; +int ReadNBlockMemMap(MemoryMappedFile &memmap, int *nnum, double *nodes, + double *node_angles, const int nnodes, + std::array d_size, const int f_size) { -int ReadNBlockMemMap( - MemoryMappedFile &memmap, - int *nnum, - double *nodes, - double *node_angles, - const int nnodes, - std::array d_size, - const int f_size) { - - int i, j, i_val, eol; - for (i = 0; i < nnodes; i++) { - // Read a line from the file - int count = memmap.current_line_length(); - - // It's possible that less nodes are written to the record than - // indicated. In this case the line starts with a -1 - if (memmap[0] == '-') { - break; - } + int i, j, i_val, eol; + for (i = 0; i < nnodes; i++) { + // Read a line from the file + int count = memmap.current_line_length(); - i_val = fast_atoi(memmap.current, d_size[0]); + // It's possible that less nodes are written to the record than + // indicated. In this case the line starts with a -1 + if (memmap[0] == '-') { + break; + } + + i_val = fast_atoi(memmap.current, d_size[0]); #ifdef DEBUG - std::cout << "Node number " << i_val << std::endl; + std::cout << "Node number " << i_val << std::endl; #endif - nnum[i] = i_val; + nnum[i] = i_val; - memmap += d_size[0]; - memmap += d_size[1]; - memmap += d_size[2]; + memmap += d_size[0]; + memmap += d_size[1]; + memmap += d_size[2]; - // read nodes - int n_read = (count - d_size[0] * 3) / f_size; + // read nodes + int n_read = (count - d_size[0] * 3) / f_size; #ifdef DEBUG - std::cout << "n_read: " << n_read << std::endl; + std::cout << "n_read: " << n_read << std::endl; #endif - int n_read_nodes = (n_read < 3) ? n_read : 3; - for (j = 0; j < n_read_nodes; j++) { - ans_strtod(memmap.current, f_size, &nodes[3 * i + j]); + int n_read_nodes = (n_read < 3) ? n_read : 3; + for (j = 0; j < n_read_nodes; j++) { + ans_strtod(memmap.current, f_size, &nodes[3 * i + j]); #ifdef DEBUG - std::cout << " " << nodes[3 * i + j] << " "; + std::cout << " " << nodes[3 * i + j] << " "; #endif - memmap.current += f_size; - } + memmap.current += f_size; + } #ifdef DEBUG - std::cout << std::endl; + std::cout << std::endl; #endif - // fill in unread nodes with zeros - for (; j < 3; j++) { - nodes[3 * i + j] = 0; - } - - // read in node angles if applicable - for (; j < n_read; j++) { - eol = ans_strtod(memmap.current, f_size, &node_angles[3 * i + (j - 3)]); - if (eol) - break; - memmap += f_size; - } - - // seek to the end of the line - memmap.seek_eol(); + // fill in unread nodes with zeros + for (; j < 3; j++) { + nodes[3 * i + j] = 0; } - return i; -} - -class Deck { - - private: - // std::string line; - bool debug; - std::string filename; - MemoryMappedFile memmap; - - public: - bool read_parameters; - - bool read_eblock; - bool eblock_is_read = false; - bool nblock_is_read = false; - - std::vector> elem_type; - std::unordered_map>> keyopt; - - // Element block - int n_elem = 0; - NDArray elem_arr; - NDArray elem_off_arr; - - // Node block - int n_nodes = 0; - int nblock_start = 0; - int nblock_end = 0; - NDArray nodes_arr; - NDArray node_angles_arr; - NDArray nnum_arr; - - // real constants - std::vector rnum; // real constant number - std::vector> rdat; // real constant data - - // components - std::unordered_map> elem_comps; - std::unordered_map> node_comps; - - Deck( - const std::string &fname, - bool readParams = false, - bool dbg = false, - bool readEblock = true) - : filename(fname), read_elements(readEblock), debug(dbg) - memmap(fname.c_str()) { - - // Likely bogus leak warnings. See: - // https://nanobind.readthedocs.io/en/latest/faq.html#why-am-i-getting-errors-about-leaked-functions-and-types - nb::set_leak_warnings(false); - + // read in node angles if applicable + for (; j < n_read; j++) { + eol = ans_strtod(memmap.current, f_size, &node_angles[3 * i + (j - 3)]); + if (eol) + break; + memmap += f_size; } - // // EXAMPLE: ET, 4, 186 - // // Element type, element number, element type - // void ReadETLine() { - // if (debug) { - // std::cout << "reading ET" << std::endl; - // } - - // // Assumes that the memory map is already at the ET line - // std::istringstream iss(memmap.line); - // std::string token; - // std::vector et_vals; - - // // Must start with ET (simply skip first token) - // std::getline(iss, token, ','); - // // if (token.compare(0, 2, "ET") != 0 && token.compare(0, 2, "et") != 0) { - // // return; - // // } - - // while (std::getline(iss, token, ',')) { - // try { - // et_vals.push_back(std::stoi(token)); - // } catch (const std::invalid_argument &) { - // return; - // } - // } - // if (et_vals.size() == 2) { - // elem_type.push_back(et_vals); - // } - // } - - // // Read ETBLOCK - // void ReadETBlock() { - // std::string token; - // std::vector values; - - // // Assumes current line is an ETBLOCK - // std::istringstream iss(memmap.line); - - // // skip first item (ETBLOCK) - // std::getline(iss, token, ','); - - // // read the number of items in the block - // std::getline(iss, token, ','); - // int n_items = 0; - // try { - // n_items = std::stoi(token); - // } catch (...) { - // std::cerr << "Failed to read ETBLOCK n_items" << std::endl; - // return; - // } - - // // Skip format line (2i9,19a9) - // memmap.seek_eol(); - - // // Read each item - // int elem_id, elem_type_id; - // for (int i = 0; i < n_items; i++) { - // // std::getline(cfile, line); - // memmap.read_line(); - // iss.clear(); - // iss.str(memmap.line); - - // // We only care about the first two integers in the block, the rest are keyopts - // if (!(iss >> elem_id >> elem_type_id)) { - // // allow a soft return - // std::cerr << "Failed to read ETBLOCK entry at line: " << memmap.line - // << std::endl; - // return; - // } - // std::vector entry; - // entry.push_back(elem_id); - // entry.push_back(elem_type_id); - // elem_type.push_back(entry); - // } - // } - - // // Read EBLOCK - // void ReadEBlock() { - // // Sometimes, DAT files contain two EBLOCKs. Read only the first block. - // if (eblock_is_read) { - // return; - // } - - // // Assumes already start of EBLOCK - // std::istringstream iss(memmap.line); - - // // Only read in SOLID eblocks - // std::transform( - // memmap.line.begin(), memmap.line.end(), memmap.line.begin(), [](unsigned char c) { - // return std::tolower(c); - // }); - // if (memmap.line.find("solid") == std::string::npos) { - // return; - // } - - // // Get size of EBLOCK from the last item in the line - // // Example: "EBLOCK,19,SOLID,,3588" - // n_elem = getSizeOfEBLOCK(memmap.line); - - // // we have to allocate memory for the maximum size since we don't know that a priori - // int *elem_data = AllocateArray(n_elem * 30); - // elem_off_arr = MakeNDArray({n_elem + 1}); - // int elem_sz = ReadEBlockMemMap(memmap, elem_off_arr.data(), elem_data, n_elem); - - // // wrap the raw data but limit it to the size of the number of elements read - // // Note: This is faster than reallocating a new array but will use more memory - // elem_arr = WrapNDarray(elem_data, {elem_sz}); - - // // must mark eblock is read since we can only read one eblock per archive file - // eblock_is_read = true; - // } - - // // Read: - // // KEYOPT, ITYPE, KNUM, VALUE - // void ReadKEYOPTLine() { - // if (debug) { - // std::cout << "reading KEYOPT" << std::endl; - // } - - // // Assumes at KEYOPT line - // std::istringstream iss(memmap.line); - // std::string token; - - // // Skip the first token (KEYOPT) - // std::getline(iss, token, ','); - - // // element number is the first item - // int etype; - // std::getline(iss, token, ','); - // try { - // etype = std::stoi(token); - // } catch (const std::invalid_argument &) { - // // soft error - // std::cerr << "Invalid format in KEYOPT line for elem_num:" << memmap.line - // << std::endl; - // return; - // } - - // // Read the rest of the values (Element Type, Key Option, Value) - // std::vector keyopt_vals; - // while (std::getline(iss, token, ',')) { - // try { - // keyopt_vals.push_back(std::stoi(token)); - // } catch (const std::invalid_argument &) { - // // soft error - // std::cerr << "Invalid format in KEYOPT line:" << memmap.line << std::endl; - // return; - // } - // } - - // if (keyopt_vals.size() == 2) { - // keyopt[etype].push_back(keyopt_vals); - // } - // } - - // // Read RLBLOCK - // void ReadRLBLOCK() { - // if (debug) { - // std::cout << "reading RLBLOCK" << std::endl; - // } - - // std::vector set_dat; - - // // Assumes line at RLBLOCK - // std::istringstream iss(memmap.line); - // std::string token; - // // while (std::getline(iss, token, ',')) { - // // set_dat.push_back(std::stoi(token)); - // // } - - // // skip starting RLBLOCK command - // std::getline(iss, token, ','); - - // int nset, maxset, maxitems, nperline; - // try { - // std::getline(iss, token, ','); - // nset = std::stoi(token); - // } catch (...) { - // std::cerr << "Error RLBLOCK line when reading nset:" << memmap.line << std::endl; - // return; - // } - // try { - // std::getline(iss, token, ','); - // maxset = std::stoi(token); - // } catch (...) { - // std::cerr << "Error RLBLOCK line when reading maxset:" << memmap.line - // << std::endl; - // return; - // } - // try { - // std::getline(iss, token, ','); - // maxitems = std::stoi(token); - // } catch (...) { - // std::cerr << "Error RLBLOCK line when reading maxitems:" << memmap.line - // << std::endl; - // return; - // } - // try { - // std::getline(iss, token, ','); - // nperline = std::stoi(token); - // } catch (...) { - // std::cerr << "Error RLBLOCK line when reading nperline:" << memmap.line - // << std::endl; - // return; - // } - - // // Skip next two lines - // memmap.seek_eol(); // (2i8,6g16.9) - // memmap.seek_eol(); // (7g16.9) - - // for (int i = 0; i < nset; i++) { - // // read real constant data in the form of: - // // 2 6 1.00000000 0.566900000E-07 0.00000000 0.000... - // memmap.read_line(); - // std::vector rcon; - - // // real constant number - // try { - // rnum.push_back(std::stoi(memmap.line.substr(0, 8))); - // } catch (...) { - // std::cerr << "Failed to read real constant ID when reading:" << memmap.line - // << std::endl; - // return; - // } - - // int ncon = std::stoi(memmap.line.substr(8, 8)); - // if (ncon > 6) { - // for (int j = 0; j < 6; j++) { - // rcon.push_back(std::stod(memmap.line.substr(16 + 16 * j, 16))); - // ncon--; - // } - - // memmap.read_line(); - - // while (ncon > 0) { - // for (int j = 0; j < 7 && ncon > 0; j++) { - // rcon.push_back(std::stod(memmap.line.substr(16 * j, 16))); - // ncon--; - // } - // if (ncon > 0) { - // memmap.read_line(); - // } - // } - // } else { - // for (int j = 0; j < ncon; j++) { - // rcon.push_back(std::stod(memmap.line.substr(16 + 16 * j, 16))); - // } - // } - - // rdat.push_back(rcon); - // } - // } + // seek to the end of the line + memmap.seek_eol(); + } - // void ReadNBlock(const int pos) { - - // // Sometimes, DAT files contains multiple node blocks. Read only the first block. - // if (nblock_is_read) { - // return; - // } - // nblock_start = pos; - - // // Get size of NBLOCK - // // Assumes line is at NBLOCK - // // std::cout << "line: " << memmap.line << std::endl; - // try { - // // Number of nodes is last item in string - // n_nodes = std::stoi(memmap.line.substr(memmap.line.rfind(',') + 1)); - // } catch (...) { - // std::cerr << "Failed to read number of nodes when reading:" << memmap.line - // << std::endl; - // return; - // } - // if (debug) { - // std::cout << "Reading " << n_nodes << " nodes" << std::endl; - // } - - // int *nnum_data = AllocateArray(n_nodes); - // double *nodes_data = AllocateArray(n_nodes * 3); - - // // often angles aren't written, so it makes sense to initialize this to - // // zero - // double *node_angles_data = AllocateArray(n_nodes * 3, true); - - // // Get format of nblock - // memmap.read_line(); - // NodeBlockFormat nfmt = GetNodeBlockFormat(memmap.line); - - // // Return actual number of nodes read and wrap the raw data - // // std::cout << memmap.tellg() << std::endl; - // n_nodes = ReadNBlockMemMap( - // memmap, - // nnum_data, - // nodes_data, - // node_angles_data, - // n_nodes, - // nfmt.d_size, - // nfmt.f_size); - // nnum_arr = WrapNDarray(nnum_data, {n_nodes}); - // nodes_arr = WrapNDarray(nodes_data, {n_nodes, 3}); - // node_angles_arr = WrapNDarray(node_angles_data, {n_nodes, 3}); - // // std::cout << memmap.tellg() << std::endl; - - // // Read final line, this is always "N,R5.3,LOC, -1," and store file - // // position. This is used for later access (or rewrite) of the node - // // block. - // memmap.seek_eol(); - // nblock_end = memmap.tellg(); - - // // Must mark nblock is read since we can only read one nblock per archive file. - // nblock_is_read = true; - // if (debug) { - // // std::cout << "Last line is: " << memmap.line << std::endl; - // std::cout << "NBLOCK complete, read " << n_nodes << " nodes." << std::endl; - // } - // } - - // // Read CMBLOCK - // void ReadCMBlock() { - // if (debug) { - // std::cout << "reading CMBLOCK" << std::endl; - // } - - // // Assumes line at CMBLOCK - // std::istringstream iss(memmap.line); - // std::string token; - // std::vector split_line; - - // while (std::getline(iss, token, ',')) { - // split_line.push_back(token); - // } - - // if (split_line.size() < 4) { - // std::cerr << "Poor formatting of CMBLOCK: " << memmap.line << std::endl; - // return; - // } - - // int ncomp; - // try { - // ncomp = std::stoi(split_line[3].substr(0, split_line[3].find('!')).c_str()); - // } catch (...) { - // std::cerr << "Poor formatting of CMBLOCK: " << memmap.line << std::endl; - // return; - // } - - // memmap.read_line(); - // int isz; - // try { - // isz = std::stoi(memmap.line.substr( - // memmap.line.find('i') + 1, - // memmap.line.find(')') - memmap.line.find('i') - 1)); - // } catch (...) { - // std::cerr << "Failed to read integer size in CMBLOCK" << std::endl; - // return; - // } - - // int nblock; - // try { - // nblock = std::stoi(memmap.line.substr( - // memmap.line.find('(') + 1, - // memmap.line.find('i') - memmap.line.find('(') - 1)); - // } catch (...) { - // std::cerr << "Failed to read number of integers per line in CMBLOCK" << std::endl; - // return; - // } - - // std::vector component(ncomp); - // std::string tempstr(isz, '\0'); - - // for (int i = 0; i < ncomp; ++i) { - // if (i % nblock == 0) { - // memmap.read_line(); - // } - // try { - // component[i] = std::stoi(memmap.line.substr(isz * (i % nblock), isz)); - // } catch (...) { - // std::cerr << "Failed to parse integer from CMBLOCK line: " << memmap.line - // << std::endl; - // return; - // } - // } - - // // get component name - // std::string comname = split_line[1]; - // comname.erase(comname.find_last_not_of(' ') + 1, std::string::npos); - - // std::string line_comp_type = split_line[2]; - - // if (line_comp_type.find("NODE") != std::string::npos) { - // node_comps[comname] = InterpretComponent(component); - // } else if (line_comp_type.find("ELEM") != std::string::npos) { - // elem_comps[comname] = InterpretComponent(component); - // } - - // if (debug) { - // std::cout << "Done reading CMBLOCK" << std::endl; - // } - // } - - // Read the LS-DYNA deck by reading the start of each line and checking if - // it matches a keyword. For example, *NODE denotes a node card. - // Note: Cards appear to start with the * char. - void Read() { - int first_char, next_char; + return i; +} - int position_start, position_end; - while (true) { - // Parse based on the first character rather than reading the entire - // line. It's faster and the parsing logic is always based on the first - // character +class Deck { - if (memmap.eof()) { +private: + // std::string line; + bool debug; + std::string filename; + MemoryMappedFile memmap; + +public: + bool read_parameters; + + bool read_eblock; + bool eblock_is_read = false; + bool nblock_is_read = false; + + std::vector> elem_type; + std::unordered_map>> keyopt; + + // Element block + int n_elem = 0; + NDArray elem_arr; + NDArray elem_off_arr; + + // Node block + int n_nodes = 0; + int nblock_start = 0; + int nblock_end = 0; + NDArray nodes_arr; + NDArray node_angles_arr; + NDArray nnum_arr; + + // real constants + std::vector rnum; // real constant number + std::vector> rdat; // real constant data + + // components + std::unordered_map> elem_comps; + std::unordered_map> node_comps; + + Deck(const std::string &fname, bool readParams = false, bool dbg = false, + bool readEblock = true) + : filename(fname), read_elements(readEblock), + debug(dbg) memmap(fname.c_str()) { + + // Likely bogus leak warnings. See: + // https://nanobind.readthedocs.io/en/latest/faq.html#why-am-i-getting-errors-about-leaked-functions-and-types + nb::set_leak_warnings(false); + } + + // // EXAMPLE: ET, 4, 186 + // // Element type, element number, element type + // void ReadETLine() { + // if (debug) { + // std::cout << "reading ET" << std::endl; + // } + + // // Assumes that the memory map is already at the ET line + // std::istringstream iss(memmap.line); + // std::string token; + // std::vector et_vals; + + // // Must start with ET (simply skip first token) + // std::getline(iss, token, ','); + // // if (token.compare(0, 2, "ET") != 0 && token.compare(0, 2, "et") != + // 0) { + // // return; + // // } + + // while (std::getline(iss, token, ',')) { + // try { + // et_vals.push_back(std::stoi(token)); + // } catch (const std::invalid_argument &) { + // return; + // } + // } + // if (et_vals.size() == 2) { + // elem_type.push_back(et_vals); + // } + // } + + // // Read ETBLOCK + // void ReadETBlock() { + // std::string token; + // std::vector values; + + // // Assumes current line is an ETBLOCK + // std::istringstream iss(memmap.line); + + // // skip first item (ETBLOCK) + // std::getline(iss, token, ','); + + // // read the number of items in the block + // std::getline(iss, token, ','); + // int n_items = 0; + // try { + // n_items = std::stoi(token); + // } catch (...) { + // std::cerr << "Failed to read ETBLOCK n_items" << std::endl; + // return; + // } + + // // Skip format line (2i9,19a9) + // memmap.seek_eol(); + + // // Read each item + // int elem_id, elem_type_id; + // for (int i = 0; i < n_items; i++) { + // // std::getline(cfile, line); + // memmap.read_line(); + // iss.clear(); + // iss.str(memmap.line); + + // // We only care about the first two integers in the block, the rest + // are keyopts if (!(iss >> elem_id >> elem_type_id)) { + // // allow a soft return + // std::cerr << "Failed to read ETBLOCK entry at line: " << + // memmap.line + // << std::endl; + // return; + // } + // std::vector entry; + // entry.push_back(elem_id); + // entry.push_back(elem_type_id); + // elem_type.push_back(entry); + // } + // } + + // // Read EBLOCK + // void ReadEBlock() { + // // Sometimes, DAT files contain two EBLOCKs. Read only the first block. + // if (eblock_is_read) { + // return; + // } + + // // Assumes already start of EBLOCK + // std::istringstream iss(memmap.line); + + // // Only read in SOLID eblocks + // std::transform( + // memmap.line.begin(), memmap.line.end(), memmap.line.begin(), + // [](unsigned char c) { + // return std::tolower(c); + // }); + // if (memmap.line.find("solid") == std::string::npos) { + // return; + // } + + // // Get size of EBLOCK from the last item in the line + // // Example: "EBLOCK,19,SOLID,,3588" + // n_elem = getSizeOfEBLOCK(memmap.line); + + // // we have to allocate memory for the maximum size since we don't know + // that a priori int *elem_data = AllocateArray(n_elem * 30); + // elem_off_arr = MakeNDArray({n_elem + 1}); + // int elem_sz = ReadEBlockMemMap(memmap, elem_off_arr.data(), elem_data, + // n_elem); + + // // wrap the raw data but limit it to the size of the number of elements + // read + // // Note: This is faster than reallocating a new array but will use more + // memory elem_arr = WrapNDarray(elem_data, {elem_sz}); + + // // must mark eblock is read since we can only read one eblock per + // archive file eblock_is_read = true; + // } + + // // Read: + // // KEYOPT, ITYPE, KNUM, VALUE + // void ReadKEYOPTLine() { + // if (debug) { + // std::cout << "reading KEYOPT" << std::endl; + // } + + // // Assumes at KEYOPT line + // std::istringstream iss(memmap.line); + // std::string token; + + // // Skip the first token (KEYOPT) + // std::getline(iss, token, ','); + + // // element number is the first item + // int etype; + // std::getline(iss, token, ','); + // try { + // etype = std::stoi(token); + // } catch (const std::invalid_argument &) { + // // soft error + // std::cerr << "Invalid format in KEYOPT line for elem_num:" << + // memmap.line + // << std::endl; + // return; + // } + + // // Read the rest of the values (Element Type, Key Option, Value) + // std::vector keyopt_vals; + // while (std::getline(iss, token, ',')) { + // try { + // keyopt_vals.push_back(std::stoi(token)); + // } catch (const std::invalid_argument &) { + // // soft error + // std::cerr << "Invalid format in KEYOPT line:" << memmap.line << + // std::endl; return; + // } + // } + + // if (keyopt_vals.size() == 2) { + // keyopt[etype].push_back(keyopt_vals); + // } + // } + + // // Read RLBLOCK + // void ReadRLBLOCK() { + // if (debug) { + // std::cout << "reading RLBLOCK" << std::endl; + // } + + // std::vector set_dat; + + // // Assumes line at RLBLOCK + // std::istringstream iss(memmap.line); + // std::string token; + // // while (std::getline(iss, token, ',')) { + // // set_dat.push_back(std::stoi(token)); + // // } + + // // skip starting RLBLOCK command + // std::getline(iss, token, ','); + + // int nset, maxset, maxitems, nperline; + // try { + // std::getline(iss, token, ','); + // nset = std::stoi(token); + // } catch (...) { + // std::cerr << "Error RLBLOCK line when reading nset:" << memmap.line + // << std::endl; return; + // } + // try { + // std::getline(iss, token, ','); + // maxset = std::stoi(token); + // } catch (...) { + // std::cerr << "Error RLBLOCK line when reading maxset:" << + // memmap.line + // << std::endl; + // return; + // } + // try { + // std::getline(iss, token, ','); + // maxitems = std::stoi(token); + // } catch (...) { + // std::cerr << "Error RLBLOCK line when reading maxitems:" << + // memmap.line + // << std::endl; + // return; + // } + // try { + // std::getline(iss, token, ','); + // nperline = std::stoi(token); + // } catch (...) { + // std::cerr << "Error RLBLOCK line when reading nperline:" << + // memmap.line + // << std::endl; + // return; + // } + + // // Skip next two lines + // memmap.seek_eol(); // (2i8,6g16.9) + // memmap.seek_eol(); // (7g16.9) + + // for (int i = 0; i < nset; i++) { + // // read real constant data in the form of: + // // 2 6 1.00000000 0.566900000E-07 0.00000000 + // 0.000... memmap.read_line(); std::vector rcon; + + // // real constant number + // try { + // rnum.push_back(std::stoi(memmap.line.substr(0, 8))); + // } catch (...) { + // std::cerr << "Failed to read real constant ID when reading:" << + // memmap.line + // << std::endl; + // return; + // } + + // int ncon = std::stoi(memmap.line.substr(8, 8)); + // if (ncon > 6) { + // for (int j = 0; j < 6; j++) { + // rcon.push_back(std::stod(memmap.line.substr(16 + 16 * j, + // 16))); ncon--; + // } + + // memmap.read_line(); + + // while (ncon > 0) { + // for (int j = 0; j < 7 && ncon > 0; j++) { + // rcon.push_back(std::stod(memmap.line.substr(16 * j, + // 16))); ncon--; + // } + // if (ncon > 0) { + // memmap.read_line(); + // } + // } + // } else { + // for (int j = 0; j < ncon; j++) { + // rcon.push_back(std::stod(memmap.line.substr(16 + 16 * j, + // 16))); + // } + // } + + // rdat.push_back(rcon); + // } + // } + + // void ReadNBlock(const int pos) { + + // // Sometimes, DAT files contains multiple node blocks. Read only the + // first block. if (nblock_is_read) { + // return; + // } + // nblock_start = pos; + + // // Get size of NBLOCK + // // Assumes line is at NBLOCK + // // std::cout << "line: " << memmap.line << std::endl; + // try { + // // Number of nodes is last item in string + // n_nodes = std::stoi(memmap.line.substr(memmap.line.rfind(',') + + // 1)); + // } catch (...) { + // std::cerr << "Failed to read number of nodes when reading:" << + // memmap.line + // << std::endl; + // return; + // } + // if (debug) { + // std::cout << "Reading " << n_nodes << " nodes" << std::endl; + // } + + // int *nnum_data = AllocateArray(n_nodes); + // double *nodes_data = AllocateArray(n_nodes * 3); + + // // often angles aren't written, so it makes sense to initialize this to + // // zero + // double *node_angles_data = AllocateArray(n_nodes * 3, true); + + // // Get format of nblock + // memmap.read_line(); + // NodeBlockFormat nfmt = GetNodeBlockFormat(memmap.line); + + // // Return actual number of nodes read and wrap the raw data + // // std::cout << memmap.tellg() << std::endl; + // n_nodes = ReadNBlockMemMap( + // memmap, + // nnum_data, + // nodes_data, + // node_angles_data, + // n_nodes, + // nfmt.d_size, + // nfmt.f_size); + // nnum_arr = WrapNDarray(nnum_data, {n_nodes}); + // nodes_arr = WrapNDarray(nodes_data, {n_nodes, 3}); + // node_angles_arr = WrapNDarray(node_angles_data, {n_nodes, + // 3}); + // // std::cout << memmap.tellg() << std::endl; + + // // Read final line, this is always "N,R5.3,LOC, -1," and store file + // // position. This is used for later access (or rewrite) of the node + // // block. + // memmap.seek_eol(); + // nblock_end = memmap.tellg(); + + // // Must mark nblock is read since we can only read one nblock per + // archive file. nblock_is_read = true; if (debug) { + // // std::cout << "Last line is: " << memmap.line << std::endl; + // std::cout << "NBLOCK complete, read " << n_nodes << " nodes." << + // std::endl; + // } + // } + + // // Read CMBLOCK + // void ReadCMBlock() { + // if (debug) { + // std::cout << "reading CMBLOCK" << std::endl; + // } + + // // Assumes line at CMBLOCK + // std::istringstream iss(memmap.line); + // std::string token; + // std::vector split_line; + + // while (std::getline(iss, token, ',')) { + // split_line.push_back(token); + // } + + // if (split_line.size() < 4) { + // std::cerr << "Poor formatting of CMBLOCK: " << memmap.line << + // std::endl; return; + // } + + // int ncomp; + // try { + // ncomp = std::stoi(split_line[3].substr(0, + // split_line[3].find('!')).c_str()); + // } catch (...) { + // std::cerr << "Poor formatting of CMBLOCK: " << memmap.line << + // std::endl; return; + // } + + // memmap.read_line(); + // int isz; + // try { + // isz = std::stoi(memmap.line.substr( + // memmap.line.find('i') + 1, + // memmap.line.find(')') - memmap.line.find('i') - 1)); + // } catch (...) { + // std::cerr << "Failed to read integer size in CMBLOCK" << std::endl; + // return; + // } + + // int nblock; + // try { + // nblock = std::stoi(memmap.line.substr( + // memmap.line.find('(') + 1, + // memmap.line.find('i') - memmap.line.find('(') - 1)); + // } catch (...) { + // std::cerr << "Failed to read number of integers per line in + // CMBLOCK" << std::endl; return; + // } + + // std::vector component(ncomp); + // std::string tempstr(isz, '\0'); + + // for (int i = 0; i < ncomp; ++i) { + // if (i % nblock == 0) { + // memmap.read_line(); + // } + // try { + // component[i] = std::stoi(memmap.line.substr(isz * (i % nblock), + // isz)); + // } catch (...) { + // std::cerr << "Failed to parse integer from CMBLOCK line: " << + // memmap.line + // << std::endl; + // return; + // } + // } + + // // get component name + // std::string comname = split_line[1]; + // comname.erase(comname.find_last_not_of(' ') + 1, std::string::npos); + + // std::string line_comp_type = split_line[2]; + + // if (line_comp_type.find("NODE") != std::string::npos) { + // node_comps[comname] = InterpretComponent(component); + // } else if (line_comp_type.find("ELEM") != std::string::npos) { + // elem_comps[comname] = InterpretComponent(component); + // } + + // if (debug) { + // std::cout << "Done reading CMBLOCK" << std::endl; + // } + // } + + // Read the LS-DYNA deck by reading the start of each line and checking if + // it matches a keyword. For example, *NODE denotes a node card. + // Note: Cards appear to start with the * char. + void Read() { + int first_char, next_char; + + int position_start, position_end; + while (true) { + // Parse based on the first character rather than reading the entire + // line. It's faster and the parsing logic is always based on the first + // character + + if (memmap.eof()) { #ifdef DEBUG - std::cout << "Reached EOF" << std::endl; + std::cout << "Reached EOF" << std::endl; #endif - break; - } + break; + } - first_char = memmap[0]; + first_char = memmap[0]; #ifdef DEBUG - std::cout << "Read character: " << static_cast(first_char) << std::endl; + std::cout << "Read character: " << static_cast(first_char) + << std::endl; #endif - // Check if starting a card - if (first_char == '*') { - memmap.read_line(); - - if (debug) { - std::cout << "Read *" << std::endl; - } - - // Record element type - if (memmap.line.compare(0, 5, "*NODE") == 0) { - ReadNodeCard(); - } - } else { - if (debug) { - std::cout << "No match, continuing..." << std::endl; - } - // Skip remainder of the line - memmap.seek_eol(); - } + // Check if starting a card + if (first_char == '*') { + memmap.read_line(); + + if (debug) { + std::cout << "Read *" << std::endl; + } + + // Record element type + if (memmap.line.compare(0, 5, "*NODE") == 0) { + ReadNodeCard(); } + } else { + if (debug) { + std::cout << "No match, continuing..." << std::endl; + } + // Skip remainder of the line + memmap.seek_eol(); + } } + } - // Read line and return the position prior to reading the line - int ReadLine() { return memmap.read_line(); } - - // // Convert ansys style connectivity to VTK connectivity - // // type_ref is a mapping between ansys element types and VTK element types - // nb::tuple ToVTK(NDArray type_ref) { - // NDArray offset_arr = MakeNDArray({n_elem + 1}); - // NDArray celltypes_arr = MakeNDArray({n_elem}); - - // // Allocate connectivity: max cell size 20 for VTK_HEXAHEDRAL - // int64_t *cell_data = new int64_t[n_elem * 20]; - - // int loc = ans_to_vtk( - // n_elem, - // elem_arr.data(), - // elem_off_arr.data(), - // type_ref.data(), - // n_nodes, - // nnum_arr.data(), - // offset_arr.data(), - // cell_data, - // celltypes_arr.data()); - - // NDArray cells_arr = WrapNDarray(cell_data, {loc}); - // return nb::make_tuple(offset_arr, celltypes_arr, cells_arr); - // } + // Read line and return the position prior to reading the line + int ReadLine() { return memmap.read_line(); } + + // // Convert ansys style connectivity to VTK connectivity + // // type_ref is a mapping between ansys element types and VTK element types + // nb::tuple ToVTK(NDArray type_ref) { + // NDArray offset_arr = MakeNDArray({n_elem + 1}); + // NDArray celltypes_arr = MakeNDArray({n_elem}); + + // // Allocate connectivity: max cell size 20 for VTK_HEXAHEDRAL + // int64_t *cell_data = new int64_t[n_elem * 20]; + + // int loc = ans_to_vtk( + // n_elem, + // elem_arr.data(), + // elem_off_arr.data(), + // type_ref.data(), + // n_nodes, + // nnum_arr.data(), + // offset_arr.data(), + // cell_data, + // celltypes_arr.data()); - ~Deck() { memmap.close_file(); } + // NDArray cells_arr = WrapNDarray(cell_data, + // {loc}); return nb::make_tuple(offset_arr, celltypes_arr, cells_arr); + // } + + ~Deck() { memmap.close_file(); } }; // Archive @@ -937,69 +942,56 @@ class Deck { // nnum: NDArray[np.int32], // ) -> Tuple[NDArray[U], NDArray[np.uint8], NDArray[U]]: ... -nb::tuple ConvertToVTK( - NDArray elem_arr, - NDArray elem_off_arr, - NDArray type_ref, - NDArray nnum_arr) { - - int n_elem = elem_off_arr.size() - 1; - int n_nodes = nnum_arr.size(); - - NDArray offset_arr = MakeNDArray({n_elem + 1}); - NDArray celltypes_arr = MakeNDArray({n_elem}); - - // Allocate connectivity: max cell size 20 for VTK_HEXAHEDRAL - int64_t *cell_data = new int64_t[n_elem * 20]; - - int loc = ans_to_vtk( - n_elem, - elem_arr.data(), - elem_off_arr.data(), - type_ref.data(), - n_nodes, - nnum_arr.data(), - offset_arr.data(), - cell_data, - celltypes_arr.data()); - - NDArray cells_arr = WrapNDarray(cell_data, {loc}); - return nb::make_tuple(offset_arr, celltypes_arr, cells_arr); +nb::tuple ConvertToVTK(NDArray elem_arr, NDArray elem_off_arr, + NDArray type_ref, NDArray nnum_arr) { + + int n_elem = elem_off_arr.size() - 1; + int n_nodes = nnum_arr.size(); + + NDArray offset_arr = MakeNDArray({n_elem + 1}); + NDArray celltypes_arr = MakeNDArray({n_elem}); + + // Allocate connectivity: max cell size 20 for VTK_HEXAHEDRAL + int64_t *cell_data = new int64_t[n_elem * 20]; + + int loc = ans_to_vtk(n_elem, elem_arr.data(), elem_off_arr.data(), + type_ref.data(), n_nodes, nnum_arr.data(), + offset_arr.data(), cell_data, celltypes_arr.data()); + + NDArray cells_arr = WrapNDarray(cell_data, {loc}); + return nb::make_tuple(offset_arr, celltypes_arr, cells_arr); } NB_MODULE(_reader, m) { - m.def("ans_to_vtk", &ConvertToVTK); - nb::class_(m, "Deck") - .def( - nb::init(), - "fname"_a, - "read_params"_a = false, - "debug"_a = false, - "read_eblock"_a = true, - "This class represents a LS-DYNA keyword deck.") - .def_ro("n_elem", &Archive::n_elem) - .def_ro("elem_type", &Archive::elem_type) - .def_ro("elem", &Archive::elem_arr, nb::rv_policy::automatic) - .def_ro("elem_off", &Archive::elem_off_arr, nb::rv_policy::automatic) - .def_ro("keyopt", &Archive::keyopt) - .def_ro("rdat", &Archive::rdat) - .def_ro("rnum", &Archive::rnum) - .def_ro("elem_comps", &Archive::elem_comps, nb::rv_policy::automatic) - .def_ro("node_comps", &Archive::node_comps, nb::rv_policy::automatic) - .def_ro("nodes", &Archive::nodes_arr, nb::rv_policy::automatic) - .def_ro("node_angles", &Archive::node_angles_arr, nb::rv_policy::automatic) - .def_ro("nnum", &Archive::nnum_arr, nb::rv_policy::automatic) - .def_ro("n_nodes", &Archive::n_nodes) - .def_ro("nblock_start", &Archive::nblock_start) - .def_ro("nblock_end", &Archive::nblock_end) - .def("to_vtk", &Archive::ToVTK) - .def("read", &Archive::Read) - .def("read_line", &Archive::ReadLine) - .def("read_nblock", &Archive::ReadNBlock) - .def("read_rlblock", &Archive::ReadRLBLOCK) - .def("read_keyopt_line", &Archive::ReadKEYOPTLine) - .def("read_et_line", &Archive::ReadETLine) - .def("read_etblock", &Archive::ReadETBlock) - .def("read_cmblock", &Archive::ReadCMBlock) - .def("read_eblock", &Archive::ReadEBlock); + m.def("ans_to_vtk", &ConvertToVTK); + nb::class_(m, "Deck") + .def(nb::init(), "fname"_a, + "read_params"_a = false, "debug"_a = false, "read_eblock"_a = true, + "This class represents a LS-DYNA keyword deck.") + .def_ro("n_elem", &Archive::n_elem) + .def_ro("elem_type", &Archive::elem_type) + .def_ro("elem", &Archive::elem_arr, nb::rv_policy::automatic) + .def_ro("elem_off", &Archive::elem_off_arr, nb::rv_policy::automatic) + .def_ro("keyopt", &Archive::keyopt) + .def_ro("rdat", &Archive::rdat) + .def_ro("rnum", &Archive::rnum) + .def_ro("elem_comps", &Archive::elem_comps, nb::rv_policy::automatic) + .def_ro("node_comps", &Archive::node_comps, nb::rv_policy::automatic) + .def_ro("nodes", &Archive::nodes_arr, nb::rv_policy::automatic) + .def_ro("node_angles", &Archive::node_angles_arr, + nb::rv_policy::automatic) + .def_ro("nnum", &Archive::nnum_arr, nb::rv_policy::automatic) + .def_ro("n_nodes", &Archive::n_nodes) + .def_ro("nblock_start", &Archive::nblock_start) + .def_ro("nblock_end", &Archive::nblock_end) + .def("to_vtk", &Archive::ToVTK) + .def("read", &Archive::Read) + .def("read_line", &Archive::ReadLine) + .def("read_nblock", &Archive::ReadNBlock) + .def("read_rlblock", &Archive::ReadRLBLOCK) + .def("read_keyopt_line", &Archive::ReadKEYOPTLine) + .def("read_et_line", &Archive::ReadETLine) + .def("read_etblock", &Archive::ReadETBlock) + .def("read_cmblock", &Archive::ReadCMBlock) + .def("read_eblock", &Archive::ReadEBlock); } diff --git a/src/vtk_support.cpp b/src/vtk_support.cpp index ad8d26b..6c408c1 100644 --- a/src/vtk_support.cpp +++ b/src/vtk_support.cpp @@ -23,24 +23,24 @@ uint8_t VTK_QUADRATIC_HEXAHEDRON = 25; // Contains data for VTK UnstructuredGrid struct VtkData { - int64_t *offset; - int64_t *cells; - uint8_t *celltypes; - int loc; // current position within cells - int *nref; // conversion between ansys and vtk numbering + int64_t *offset; + int64_t *cells; + uint8_t *celltypes; + int loc; // current position within cells + int *nref; // conversion between ansys and vtk numbering }; struct VtkData vtk_data; // Populate offset, cell type, and prepare the cell array for a cell static inline void add_cell(int n_points, uint8_t celltype) { - vtk_data.offset[0] = vtk_data.loc; - vtk_data.offset++; + vtk_data.offset[0] = vtk_data.loc; + vtk_data.offset++; - vtk_data.celltypes[0] = celltype; - vtk_data.celltypes++; + vtk_data.celltypes[0] = celltype; + vtk_data.celltypes++; - /* printf("Finished adding cell\n"); */ - return; + /* printf("Finished adding cell\n"); */ + return; } /* ============================================================================ @@ -74,34 +74,34 @@ static inline void add_cell(int n_points, uint8_t celltype) { * 19 (3, 7) */ static inline void add_hex(const int *elem, int nnode) { - int i; - bool quad = nnode > 8; - if (quad) { - add_cell(20, VTK_QUADRATIC_HEXAHEDRON); - } else { - add_cell(8, VTK_HEXAHEDRON); - } - - // Always add linear - for (i = 0; i < 8; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + int i; + bool quad = nnode > 8; + if (quad) { + add_cell(20, VTK_QUADRATIC_HEXAHEDRON); + } else { + add_cell(8, VTK_HEXAHEDRON); + } + + // Always add linear + for (i = 0; i < 8; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + } + + // translate connectivity + if (quad) { + for (i = 8; i < nnode; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + /* printf("added %d at %d using node %d\n", vtk_data.nref[elem[i]], + * vtk_data.loc - 1, elem[i]); */ } - - // translate connectivity - if (quad) { - for (i = 8; i < nnode; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; - /* printf("added %d at %d using node %d\n", vtk_data.nref[elem[i]], - * vtk_data.loc - 1, elem[i]); */ - } - // ANSYS sometimes doesn't write 0 at the end of the element block - // and quadratic cells always contain 10 nodes - for (i = nnode; i < 20; i++) { - vtk_data.cells[vtk_data.loc++] = -1; - } + // ANSYS sometimes doesn't write 0 at the end of the element block + // and quadratic cells always contain 10 nodes + for (i = nnode; i < 20; i++) { + vtk_data.cells[vtk_data.loc++] = -1; } + } - return; + return; } /* ============================================================================ @@ -128,60 +128,60 @@ on the edges defined by : (0,1), (1,2), (2,0), (3,4), (4,5), (5,3), (0,3), (1,4), (2,5) */ static inline void add_wedge(const int *elem, int nnode) { - bool quad = nnode > 8; - if (quad) { - add_cell(15, VTK_QUADRATIC_WEDGE); - } else { - add_cell(6, VTK_WEDGE); - } - - // [0, 1, 2, 2, 3, 4, 5, 5] - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[6]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[5]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; - - if (quad) { // todo: check if index > nnode - 1 - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[9]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[8]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[11]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[13]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[12]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[15]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[18]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[17]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[16]]; - } - - return; + bool quad = nnode > 8; + if (quad) { + add_cell(15, VTK_QUADRATIC_WEDGE); + } else { + add_cell(6, VTK_WEDGE); + } + + // [0, 1, 2, 2, 3, 4, 5, 5] + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[6]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[5]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; + + if (quad) { // todo: check if index > nnode - 1 + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[9]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[8]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[11]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[13]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[12]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[15]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[18]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[17]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[16]]; + } + + return; } static inline void add_pyr(const int *elem, int nnode) { - int i; // counter - bool quad = nnode > 8; - if (quad) { - add_cell(13, VTK_QUADRATIC_PYRAMID); - } else { - add_cell(5, VTK_PYRAMID); + int i; // counter + bool quad = nnode > 8; + if (quad) { + add_cell(13, VTK_QUADRATIC_PYRAMID); + } else { + add_cell(5, VTK_PYRAMID); + } + + // [0, 1, 2, 3, 4, X, X, X] + for (i = 0; i < 5; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + } + + if (quad) { // todo: check if index > nnode - 1 + for (i = 8; i < 12; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } - // [0, 1, 2, 3, 4, X, X, X] - for (i = 0; i < 5; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + for (i = 16; i < 20; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } - - if (quad) { // todo: check if index > nnode - 1 - for (i = 8; i < 12; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; - } - - for (i = 16; i < 20; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; - } - } - return; + } + return; } /* ============================================================================ @@ -203,164 +203,164 @@ static inline void add_pyr(const int *elem, int nnode) { * (0,1), (1,2), (2,0), (0,3), (1,3), and (2,3) ============================================================================ */ static inline void add_tet(const int *elem, int nnode) { - bool quad = nnode > 8; - if (quad) { - add_cell(10, VTK_QUADRATIC_TETRA); - } else { - add_cell(4, VTK_TETRA); - } - - // edge nodes - // [0, 1, 2, 2, 3, 3, 3, 3] - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; - - if (quad) { // todo: check if index > nnode - 1 - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[8]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[9]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[11]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[16]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[17]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[18]]; - } - - return; + bool quad = nnode > 8; + if (quad) { + add_cell(10, VTK_QUADRATIC_TETRA); + } else { + add_cell(4, VTK_TETRA); + } + + // edge nodes + // [0, 1, 2, 2, 3, 3, 3, 3] + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; + + if (quad) { // todo: check if index > nnode - 1 + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[8]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[9]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[11]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[16]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[17]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[18]]; + } + + return; } // ANSYS Tetrahedral style [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] static inline void add_tet10(const int *elem, int nnode) { - int i; // counter - bool quad = nnode > 4; - if (quad) { - add_cell(10, VTK_QUADRATIC_TETRA); - } else { - add_cell(4, VTK_TETRA); + int i; // counter + bool quad = nnode > 4; + if (quad) { + add_cell(10, VTK_QUADRATIC_TETRA); + } else { + add_cell(4, VTK_TETRA); + } + + // edge nodes + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[3]]; + + if (quad) { + for (i = 4; i < nnode; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } - - // edge nodes - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[3]]; - - if (quad) { - for (i = 4; i < nnode; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; - } - // ANSYS sometimes doesn't write 0 at the end of the element block - // and quadratic cells always contain 10 nodes - for (i = nnode; i < 10; i++) { - vtk_data.cells[vtk_data.loc++] = -1; - } + // ANSYS sometimes doesn't write 0 at the end of the element block + // and quadratic cells always contain 10 nodes + for (i = nnode; i < 10; i++) { + vtk_data.cells[vtk_data.loc++] = -1; } + } - return; + return; } static inline void add_quad(const int *elem, bool is_quad) { - int i; - int n_points; - if (is_quad) { - n_points = 8; - - if (elem[6] == elem[7]) { // edge case: check if repeated - n_points = 4; - /* printf("Duplicate midside points found...\nCelltype - Linear\n"); */ - add_cell(n_points, VTK_QUAD); - } else { - /* printf("Celltype - Quadratic\n"); */ - add_cell(n_points, VTK_QUADRATIC_QUAD); - } + int i; + int n_points; + if (is_quad) { + n_points = 8; + + if (elem[6] == elem[7]) { // edge case: check if repeated + n_points = 4; + /* printf("Duplicate midside points found...\nCelltype - Linear\n"); */ + add_cell(n_points, VTK_QUAD); } else { - n_points = 4; - /* printf("Celltype - Linear\n"); */ - add_cell(n_points, VTK_QUAD); - } - - for (i = 0; i < n_points; i++) { - /* printf("(%i) %i --> ", i, elem[i]); */ - /* printf("%i, \n", vtk_data.nref[elem[i]]); */ - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + /* printf("Celltype - Quadratic\n"); */ + add_cell(n_points, VTK_QUADRATIC_QUAD); } - /* printf("\n"); */ - - return; + } else { + n_points = 4; + /* printf("Celltype - Linear\n"); */ + add_cell(n_points, VTK_QUAD); + } + + for (i = 0; i < n_points; i++) { + /* printf("(%i) %i --> ", i, elem[i]); */ + /* printf("%i, \n", vtk_data.nref[elem[i]]); */ + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + } + /* printf("\n"); */ + + return; } void add_tri(const int *elem, bool quad) { - if (quad) { - add_cell(6, VTK_QUADRATIC_TRIANGLE); - } else { - add_cell(3, VTK_TRIANGLE); - } + if (quad) { + add_cell(6, VTK_QUADRATIC_TRIANGLE); + } else { + add_cell(3, VTK_TRIANGLE); + } - // edge nodes - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + // edge nodes + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - if (quad) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[5]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[7]]; - } + if (quad) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[5]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[7]]; + } - return; + return; } void add_tri_missing_midside(const int *elem, const int n_mid) { - add_cell(6, VTK_QUADRATIC_TRIANGLE); + add_cell(6, VTK_QUADRATIC_TRIANGLE); - int i; - int elem_indices[] = {4, 5, 7}; + int i; + int elem_indices[] = {4, 5, 7}; - // edge nodes - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + // edge nodes + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - // add midside nodes and populate any missing midside with -1 - for (i = 0; i < n_mid && i < 3; i++) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[elem_indices[i]]]; - } - for (i = 3; i > n_mid; i--) { - vtk_data.cells[vtk_data.loc++] = -1; - } + // add midside nodes and populate any missing midside with -1 + for (i = 0; i < n_mid && i < 3; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[elem_indices[i]]]; + } + for (i = 3; i > n_mid; i--) { + vtk_data.cells[vtk_data.loc++] = -1; + } - return; + return; } static inline void add_line(const int *elem, int nnode) { - bool is_quad; - if (nnode > 2) { - is_quad = elem[2] > 0; - } else { - is_quad = false; - } - - /* printf("is_quad, %i\n", is_quad); */ - if (is_quad) { - add_cell(3, VTK_QUADRATIC_EDGE); - } else { - add_cell(2, VTK_LINE); - } - - // edge nodes - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - if (is_quad) { - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - } + bool is_quad; + if (nnode > 2) { + is_quad = elem[2] > 0; + } else { + is_quad = false; + } + + /* printf("is_quad, %i\n", is_quad); */ + if (is_quad) { + add_cell(3, VTK_QUADRATIC_EDGE); + } else { + add_cell(2, VTK_LINE); + } + + // edge nodes + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + if (is_quad) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + } - return; + return; } static inline void add_point(const int *elem) { - add_cell(1, VTK_VERTEX); - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; - return; + add_cell(1, VTK_VERTEX); + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + return; } /* @@ -435,139 +435,132 @@ on the edges defined by : * * ==========================================================================*/ -int ans_to_vtk( - const int nelem, - const int *elem, - const int *elem_off, - const int *type_ref, - const int nnode, - const int *nnum, - int64_t *offset, - int64_t *cells, - uint8_t *celltypes) { - bool is_quad; - int i; // counter - int nnode_elem; // number of nodes belonging to the element - int off; // location of the element nodes - int etype; // ANSYS element type number - - // index ansys node number to VTK C based compatible indexing - // max node number should be last node - // Consider using a hash table here instead - int *nref = (int *)malloc((nnum[nnode - 1] + 1) * sizeof(int)); - nref[0] = -1; // for missing midside nodes ANSYS uses a node number of 0 - for (i = 0; i < nnode; i++) { - nref[nnum[i]] = i; - } +int ans_to_vtk(const int nelem, const int *elem, const int *elem_off, + const int *type_ref, const int nnode, const int *nnum, + int64_t *offset, int64_t *cells, uint8_t *celltypes) { + bool is_quad; + int i; // counter + int nnode_elem; // number of nodes belonging to the element + int off; // location of the element nodes + int etype; // ANSYS element type number + + // index ansys node number to VTK C based compatible indexing + // max node number should be last node + // Consider using a hash table here instead + int *nref = (int *)malloc((nnum[nnode - 1] + 1) * sizeof(int)); + nref[0] = -1; // for missing midside nodes ANSYS uses a node number of 0 + for (i = 0; i < nnode; i++) { + nref[nnum[i]] = i; + } + + // populate global vtk data + vtk_data.offset = offset; + vtk_data.cells = cells; + vtk_data.celltypes = celltypes; + vtk_data.nref = nref; + vtk_data.loc = 0; + + // Convert each element from ANSYS connectivity to VTK connectivity + for (i = 0; i < nelem; i++) { + // etype + etype = elem[elem_off[i] + 1]; + off = elem_off[i] + 10; // to the start of the nodes + nnode_elem = elem_off[i + 1] - off; // number of nodes in element + + switch (type_ref[etype]) { + case 0: // not supported or not set + add_cell(0, VTK_EMPTY_CELL); + break; + case 1: // point + add_point(&elem[off]); + break; + case 2: // line + add_line(&elem[off], nnode_elem); + break; + case 3: // shell + if (nnode_elem == 3) { + /*General triangular elements*/ + add_tri(&elem[off], false); + + } else if (nnode_elem == 4 && elem[off + 2] == elem[off + 3]) { + /* Degenerated linear quad elements */ + add_tri(&elem[off], false); + + } else if (nnode_elem == 4) { + /* General linear quadrangle */ + add_quad(&elem[off], false); + + } else if (nnode_elem == 6) { + /* Quadratic tri elements */ + add_tri(&elem[off], true); + + } else if (nnode_elem == 8 && elem[off + 2] == elem[off + 3]) { + /* Degenerated quadratic quadrangle elements */ + add_tri(&elem[off], true); + + } else if (nnode_elem == 8) { + /* General quadratic quadrangle */ + add_quad(&elem[off], true); + + } else if (nnode_elem == 5) { + // For shell/plane elements with one extra node. + // Check element SURF 152, with KEYOPT 5,1. + add_quad(&elem[off], false); + + } else if (nnode_elem == 10) { + // For quadratic shell/plane elements with two extra nodes. + // Check element TARGET170. + add_quad(&elem[off], true); + + } else { + // Any other case. Possible when missing midside nodes. + is_quad = nnode_elem > 5; + + // Check degenerate triangle + if (elem[off + 2] == elem[off + 3]) { + if (is_quad) { + add_tri_missing_midside(&elem[off], 3 - (8 - nnode_elem)); + } else { + add_tri(&elem[off], false); + } + } else { + // printf(" The type could not be identified. Check vtk_support.c + // file"); printf("Number of elements is %d\n" , nnode_elem); - // populate global vtk data - vtk_data.offset = offset; - vtk_data.cells = cells; - vtk_data.celltypes = celltypes; - vtk_data.nref = nref; - vtk_data.loc = 0; - - // Convert each element from ANSYS connectivity to VTK connectivity - for (i = 0; i < nelem; i++) { - // etype - etype = elem[elem_off[i] + 1]; - off = elem_off[i] + 10; // to the start of the nodes - nnode_elem = elem_off[i + 1] - off; // number of nodes in element - - switch (type_ref[etype]) { - case 0: // not supported or not set - add_cell(0, VTK_EMPTY_CELL); - break; - case 1: // point - add_point(&elem[off]); - break; - case 2: // line - add_line(&elem[off], nnode_elem); - break; - case 3: // shell - if (nnode_elem == 3) { - /*General triangular elements*/ - add_tri(&elem[off], false); - - } else if (nnode_elem == 4 && elem[off + 2] == elem[off + 3]) { - /* Degenerated linear quad elements */ - add_tri(&elem[off], false); - - } else if (nnode_elem == 4) { - /* General linear quadrangle */ - add_quad(&elem[off], false); - - } else if (nnode_elem == 6) { - /* Quadratic tri elements */ - add_tri(&elem[off], true); - - } else if (nnode_elem == 8 && elem[off + 2] == elem[off + 3]) { - /* Degenerated quadratic quadrangle elements */ - add_tri(&elem[off], true); - - } else if (nnode_elem == 8) { - /* General quadratic quadrangle */ - add_quad(&elem[off], true); - - } else if (nnode_elem == 5) { - // For shell/plane elements with one extra node. - // Check element SURF 152, with KEYOPT 5,1. - add_quad(&elem[off], false); - - } else if (nnode_elem == 10) { - // For quadratic shell/plane elements with two extra nodes. - // Check element TARGET170. - add_quad(&elem[off], true); - - } else { - // Any other case. Possible when missing midside nodes. - is_quad = nnode_elem > 5; - - // Check degenerate triangle - if (elem[off + 2] == elem[off + 3]) { - if (is_quad) { - add_tri_missing_midside(&elem[off], 3 - (8 - nnode_elem)); - } else { - add_tri(&elem[off], false); - } - } else { - // printf(" The type could not be identified. Check vtk_support.c - // file"); printf("Number of elements is %d\n" , nnode_elem); - - // Assume quad - add_quad(&elem[off], is_quad); - } - } - break; - case 4: // solid - /* printf("Adding solid "); */ - if (elem[off + 6] != elem[off + 7]) { // hexahedral - /* printf(" subtype hexahedral\n"); */ - add_hex(&elem[off], nnode_elem); - } else if (elem[off + 5] != elem[off + 6]) { // wedge - /* printf(" subtype wedge\n"); */ - add_wedge(&elem[off], nnode_elem); - } else if (elem[off + 2] != elem[off + 3]) { // pyramid - /* printf(" subtype pyramid\n"); */ - add_pyr(&elem[off], nnode_elem); - } else { // tetrahedral - /* printf(" subtype tetrahedral\n"); */ - add_tet(&elem[off], nnode_elem); - } - break; - case 5: // tetrahedral - /* printf("Adding tetrahedral\n"); */ - add_tet10(&elem[off], nnode_elem); - break; - case 6: // linear line - add_line(&elem[off], false); - // should never reach here - } // end of switch - } // end of loop - - /* printf("Done\n"); */ - offset[nelem] = vtk_data.loc; - - free(nref); - return vtk_data.loc; + // Assume quad + add_quad(&elem[off], is_quad); + } + } + break; + case 4: // solid + /* printf("Adding solid "); */ + if (elem[off + 6] != elem[off + 7]) { // hexahedral + /* printf(" subtype hexahedral\n"); */ + add_hex(&elem[off], nnode_elem); + } else if (elem[off + 5] != elem[off + 6]) { // wedge + /* printf(" subtype wedge\n"); */ + add_wedge(&elem[off], nnode_elem); + } else if (elem[off + 2] != elem[off + 3]) { // pyramid + /* printf(" subtype pyramid\n"); */ + add_pyr(&elem[off], nnode_elem); + } else { // tetrahedral + /* printf(" subtype tetrahedral\n"); */ + add_tet(&elem[off], nnode_elem); + } + break; + case 5: // tetrahedral + /* printf("Adding tetrahedral\n"); */ + add_tet10(&elem[off], nnode_elem); + break; + case 6: // linear line + add_line(&elem[off], false); + // should never reach here + } // end of switch + } // end of loop + + /* printf("Done\n"); */ + offset[nelem] = vtk_data.loc; + + free(nref); + return vtk_data.loc; } diff --git a/src/vtk_support.h b/src/vtk_support.h index b98825e..e314bcb 100644 --- a/src/vtk_support.h +++ b/src/vtk_support.h @@ -3,15 +3,7 @@ #include -int ans_to_vtk( - const int, - const int *, - const int *, - const int *, - const int, - const int *, - int64_t *, - int64_t *, - uint8_t *); +int ans_to_vtk(const int, const int *, const int *, const int *, const int, + const int *, int64_t *, int64_t *, uint8_t *); #endif