3D triangulation computes edge-to-cells bindings, Tetrahedron::FaceTy…

…pe exposes the ids of its 2 adjacent cells, Tetrahedron::EdgeType exposes the set of cells ids insisting on it (currently without any ordering guarantee)
fdaPDE · Jun 14, 2024 · ec149a3 · ec149a3
1 parent d3e3e4f
commit ec149a3
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 5 deletions.
diff --git a/fdaPDE/geometry/tetrahedron.h b/fdaPDE/geometry/tetrahedron.h
@@ -62,6 +62,7 @@ class Tetrahedron : public Simplex<Triangulation::local_dim, Triangulation::embe
       bool on_boundary() const { return mesh_->is_edge_on_boundary(edge_id_); }
       DVector<int> node_ids() const { return mesh_->edges().row(edge_id_); }
       int id() const { return edge_id_; }
+      const std::unordered_set<int>& adjacent_cells() const { return mesh_->edge_to_cells().at(edge_id_); }
     };
     // a triangulation-aware view of a tetrahedron face
     class FaceType : public Simplex<2, Triangulation::embed_dim> {
@@ -80,6 +81,7 @@ class Tetrahedron : public Simplex<Triangulation::local_dim, Triangulation::embe
       DVector<int> edge_ids() const { return mesh_->face_to_edges().row(face_id_); }
       int id() const { return face_id_; }
       EdgeType edge(int n) const { return EdgeType(mesh_->face_to_edges()(face_id_, n), mesh_); }
+      DVector<int> adjacent_cells() const { return mesh_->face_to_cells().row(face_id_); }
     };
 
     // getters
@@ -128,7 +130,7 @@ class Tetrahedron : public Simplex<Triangulation::local_dim, Triangulation::embe
     face_iterator faces_begin() const { return face_iterator(0, this); }
     face_iterator faces_end() const { return face_iterator(this->n_faces, this); }
    private:
-    int id_ = 0;                     // tetrahedron identifier in the physical mesh
+    int id_ = 0;                    // tetrahedron identifier in the physical mesh
     std::array<int, 6> edge_ids_;   // edges identifiers int the physical mesh
     const Triangulation* mesh_ = nullptr;
     bool boundary_ = false;   // true if the element has at least one vertex on the boundary

diff --git a/fdaPDE/geometry/triangulation.h b/fdaPDE/geometry/triangulation.h
@@ -366,18 +366,24 @@ template <> class Triangulation<3, 3> : public TriangulationBase<3, 3, Triangula
                         if (it == edges_map.end()) {
                             edges_.insert(edges_.end(), edge.begin(), edge.end());
                             face_to_edges_.push_back(edge_id);
+                            edge_to_cells_[edge_id].insert(i);   // store (edge, cell) binding
                             edges_map.emplace(edge, edge_id);
 			    edges_markers.push_back(nodes_markers_[edge[0]] && nodes_markers_[edge[1]]);
                             edge_id++;
                         } else {
                             face_to_edges_.push_back(edges_map.at(edge));
+			    edge_to_cells_[edges_map.at(edge)].insert(i);   // store (edge, cell) binding
                         }
                     }
                 } else {
                     const auto& [h, k] = it->second;
                     // elements k and i are neighgbors (they share a face)
                     neighbors_(k, node_opposite_to_face(h, k)) = i;
                     neighbors_(i, node_opposite_to_face(h, i)) = k;
+		    // store (edge, cell) binding for each edge of this face
+                    for (int edge = 0; edge < n_edges_per_face; edge++) {
+                        edge_to_cells_.at(face_to_edges_[h * n_edges_per_face + edge]).insert(i);
+                    }
                     cell_to_faces_(i, j) = h;
 		    face_to_cells_[2 * h + 1] = i;
 		    faces_markers[h] = false;
@@ -404,6 +410,10 @@ template <> class Triangulation<3, 3> : public TriangulationBase<3, 3, Triangula
     Eigen::Map<const DMatrix<int, Eigen::RowMajor>> face_to_edges() const {
         return Eigen::Map<const DMatrix<int, Eigen::RowMajor>>(face_to_edges_.data(), n_faces_, n_edges_per_face);
     }
+    Eigen::Map<const DMatrix<int, Eigen::RowMajor>> face_to_cells() const {
+        return Eigen::Map<const DMatrix<int, Eigen::RowMajor>>(face_to_cells_.data(), n_faces_, 2);
+    }
+    const std::unordered_map<int, std::unordered_set<int>>& edge_to_cells() const { return edge_to_cells_; }
     const BinaryVector<Dynamic>& boundary_faces() const { return faces_markers_; }
     int n_faces() const { return n_faces_; }
     int n_edges() const { return n_edges_; }
@@ -473,10 +483,11 @@ template <> class Triangulation<3, 3> : public TriangulationBase<3, 3, Triangula
    protected:
     std::vector<int> faces_, edges_;   // nodes (as row indexes in nodes_ matrix) composing each face and edge
     std::vector<int> face_to_cells_;   // for each face, the ids of adjacent cells
-    DMatrix<int, Eigen::RowMajor> cell_to_faces_ {};   // ids of faces composing each cell
-    std::vector<int> face_to_edges_;                   // ids of edges composing each face
-    BinaryVector<fdapde::Dynamic> faces_markers_ {};   // j-th element is 1 \iff face j is on boundary
-    BinaryVector<fdapde::Dynamic> edges_markers_ {};   // j-th element is 1 \iff edge j is on boundary
+    std::unordered_map<int, std::unordered_set<int>> edge_to_cells_;   // for each edge, the ids of insisting cells
+    DMatrix<int, Eigen::RowMajor> cell_to_faces_ {};            // ids of faces composing each cell
+    std::vector<int> face_to_edges_;                            // ids of edges composing each face
+    BinaryVector<fdapde::Dynamic> faces_markers_ {};            // j-th element is 1 \iff face j is on boundary
+    BinaryVector<fdapde::Dynamic> edges_markers_ {};            // j-th element is 1 \iff edge j is on boundary
     int n_faces_ = 0, n_edges_ = 0;
     mutable std::optional<LocationPolicy> location_policy_ {};
 };