Update Local 3rd-Order Tensor Jacobian Layout

cmake/recipes/finite_diff.cmake

@@ -7,7 +7,7 @@ endif()
 message(STATUS "Third-party: creating target 'finitediff::finitediff'")
 
 include(CPM)
-CPMAddPackage("gh:zfergus/finite-diff@1.0.3")
+CPMAddPackage("gh:zfergus/finite-diff@1.0.4")
 
 # Folder name for IDE
 set_target_properties(finitediff_finitediff PROPERTIES FOLDER "ThirdParty")

docs/source/cpp-api/tangent.rst

@@ -28,18 +28,18 @@ Relative Velocity
 Relative Velocity as Multiplier Matricies
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. doxygenfunction:: ipc::point_point_relative_velocity_matrix
-.. doxygenfunction:: ipc::point_edge_relative_velocity_matrix
-.. doxygenfunction:: ipc::edge_edge_relative_velocity_matrix
-.. doxygenfunction:: ipc::point_triangle_relative_velocity_matrix
+.. doxygenfunction:: ipc::point_point_relative_velocity_jacobian
+.. doxygenfunction:: ipc::point_edge_relative_velocity_jacobian
+.. doxygenfunction:: ipc::edge_edge_relative_velocity_jacobian
+.. doxygenfunction:: ipc::point_triangle_relative_velocity_jacobian
 
 Relative Velocity Matrix Jacobians
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. doxygenfunction:: ipc::point_point_relative_velocity_matrix_jacobian
-.. doxygenfunction:: ipc::point_edge_relative_velocity_matrix_jacobian
-.. doxygenfunction:: ipc::edge_edge_relative_velocity_matrix_jacobian
-.. doxygenfunction:: ipc::point_triangle_relative_velocity_matrix_jacobian
+.. doxygenfunction:: ipc::point_point_relative_velocity_dx_dbeta
+.. doxygenfunction:: ipc::point_edge_relative_velocity_dx_dbeta
+.. doxygenfunction:: ipc::edge_edge_relative_velocity_dx_dbeta
+.. doxygenfunction:: ipc::point_triangle_relative_velocity_dx_dbeta
 
 Closet Points
 -------------

docs/source/python-api/tangent.rst

@@ -28,18 +28,18 @@ Relative Velocity
 Relative Velocity as Multiplier Matricies
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. autofunction:: ipctk.point_point_relative_velocity_matrix
-.. autofunction:: ipctk.point_edge_relative_velocity_matrix
-.. autofunction:: ipctk.edge_edge_relative_velocity_matrix
-.. autofunction:: ipctk.point_triangle_relative_velocity_matrix
+.. autofunction:: ipctk.point_point_relative_velocity_jacobian
+.. autofunction:: ipctk.point_edge_relative_velocity_jacobian
+.. autofunction:: ipctk.edge_edge_relative_velocity_jacobian
+.. autofunction:: ipctk.point_triangle_relative_velocity_jacobian
 
 Relative Velocity Matrix Jacobians
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. autofunction:: ipctk.point_point_relative_velocity_matrix_jacobian
-.. autofunction:: ipctk.point_edge_relative_velocity_matrix_jacobian
-.. autofunction:: ipctk.edge_edge_relative_velocity_matrix_jacobian
-.. autofunction:: ipctk.point_triangle_relative_velocity_matrix_jacobian
+.. autofunction:: ipctk.point_point_relative_velocity_dx_dbeta
+.. autofunction:: ipctk.point_edge_relative_velocity_dx_dbeta
+.. autofunction:: ipctk.edge_edge_relative_velocity_dx_dbeta
+.. autofunction:: ipctk.point_triangle_relative_velocity_dx_dbeta
 
 Closet Points
 -------------

notebooks/generate_cpp_code.py

@@ -1,23 +1,36 @@
 import numpy as np
 import sympy
 from sympy import Matrix, MatrixSymbol
+from sympy.core.relational import Relational
+from sympy.logic.boolalg import Boolean
 from sympy.printing import ccode
 import subprocess
+import re
 
 from utils import jacobian
 
 
 def generate_code(expr, out_var_name=None):
-    CSE_results = sympy.cse(
-        expr, sympy.numbered_symbols("t"), optimizations='basic')
+    CSE_results = sympy.cse(expr, sympy.numbered_symbols("t"), optimizations="basic")
     lines = []
     for helper in CSE_results[0]:
+        var_name = sympy.cxxcode(helper[0])
+        value = helper[1]
+
+        # 1. Determine the Correct C++ Type
         if isinstance(helper[1], MatrixSymbol):
-            lines.append(f'const auto {helper[0]}[{helper[1].size}];')
-            lines.append(sympy.cxxcode(helper[1], helper[0]))
+            cpp_type = "const double"
+            lines.append(f"{cpp_type} {var_name}[{value.size}];")
+            lines.append(sympy.cxxcode(value, var_name))
+            continue  # Skip the standard assignment line below
+
+        if isinstance(value, (Relational, Boolean)) or value.is_Boolean:
+            cpp_type = "const bool"
         else:
-            lines.append(
-                f'const auto {sympy.cxxcode(helper[0])} = {sympy.cxxcode(helper[1])};')
+            cpp_type = "const double"
+
+        # 2. Append the formatted assignment
+        lines.append(f"{cpp_type} {var_name} = {sympy.cxxcode(value)};")
 
     if out_var_name != None:
         for i, result in enumerate(CSE_results[1]):
@@ -26,7 +39,18 @@ def generate_code(expr, out_var_name=None):
         for i, result in enumerate(CSE_results[1]):
             lines.append(f"return {sympy.cxxcode(result)};")
 
-    return '\n'.join(lines)
+    code = "\n".join(lines)
+
+    # Regex replacements for better C++ code
+    code = re.sub(r"std::pow\(\s*([^,]+),\s*2\s*\)", r"((\1) * (\1))", code)
+    code = re.sub(
+        r"std::pow\(\s*([^,]+),\s*3\.0\s*/\s*2\.0\s*\)", r"((\1) * std::sqrt(\1))", code
+    )
+    code = re.sub(
+        r"std::pow\(\s*([^,]+),\s*-1\.0\s*/\s*2\.0\s*\)", r"(1.0 / std::sqrt(\1))", code
+    )
+
+    return code
 
 
 class CXXFunctionGenerator:
@@ -79,26 +103,30 @@ def __init__(self, expr, params, name, out_param_name="hess"):
 
 def generate_hpp_file(code_generators, file_name, transformer=None):
     newline = "\n"
-    with open(file_name, 'w') as f:
-        f.write(f"""\
+    with open(file_name, "w") as f:
+        f.write(
+            f"""\
 #pragma once
 
 namespace ipc::autogen{{
     {newline.join(code_generator.signature()
                   for code_generator in code_generators)}
 }}
-""")
+"""
+        )
     subprocess.run(["clang-format", str(file_name), "-i"])
 
 
 def generate_cpp_file(code_generators, file_name):
     newline = "\n"
-    with open(file_name, 'w') as f:
-        f.write(f"""\
+    with open(file_name, "w") as f:
+        f.write(
+            f"""\
 #include <{file_name[:-4]}.hpp>
 
 namespace ipc::autogen{{
     {newline.join(code_generator() for code_generator in code_generators)}
 }}
-""")
+"""
+        )
     subprocess.run(["clang-format", str(file_name), "-i"])

notebooks/tangent_basis_grad.ipynb

@@ -112,14 +112,17 @@
     "cross_x = numpy.cross(numpy.array([1, 0, 0]), p0_to_p1)\n",
     "cross_y = numpy.cross(numpy.array([0, 1, 0]), p0_to_p1)\n",
     "\n",
-    "cross = numpy.array([\n",
-    "    sympy.Piecewise(\n",
-    "        (cross_x[i], sq_norm(cross_x) > sq_norm(cross_y)), \n",
-    "        (cross_y[i], True))\n",
-    "    for i in range(3)])\n",
+    "cross = numpy.array(\n",
+    "    [\n",
+    "        sympy.Piecewise(\n",
+    "            (cross_x[i], sq_norm(cross_x) > sq_norm(cross_y)), (cross_y[i], True)\n",
+    "        )\n",
+    "        for i in range(3)\n",
+    "    ]\n",
+    ")\n",
     "\n",
-    "T[:,0] = normalize(cross)\n",
-    "T[:,1] = normalize(numpy.cross(p0_to_p1, cross))\n",
+    "T[:, 0] = normalize(cross)\n",
+    "T[:, 1] = normalize(numpy.cross(p0_to_p1, cross))\n",
     "\n",
     "sympy.Matrix(T).diff(x).shape"
    ]
@@ -250,8 +253,8 @@
     "ea = ea1 - ea0\n",
     "T[:, 0] = normalize(ea)\n",
     "\n",
-    "normal = numpy.cross(ea, eb1 - eb0);\n",
-    "T[:, 1] = normalize(numpy.cross(normal, ea));\n",
+    "normal = numpy.cross(ea, eb1 - eb0)\n",
+    "T[:, 1] = normalize(numpy.cross(normal, ea))\n",
     "\n",
     "generators.append(CXXJacobianGenerator(T, x, \"edge_edge_tangent_basis_jacobian\"))"
    ]
@@ -328,7 +331,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.1"
+   "version": "3.13.9"
   },
   "varInspector": {
    "cols": {

notebooks/utils.py

@@ -18,9 +18,9 @@ def normalize(x):
 
 
 def jacobian(F, x):
-    J = np.empty((x.size * F.shape[0], F.shape[1]), dtype=object)
+    """Compute the Jacobian of a matrix-valued function F with respect to vector x."""
+    J = np.empty((F.size, x.size), dtype=object)
     for xi in range(x.size):
-        for Fi in range(F.shape[0]):
-            for Fj in range(F.shape[1]):
-                J[xi * F.shape[0] + Fi, Fj] = F[Fi, Fj].diff(x[xi])
+        # Flatten column-major order
+        J[:, xi] = np.array(sympy.Matrix(F).diff(x[xi])).reshape(F.size, order="F")
     return J

python/src/collisions/tangential/tangential_collision.cpp

@@ -78,9 +78,9 @@ void define_tangential_collision(py::module_& m)
             )ipc_Qu8mg5v7",
             "velocities"_a)
         .def(
-            "relative_velocity_matrix",
+            "relative_velocity_jacobian",
             py::overload_cast<>(
-                &TangentialCollision::relative_velocity_matrix, py::const_),
+                &TangentialCollision::relative_velocity_jacobian, py::const_),
             R"ipc_Qu8mg5v7(
             Construct the premultiplier matrix for the relative velocity.
 
@@ -91,9 +91,9 @@ void define_tangential_collision(py::module_& m)
                 A matrix M such that `relative_velocity = M * velocities`.
             )ipc_Qu8mg5v7")
         .def(
-            "relative_velocity_matrix",
+            "relative_velocity_jacobian",
             py::overload_cast<Eigen::ConstRef<VectorMax2d>>(
-                &TangentialCollision::relative_velocity_matrix, py::const_),
+                &TangentialCollision::relative_velocity_jacobian, py::const_),
             R"ipc_Qu8mg5v7(
             Construct the premultiplier matrix for the relative velocity.
 
@@ -105,8 +105,8 @@ void define_tangential_collision(py::module_& m)
             )ipc_Qu8mg5v7",
             "closest_point"_a)
         .def(
-            "relative_velocity_matrix_jacobian",
-            &TangentialCollision::relative_velocity_matrix_jacobian,
+            "relative_velocity_dx_dbeta",
+            &TangentialCollision::relative_velocity_dx_dbeta,
             R"ipc_Qu8mg5v7(
             Construct the Jacobian of the relative velocity premultiplier wrt the closest points.
 

python/src/tangent/relative_velocity.cpp

@@ -21,8 +21,8 @@ void define_relative_velocity(py::module_& m)
         "dp0"_a, "dp1"_a);
 
     m.def(
-        "point_point_relative_velocity_matrix",
-        &point_point_relative_velocity_matrix,
+        "point_point_relative_velocity_jacobian",
+        &point_point_relative_velocity_jacobian,
         R"ipc_Qu8mg5v7(
         Compute the point-point relative velocity premultiplier matrix
 
@@ -35,8 +35,8 @@ void define_relative_velocity(py::module_& m)
         "dim"_a);
 
     m.def(
-        "point_point_relative_velocity_matrix_jacobian",
-        &point_point_relative_velocity_matrix_jacobian,
+        "point_point_relative_velocity_dx_dbeta",
+        &point_point_relative_velocity_dx_dbeta,
         R"ipc_Qu8mg5v7(
         Compute the Jacobian of the relative velocity premultiplier matrix
 
@@ -65,8 +65,8 @@ void define_relative_velocity(py::module_& m)
         "dp"_a, "de0"_a, "de1"_a, "alpha"_a);
 
     m.def(
-        "point_edge_relative_velocity_matrix",
-        &point_edge_relative_velocity_matrix,
+        "point_edge_relative_velocity_jacobian",
+        &point_edge_relative_velocity_jacobian,
         R"ipc_Qu8mg5v7(
         Compute the point-edge relative velocity premultiplier matrix
 
@@ -80,8 +80,8 @@ void define_relative_velocity(py::module_& m)
         "dim"_a, "alpha"_a);
 
     m.def(
-        "point_edge_relative_velocity_matrix_jacobian",
-        &point_edge_relative_velocity_matrix_jacobian,
+        "point_edge_relative_velocity_dx_dbeta",
+        &point_edge_relative_velocity_dx_dbeta,
         R"ipc_Qu8mg5v7(
         Compute the Jacobian of the relative velocity premultiplier matrix
 
@@ -112,34 +112,32 @@ void define_relative_velocity(py::module_& m)
         "dea0"_a, "dea1"_a, "deb0"_a, "deb1"_a, "coords"_a);
 
     m.def(
-        "edge_edge_relative_velocity_matrix",
-        &edge_edge_relative_velocity_matrix,
+        "edge_edge_relative_velocity_jacobian",
+        &edge_edge_relative_velocity_jacobian,
         R"ipc_Qu8mg5v7(
         Compute the edge-edge relative velocity matrix.
 
         Parameters:
-            dim: Dimension (2 or 3)
             coords: Two parametric coordinates of the closest points on the edges
 
         Returns:
             The relative velocity matrix
         )ipc_Qu8mg5v7",
-        "dim"_a, "coords"_a);
+        "coords"_a);
 
     m.def(
-        "edge_edge_relative_velocity_matrix_jacobian",
-        &edge_edge_relative_velocity_matrix_jacobian,
+        "edge_edge_relative_velocity_dx_dbeta",
+        &edge_edge_relative_velocity_dx_dbeta,
         R"ipc_Qu8mg5v7(
         Compute the Jacobian of the edge-edge relative velocity matrix.
 
         Parameters:
-            dim: Dimension (2 or 3)
             coords: Two parametric coordinates of the closest points on the edges
 
         Returns:
             The Jacobian of the relative velocity matrix
         )ipc_Qu8mg5v7",
-        "dim"_a, "coords"_a);
+        "coords"_a);
 
     m.def(
         "point_triangle_relative_velocity", &point_triangle_relative_velocity,
@@ -159,32 +157,30 @@ void define_relative_velocity(py::module_& m)
         "dp"_a, "dt0"_a, "dt1"_a, "dt2"_a, "coords"_a);
 
     m.def(
-        "point_triangle_relative_velocity_matrix",
-        &point_triangle_relative_velocity_matrix,
+        "point_triangle_relative_velocity_jacobian",
+        &point_triangle_relative_velocity_jacobian,
         R"ipc_Qu8mg5v7(
         Compute the point-triangle relative velocity matrix.
 
         Parameters:
-            dim: Dimension (2 or 3)
             coords: Barycentric coordinates of the closest point on the triangle
 
         Returns:
             The relative velocity matrix
         )ipc_Qu8mg5v7",
-        "dim"_a, "coords"_a);
+        "coords"_a);
 
     m.def(
-        "point_triangle_relative_velocity_matrix_jacobian",
-        &point_triangle_relative_velocity_matrix_jacobian,
+        "point_triangle_relative_velocity_dx_dbeta",
+        &point_triangle_relative_velocity_dx_dbeta,
         R"ipc_Qu8mg5v7(
         Compute the Jacobian of the point-triangle relative velocity matrix.
 
         Parameters:
-            dim: Dimension (2 or 3)
             coords: Barycentric coordinates of the closest point on the triangle
 
         Returns:
             The Jacobian of the relative velocity matrix
         )ipc_Qu8mg5v7",
-        "dim"_a, "coords"_a);
+        "coords"_a);
 }