Feature: OpenFOAM reader [2/n] - Hex meshes without connectivity

src/t8_openfoam/t8_openfoam_reader.hxx

@@ -464,15 +464,325 @@ 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<u_int8_t, u_int8_t> 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 <typename ListList, typename Forbidden>
+  static std::optional<size_t>
+  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 <typename ListList, typename Required>
+  static std::optional<size_t>
+  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
    */
   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<t8_geometry_linear> (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<size_t> current_cell_face_ids;
+      current_cell_face_ids.reserve (cell.size ());
+      std::vector<std::span<size_t>> current_cell_faces;
+      current_cell_faces.reserve (cell.size ());
+      std::vector<char> 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<size_t> face_ids, std::vector<std::span<size_t>> face_point_ids,
+                        std::vector<char> 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<std::optional<size_t>, t8_eclass_num_faces[T8_ECLASS_HEX]> faces {};
+    std::array<std::optional<size_t>, t8_eclass_num_vertices[T8_ECLASS_HEX]> points {};
+    std::array<int8_t, t8_eclass_num_faces[T8_ECLASS_HEX]> 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 <orientation.first>
+     * same normal direction <orientation.second> = 1, different normal direction<orientation.second> = -1
+     */
+    std::pair<u_int8_t, bool> 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<size_t> 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<size_t>::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<size_t> 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<size_t, t8_eclass_num_vertices[T8_ECLASS_QUAD]> 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<double, t8_eclass_num_vertices[T8_ECLASS_HEX] * 3> 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. */
@@ -482,6 +792,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<int, 4> quad_conversion = { { 0, 1, 3, 2 } };
+
   /** Holds all points of the mesh. point_id -> x, y, z */
   std::vector<t8_3D_vec> m_points;
   /** Holds all faces of the mesh. face_id -> point_id 0, ..., point_id n */