Fix shape derivative of mollified edge-edge contact

notebooks/ee-mollifier-shape-derivative.ipynb

@@ -0,0 +1,268 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sympy\n",
+    "\n",
+    "from generate_cpp_code import *"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Mollifier"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\begin{cases} \\frac{x \\left(2 eps_{x} - x\\right)}{eps_{x}^{2}} & \\text{for}\\: eps_{x} > x \\\\1.0 & \\text{otherwise} \\end{cases}$"
+      ],
+      "text/plain": [
+       "Piecewise((x*(2*eps_x - x)/eps_x**2, eps_x > x), (1.0, True))"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x = sympy.Symbol('x')\n",
+    "eps_x = sympy.Symbol('eps_x')\n",
+    "\n",
+    "m = sympy.Piecewise(\n",
+    "    ((-x / eps_x + 2) * x / eps_x, x < eps_x),\n",
+    "    (1.0, True)\n",
+    ")\n",
+    "\n",
+    "m.simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\begin{cases} \\frac{2 \\left(eps_{x} - x\\right)}{eps_{x}^{2}} & \\text{for}\\: eps_{x} > x \\\\0 & \\text{otherwise} \\end{cases}$"
+      ],
+      "text/plain": [
+       "Piecewise((2*(eps_x - x)/eps_x**2, eps_x > x), (0, True))"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "m.diff(x).simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\begin{cases} - \\frac{2}{eps_{x}^{2}} & \\text{for}\\: eps_{x} > x \\\\0 & \\text{otherwise} \\end{cases}$"
+      ],
+      "text/plain": [
+       "Piecewise((-2/eps_x**2, eps_x > x), (0, True))"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "m.diff(x).diff(x).simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "return ((eps_x > x) ? (\n",
+      "   -2*x*(eps_x - x)/std::pow(eps_x, 3)\n",
+      ")\n",
+      ": (\n",
+      "   0\n",
+      "));\n"
+     ]
+    },
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\begin{cases} \\frac{2 x \\left(- eps_{x} + x\\right)}{eps_{x}^{3}} & \\text{for}\\: eps_{x} > x \\\\0 & \\text{otherwise} \\end{cases}$"
+      ],
+      "text/plain": [
+       "Piecewise((2*x*(-eps_x + x)/eps_x**3, eps_x > x), (0, True))"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "print(generate_code(m.diff(eps_x).simplify()))\n",
+    "m.diff(eps_x).simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "return ((eps_x > x) ? (\n",
+      "   -2*(eps_x - 2*x)/std::pow(eps_x, 3)\n",
+      ")\n",
+      ": (\n",
+      "   0\n",
+      "));\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(generate_code(m.diff(x).diff(eps_x).simplify()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# $\\nabla_x \\varepsilon(x)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle scale \\left(\\left(ea0x - ea1x\\right)^{2} + \\left(ea0y - ea1y\\right)^{2} + \\left(ea0z - ea1z\\right)^{2}\\right) \\left(\\left(eb0x - eb1x\\right)^{2} + \\left(eb0y - eb1y\\right)^{2} + \\left(eb0z - eb1z\\right)^{2}\\right)$"
+      ],
+      "text/plain": [
+       "scale*((ea0x - ea1x)**2 + (ea0y - ea1y)**2 + (ea0z - ea1z)**2)*((eb0x - eb1x)**2 + (eb0y - eb1y)**2 + (eb0z - eb1z)**2)"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ea0 = sympy.Matrix(sympy.symbols('ea0x ea0y ea0z'))\n",
+    "ea1 = sympy.Matrix(sympy.symbols('ea1x ea1y ea1z'))\n",
+    "eb0 = sympy.Matrix(sympy.symbols('eb0x eb0y eb0z'))\n",
+    "eb1 = sympy.Matrix(sympy.symbols('eb1x eb1y eb1z'))\n",
+    "\n",
+    "scale = sympy.Symbol('scale')\n",
+    "\n",
+    "eps_x = (scale * (ea0 - ea1).dot(ea0 - ea1) * (eb0 - eb1).dot(eb0 - eb1)).simplify()\n",
+    "eps_x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "\n",
+      "void edge_edge_mollifier_threshold_gradient(double ea0x, double ea0y, double ea0z, double ea1x, double ea1y, double ea1z, double eb0x, double eb0y, double eb0z, double eb1x, double eb1y, double eb1z, double grad[12]){\n",
+      "const auto t0 = ea0x - ea1x;\n",
+      "const auto t1 = eb0x - eb1x;\n",
+      "const auto t2 = eb0y - eb1y;\n",
+      "const auto t3 = eb0z - eb1z;\n",
+      "const auto t4 = 2*scale;\n",
+      "const auto t5 = t4*(std::pow(t1, 2) + std::pow(t2, 2) + std::pow(t3, 2));\n",
+      "const auto t6 = t0*t5;\n",
+      "const auto t7 = ea0y - ea1y;\n",
+      "const auto t8 = t5*t7;\n",
+      "const auto t9 = ea0z - ea1z;\n",
+      "const auto t10 = t5*t9;\n",
+      "const auto t11 = t4*(std::pow(t0, 2) + std::pow(t7, 2) + std::pow(t9, 2));\n",
+      "const auto t12 = t1*t11;\n",
+      "const auto t13 = t11*t2;\n",
+      "const auto t14 = t11*t3;\n",
+      "grad[0] = t6;\n",
+      "grad[1] = t8;\n",
+      "grad[2] = t10;\n",
+      "grad[3] = -t6;\n",
+      "grad[4] = -t8;\n",
+      "grad[5] = -t10;\n",
+      "grad[6] = t12;\n",
+      "grad[7] = t13;\n",
+      "grad[8] = t14;\n",
+      "grad[9] = -t12;\n",
+      "grad[10] = -t13;\n",
+      "grad[11] = -t14;\n",
+      "}\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "x = sympy.Matrix([ea0, ea1, eb0, eb1])\n",
+    "\n",
+    "g = CXXGradientGenerator(eps_x, x, 'edge_edge_mollifier_threshold_gradient')\n",
+    "\n",
+    "print(g())"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

notebooks/generate_cpp_code.py

@@ -7,7 +7,7 @@
 from utils import jacobian
 
 
-def generate_code(expr, out_var_name="J"):
+def generate_code(expr, out_var_name=None):
     CSE_results = sympy.cse(
         expr, sympy.numbered_symbols("t"), optimizations='basic')
     lines = []

python/src/collisions/collision_constraints.cpp

@@ -27,7 +27,7 @@ void define_collision_constraints(py::module_& m)
             )ipc_Qu8mg5v7",
             py::arg("mesh"), py::arg("vertices"), py::arg("dhat"),
             py::arg("dmin") = 0,
-            py::arg("broad_phase_method") = BroadPhaseMethod::HASH_GRID)
+            py::arg("broad_phase_method") = DEFAULT_BROAD_PHASE_METHOD)
         .def(
             "build",
             py::overload_cast<
@@ -147,6 +147,69 @@ void define_collision_constraints(py::module_& m)
                 A reference to the constraint.
             )ipc_Qu8mg5v7",
             py::arg("idx"))
+        .def(
+            "is_vertex_vertex", &CollisionConstraints::is_vertex_vertex,
+            R"ipc_Qu8mg5v7(
+            Get if the constraint at idx is a vertex-vertex constraint.
+
+            Parameters:
+                idx: The index of the constraint.
+
+            Returns:
+                If the constraint at idx is a vertex-vertex constraint.
+            )ipc_Qu8mg5v7",
+            py::arg("idx"))
+        .def(
+            "is_edge_vertex", &CollisionConstraints::is_edge_vertex,
+            R"ipc_Qu8mg5v7(
+            Get if the constraint at idx is an edge-vertex constraint.
+
+            Parameters:
+                idx: The index of the constraint.
+
+            Returns:
+                If the constraint at idx is an edge-vertex constraint.
+            )ipc_Qu8mg5v7",
+            py::arg("idx"))
+        .def(
+            "is_edge_edge", &CollisionConstraints::is_edge_edge,
+            R"ipc_Qu8mg5v7(
+            Get if the constraint at idx is an edge-edge constraint.
+
+            Parameters:
+                idx: The index of the constraint.
+
+            Returns:
+                If the constraint at idx is an edge-edge constraint.
+            )ipc_Qu8mg5v7",
+            py::arg("idx"))
+        .def(
+            "is_face_vertex", &CollisionConstraints::is_face_vertex,
+            R"ipc_Qu8mg5v7(
+            Get if the constraint at idx is an face-vertex constraint.
+
+            Parameters:
+                idx: The index of the constraint.
+
+            Returns:
+                If the constraint at idx is an face-vertex constraint.
+            )ipc_Qu8mg5v7",
+            py::arg("idx"))
+        .def(
+            "is_plane_vertex", &CollisionConstraints::is_plane_vertex,
+            R"ipc_Qu8mg5v7(
+            Get if the constraint at idx is an plane-vertex constraint.
+
+            Parameters:
+                idx: The index of the constraint.
+
+            Returns:
+                If the constraint at idx is an plane-vertex constraint.
+            )ipc_Qu8mg5v7",
+            py::arg("idx"))
+        .def(
+            "to_string", &CollisionConstraints::to_string, "", py::arg("mesh"),
+            py::arg("vertices"))
         .def_property(
             "use_convergent_formulation",
             &CollisionConstraints::use_convergent_formulation,