Skip to content
This repository was archived by the owner on Oct 22, 2025. It is now read-only.

Commit ed0a19b

Browse files
HuangzizhouHuangzizhou
committed
fix python binding
1 parent 7185ecd commit ed0a19b

18 files changed

+867
-73
lines changed

‎python/new-formulation-high-order.ipynb‎

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import numpy as np\n",
10+
"import plotly.figure_factory as ff\n",
11+
"import plotly.graph_objects as go\n",
12+
"from plotly.subplots import make_subplots\n",
13+
"import meshio, numpy, copy\n",
14+
"import meshplot as mp\n",
15+
"# import sys\n",
16+
"# import pathlib\n",
17+
"# sys.path.append(str(pathlib.Path(\"\").resolve().parent / \"build\" / \"release\" / \"python\")) # noqa\n",
18+
"\n",
19+
"from ipctk import *"
20+
]
21+
},
22+
{
23+
"cell_type": "code",
24+
"execution_count": null,
25+
"metadata": {},
26+
"outputs": [],
27+
"source": [
28+
"V = np.array([[0.5,0],[1.0,0],[1.5,0],[0.9,0.3],[1,0.1],[1.1,0.3],[1,-0.5]])\n",
29+
"E = np.array([[1,0],[2,1],[3,4],[4,5],[5,3],[0,6],[6,2]], dtype=int)\n",
30+
"F = np.array([])\n",
31+
"\n",
32+
"p = mp.plot(V, F, shading={\"width\": 200, \"height\": 200})\n",
33+
"p.add_points(V, shading={\"point_color\": \"green\", \"point_size\": 0.1})\n",
34+
"p.add_lines(V[E[:,0],:],V[E[:,1],:], shading={\"line_width\": 2})\n",
35+
"\n",
36+
"dhat=0.12\n",
37+
"param = ParameterType(dhat, 1, 0, 1, 0, 1, 200)\n",
38+
"mesh = CollisionMesh(V, E, F)"
39+
]
40+
},
41+
{
42+
"cell_type": "code",
43+
"execution_count": null,
44+
"metadata": {},
45+
"outputs": [],
46+
"source": [
47+
"# high order\n",
48+
"\n",
49+
"def potential_high_order(points):\n",
50+
" cs1 = SmoothCollisions2()\n",
51+
" cs1.use_high_order_quadrature = True\n",
52+
" cs1.build(mesh, points, param)\n",
53+
"\n",
54+
" B = SmoothPotential(param)\n",
55+
" # g = B.gradient(cs1, mesh, points)\n",
56+
"\n",
57+
" return B(cs1, mesh, points)"
58+
]
59+
},
60+
{
61+
"cell_type": "code",
62+
"execution_count": null,
63+
"metadata": {},
64+
"outputs": [],
65+
"source": [
66+
"cs1 = SmoothCollisions2()\n",
67+
"cs1.use_high_order_quadrature = True\n",
68+
"cs1.build(mesh, V, param)\n",
69+
"print(cs1.to_string(mesh, V, param))\n",
70+
"print(mesh.edges)"
71+
]
72+
},
73+
{
74+
"cell_type": "code",
75+
"execution_count": null,
76+
"metadata": {},
77+
"outputs": [],
78+
"source": [
79+
"# single point\n",
80+
"\n",
81+
"def potential_single_point(points):\n",
82+
" cs2 = SmoothCollisions2()\n",
83+
" cs2.use_high_order_quadrature = False\n",
84+
" cs2.build(mesh, points, param)\n",
85+
" \n",
86+
" B = SmoothPotential(param)\n",
87+
" # g = B.gradient(cs1, mesh, points)\n",
88+
"\n",
89+
" return B(cs2, mesh, points)"
90+
]
91+
},
92+
{
93+
"cell_type": "code",
94+
"execution_count": null,
95+
"metadata": {},
96+
"outputs": [],
97+
"source": [
98+
"# IPC\n",
99+
"\n",
100+
"def potential_IPC(points):\n",
101+
" cs3 = Collisions()\n",
102+
" cs3.build(mesh, points, param.dhat)\n",
103+
"\n",
104+
" B = BarrierPotential(param.dhat)\n",
105+
" # g = B.gradient(cs3, mesh, points)\n",
106+
"\n",
107+
" return B(cs3, mesh, points)"
108+
]
109+
},
110+
{
111+
"cell_type": "code",
112+
"execution_count": null,
113+
"metadata": {},
114+
"outputs": [],
115+
"source": [
116+
"print(potential_IPC(V), potential_single_point(V), potential_high_order(V))"
117+
]
118+
},
119+
{
120+
"cell_type": "code",
121+
"execution_count": null,
122+
"metadata": {},
123+
"outputs": [],
124+
"source": [
125+
"displacements = np.linspace(-0.1, 0.1, 200)\n",
126+
"pA = []\n",
127+
"pB = []\n",
128+
"pC = []\n",
129+
"\n",
130+
"for disp in displacements:\n",
131+
" V_disp = copy.deepcopy(V)\n",
132+
" V_disp[[3,4,5], 0] += disp\n",
133+
" \n",
134+
" pA.append(potential_IPC(V_disp))\n",
135+
" pB.append(potential_single_point(V_disp))\n",
136+
" pC.append(potential_high_order(V_disp))\n"
137+
]
138+
},
139+
{
140+
"cell_type": "code",
141+
"execution_count": null,
142+
"metadata": {},
143+
"outputs": [],
144+
"source": [
145+
"fig = go.Figure(data=[\n",
146+
" go.Scatter(x=displacements, y=pA, name=\"IPC\"),\n",
147+
" go.Scatter(x=displacements, y=pB, name=\"Single Point\"),\n",
148+
" go.Scatter(x=displacements, y=pC, name=\"Exact Quadrature\"),\n",
149+
"], layout=go.Layout(width=800, height=400))\n",
150+
"\n",
151+
"# fig.update_layout(\n",
152+
"# yaxis_type=\"log\"\n",
153+
"# )\n",
154+
"\n",
155+
"fig.show()"
156+
]
157+
},
158+
{
159+
"cell_type": "code",
160+
"execution_count": null,
161+
"metadata": {},
162+
"outputs": [],
163+
"source": []
164+
}
165+
],
166+
"metadata": {
167+
"kernelspec": {
168+
"display_name": "base",
169+
"language": "python",
170+
"name": "python3"
171+
},
172+
"language_info": {
173+
"codemirror_mode": {
174+
"name": "ipython",
175+
"version": 3
176+
},
177+
"file_extension": ".py",
178+
"mimetype": "text/x-python",
179+
"name": "python",
180+
"nbconvert_exporter": "python",
181+
"pygments_lexer": "ipython3",
182+
"version": "3.11.8"
183+
}
184+
},
185+
"nbformat": 4,
186+
"nbformat_minor": 2
187+
}

‎python/src/potentials/barrier_potential.cpp‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -45,15 +45,16 @@ void define_smooth_potential(py::module_& m)
4545
.def(
4646
py::init<
4747
const double&, const double&, const double&, const double&,
48-
const double&, const int&>(),
48+
const double&, const int&, const int&>(),
4949
R"ipc_Qu8mg5v7(
5050
Construct parameter set for smooth contact.
5151
5252
Parameters:
5353
dhat, alpha_t, beta_t, alpha_n, beta_n, r
5454
)ipc_Qu8mg5v7",
5555
py::arg("dhat"), py::arg("alpha_t"), py::arg("beta_t"),
56-
py::arg("alpha_n"), py::arg("beta_n"), py::arg("r"))
56+
py::arg("alpha_n"), py::arg("beta_n"), py::arg("r"),
57+
py::arg("n_quadrature_samples"))
5758
.def(
5859
py::init<const double&, const double&, const double&, const int&>(),
5960
R"ipc_Qu8mg5v7(
@@ -69,7 +70,9 @@ void define_smooth_potential(py::module_& m)
6970
.def_readonly("beta_t", &ParameterType::beta_t)
7071
.def_readonly("alpha_n", &ParameterType::alpha_n)
7172
.def_readonly("beta_n", &ParameterType::beta_n)
72-
.def_readonly("r", &ParameterType::r);
73+
.def_readonly("r", &ParameterType::r)
74+
.def_readonly(
75+
"n_quadrature_samples", &ParameterType::n_quadrature_samples);
7376

7477
py::class_<SmoothContactPotential>(m, "SmoothPotential")
7578
.def(

‎src/ipc/smooth_contact/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,3 +18,4 @@ target_sources(ipc_toolkit PRIVATE ${SOURCES})
1818
add_subdirectory(collisions)
1919
add_subdirectory(distance)
2020
add_subdirectory(primitives)
21+
add_subdirectory(high_order)

‎src/ipc/smooth_contact/common.hpp‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -29,28 +29,28 @@ struct ParameterType {
2929
const double _dhat,
3030
const double _alpha_t,
3131
const double _beta_t,
32-
const int _r)
33-
: ParameterType(_dhat, _alpha_t, _beta_t, 0, 0.1, _r)
34-
{
35-
}
32+
const int _r):
33+
ParameterType(_dhat, _alpha_t, _beta_t, 0, 0.1, _r, 1000) {}
3634
ParameterType(
3735
const double _dhat,
3836
const double _alpha_t,
3937
const double _beta_t,
4038
const double _alpha_n,
4139
const double _beta_n,
42-
const int _r)
40+
const int _r,
41+
const int _n_quadrature_samples = 1)
4342
: dhat(_dhat)
4443
, alpha_t(_alpha_t)
4544
, beta_t(_beta_t)
4645
, alpha_n(_alpha_n)
4746
, beta_n(_beta_n)
4847
, r(_r)
48+
, n_quadrature_samples(_n_quadrature_samples)
4949
{
50-
if (!(r > 0) || !(dhat > 0) || !(abs(alpha_t) <= 1)
51-
|| !(abs(alpha_n) <= 1) || !(abs(beta_t) <= 1)
52-
|| !(abs(beta_n) <= 1) || !(beta_t + alpha_t > 1e-6)
53-
|| !(beta_n + alpha_n > 1e-6))
50+
if (!(r > 0) || !(dhat > 0) ||
51+
!(abs(alpha_t) <= 1) || !(abs(alpha_n) <= 1) ||
52+
!(abs(beta_t) <= 1) || !(abs(beta_n) <= 1) ||
53+
!(beta_t + alpha_t > 1e-6) || !(beta_n + alpha_n > 1e-6) || (n_quadrature_samples < 1))
5454
logger().error(
5555
"Wrong parameters for smooth contact! dhat {} alpha_t {} beta_t {} alpha_n {} beta_n {} r {}",
5656
dhat, alpha_t, beta_t, alpha_n, beta_n, r);
@@ -67,6 +67,7 @@ struct ParameterType {
6767
double alpha_t = 1, beta_t = 0;
6868
double alpha_n = 0.1, beta_n = 0;
6969
int r = 2;
70+
int n_quadrature_samples = 1000;
7071

7172
private:
7273
double adaptive_dhat_ratio = 0.5;

‎src/ipc/smooth_contact/high_order/CMakeLists.txt‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
set(SOURCES
2+
edge_edge_2d.hpp
3+
edge_edge_2d.cpp
4+
edge_vert_2d.hpp
5+
edge_vert_2d.cpp
6+
vert_vert_2d.hpp
7+
vert_vert_2d.cpp
8+
)
9+
10+
source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES})
11+
target_sources(ipc_toolkit PRIVATE ${SOURCES})

‎src/ipc/smooth_contact/high_order/edge_edge_2d.cpp‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
#include "edge_edge_2d.hpp"
2+
#include <ipc/utils/math.hpp>
3+
#include <ipc/distance/point_edge.hpp>
4+
#include <ipc/distance/edge_edge.hpp>
5+
6+
namespace ipc {
7+
EdgeEdge2D::EdgeEdge2D(
8+
const long& id0,
9+
const long& id1,
10+
const CollisionMesh& mesh,
11+
const double& dhat,
12+
const Eigen::MatrixXd& V)
13+
: SmoothCollision(id0, id1, dhat, mesh)
14+
{
15+
vertex_ids_[0] = mesh.edges()(id0, 0);
16+
vertex_ids_[1] = mesh.edges()(id0, 1);
17+
vertex_ids_[2] = mesh.edges()(id1, 0);
18+
vertex_ids_[3] = mesh.edges()(id1, 1);
19+
20+
Eigen::Matrix<double, 4, 3> P;
21+
P.setZero();
22+
P.leftCols<2>() << V.row(vertex_ids_[0]), V.row(vertex_ids_[1]),
23+
V.row(vertex_ids_[2]), V.row(vertex_ids_[3]);
24+
25+
const double dist =
26+
edge_edge_distance(P.row(0), P.row(1), P.row(2), P.row(3));
27+
28+
is_active_ = dist < dhat * dhat;
29+
}
30+
31+
double EdgeEdge2D::operator()(
32+
Eigen::ConstRef<Vector<double, -1, SmoothCollision::element_size>> positions,
33+
const ParameterType& params) const
34+
{
35+
assert(positions.size() == n_dofs());
36+
const Eigen::Vector2d x0 = positions.segment<2>(0);
37+
const Eigen::Vector2d x1 = positions.segment<2>(4);
38+
const Eigen::Vector2d e0 = positions.segment<2>(2) - x0;
39+
const Eigen::Vector2d e1 = positions.segment<2>(6) - x1;
40+
const double l0 = e0.norm();
41+
const double l1 = e1.norm();
42+
43+
Eigen::Vector2d n0, n1;
44+
n0 << e0(1), -e0(0);
45+
n1 << e1(1), -e1(0);
46+
n0.normalize();
47+
n1.normalize();
48+
49+
auto integrand = [&](double theta0, double theta1) {
50+
const Eigen::Vector2d p0 = x0 + theta0 * e0;
51+
const Eigen::Vector2d p1 = x1 + theta1 * e1;
52+
const Eigen::Vector2d d = p1 - p0;
53+
const double dist = d.norm();
54+
55+
return Math<double>::smooth_heaviside(
56+
d.dot(n0) / dist - 1., params.alpha_n, params.beta_n)
57+
* Math<double>::smooth_heaviside(
58+
-d.dot(n1) / dist - 1., params.alpha_n, params.beta_n)
59+
* Math<double>::inv_barrier(dist / params.dhat, params.r);
60+
};
61+
62+
// uniform quadrature
63+
const int n_samples = params.n_quadrature_samples;
64+
double result = 0.;
65+
for (int i = 0; i <= n_samples; i++) {
66+
for (int j = 0; j <= n_samples; j++) {
67+
double f = integrand(static_cast<double>(i) / n_samples, static_cast<double>(j) / n_samples);
68+
result += f;
69+
}
70+
}
71+
result *= (l0 * l1) / (n_samples * n_samples);
72+
73+
return result;
74+
}
75+
76+
Vector<double, -1, SmoothCollision::element_size> EdgeEdge2D::gradient(
77+
Eigen::ConstRef<Vector<double, -1, SmoothCollision::element_size>> positions,
78+
const ParameterType& params) const
79+
{
80+
return Eigen::VectorXd::Zero(positions.size());
81+
}
82+
83+
MatrixMax<double, SmoothCollision::element_size, SmoothCollision::element_size>
84+
EdgeEdge2D::hessian(
85+
Eigen::ConstRef<Vector<double, -1, SmoothCollision::element_size>> positions,
86+
const ParameterType& params) const
87+
{
88+
return Eigen::MatrixXd::Zero(positions.size(), positions.size());
89+
}
90+
} // namespace ipc

0 commit comments

Comments
 (0)