Skip to content

Commit bad4c16

Browse files
zferguszfergus
authored
Fix shape derivative of mollified edge-edge contact (#76)
* Fix shape derivative of mollified edge-edge contact
1 parent 8768a78 commit bad4c16

15 files changed

+1210
-51
lines changed

‎notebooks/ee-mollifier-shape-derivative.ipynb‎

Lines changed: 268 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,268 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 1,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import sympy\n",
10+
"\n",
11+
"from generate_cpp_code import *"
12+
]
13+
},
14+
{
15+
"cell_type": "markdown",
16+
"metadata": {},
17+
"source": [
18+
"# Mollifier"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": 2,
24+
"metadata": {},
25+
"outputs": [
26+
{
27+
"data": {
28+
"text/latex": [
29+
"$\\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}$"
30+
],
31+
"text/plain": [
32+
"Piecewise((x*(2*eps_x - x)/eps_x**2, eps_x > x), (1.0, True))"
33+
]
34+
},
35+
"execution_count": 2,
36+
"metadata": {},
37+
"output_type": "execute_result"
38+
}
39+
],
40+
"source": [
41+
"x = sympy.Symbol('x')\n",
42+
"eps_x = sympy.Symbol('eps_x')\n",
43+
"\n",
44+
"m = sympy.Piecewise(\n",
45+
" ((-x / eps_x + 2) * x / eps_x, x < eps_x),\n",
46+
" (1.0, True)\n",
47+
")\n",
48+
"\n",
49+
"m.simplify()"
50+
]
51+
},
52+
{
53+
"cell_type": "code",
54+
"execution_count": 3,
55+
"metadata": {},
56+
"outputs": [
57+
{
58+
"data": {
59+
"text/latex": [
60+
"$\\displaystyle \\begin{cases} \\frac{2 \\left(eps_{x} - x\\right)}{eps_{x}^{2}} & \\text{for}\\: eps_{x} > x \\\\0 & \\text{otherwise} \\end{cases}$"
61+
],
62+
"text/plain": [
63+
"Piecewise((2*(eps_x - x)/eps_x**2, eps_x > x), (0, True))"
64+
]
65+
},
66+
"execution_count": 3,
67+
"metadata": {},
68+
"output_type": "execute_result"
69+
}
70+
],
71+
"source": [
72+
"m.diff(x).simplify()"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": 4,
78+
"metadata": {},
79+
"outputs": [
80+
{
81+
"data": {
82+
"text/latex": [
83+
"$\\displaystyle \\begin{cases} - \\frac{2}{eps_{x}^{2}} & \\text{for}\\: eps_{x} > x \\\\0 & \\text{otherwise} \\end{cases}$"
84+
],
85+
"text/plain": [
86+
"Piecewise((-2/eps_x**2, eps_x > x), (0, True))"
87+
]
88+
},
89+
"execution_count": 4,
90+
"metadata": {},
91+
"output_type": "execute_result"
92+
}
93+
],
94+
"source": [
95+
"m.diff(x).diff(x).simplify()"
96+
]
97+
},
98+
{
99+
"cell_type": "code",
100+
"execution_count": 5,
101+
"metadata": {},
102+
"outputs": [
103+
{
104+
"name": "stdout",
105+
"output_type": "stream",
106+
"text": [
107+
"return ((eps_x > x) ? (\n",
108+
" -2*x*(eps_x - x)/std::pow(eps_x, 3)\n",
109+
")\n",
110+
": (\n",
111+
" 0\n",
112+
"));\n"
113+
]
114+
},
115+
{
116+
"data": {
117+
"text/latex": [
118+
"$\\displaystyle \\begin{cases} \\frac{2 x \\left(- eps_{x} + x\\right)}{eps_{x}^{3}} & \\text{for}\\: eps_{x} > x \\\\0 & \\text{otherwise} \\end{cases}$"
119+
],
120+
"text/plain": [
121+
"Piecewise((2*x*(-eps_x + x)/eps_x**3, eps_x > x), (0, True))"
122+
]
123+
},
124+
"execution_count": 5,
125+
"metadata": {},
126+
"output_type": "execute_result"
127+
}
128+
],
129+
"source": [
130+
"print(generate_code(m.diff(eps_x).simplify()))\n",
131+
"m.diff(eps_x).simplify()"
132+
]
133+
},
134+
{
135+
"cell_type": "code",
136+
"execution_count": 6,
137+
"metadata": {},
138+
"outputs": [
139+
{
140+
"name": "stdout",
141+
"output_type": "stream",
142+
"text": [
143+
"return ((eps_x > x) ? (\n",
144+
" -2*(eps_x - 2*x)/std::pow(eps_x, 3)\n",
145+
")\n",
146+
": (\n",
147+
" 0\n",
148+
"));\n"
149+
]
150+
}
151+
],
152+
"source": [
153+
"print(generate_code(m.diff(x).diff(eps_x).simplify()))"
154+
]
155+
},
156+
{
157+
"cell_type": "markdown",
158+
"metadata": {},
159+
"source": [
160+
"# $\\nabla_x \\varepsilon(x)$"
161+
]
162+
},
163+
{
164+
"cell_type": "code",
165+
"execution_count": 7,
166+
"metadata": {},
167+
"outputs": [
168+
{
169+
"data": {
170+
"text/latex": [
171+
"$\\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)$"
172+
],
173+
"text/plain": [
174+
"scale*((ea0x - ea1x)**2 + (ea0y - ea1y)**2 + (ea0z - ea1z)**2)*((eb0x - eb1x)**2 + (eb0y - eb1y)**2 + (eb0z - eb1z)**2)"
175+
]
176+
},
177+
"execution_count": 7,
178+
"metadata": {},
179+
"output_type": "execute_result"
180+
}
181+
],
182+
"source": [
183+
"ea0 = sympy.Matrix(sympy.symbols('ea0x ea0y ea0z'))\n",
184+
"ea1 = sympy.Matrix(sympy.symbols('ea1x ea1y ea1z'))\n",
185+
"eb0 = sympy.Matrix(sympy.symbols('eb0x eb0y eb0z'))\n",
186+
"eb1 = sympy.Matrix(sympy.symbols('eb1x eb1y eb1z'))\n",
187+
"\n",
188+
"scale = sympy.Symbol('scale')\n",
189+
"\n",
190+
"eps_x = (scale * (ea0 - ea1).dot(ea0 - ea1) * (eb0 - eb1).dot(eb0 - eb1)).simplify()\n",
191+
"eps_x"
192+
]
193+
},
194+
{
195+
"cell_type": "code",
196+
"execution_count": 8,
197+
"metadata": {},
198+
"outputs": [
199+
{
200+
"name": "stdout",
201+
"output_type": "stream",
202+
"text": [
203+
"\n",
204+
"\n",
205+
"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",
206+
"const auto t0 = ea0x - ea1x;\n",
207+
"const auto t1 = eb0x - eb1x;\n",
208+
"const auto t2 = eb0y - eb1y;\n",
209+
"const auto t3 = eb0z - eb1z;\n",
210+
"const auto t4 = 2*scale;\n",
211+
"const auto t5 = t4*(std::pow(t1, 2) + std::pow(t2, 2) + std::pow(t3, 2));\n",
212+
"const auto t6 = t0*t5;\n",
213+
"const auto t7 = ea0y - ea1y;\n",
214+
"const auto t8 = t5*t7;\n",
215+
"const auto t9 = ea0z - ea1z;\n",
216+
"const auto t10 = t5*t9;\n",
217+
"const auto t11 = t4*(std::pow(t0, 2) + std::pow(t7, 2) + std::pow(t9, 2));\n",
218+
"const auto t12 = t1*t11;\n",
219+
"const auto t13 = t11*t2;\n",
220+
"const auto t14 = t11*t3;\n",
221+
"grad[0] = t6;\n",
222+
"grad[1] = t8;\n",
223+
"grad[2] = t10;\n",
224+
"grad[3] = -t6;\n",
225+
"grad[4] = -t8;\n",
226+
"grad[5] = -t10;\n",
227+
"grad[6] = t12;\n",
228+
"grad[7] = t13;\n",
229+
"grad[8] = t14;\n",
230+
"grad[9] = -t12;\n",
231+
"grad[10] = -t13;\n",
232+
"grad[11] = -t14;\n",
233+
"}\n",
234+
"\n"
235+
]
236+
}
237+
],
238+
"source": [
239+
"x = sympy.Matrix([ea0, ea1, eb0, eb1])\n",
240+
"\n",
241+
"g = CXXGradientGenerator(eps_x, x, 'edge_edge_mollifier_threshold_gradient')\n",
242+
"\n",
243+
"print(g())"
244+
]
245+
}
246+
],
247+
"metadata": {
248+
"kernelspec": {
249+
"display_name": "Python 3",
250+
"language": "python",
251+
"name": "python3"
252+
},
253+
"language_info": {
254+
"codemirror_mode": {
255+
"name": "ipython",
256+
"version": 3
257+
},
258+
"file_extension": ".py",
259+
"mimetype": "text/x-python",
260+
"name": "python",
261+
"nbconvert_exporter": "python",
262+
"pygments_lexer": "ipython3",
263+
"version": "3.11.6"
264+
}
265+
},
266+
"nbformat": 4,
267+
"nbformat_minor": 2
268+
}

‎notebooks/generate_cpp_code.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
from utils import jacobian
88

99

10-
def generate_code(expr, out_var_name="J"):
10+
def generate_code(expr, out_var_name=None):
1111
CSE_results = sympy.cse(
1212
expr, sympy.numbered_symbols("t"), optimizations='basic')
1313
lines = []

‎python/src/collisions/collision_constraints.cpp‎

Lines changed: 64 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ void define_collision_constraints(py::module_& m)
2727
)ipc_Qu8mg5v7",
2828
py::arg("mesh"), py::arg("vertices"), py::arg("dhat"),
2929
py::arg("dmin") = 0,
30-
py::arg("broad_phase_method") = BroadPhaseMethod::HASH_GRID)
30+
py::arg("broad_phase_method") = DEFAULT_BROAD_PHASE_METHOD)
3131
.def(
3232
"build",
3333
py::overload_cast<
@@ -147,6 +147,69 @@ void define_collision_constraints(py::module_& m)
147147
A reference to the constraint.
148148
)ipc_Qu8mg5v7",
149149
py::arg("idx"))
150+
.def(
151+
"is_vertex_vertex", &CollisionConstraints::is_vertex_vertex,
152+
R"ipc_Qu8mg5v7(
153+
Get if the constraint at idx is a vertex-vertex constraint.
154+
155+
Parameters:
156+
idx: The index of the constraint.
157+
158+
Returns:
159+
If the constraint at idx is a vertex-vertex constraint.
160+
)ipc_Qu8mg5v7",
161+
py::arg("idx"))
162+
.def(
163+
"is_edge_vertex", &CollisionConstraints::is_edge_vertex,
164+
R"ipc_Qu8mg5v7(
165+
Get if the constraint at idx is an edge-vertex constraint.
166+
167+
Parameters:
168+
idx: The index of the constraint.
169+
170+
Returns:
171+
If the constraint at idx is an edge-vertex constraint.
172+
)ipc_Qu8mg5v7",
173+
py::arg("idx"))
174+
.def(
175+
"is_edge_edge", &CollisionConstraints::is_edge_edge,
176+
R"ipc_Qu8mg5v7(
177+
Get if the constraint at idx is an edge-edge constraint.
178+
179+
Parameters:
180+
idx: The index of the constraint.
181+
182+
Returns:
183+
If the constraint at idx is an edge-edge constraint.
184+
)ipc_Qu8mg5v7",
185+
py::arg("idx"))
186+
.def(
187+
"is_face_vertex", &CollisionConstraints::is_face_vertex,
188+
R"ipc_Qu8mg5v7(
189+
Get if the constraint at idx is an face-vertex constraint.
190+
191+
Parameters:
192+
idx: The index of the constraint.
193+
194+
Returns:
195+
If the constraint at idx is an face-vertex constraint.
196+
)ipc_Qu8mg5v7",
197+
py::arg("idx"))
198+
.def(
199+
"is_plane_vertex", &CollisionConstraints::is_plane_vertex,
200+
R"ipc_Qu8mg5v7(
201+
Get if the constraint at idx is an plane-vertex constraint.
202+
203+
Parameters:
204+
idx: The index of the constraint.
205+
206+
Returns:
207+
If the constraint at idx is an plane-vertex constraint.
208+
)ipc_Qu8mg5v7",
209+
py::arg("idx"))
210+
.def(
211+
"to_string", &CollisionConstraints::to_string, "", py::arg("mesh"),
212+
py::arg("vertices"))
150213
.def_property(
151214
"use_convergent_formulation",
152215
&CollisionConstraints::use_convergent_formulation,

0 commit comments

Comments
 (0)