diff --git a/src/archive.cpp b/src/archive.cpp index 6214d16..9e94221 100644 --- a/src/archive.cpp +++ b/src/archive.cpp @@ -16,13 +16,6 @@ #include #endif -// Standard VTK int size -#if defined(__linux__) || defined(__APPLE__) -typedef int64_t vtk_int; -#else -typedef int32_t vtk_int; -#endif - // VTK cell types #define VTK_TRIANGLE 5 #define VTK_QUAD 9 @@ -137,8 +130,8 @@ void WriteEblock( 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 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 nodenum_arr, // ANSYS node numbering std::string &mode) { @@ -151,8 +144,8 @@ void WriteEblock( const int *rcon = rcon_arr.data(); const int *elem_nnodes = elem_nnodes_arr.data(); const uint8_t *celltypes = celltypes_arr.data(); - const vtk_int *offset = offset_arr.data(); - const vtk_int *cells = cells_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(); @@ -431,17 +424,17 @@ NDArray CmblockItems(const NDArray array) { // 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 vtk_int *cells, T *points) { - vtk_int ind0 = cells[0] * 3; - vtk_int ind1 = cells[1] * 3; - vtk_int ind2 = cells[2] * 3; - vtk_int ind3 = cells[3] * 3; - vtk_int ind4 = cells[4] * 3; - vtk_int ind5 = cells[5] * 3; - vtk_int ind6 = cells[6] * 3; - vtk_int ind7 = cells[7] * 3; - vtk_int ind8 = cells[8] * 3; - vtk_int ind9 = cells[9] * 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; @@ -462,20 +455,20 @@ template inline void ResetMidTet(const vtk_int *cells, T *points) { // 10(1, 4) // 11(2, 4) // 12(3, 4) -template inline void ResetMidPyr(const vtk_int *cells, T *points) { - vtk_int ind0 = cells[0] * 3; - vtk_int ind1 = cells[1] * 3; - vtk_int ind2 = cells[2] * 3; - vtk_int ind3 = cells[3] * 3; - vtk_int ind4 = cells[4] * 3; - vtk_int ind5 = cells[5] * 3; - vtk_int ind6 = cells[6] * 3; - vtk_int ind7 = cells[7] * 3; - vtk_int ind8 = cells[8] * 3; - vtk_int ind9 = cells[9] * 3; - vtk_int ind10 = cells[10] * 3; - vtk_int ind11 = cells[11] * 3; - vtk_int ind12 = cells[12] * 3; +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; @@ -489,7 +482,7 @@ template inline void ResetMidPyr(const vtk_int *cells, T *points) { } } -template inline void ResetMidWeg(const vtk_int *cells, T *points) { +template inline void ResetMidWeg(const int64_t *cells, T *points) { // Reset midside nodes of a wedge cell: // 6 (0,1) // 7 (1,2) @@ -500,21 +493,21 @@ template inline void ResetMidWeg(const vtk_int *cells, T *points) { // 12 (0,3) // 13 (1,4) // 14 (2,5) - vtk_int ind0 = cells[0] * 3; - vtk_int ind1 = cells[1] * 3; - vtk_int ind2 = cells[2] * 3; - vtk_int ind3 = cells[3] * 3; - vtk_int ind4 = cells[4] * 3; - vtk_int ind5 = cells[5] * 3; - vtk_int ind6 = cells[6] * 3; - vtk_int ind7 = cells[7] * 3; - vtk_int ind8 = cells[8] * 3; - vtk_int ind9 = cells[9] * 3; - vtk_int ind10 = cells[10] * 3; - vtk_int ind11 = cells[11] * 3; - vtk_int ind12 = cells[12] * 3; - vtk_int ind13 = cells[13] * 3; - vtk_int ind14 = cells[14] * 3; + 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; @@ -542,27 +535,27 @@ template inline void ResetMidWeg(const vtk_int *cells, T *points) { // 17 (1,5) // 18 (2,6) // 19 (3,7) -template inline void ResetMidHex(const vtk_int *cells, T *points) { - vtk_int ind0 = cells[0] * 3; - vtk_int ind1 = cells[1] * 3; - vtk_int ind2 = cells[2] * 3; - vtk_int ind3 = cells[3] * 3; - vtk_int ind4 = cells[4] * 3; - vtk_int ind5 = cells[5] * 3; - vtk_int ind6 = cells[6] * 3; - vtk_int ind7 = cells[7] * 3; - vtk_int ind8 = cells[8] * 3; - vtk_int ind9 = cells[9] * 3; - vtk_int ind10 = cells[10] * 3; - vtk_int ind11 = cells[11] * 3; - vtk_int ind12 = cells[12] * 3; - vtk_int ind13 = cells[13] * 3; - vtk_int ind14 = cells[14] * 3; - vtk_int ind15 = cells[15] * 3; - vtk_int ind16 = cells[16] * 3; - vtk_int ind17 = cells[17] * 3; - vtk_int ind18 = cells[18] * 3; - vtk_int ind19 = cells[19] * 3; +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; @@ -583,18 +576,18 @@ template inline void ResetMidHex(const vtk_int *cells, T *points) { template void ResetMidside( NDArray celltypes_arr, - NDArray cells_arr, - NDArray offset_arr, + NDArray cells_arr, + NDArray offset_arr, NDArray points_arr) { int n_cells = celltypes_arr.size(); - const vtk_int *cells = cells_arr.data(); - const vtk_int *offset = offset_arr.data(); + 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++) { - vtk_int c = offset[i]; + int64_t c = offset[i]; switch (celltypes[i]) { case VTK_QUADRATIC_TETRA: diff --git a/src/reader.cpp b/src/reader.cpp index d932893..d4de5de 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -489,10 +489,14 @@ NDArray InterpretComponent(const std::vector &component) { class Archive { - public: + private: + std::string line; + bool debug; std::string filename; + + public: bool read_parameters; - bool debug; + bool read_eblock; bool eblock_is_read = false; bool nblock_is_read = false; @@ -534,8 +538,6 @@ class Archive { // https://nanobind.readthedocs.io/en/latest/faq.html#why-am-i-getting-errors-about-leaked-functions-and-types nb::set_leak_warnings(false); - std::ifstream cfile(filename); - if (!cfile.is_open()) { throw std::runtime_error("No such file or directory: '" + filename + "'"); } @@ -552,8 +554,6 @@ class Archive { std::cout << "reading ET" << std::endl; } - std::string line; - std::getline(cfile, line); std::istringstream iss(line); std::string token; std::vector et_vals; @@ -576,14 +576,12 @@ class Archive { } } - // Read ETBLOCK line + // Read ETBLOCK void ReadETBlock() { - std::string line; std::string token; std::vector values; - // start by reading in the first line, which must be an ETBLOCK - std::getline(cfile, line); + // Assumes current line is an ETBLOCK std::istringstream iss(line); // skip first item (ETBLOCK) @@ -622,17 +620,14 @@ class Archive { } } + // Read EBLOCK void ReadEBlock() { - - // start by reading in the first line, which must be an EBLOCK - std::string line; - std::getline(cfile, line); - // Sometimes, DAT files contain two EBLOCKs. Read only the first block. if (eblock_is_read) { return; } + // Assumes already start of EBLOCK std::istringstream iss(line); // Only read in SOLID eblocks @@ -667,8 +662,7 @@ class Archive { std::cout << "reading KEYOPT" << std::endl; } - std::string line; - std::getline(cfile, line); + // Assumes at KEYOPT line std::istringstream iss(line); std::string token; @@ -703,14 +697,15 @@ class Archive { } } + // Read RLBLOCK void ReadRLBLOCK() { - std::string line; - std::getline(cfile, line); if (debug) { std::cout << "reading RLBLOCK" << std::endl; } std::vector set_dat; + + // Assumes line at RLBLOCK std::istringstream iss(line); std::string token; // while (std::getline(iss, token, ',')) { @@ -797,22 +792,19 @@ class Archive { } } - void ReadNBlock() { - - // Before reading NBLOCK, save where the nblock started - nblock_start = cfile.tellg(); - - // start by reading in the first line, which must be a NBLOCK - std::string line; - std::getline(cfile, line); + 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: " << line << std::endl; try { + // Number of nodes is last item in string n_nodes = std::stoi(line.substr(line.rfind(',') + 1)); } catch (...) { std::cerr << "Failed to read number of nodes when reading:" << line << std::endl; @@ -859,18 +851,13 @@ class Archive { } } + // Read CMBLOCK void ReadCMBlock() { - std::string line; - std::getline(cfile, line); - - if (line.compare(0, 8, "CMBLOCK,") != 0 && line.compare(0, 8, "cmblock,") != 0) { - return; - } - if (debug) { std::cout << "reading CMBLOCK" << std::endl; } + // Assumes line at CMBLOCK std::istringstream iss(line); std::string token; std::vector split_line; @@ -938,11 +925,14 @@ class Archive { } else if (line_comp_type.find("ELEM") != std::string::npos) { elem_comps[comname] = InterpretComponent(component); } + + if (debug) { + std::cout << "Done reading CMBLOCK" << std::endl; + } } void Read() { int first_char, next_char; - std::string line; int position_start, position_end; while (true) { @@ -950,112 +940,84 @@ class Archive { // line. It's faster and the parsing logic is always based on the first // character first_char = cfile.peek(); - // if (debug) { - // std::cout << "Read character: " << static_cast(first_char) << - // std::endl; - // } +#ifdef DEBUG + std::cout << "Read character: " << static_cast(first_char) << std::endl; +#endif if (cfile.eof()) { +#ifdef DEBUG + std::cout << "Reached EOF" << std::endl; +#endif break; } else if (first_char == 'E' || first_char == 'e') { + std::getline(cfile, line); + // E commands (ET or ETBLOCK) if (debug) { std::cout << "Read E" << std::endl; } - // get line but do not advance - position_start = cfile.tellg(); - std::getline(cfile, line); - position_end = cfile.tellg(); - cfile.seekg(position_start); - // Record element type if (line.compare(0, 3, "ET,") == 0 || line.compare(0, 3, "et,") == 0) { ReadETLine(); } else if ( - line.compare(0, 7, "ETBLOCK") == 0 || - line.compare(0, 7, "etblock") == 0) { + line.compare(0, 6, "ETBLOC") == 0 || line.compare(0, 6, "etbloc") == 0) { ReadETBlock(); } else if ( line.compare(0, 6, "EBLOCK") == 0 || line.compare(0, 6, "eblock") == 0 && read_eblock) { ReadEBlock(); - } else { - cfile.seekg(position_end); } } else if (first_char == 'K' || first_char == 'k') { + std::getline(cfile, line); if (debug) { std::cout << "Read K" << std::endl; } - // get line but do not advance - position_start = cfile.tellg(); - std::getline(cfile, line); - position_end = cfile.tellg(); - cfile.seekg(position_start); - // Record keyopt if (line.compare(0, 5, "KEYOP") == 0 || line.compare(0, 5, "keyop") == 0) { ReadKEYOPTLine(); - } else { - cfile.seekg(position_end); } } else if (first_char == 'R' || first_char == 'r') { + std::getline(cfile, line); // test for RLBLOCK if (debug) { std::cout << "Read R" << std::endl; } - // get line but do not advance - position_start = cfile.tellg(); - std::getline(cfile, line); - position_end = cfile.tellg(); - cfile.seekg(position_start); - // Record keyopt if (line.compare(0, 5, "RLBLO") == 0 || line.compare(0, 5, "rlblo") == 0) { ReadRLBLOCK(); - } else { - cfile.seekg(position_end); } } else if (first_char == 'N' || first_char == 'n') { + // store current position if we read in the node block + const int pos = cfile.tellg(); + std::getline(cfile, line); // test for NBLOCK if (debug) { std::cout << "Read N" << std::endl; } - // get line but do not advance - position_start = cfile.tellg(); - std::getline(cfile, line); - position_end = cfile.tellg(); - cfile.seekg(position_start); - // Record node block if (line.compare(0, 5, "NBLOC") == 0 || line.compare(0, 5, "nbloc") == 0) { - ReadNBlock(); - } else { - cfile.seekg(position_end); + ReadNBlock(pos); } + } else if (first_char == 'C' || first_char == 'c') { + std::getline(cfile, line); + if (debug) { std::cout << "Read C" << std::endl; } - // get line but do not advance - position_start = cfile.tellg(); - std::getline(cfile, line); - position_end = cfile.tellg(); - cfile.seekg(position_start); - // Record component block if (line.compare(0, 5, "CMBLO") == 0 || line.compare(0, 5, "cmblo") == 0) { + // std::cout << "Reading CMBLOCK" << std::endl; ReadCMBlock(); - } else { - cfile.seekg(position_end); } } else { if (debug) { @@ -1067,6 +1029,13 @@ class Archive { } } + // Read line and return the position prior to reading the line + int ReadLine() { + int pos = cfile.tellg(); + std::getline(cfile, line); + return pos; + } + // 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) { @@ -1159,9 +1128,9 @@ NB_MODULE(_reader, m) { .def_ro("n_nodes", &Archive::n_nodes) .def_ro("nblock_start", &Archive::nblock_start) .def_ro("nblock_end", &Archive::nblock_end) - .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) diff --git a/tests/test__reader.py b/tests/test__reader.py index 30934b3..e4924f3 100644 --- a/tests/test__reader.py +++ b/tests/test__reader.py @@ -139,6 +139,7 @@ def test_read_et(tmp_path: Path) -> None: fid.write(et_line) archive = _reader.Archive(filename) + archive.read_line() archive.read_et_line() assert np.array_equal([[4, 186]], archive.elem_type) @@ -149,6 +150,7 @@ def test_read_etblock(tmp_path: Path) -> None: fid.write(ETBLOCK_STR) archive = _reader.Archive(filename) + archive.read_line() archive.read_etblock() assert np.array_equal([[1, 181]], archive.elem_type) @@ -159,6 +161,7 @@ def test_read_eblock_not_solid(tmp_path: Path) -> None: fid.write(EBLOCK_STR_NOT_SOLID) archive = _reader.Archive(filename) + archive.read_line() archive.read_eblock() assert archive.n_elem == 0 @@ -169,6 +172,7 @@ def test_read_eblock(tmp_path: Path) -> None: fid.write(EBLOCK_STR) archive = _reader.Archive(filename) + archive.read_line() archive.read_eblock() assert archive.n_elem == 4 @@ -190,7 +194,9 @@ def test_read_keyopt(tmp_path: Path) -> None: fid.write(KEYOPT_STR) archive = _reader.Archive(filename) + archive.read_line() archive.read_keyopt_line() + archive.read_line() archive.read_keyopt_line() assert archive.keyopt[1][0] == [1, 0] assert archive.keyopt[1][1] == [2, 0] @@ -202,6 +208,7 @@ def test_read_rlblock(tmp_path: Path) -> None: fid.write(RLBLOCK_STR) archive = _reader.Archive(filename) + archive.read_line() archive.read_rlblock() assert archive.rnum == [2] # set number assert len(archive.rdat) == 1 @@ -216,10 +223,12 @@ def test_read_nblock(tmp_path: Path) -> None: fid.write(NBLOCK_STR) archive = _reader.Archive(filename, debug=True) - archive.read_nblock() + pos = archive.read_line() + archive.read_nblock(pos) archive.n_nodes == 9 assert np.array_equal(archive.nnum, NBLOCK_NODE_ID_ARR) assert np.allclose(archive.nodes, NBLOCK_POS_ARRAY) + assert archive.nblock_start == pos def test_read_nblock_incomplete(tmp_path: Path) -> None: @@ -228,7 +237,8 @@ def test_read_nblock_incomplete(tmp_path: Path) -> None: fid.write(NBLOCK_INCOMPLETE) archive = _reader.Archive(filename, debug=True) - archive.read_nblock() + pos = archive.read_line() + archive.read_nblock(pos) archive.n_nodes == 3 assert np.array_equal(archive.nnum, [1, 2, 3]) @@ -243,6 +253,7 @@ def test_read_cmblock_node(tmp_path: Path) -> None: fid.write(CMBLOCK_NODE_STR) archive = _reader.Archive(filename, debug=True) + archive.read_line() archive.read_cmblock() assert "INTERFACE" in archive.node_comps assert np.array_equal(archive.node_comps["INTERFACE"], NCOMP_INTERFACE) @@ -254,6 +265,7 @@ def test_read_cmblock_elem(tmp_path: Path) -> None: fid.write(CMBLOCK_ELEM_STR) archive = _reader.Archive(filename, debug=True) + archive.read_line() archive.read_cmblock() assert "ELMISC" in archive.elem_comps assert np.array_equal(archive.elem_comps["ELMISC"], ECOMP_ELMISC) @@ -263,8 +275,8 @@ def test_read_mesh200() -> None: filename = os.path.join(TESTFILES_PATH, "mesh200.cdb") archive = _reader.Archive(filename, debug=True) archive.read() - assert archive.n_elem == 1000 assert archive.n_nodes == 4961 + assert archive.n_elem == 1000 # spot check # 1290 0 0 4.0000000000000E-001 6.0000000000000E-001