Made area jacobians computation on demand

Author: zfergus

python/src/collision_mesh.cpp

@@ -52,19 +52,39 @@ void define_collision_mesh(py::module_& m)
             py::overload_cast<const Eigen::SparseMatrix<double>&>(
                 &CollisionMesh::to_full_dof, py::const_),
             "", py::arg("X"))
+        .def(
+            "init_adjacencies", &CollisionMesh::init_adjacencies,
+            "Initialize vertex-vertex and edge-vertex adjacencies.")
         .def_property_readonly(
             "vertex_vertex_adjacencies",
             &CollisionMesh::vertex_vertex_adjacencies, "")
         .def_property_readonly(
             "edge_vertex_adjacencies", &CollisionMesh::edge_vertex_adjacencies,
             "")
+        .def(
+            "are_adjacencies_initialized",
+            &CollisionMesh::are_adjacencies_initialized,
+            "Determine if the adjacencies have been initialized by calling init_adjacencies().")
         .def(
             "is_vertex_on_boundary", &CollisionMesh::is_vertex_on_boundary, "",
             py::arg("i"))
-        .def("vertex_area", &CollisionMesh::vertex_area, "", py::arg("pi"))
+        .def("vertex_area", &CollisionMesh::vertex_area, "", py::arg("vi"))
         .def_property_readonly("vertex_areas", &CollisionMesh::vertex_areas, "")
+        // TODO: enable these functions when Pybind11 supports
+        // Eigen::SparseVector
+        // .def(
+        //     "vertex_area_gradient", &CollisionMesh::vertex_area_gradient, "",
+        //     py::arg("vi"))
         .def("edge_area", &CollisionMesh::edge_area, "", py::arg("ei"))
         .def_property_readonly("edge_areas", &CollisionMesh::edge_areas, "")
+        // TODO: enable these functions when Pybind11 supports
+        // Eigen::SparseVector
+        // .def(
+        //     "edge_area_gradient", &CollisionMesh::edge_area_gradient, "",
+        //     py::arg("ei"))
+        .def(
+            "are_area_jacobians_initialized",
+            &CollisionMesh::are_area_jacobians_initialized, "")
         .def_static(
             "construct_is_on_surface", &CollisionMesh::construct_is_on_surface,
             "", py::arg("num_vertices"), py::arg("edges"))

src/ipc/collision_mesh.cpp

@@ -104,8 +104,10 @@ CollisionMesh::CollisionMesh(
 
     m_faces_to_edges = construct_faces_to_edges(m_faces, m_edges);
 
-    init_adjacencies();
     init_areas();
+    // Compute these manually if needed.
+    // init_adjacencies();
+    // init_area_jacobian();
 }
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -203,66 +205,40 @@ void CollisionMesh::init_areas()
     // Compute vertex areas as the sum of ½ the length of connected edges
     Eigen::VectorXd vertex_edge_areas =
         Eigen::VectorXd::Constant(num_vertices(), -1);
-    m_vertex_area_jacobian.resize(
-        num_vertices(), Eigen::SparseVector<double>(ndof()));
     for (int i = 0; i < m_edges.rows(); i++) {
         const auto& e0 = m_vertices_at_rest.row(m_edges(i, 0));
         const auto& e1 = m_vertices_at_rest.row(m_edges(i, 1));
         double edge_len = (e1 - e0).norm();
 
-        VectorMax6d edge_len_gradient;
-        edge_length_gradient(e0, e1, edge_len_gradient);
-
         for (int j = 0; j < m_edges.cols(); j++) {
             if (vertex_edge_areas[m_edges(i, j)] < 0) {
                 vertex_edge_areas[m_edges(i, j)] = 0;
             }
             vertex_edge_areas[m_edges(i, j)] += edge_len / 2;
-
-            local_gradient_to_global_gradient(
-                edge_len_gradient / 2, m_edges.row(i), dim(),
-                m_vertex_area_jacobian[m_edges(i, j)]);
         }
     }
 
     // Compute vertex/edge areas as the sum of ⅓ the area of connected face
     Eigen::VectorXd vertex_face_areas =
         Eigen::VectorXd::Constant(num_vertices(), -1);
     m_edge_areas.setConstant(m_edges.rows(), -1);
-    m_edge_area_jacobian.resize(
-        m_edges.rows(), Eigen::SparseVector<double>(ndof()));
     if (dim() == 3) {
         for (int i = 0; i < m_faces.rows(); i++) {
             const auto& f0 = m_vertices_at_rest.row(m_faces(i, 0));
             const auto& f1 = m_vertices_at_rest.row(m_faces(i, 1));
             const auto& f2 = m_vertices_at_rest.row(m_faces(i, 2));
             double face_area = cross(f1 - f0, f2 - f0).norm() / 2;
 
-            VectorMax9d face_area_gradient;
-            triangle_area_gradient(f0, f1, f2, face_area_gradient);
-
             for (int j = 0; j < m_faces.cols(); ++j) {
                 if (vertex_face_areas[m_faces(i, j)] < 0) {
                     vertex_face_areas[m_faces(i, j)] = 0;
-                    // remove the computed value from vertex_edge_areas
-                    m_vertex_area_jacobian[m_faces(i, j)].setZero();
                 }
                 vertex_face_areas[m_faces(i, j)] += face_area / 3;
 
                 if (m_edge_areas[m_faces_to_edges(i, j)] < 0) {
                     m_edge_areas[m_faces_to_edges(i, j)] = 0;
                 }
                 m_edge_areas[m_faces_to_edges(i, j)] += face_area / 3;
-
-                // compute gradient of area
-
-                local_gradient_to_global_gradient(
-                    face_area_gradient / 3, m_faces.row(i), dim(),
-                    m_vertex_area_jacobian[m_faces(i, j)]);
-
-                local_gradient_to_global_gradient(
-                    face_area_gradient / 3, m_faces.row(i), dim(),
-                    m_edge_area_jacobian[m_faces_to_edges(i, j)]);
             }
         }
     }
@@ -278,6 +254,59 @@ void CollisionMesh::init_areas()
     m_edge_areas = (m_edge_areas.array() < 0).select(1, m_edge_areas);
 }
 
+void CollisionMesh::init_area_jacobians()
+{
+    // Compute vertex areas as the sum of ½ the length of connected edges
+    m_vertex_area_jacobian.resize(
+        num_vertices(), Eigen::SparseVector<double>(ndof()));
+    for (int i = 0; i < m_edges.rows(); i++) {
+        const auto& e0 = m_vertices_at_rest.row(m_edges(i, 0));
+        const auto& e1 = m_vertices_at_rest.row(m_edges(i, 1));
+
+        VectorMax6d edge_len_gradient;
+        edge_length_gradient(e0, e1, edge_len_gradient);
+
+        for (int j = 0; j < m_edges.cols(); j++) {
+            local_gradient_to_global_gradient(
+                edge_len_gradient / 2, m_edges.row(i), dim(),
+                m_vertex_area_jacobian[m_edges(i, j)]);
+        }
+    }
+
+    // Compute vertex/edge areas as the sum of ⅓ the area of connected face
+    m_edge_area_jacobian.resize(
+        m_edges.rows(), Eigen::SparseVector<double>(ndof()));
+    if (dim() == 3) {
+        std::vector<bool> visited_vertex_before(num_vertices(), false);
+        for (int i = 0; i < m_faces.rows(); i++) {
+            const auto& f0 = m_vertices_at_rest.row(m_faces(i, 0));
+            const auto& f1 = m_vertices_at_rest.row(m_faces(i, 1));
+            const auto& f2 = m_vertices_at_rest.row(m_faces(i, 2));
+
+            VectorMax9d face_area_gradient;
+            triangle_area_gradient(f0, f1, f2, face_area_gradient);
+
+            for (int j = 0; j < m_faces.cols(); ++j) {
+                if (!visited_vertex_before[m_faces(i, j)]) {
+                    // remove the computed value from vertex_edge_areas
+                    m_vertex_area_jacobian[m_faces(i, j)].setZero();
+                    visited_vertex_before[m_faces(i, j)] = true;
+                }
+
+                // compute gradient of area
+
+                local_gradient_to_global_gradient(
+                    face_area_gradient / 3, m_faces.row(i), dim(),
+                    m_vertex_area_jacobian[m_faces(i, j)]);
+
+                local_gradient_to_global_gradient(
+                    face_area_gradient / 3, m_faces.row(i), dim(),
+                    m_edge_area_jacobian[m_faces_to_edges(i, j)]);
+            }
+        }
+    }
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 
 Eigen::MatrixXd

src/ipc/collision_mesh.hpp

@@ -55,6 +55,14 @@ class CollisionMesh {
             full_vertices_at_rest, edges, faces);
     }
 
+    // The following functions are used to initialize optional data structures.
+
+    /// @brief Initialize vertex-vertex and edge-vertex adjacencies.
+    void init_adjacencies();
+
+    /// @brief Initialize vertex and edge areas.
+    void init_area_jacobians();
+
     /// @brief Destroy the Collision Mesh object
     ~CollisionMesh() { }
 
@@ -152,15 +160,30 @@ class CollisionMesh {
     /// @brief Get the vertex-vertex adjacency matrix.
     const std::vector<unordered_set<int>>& vertex_vertex_adjacencies() const
     {
+        if (!are_adjacencies_initialized()) {
+            throw std::runtime_error(
+                "Vertex-vertex adjacencies not initialized. Call init_adjacencies() first.");
+        }
         return m_vertex_vertex_adjacencies;
     }
 
     /// @brief Get the edge-vertex adjacency matrix.
     const std::vector<unordered_set<int>>& edge_vertex_adjacencies() const
     {
+        if (!are_adjacencies_initialized()) {
+            throw std::runtime_error(
+                "Edge-vertex adjacencies not initialized. Call init_area_jacobians() first.");
+        }
         return m_edge_vertex_adjacencies;
     }
 
+    /// @brief Determine if the adjacencies have been initialized by calling init_adjacencies().
+    bool are_adjacencies_initialized() const
+    {
+        return !m_vertex_vertex_adjacencies.empty()
+            && !m_edge_vertex_adjacencies.empty();
+    }
+
     /// @brief Is a vertex on the boundary of the collision mesh?
     /// @param vi Vertex ID.
     /// @return True if the vertex is on the boundary of the collision mesh.
@@ -183,6 +206,10 @@ class CollisionMesh {
     const Eigen::SparseVector<double>&
     vertex_area_gradient(const size_t vi) const
     {
+        if (!are_area_jacobians_initialized()) {
+            throw std::runtime_error(
+                "Vertex area Jacobian not initialized. Call init_area_jacobians() first.");
+        }
         return m_vertex_area_jacobian[vi];
     }
 
@@ -199,9 +226,19 @@ class CollisionMesh {
     /// @return Gradient of the barycentric area of edge ei wrt the rest positions of all points.
     const Eigen::SparseVector<double>& edge_area_gradient(const size_t ei) const
     {
+        if (!are_area_jacobians_initialized()) {
+            throw std::runtime_error(
+                "Edge area Jacobian not initialized. Call init_area_jacobians() first.");
+        }
         return m_edge_area_jacobian[ei];
     }
 
+    /// @brief Determine if the area Jacobians have been initialized by calling init_area_jacobians().
+    bool are_area_jacobians_initialized() const
+    {
+        return !m_vertex_area_jacobian.empty() && !m_edge_area_jacobian.empty();
+    }
+
     // -----------------------------------------------------------------------
 
     /// @brief Construct a vector of bools indicating whether each vertex is on the surface.
@@ -231,9 +268,6 @@ class CollisionMesh {
     /// @brief Initialize the selection matrix from full vertices/DOF to collision vertices/DOF.
     void init_selection_matrices(const int dim);
 
-    /// @brief Initialize vertex-vertex and edge-vertex adjacencies.
-    void init_adjacencies();
-
     /// @brief Initialize vertex and edge areas.
     void init_areas();
 

src/ipc/collisions/constraints.cpp

@@ -261,6 +261,7 @@ void Constraints::edge_vertex_candiates_to_constraints(
 
         Eigen::SparseVector<double> weight_gradient;
         if (compute_shape_derivatives) {
+            assert(mesh.are_area_jacobians_initialized());
             weight_gradient = use_convergent_formulation
                 ? (mesh.vertex_area_gradient(vi) / 2)
                 : Eigen::SparseVector<double>(V.size());
@@ -330,6 +331,7 @@ void Constraints::edge_edge_candiates_to_constraints(
 
         Eigen::SparseVector<double> weight_gradient;
         if (compute_shape_derivatives) {
+            assert(mesh.are_area_jacobians_initialized());
             weight_gradient = use_convergent_formulation
                 ? ((mesh.edge_area_gradient(eai) + mesh.edge_area_gradient(ebi))
                    / 4)
@@ -437,6 +439,7 @@ void Constraints::face_vertex_candiates_to_constraints(
 
         Eigen::SparseVector<double> weight_gradient;
         if (compute_shape_derivatives) {
+            assert(mesh.are_area_jacobians_initialized());
             weight_gradient = use_convergent_formulation
                 ? (mesh.vertex_area_gradient(vi) / 4)
                 : Eigen::SparseVector<double>(V.size());