Repository: ipc-sim/ipc-toolkit (language: python)

Issue:
CollisionMesh Python constructor requires faces array in 2D
**Describe the bug**
Faces array is required in the CollisionMesh Python constructor even in 2D.

**To Reproduce**
The short script below raises the following TypeError.  Modifiying the constructor to take `array([[0, 1, 2]])` as the face parameter runs normally, even without modifying the points to have an additional z coordinate.
```python
from numpy import array
import ipctk

points = array([[-0.5, 0.5], [0.5, 0.5], [0.5, -0.5]])
edges = array([[0, 1], [1, 2], [2, 0]])

mesh = ipctk.CollisionMesh(points, edges)
```
```
Traceback (most recent call last):
  File "./test.py", line 7, in <module>
    mesh = ipctk.CollisionMesh(points, edges)
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. ipctk.CollisionMesh(rest_positions: numpy.ndarray[numpy.float64[m, n]], edges: numpy.ndarray[numpy.int32[m, n]], faces: numpy.ndarray[numpy.int32[m, n]], displacement_map: scipy.sparse.csc_matrix[numpy.float64] = <0x0 sparse matrix of type '<class 'numpy.float64'>'
        with 0 stored elements in Compressed Sparse Column format>)
    2. ipctk.CollisionMesh(include_vertex: List[bool], full_rest_positions: numpy.ndarray[numpy.float64[m, n]], edges: numpy.ndarray[numpy.int32[m, n]], faces: numpy.ndarray[numpy.int32[m, n]], displacement_map: scipy.sparse.csc_matrix[numpy.float64] = <0x0 sparse matrix of type '<class 'numpy.float64'>'
        with 0 stored elements in Compressed Sparse Column format>)

Invoked with: array([[-0.5,  0.5],
       [ 0.5,  0.5],
       [ 0.5, -0.5]]), array([[0, 1],
       [1, 2],
       [2, 0]])
```
**Expected behavior**
The Getting Started tutorial has a note that mentions that `In 2D only the edges matrix is required. In 3D both edges and faces are required.`, however, this is not reflected in the behaviour.

**Environment (please complete the following information):**
 - OS and Version: Debian Sid
 - Compiler and Version: Python 3.8.19


Discussion / hints:
Thank you for bringing this to my attention. You are correct the face matrix is currently required even for 2D. I will try to fix this. For now, you can pass an empty array (`array([])`) as input for the face matrix and it will work.

Fix patch:
diff --git a/python/src/collision_mesh.cpp b/python/src/collision_mesh.cpp
index 866647823..4467b3562 100644
--- a/python/src/collision_mesh.cpp
+++ b/python/src/collision_mesh.cpp
@@ -21,7 +21,8 @@ void define_collision_mesh(py::module_& m)
                 faces: The faces of the collision mesh (#F × 3).
                 displacement_map: The displacement mapping from displacements on the full mesh to the collision mesh.
             )ipc_Qu8mg5v7",
-            py::arg("rest_positions"), py::arg("edges"), py::arg("faces"),
+            py::arg("rest_positions"), py::arg("edges") = Eigen::MatrixXi(),
+            py::arg("faces") = Eigen::MatrixXi(),
             py::arg("displacement_map") = Eigen::SparseMatrix<double>())
         .def(
             py::init<
@@ -39,7 +40,8 @@ void define_collision_mesh(py::module_& m)
                 displacement_map: The displacement mapping from displacements on the full mesh to the collision mesh.
             )ipc_Qu8mg5v7",
             py::arg("include_vertex"), py::arg("full_rest_positions"),
-            py::arg("edges"), py::arg("faces"),
+            py::arg("edges") = Eigen::MatrixXi(),
+            py::arg("faces") = Eigen::MatrixXi(),
             py::arg("displacement_map") = Eigen::SparseMatrix<double>())
         .def_static(
             "build_from_full_mesh", &CollisionMesh::build_from_full_mesh,
diff --git a/src/ipc/collision_mesh.hpp b/src/ipc/collision_mesh.hpp
index a33837976..f373d5e44 100644
--- a/src/ipc/collision_mesh.hpp
+++ b/src/ipc/collision_mesh.hpp
@@ -21,8 +21,8 @@ class CollisionMesh {
     /// @param displacement_map The displacement mapping from displacements on the full mesh to the collision mesh.
     CollisionMesh(
         const Eigen::MatrixXd& rest_positions,
-        const Eigen::MatrixXi& edges,
-        const Eigen::MatrixXi& faces,
+        const Eigen::MatrixXi& edges = Eigen::MatrixXi(),
+        const Eigen::MatrixXi& faces = Eigen::MatrixXi(),
         const Eigen::SparseMatrix<double>& displacement_map =
             Eigen::SparseMatrix<double>());
 
@@ -35,8 +35,8 @@ class CollisionMesh {
     CollisionMesh(
         const std::vector<bool>& include_vertex,
         const Eigen::MatrixXd& full_rest_positions,
-        const Eigen::MatrixXi& edges,
-        const Eigen::MatrixXi& faces,
+        const Eigen::MatrixXi& edges = Eigen::MatrixXi(),
+        const Eigen::MatrixXi& faces = Eigen::MatrixXi(),
         const Eigen::SparseMatrix<double>& displacement_map =
             Eigen::SparseMatrix<double>());
 
@@ -48,7 +48,7 @@ class CollisionMesh {
     static CollisionMesh build_from_full_mesh(
         const Eigen::MatrixXd& full_rest_positions,
         const Eigen::MatrixXi& edges,
-        const Eigen::MatrixXi& faces)
+        const Eigen::MatrixXi& faces = Eigen::MatrixXi())
     {
         return CollisionMesh(
             construct_is_on_surface(full_rest_positions.rows(), edges),


Test patch:
diff --git a/python/tests/test_collision_mesh.py b/python/tests/test_collision_mesh.py
new file mode 100644
index 000000000..c0d55085b
--- /dev/null
+++ b/python/tests/test_collision_mesh.py
@@ -0,0 +1,87 @@
+import numpy as np
+import scipy
+
+import find_ipctk
+from ipctk import CollisionMesh
+
+
+def test_collision_mesh():
+    V = np.array([[0, 0], [1, 0], [0, 1], [1, 1]], dtype=float)
+    E = np.array([[0, 1], [1, 3], [3, 2], [2, 0]], dtype=int)
+
+    W = scipy.sparse.lil_matrix((4, 3), dtype=float)
+    W[0, 0] = W[1, 1] = W[2, 2] = W[3, 1] = W[3, 2] = 1
+
+    mesh = CollisionMesh(V, E, displacement_map=W)
+
+    U = np.array([[0, 0], [1, 1], [0, 0]], dtype=float)
+
+    Uc = mesh.map_displacements(U)
+    expected_Uc = np.array([[0, 0], [1, 1], [0, 0], [1, 1]], dtype=float)
+    assert (Uc == expected_Uc).all()
+
+    Vc = mesh.displace_vertices(U)
+    expected_Vc = np.array([[0, 0], [2, 1], [0, 1], [2, 2]], dtype=float)
+    assert (Vc == expected_Vc).all()
+
+    g = np.array([1, 1, -1, 1, 1, -1, -1, -1], dtype=float)
+    gf = mesh.to_full_dof(g)
+    expected_gf = np.array([1, 1, -2, 0, 0, -2], dtype=float)
+    assert (gf == expected_gf).all()
+
+    H = np.eye(8)
+    Hf = mesh.to_full_dof(scipy.sparse.csc_matrix(H)).toarray()
+    expected_Hf = np.array([
+        [1, 0, 0, 0, 0, 0],
+        [0, 1, 0, 0, 0, 0],
+        [0, 0, 2, 0, 1, 0],
+        [0, 0, 0, 2, 0, 1],
+        [0, 0, 1, 0, 2, 0],
+        [0, 0, 0, 1, 0, 2],
+    ], dtype=float)
+    assert (Hf == expected_Hf).all()
+
+
+def check_faces_to_edges(E, expected_F2E):
+    F = np.array([[0, 1, 2]])
+    assert (CollisionMesh.construct_faces_to_edges(F, E) == expected_F2E).all()
+
+
+def test_faces_to_edges():
+    yield check_faces_to_edges, np.array([[0, 1], [1, 2], [2, 0]]), np.array([0, 1, 2])
+    yield check_faces_to_edges, np.array([[2, 0], [2, 1], [1, 0]]), np.array([2, 1, 0])
+    yield check_faces_to_edges, np.array([[0, 1], [2, 0], [2, 1]]), np.array([0, 2, 1])
+    # Shouldnt work
+    try:
+        check_faces_to_edges(np.array([[0, 1], [1, 2], [0, 3]]), None)
+        assert False
+    except RuntimeError as e:
+        assert str(e) == "Unable to find edge!"
+    except:
+        assert False
+
+
+def test_codim_points_collision_mesh():
+    V = np.array([[0, 0], [1, 0], [0, 1], [1, 1]], dtype=float)
+    mesh = CollisionMesh(V)
+    expected_codim_vertices = np.array([0, 1, 2, 3], dtype=int)
+    assert (mesh.codim_vertices == expected_codim_vertices).all()
+
+
+def test_collision_mesh_no_faces():
+    # Based on https://github.com/ipc-sim/ipc-toolkit/issues/103
+    points = np.array([[-0.5, 0.5], [0.5, 0.5], [0.5, -0.5]])
+    edges = np.array([[0, 1], [1, 2], [2, 0]])
+    mesh = CollisionMesh(points, edges)
+    assert mesh.faces.size == 0
+
+
+def test_collision_mesh_does_not_segfault():
+    # Based on https://github.com/ipc-sim/ipc-toolkit/issues/102
+    points = np.array([[-0.5, 0.5], [0.5, 0.5], [0.5, -0.5]])
+    edges = np.array([[0, 1], [1, 2], [2, 0]])
+    faces = np.array([[0, 1, 2]])
+    mesh = CollisionMesh(points, edges, faces)
+    assert (mesh.rest_positions == points).all()
+    assert (mesh.edges == edges).all()
+    assert (mesh.faces == faces).all()
diff --git a/tests/src/tests/broad_phase/test_broad_phase.cpp b/tests/src/tests/broad_phase/test_broad_phase.cpp
index ef7d6a145..483671bc2 100644
--- a/tests/src/tests/broad_phase/test_broad_phase.cpp
+++ b/tests/src/tests/broad_phase/test_broad_phase.cpp
@@ -120,7 +120,7 @@ TEST_CASE("Vertex-Vertex Broad Phase", "[ccd][broad_phase][2D]")
     E.row(0) << 1, 0;
     E.row(1) << 2, 3;
 
-    CollisionMesh mesh(V0, E, /*F=*/Eigen::MatrixXi());
+    CollisionMesh mesh(V0, E);
 
     const BroadPhaseMethod method = GENERATE_BROAD_PHASE_METHODS();
 
@@ -144,8 +144,7 @@ TEST_CASE("Broad Phase: 2D Mesh", "[ccd][broad_phase][2D][.]")
     REQUIRE(igl::readCSV((tests::DATA_DIR / "mesh-2D/E.csv").string(), E));
     E.array() -= 1; // NOTE: Convert from OBJ format to index
 
-    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(
-        V0_full, E, /*F=*/Eigen::MatrixXi());
+    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(V0_full, E);
 
     const Eigen::MatrixXd V0 = mesh.vertices(V0_full);
     const Eigen::MatrixXd V1 = mesh.vertices(V1_full);
diff --git a/tests/src/tests/collisions/test_collisions.cpp b/tests/src/tests/collisions/test_collisions.cpp
index 775a00961..f996fa634 100644
--- a/tests/src/tests/collisions/test_collisions.cpp
+++ b/tests/src/tests/collisions/test_collisions.cpp
@@ -24,7 +24,7 @@ TEST_CASE("Codim. vertex-vertex collisions", "[collisions][codim]")
         1, 1, 1;
     vertices.rowwise() -= vertices.colwise().mean();
 
-    CollisionMesh mesh(vertices, Eigen::MatrixXi(), Eigen::MatrixXi());
+    CollisionMesh mesh(vertices);
     mesh.init_area_jacobians();
 
     CHECK(mesh.num_vertices() == 8);
@@ -112,7 +112,7 @@ TEST_CASE("Codim. edge-vertex collisions", "[collisions][codim]")
         0, 3,      //
         0, 4;
 
-    CollisionMesh mesh(vertices, edges, Eigen::MatrixXi());
+    CollisionMesh mesh(vertices, edges);
     mesh.init_area_jacobians();
 
     CHECK(mesh.num_vertices() == 8);
diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp
index 846c0705f..35a5909aa 100644
--- a/tests/src/tests/potential/test_barrier_potential.cpp
+++ b/tests/src/tests/potential/test_barrier_potential.cpp
@@ -416,8 +416,7 @@ TEST_CASE(
     Eigen::MatrixXd X = data["boundary_nodes_pos"];
     Eigen::MatrixXd vertices = data["displaced"];
 
-    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(
-        X, edges, /*faces=*/Eigen::MatrixXi());
+    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(X, edges);
     mesh.init_area_jacobians();
     REQUIRE(mesh.are_area_jacobians_initialized());
 
@@ -537,8 +536,7 @@ TEST_CASE(
     Eigen::MatrixXd X = data["boundary_nodes_pos"];
     Eigen::MatrixXd vertices = data["displaced"];
 
-    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(
-        X, edges, /*faces=*/Eigen::MatrixXi());
+    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(X, edges);
     mesh.init_area_jacobians();
     REQUIRE(mesh.are_area_jacobians_initialized());
 
diff --git a/tests/src/tests/test_collision_mesh.cpp b/tests/src/tests/test_collision_mesh.cpp
index 28f6bcec6..870f3994f 100644
--- a/tests/src/tests/test_collision_mesh.cpp
+++ b/tests/src/tests/test_collision_mesh.cpp
@@ -21,7 +21,7 @@ TEST_CASE("Collision mesh", "[collision_mesh]")
     Eigen::SparseMatrix<double> W(4, 3);
     W.setFromTriplets(weights.begin(), weights.end());
 
-    CollisionMesh mesh(V, E, Eigen::MatrixXi(), W);
+    CollisionMesh mesh(V, E, /*F=*/Eigen::MatrixXi(), W);
 
     Eigen::MatrixXd U(3, 2);
     U << 0, 0, 1, 1, 0, 0;
@@ -43,15 +43,7 @@ TEST_CASE("Collision mesh", "[collision_mesh]")
     expected_gf << 1, 1, -2, 0, 0, -2;
     CHECK(gf == expected_gf);
 
-    Eigen::MatrixXd H(8, 8);
-    H << 1, 0, 0, 0, 0, 0, 0, 0, //
-        0, 1, 0, 0, 0, 0, 0, 0,  //
-        0, 0, 1, 0, 0, 0, 0, 0,  //
-        0, 0, 0, 1, 0, 0, 0, 0,  //
-        0, 0, 0, 0, 1, 0, 0, 0,  //
-        0, 0, 0, 0, 0, 1, 0, 0,  //
-        0, 0, 0, 0, 0, 0, 1, 0,  //
-        0, 0, 0, 0, 0, 0, 0, 1;  //
+    Eigen::MatrixXd H = Eigen::MatrixXd::Identity(8, 8);
     Eigen::MatrixXd Hf =
         mesh.to_full_dof(Eigen::SparseMatrix<double>(H.sparseView()));
 
@@ -113,7 +105,7 @@ TEST_CASE("Codim points collision mesh", "[collision_mesh]")
     Eigen::MatrixXd V(4, 2);
     V << 0, 0, 1, 0, 0, 1, 1, 1;
 
-    CollisionMesh mesh(V, Eigen::MatrixXi(), Eigen::MatrixXi());
+    CollisionMesh mesh(V);
 
     Eigen::VectorXi expected_codim_vertices(4);
     expected_codim_vertices << 0, 1, 2, 3;