Improved Max Approximator

.github/workflows/continuous.yml

@@ -7,7 +7,7 @@ on:
     paths:
       - '.github/workflows/continuous.yml'
       - 'cmake/**'
-      - 'src/**' 
+      - 'src/**'
       - 'tests/**'
       - 'CMakeLists.txt'
 
@@ -70,7 +70,8 @@ jobs:
       - name: Prepare ccache
         run: |
           ccache --max-size=1.0G
-          ccache -V && ccache --show-stats && ccache --zero-stats
+          ccache -V && ccache --show-config
+          ccache --show-stats && ccache --zero-stats
 
       - name: Configure (Linux/macOS)
         if: runner.os != 'Windows'

python/src/bindings.cpp

@@ -40,6 +40,7 @@ PYBIND11_MODULE(ipctk, m)
     define_point_static_plane(m);
 
     // collisions
+    define_distance_type(m); // define early because it is used next
     define_collision_constraint(m);
     define_collision_constraints(m);
     define_edge_edge_constraint(m);
@@ -49,7 +50,6 @@ PYBIND11_MODULE(ipctk, m)
     define_vertex_vertex_constraint(m);
 
     // distance
-    define_distance_type(m);
     define_edge_edge_mollifier(m);
     define_edge_edge_distance(m);
     define_line_line_distance(m);

python/src/collisions/edge_edge.cpp

@@ -10,11 +10,14 @@ void define_edge_edge_constraint(py::module_& m)
     py::class_<EdgeEdgeConstraint, EdgeEdgeCandidate, CollisionConstraint>(
         m, "EdgeEdgeConstraint")
         .def(
-            py::init<long, long, double>(), "", py::arg("edge0_id"),
-            py::arg("edge1_id"), py::arg("eps_x"))
+            py::init<long, long, double, ipc::EdgeEdgeDistanceType>(), "",
+            py::arg("edge0_id"), py::arg("edge1_id"), py::arg("eps_x"),
+            py::arg("dtype") = ipc::EdgeEdgeDistanceType::AUTO)
         .def(
-            py::init<const EdgeEdgeCandidate&, double>(), "",
-            py::arg("candidate"), py::arg("eps_x"))
+            py::init<
+                const EdgeEdgeCandidate&, double, ipc::EdgeEdgeDistanceType>(),
+            "", py::arg("candidate"), py::arg("eps_x"),
+            py::arg("dtype") = ipc::EdgeEdgeDistanceType::AUTO)
         .def(
             "compute_potential", &EdgeEdgeConstraint::compute_potential, "",
             py::arg("vertices"), py::arg("edges"), py::arg("faces"),

src/ipc/collision_mesh.cpp

@@ -105,8 +105,8 @@ CollisionMesh::CollisionMesh(
     m_faces_to_edges = construct_faces_to_edges(m_faces, m_edges);
 
     init_areas();
+    init_adjacencies();
     // Compute these manually if needed.
-    // init_adjacencies();
     // init_area_jacobian();
 }
 
@@ -152,10 +152,13 @@ Eigen::SparseMatrix<double> CollisionMesh::vertex_matrix_to_dof_matrix(
 void CollisionMesh::init_adjacencies()
 {
     m_vertex_vertex_adjacencies.resize(num_vertices());
+    m_vertex_edge_adjacencies.resize(num_vertices());
     // Edges includes the edges of the faces
     for (int i = 0; i < m_edges.rows(); i++) {
         m_vertex_vertex_adjacencies[m_edges(i, 0)].insert(m_edges(i, 1));
         m_vertex_vertex_adjacencies[m_edges(i, 1)].insert(m_edges(i, 0));
+        m_vertex_edge_adjacencies[m_edges(i, 0)].insert(i);
+        m_vertex_edge_adjacencies[m_edges(i, 1)].insert(i);
     }
 
     m_edge_vertex_adjacencies.resize(m_edges.rows());
@@ -250,59 +253,78 @@ void CollisionMesh::init_areas()
                 (vertex_edge_areas.array() < 0).select(1, vertex_edge_areas),
                 vertex_face_areas);
 
-    // Select the area based on the order face, codim
-    m_edge_areas = (m_edge_areas.array() < 0).select(1, m_edge_areas);
+    for (int i = 0; i < m_edge_areas.size(); i++) {
+        if (m_edge_areas[i] < 0) {
+            // Use the edge length for codim edges
+            const VectorMax3d e0 = m_rest_positions.row(m_edges(i, 0));
+            const VectorMax3d e1 = m_rest_positions.row(m_edges(i, 1));
+            m_edge_areas[i] = (e1 - e0).norm();
+        }
+    }
 }
 
 void CollisionMesh::init_area_jacobians()
 {
-    // Compute vertex areas as the sum of ½ the length of connected edges
+    std::vector<bool> was_vertex_visited(num_vertices(), false);
+    std::vector<bool> was_edge_visited(num_edges(), false);
+
     m_vertex_area_jacobian.resize(
         num_vertices(), Eigen::SparseVector<double>(ndof()));
-    for (int i = 0; i < m_edges.rows(); i++) {
-        const VectorMax3d e0 = m_rest_positions.row(m_edges(i, 0));
-        const VectorMax3d e1 = m_rest_positions.row(m_edges(i, 1));
+    m_edge_area_jacobian.resize(
+        num_edges(), Eigen::SparseVector<double>(ndof()));
+
+    // Compute vertex/edge areas as the sum of ⅓ the area of connected face
+    for (int i = 0; i < m_faces.rows(); i++) {
+        assert(dim() == 3);
+        const Eigen::Vector3d f0 = m_rest_positions.row(m_faces(i, 0));
+        const Eigen::Vector3d f1 = m_rest_positions.row(m_faces(i, 1));
+        const Eigen::Vector3d f2 = m_rest_positions.row(m_faces(i, 2));
 
-        const VectorMax6d edge_len_gradient = edge_length_gradient(e0, e1) / 2;
+        const Vector9d face_area_gradient =
+            triangle_area_gradient(f0, f1, f2) / 3.0;
 
-        for (int j = 0; j < m_edges.cols(); j++) {
+        for (int j = 0; j < m_faces.cols(); ++j) {
+            // compute gradient of area
+
+            was_vertex_visited[m_faces(i, j)] = true;
             local_gradient_to_global_gradient(
-                edge_len_gradient, m_edges.row(i), dim(),
-                m_vertex_area_jacobian[m_edges(i, j)]);
+                face_area_gradient, m_faces.row(i), dim(),
+                m_vertex_area_jacobian[m_faces(i, j)]);
+
+            was_edge_visited[m_faces_to_edges(i, j)] = true;
+            local_gradient_to_global_gradient(
+                face_area_gradient, m_faces.row(i), dim(),
+                m_edge_area_jacobian[m_faces_to_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 Eigen::Vector3d f0 = m_rest_positions.row(m_faces(i, 0));
-            const Eigen::Vector3d f1 = m_rest_positions.row(m_faces(i, 1));
-            const Eigen::Vector3d f2 = m_rest_positions.row(m_faces(i, 2));
-
-            const Vector9d face_area_gradient =
-                triangle_area_gradient(f0, f1, f2) / 3.0;
-
-            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 unvisited vertex areas as the sum of ½ the length of connected
+    // edges
+    for (int i = 0; i < m_edges.rows(); i++) {
+        const int e0i = m_edges(i, 0), e1i = m_edges(i, 1);
+        const VectorMax3d e0 = m_rest_positions.row(e0i);
+        const VectorMax3d e1 = m_rest_positions.row(e1i);
 
-                // compute gradient of area
+        assert(was_vertex_visited[e0i] == was_vertex_visited[e1i]);
+        if (was_vertex_visited[e0i] && was_edge_visited[i]) {
+            continue;
+        }
 
-                local_gradient_to_global_gradient(
-                    face_area_gradient, m_faces.row(i), dim(),
-                    m_vertex_area_jacobian[m_faces(i, j)]);
+        const VectorMax6d edge_len_gradient = edge_length_gradient(e0, e1);
 
+        if (!was_vertex_visited[e0i]) {
+            for (int j = 0; j < m_edges.cols(); j++) {
                 local_gradient_to_global_gradient(
-                    face_area_gradient, m_faces.row(i), dim(),
-                    m_edge_area_jacobian[m_faces_to_edges(i, j)]);
+                    edge_len_gradient / 2, m_edges.row(i), dim(),
+                    m_vertex_area_jacobian[m_edges(i, j)]);
             }
         }
+
+        if (!was_edge_visited[i]) {
+            local_gradient_to_global_gradient(
+                edge_len_gradient, m_edges.row(i), dim(),
+                m_edge_area_jacobian[i]);
+        }
     }
 }
 

src/ipc/collision_mesh.hpp

@@ -164,12 +164,22 @@ class CollisionMesh {
         return m_vertex_vertex_adjacencies;
     }
 
+    /// @brief Get the vertex-edge adjacency matrix.
+    const std::vector<unordered_set<int>>& vertex_edge_adjacencies() const
+    {
+        if (!are_adjacencies_initialized()) {
+            throw std::runtime_error(
+                "Vertex-edge adjacencies not initialized. Call init_adjacencies() first.");
+        }
+        return m_vertex_edge_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.");
+                "Edge-vertex adjacencies not initialized. Call init_adjacencies() first.");
         }
         return m_edge_vertex_adjacencies;
     }
@@ -178,6 +188,7 @@ class CollisionMesh {
     bool are_adjacencies_initialized() const
     {
         return !m_vertex_vertex_adjacencies.empty()
+            && !m_vertex_edge_adjacencies.empty()
             && !m_edge_vertex_adjacencies.empty();
     }
 
@@ -305,7 +316,9 @@ class CollisionMesh {
 
     /// @brief Vertices adjacent to vertices
     std::vector<unordered_set<int>> m_vertex_vertex_adjacencies;
-    /// @brief Edges adjacent to edges
+    /// @brief Edges adjacent to vertices
+    std::vector<unordered_set<int>> m_vertex_edge_adjacencies;
+    /// @brief Vertices adjacent to edges
     std::vector<unordered_set<int>> m_edge_vertex_adjacencies;
 
     // std::vector<std::vector<int>> m_vertices_to_faces;

src/ipc/collisions/collision_constraint.cpp

@@ -4,6 +4,13 @@
 
 namespace ipc {
 
+CollisionConstraint::CollisionConstraint(
+    const double weight, const Eigen::SparseVector<double>& weight_gradient)
+    : weight(weight)
+    , weight_gradient(weight_gradient)
+{
+}
+
 double CollisionConstraint::compute_potential(
     const Eigen::MatrixXd& vertices,
     const Eigen::MatrixXi& edges,

src/ipc/collisions/collision_constraint.hpp

@@ -11,6 +11,12 @@ namespace ipc {
 
 class CollisionConstraint : virtual public CollisionStencil {
 public:
+    CollisionConstraint() = default;
+
+    CollisionConstraint(
+        const double weight,
+        const Eigen::SparseVector<double>& weight_gradient);
+
     virtual ~CollisionConstraint() { }
 
     virtual double compute_potential(