From 7b4931f34272b23d510af6020a1868d4f37352c3 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Wed, 1 Apr 2026 15:37:27 +0200 Subject: [PATCH] add ability to read OpenFOAM hex meshes without connectivity --- src/t8_openfoam/t8_openfoam_reader.hxx | 319 ++++++++++++++++++++++++- 1 file changed, 317 insertions(+), 2 deletions(-) diff --git a/src/t8_openfoam/t8_openfoam_reader.hxx b/src/t8_openfoam/t8_openfoam_reader.hxx index 6d126b1503..db63bc0213 100644 --- a/src/t8_openfoam/t8_openfoam_reader.hxx +++ b/src/t8_openfoam/t8_openfoam_reader.hxx @@ -458,6 +458,126 @@ struct t8_openfoam_reader return true; } + /** + * Compute the cell shape of the OpenFOAM cell + * \param [in] cell_id The id of the cell + * \return The eclass of the cell + */ + t8_eclass_t + get_cell_eclass (const size_t cell_id) + { + /* Count points of the cell. The points are counted for every face, + so a hex for example has 6 (faces) * 4 (vertices per face) = 24 points. */ + int num_points = 0; + for (const auto& face : m_cell_faces[cell_id]) { + num_points += m_face_points[face.first].size (); + } + + /* OpenFOAM only has 3D cells, so we do not check 0-2 dimensional cells. */ + switch (m_cell_faces[cell_id].size ()) { + case 4: + /* This cell is hopefully a tet with 4*3=12 points. */ + if (num_points == 12) + return T8_ECLASS_TET; + break; + case 5: + /* This cell could either be a prism (3 * 4 + 2 * 3 = 18) or pyramid (4 * 3 + 1 * 4 = 16). */ + if (num_points == 18) + return T8_ECLASS_PRISM; + if (num_points == 16) + return T8_ECLASS_PYRAMID; + break; + case 6: + /* This cell is hopefully a hex with 6 * 4 = 24 points. */ + if (num_points == 24) + return T8_ECLASS_HEX; + break; + default: + break; + } + return T8_ECLASS_INVALID; + } + + /** + * Given a face eclass and orientation converts the face id of a t8code face vertex id to the corresponding + * OpenFOAM face id. + * \param [in] face_class The eclass of the face (triangle or quad) + * \param [in] orientation The orientation of the face (int and bool, the first is the OF vertex, + * which corresponds to the t8code vertex 0 and the second is 1, if both face + * normals point in the same direction and 0 otherwise) + * \param [in] t8_face_vertex_id The t8 face vertex id to convert + * \return The converted OF face vertex id + */ + constexpr static u_int8_t + t8_face_vertex_to_of_point (t8_eclass_t face_class, std::pair orientation, + u_int8_t t8_face_vertex_id) + { + /* While t8code uses the z-order for quads, OpenFOAM uses an anti-clockwise numeration. + * So if the face is a quad, we have to switch indices in accordance with this LUT. + * This is not necessary for triangles. */ + if (face_class == T8_ECLASS_QUAD) { + t8_face_vertex_id = quad_conversion[t8_face_vertex_id]; + } + + /* The OpenFOAM face point belonging to the current t8_face_vertex is calculated via the orientation + * (orientation.second tells us if we count up +1 or down -1 depending on the fact that the face normals point + * in the same direction and orientation.first tells us which OF point is located at t8 vertex 0) and the current t8_face_vertex. + * We add t8_eclass_num_vertices[face_class] once, so that we never take a mod of a negative value. */ + const int order = orientation.second ? 1 : -1; + return (order * t8_face_vertex_id + orientation.first + t8_eclass_num_vertices[face_class]) + % t8_eclass_num_vertices[face_class]; + } + + /** + * Finds a list inside a list which does not contain given values + * \param [in] listlist A list of lists (range of ranges) + * \param [in] forbidden A range of forbidden values + * \return The location of the first sublist, which does not contain any of \a forbidden. Empty on failure. + */ + template + static std::optional + find_list_not_containing (const ListList& listlist, const Forbidden& forbidden) + { + /* Find the first index where the corresponding list contains none of the forbidden values. */ + auto it = std::ranges::find_if (listlist, [&] (const auto& vec) { + /* Return true if the range contains none of the forbidden values. */ + for (const auto& forbidden_val : forbidden) { + if (std::ranges::find (vec, forbidden_val) != std::ranges::end (vec)) { + return false; + } + } + return true; + }); + + if (it != std::ranges::end (listlist)) { + return std::distance (std::ranges::begin (listlist), it); + } + return std::nullopt; + } + + /** + * Finds a list inside a list which contains given values + * \param [in] listlist A list of lists (range of ranges) + * \param [in] required A range of required values + * \return The location of the first sublist, which contains all values of \a required. Empty on failure. + */ + template + static std::optional + find_list_containing (const ListList& listlist, const Required& required) + { + /* Find the first index where the corresponding list contains all of the required values. */ + auto it = std::ranges::find_if (listlist, [&] (const auto& vec) { + /* Return true if the range contains all of the required values. */ + return std::ranges::all_of ( + required, [&] (const auto& val) { return std::ranges::find (vec, val) != std::ranges::end (vec); }); + }); + + if (it != std::ranges::end (listlist)) { + return std::distance (std::ranges::begin (listlist), it); + } + return std::nullopt; + } + /** * Build the cmesh from the raw mesh data. * \return True on success @@ -465,8 +585,198 @@ struct t8_openfoam_reader bool build_cmesh () { - t8_global_errorf ("ERROR: Not implemented yet. \n"); - return 0; + T8_ASSERT (m_cell_faces.size () != 0); + T8_ASSERT (m_face_points.size () != 0); + T8_ASSERT (m_points.size () != 0); + + t8_cmesh_init (&m_cmesh); + t8_cmesh_register_geometry (m_cmesh); + /* Reconstruct cells */ + size_t cell_id = 0; + for (auto& cell : m_cell_faces) { + const t8_eclass_t eclass = get_cell_eclass (cell_id); + if (eclass == T8_ECLASS_INVALID) { + t8_errorf ("ERROR: Encountered invalid polyhedral cell with id %li.\n", cell_id); + return false; + } + + /* To make this code more readable, we gather all faces with points for the current cell. */ + std::vector current_cell_face_ids; + current_cell_face_ids.reserve (cell.size ()); + std::vector> current_cell_faces; + current_cell_faces.reserve (cell.size ()); + std::vector current_cell_face_normals; + current_cell_face_normals.reserve (cell.size ()); + for (const auto& i_face : cell) { + current_cell_face_ids.emplace_back (i_face.first); + current_cell_faces.emplace_back (m_face_points[i_face.first]); + current_cell_face_normals.emplace_back (i_face.second); + } + + switch (eclass) { + case T8_ECLASS_HEX: + reconstruct_hex_cell (cell_id, current_cell_face_ids, current_cell_faces, current_cell_face_normals); + break; + default: + SC_ABORT_NOT_REACHED (); + } + ++cell_id; + } + t8_cmesh_commit (m_cmesh, m_comm); + return 1; + } + + /** + * Use the gathered cell information to build a hex cell and add it to the cmesh. + * \param [in] cell_id The id of the cell. + * \param [in] face_ids The ids of the cell faces. + * \param [in] face_point_ids The ids of the points of the cell faces. + * \param [in] face_normals The normals of the cell faces. + */ + void + reconstruct_hex_cell (size_t cell_id, std::vector face_ids, std::vector> face_point_ids, + std::vector face_normals) + { + /* To reconstruct the Hexahedron we will follow these stages: + * 0. Initialize some needed variables. + * 1. Assign the first OpenFOAM face to face 0 of the t8 hexahedron. + * 2. Search for the opposite OpenFOAM face (does not share any vertices with face 0) and make it face 1. + * 3. Search for a third face to find out the orientation of face 1. + * 4. We now have all 8 vertices and can use them to find and orient the last four faces. + * 5. Use the vertices and eclass to set the cmesh cell. + * 6. Use the face information to link the cell to its neighbors and apply the boundary conditions. */ + + /* ------------------------- 0. Initialization ------------------------- */ + + std::array, t8_eclass_num_faces[T8_ECLASS_HEX]> faces {}; + std::array, t8_eclass_num_vertices[T8_ECLASS_HEX]> points {}; + std::array face_orientations {}; + constexpr t8_eclass_t eclass = T8_ECLASS_HEX; + + /** Orientation of the OpenFOAM face with regard to the t8code face. + * t8 corner 0 -> OF corner + * same normal direction = 1, different normal direction = -1 + */ + std::pair orientation; + /** Normal of the t8 face. 1 points outward, 0 inward. */ + bool t8_face_normal; + + /* ------------------------- 1. Assign face 0 ------------------------- */ + t8_face_normal = t8_eclass_face_orientation[eclass][0]; + /* The orientation is 1, if both face normals point in the same direction */ + orientation.second = t8_face_normal == face_normals[0]; + orientation.first = 0; + + /* Now, we assign the first face to be face 0. */ + faces[0].emplace (face_ids[0]); + face_orientations[0] = 0; + const std::span face_0_points = face_point_ids[0]; + for (u_int8_t i_face_point = 0; i_face_point < t8_eclass_num_vertices[T8_ECLASS_QUAD]; ++i_face_point) { + const u_int8_t of_point_id = t8_face_vertex_to_of_point (T8_ECLASS_QUAD, orientation, i_face_point); + points[t8_face_vertex_to_tree_vertex[T8_ECLASS_HEX][0][i_face_point]].emplace (face_0_points[of_point_id]); + } + + /* ------------------------- 2. Find face 1 ------------------------- */ + + /* After assigning face 0, we have to find the opposite face, face 1. + * It should not share any points with face 0. */ + const auto opposite_face = find_list_not_containing (face_point_ids, face_point_ids[0]); + T8_ASSERTF ( + opposite_face.has_value (), + "ERROR: Could not find the opposite face for hex reconstruction. Maybe this cell was falsely classified " + "as a hex\n"); + t8_face_normal = t8_eclass_face_orientation[eclass][1]; + orientation.second = t8_face_normal == face_normals[*opposite_face]; + faces[1].emplace (face_ids[*opposite_face]); + const auto face_1_points = face_point_ids[*opposite_face]; + + /* Now we can drop face 0 and the opposite face, since we do not need it anymore. + * We delete the opposite face first, since deleting the 0th face would shift the ids. */ + face_ids.erase (face_ids.begin () + *opposite_face); + face_point_ids.erase (face_point_ids.begin () + *opposite_face); + face_normals.erase (face_normals.begin () + *opposite_face); + face_ids.erase (face_ids.begin ()); + face_point_ids.erase (face_point_ids.begin ()); + face_normals.erase (face_normals.begin ()); + + /* ------------------------- 3. Compute orientation of face 1 ------------------------- */ + + /* We now have to find the orientation of the face. For this we can take any OpenFOAM face which contains the 0th vertex. + * The 1st vertex will be saved right next to it. */ + const size_t face_0_vertex_0 = points[0].value (); + std::span::iterator found_face_0_vertex_0; + const auto check_face = std::find_if ( + face_point_ids.begin (), face_point_ids.end (), + [&face_0_vertex_0, &found_face_0_vertex_0] (std::span current_face_points) { + found_face_0_vertex_0 = std::find (current_face_points.begin (), current_face_points.end (), face_0_vertex_0); + return found_face_0_vertex_0 != current_face_points.end (); + }); + T8_ASSERT (check_face != face_point_ids.end ()); + + /* We found a face containing the 0th vertex. Now we can check if the next vertex is not part of face 0. If true, + * the next vertex is vertex 1 and therefore vertex 0 of face 1. Otherwise, the previous point is vertex 0 of face 1. + * We iterate cyclic over the face vertices. So if we arrive at the end, we begin at the start. */ + auto next_vertex + = (found_face_0_vertex_0 + 1) == check_face->end () ? check_face->begin () : found_face_0_vertex_0 + 1; + /* Check if face 0 contains this vertex */ + size_t face_1_vertex_0; + if (std::find (face_1_points.begin (), face_1_points.end (), *next_vertex) != face_1_points.end ()) { + face_1_vertex_0 = *next_vertex; + } + else { + const auto previous_vertex + = found_face_0_vertex_0 == check_face->begin () ? check_face->end () : found_face_0_vertex_0 - 1; + face_1_vertex_0 = *previous_vertex; + } + + /* Now we can compute the complete orientation of face 1. */ + const auto face_1_vertex_0_iterator = std::find (face_1_points.begin (), face_1_points.end (), face_1_vertex_0); + T8_ASSERT (face_1_vertex_0_iterator != face_1_points.end ()); + orientation.first = face_1_vertex_0_iterator - face_1_points.begin (); + + /* After that we add the remaining four points, completing the vertices of the cell */ + face_orientations[1] = orientation.first; + for (u_int8_t i_face_point = 0; i_face_point < t8_eclass_num_vertices[T8_ECLASS_QUAD]; ++i_face_point) { + const u_int8_t of_point_id = t8_face_vertex_to_of_point (T8_ECLASS_QUAD, orientation, i_face_point); + points[t8_face_vertex_to_tree_vertex[T8_ECLASS_HEX][1][i_face_point]].emplace (face_1_points[of_point_id]); + } + + /* ------------------------- 4. Add remaining four faces ------------------------- */ + + /* Lastly, we have to match the rest of the faces so we can add the boundaries to the cmesh and + * to get the orientations and faces for neighbor linkage. */ + for (size_t i_t8_face = 2; i_t8_face < t8_eclass_num_faces[T8_ECLASS_HEX]; ++i_t8_face) { + std::array face_vertices; + for (size_t i_face_vertex = 0; i_face_vertex < t8_eclass_num_vertices[T8_ECLASS_QUAD]; ++i_face_vertex) { + face_vertices[i_face_vertex] + = points[t8_face_vertex_to_tree_vertex[T8_ECLASS_HEX][i_t8_face][i_face_vertex]].value (); + } + const auto face_id = find_list_containing (face_point_ids, face_vertices); + T8_ASSERTF ( + face_id.has_value (), + "ERROR: Could not find a specific face for hex reconstruction. Maybe this cell was falsely classified " + "as a hex\n"); + //t8_face_normal = t8_eclass_face_orientation[eclass][1]; + //orientation.second = (face_normals[*face_id] == t8_face_normal); + const size_t vertex_0 = *points[t8_face_vertex_to_tree_vertex[T8_ECLASS_HEX][i_t8_face][0]]; + face_orientations[i_t8_face] + = std::distance (face_point_ids[*face_id].begin (), + std::find (face_point_ids[*face_id].begin (), face_point_ids[*face_id].end (), vertex_0)); + } + + /* ------------------------- 5. Set cmesh cell ------------------------- */ + + std::array vertices {}; + for (size_t i_vertex = 0; i_vertex < t8_eclass_num_vertices[T8_ECLASS_HEX]; ++i_vertex) { + for (size_t i_coord = 0; i_coord < 3; ++i_coord) { + vertices[i_vertex * 3 + i_coord] = m_points[points[i_vertex].value ()][i_coord]; + } + } + t8_cmesh_set_tree_class (m_cmesh, cell_id, T8_ECLASS_HEX); + t8_cmesh_set_tree_vertices (m_cmesh, cell_id, vertices.data (), t8_eclass_num_vertices[T8_ECLASS_HEX]); + + /* ------------------------- 6. Set neighbors and boundary conditions ------------------------- */ + //TODO } /** Path to the OpenFOAM case. */ @@ -476,6 +786,11 @@ struct t8_openfoam_reader /** The assigned communicator. */ sc_MPI_Comm m_comm; + /** Conversion table between OpenFOAM and t8code quad numeration. + * We only need this for quads, since triangles numerated equally. + */ + static constexpr std::array quad_conversion = { { 0, 1, 3, 2 } }; + /** Holds all points of the mesh. point_id -> x, y, z */ std::vector m_points; /** Holds all faces of the mesh. face_id -> point_id 0, ..., point_id n */