Skip to content

Commit 112712c

Browse files
zferguszfergus
authored
Refactor to use Eigen::ConstRef for vertices (#200)
- Remove Vector<T, dim, max_dim> and replace it with either Eigen::Vector<T, dim> or VectorMax<T, max_dim> - Updated constructors in Edge, Edge2, Edge3, Face, Point2, and Point3 classes to accept Eigen::ConstRef<Eigen::MatrixXd> instead of Eigen::MatrixXd for vertices. - Renamed GradType to GradientType in edge.hpp, edge3.hpp, and point3.hpp for consistency. - Adjusted gradient and hessian function signatures in Point2 and Point3 to use Eigen::Vector and VectorMax types. - Modified utility functions in smooth_collisions_builder to accept Eigen::ConstRef for vertices.
1 parent e62964b commit 112712c

36 files changed

+425
-391
lines changed

‎src/ipc/collisions/tangential/tangential_collisions.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ void TangentialCollisions::build(
116116

117117
void TangentialCollisions::build(
118118
const CollisionMesh& mesh,
119-
const Eigen::MatrixXd& vertices,
119+
Eigen::ConstRef<Eigen::MatrixXd> vertices,
120120
const SmoothCollisions& collisions,
121121
const SmoothContactParameters& params,
122122
const double normal_stiffness,

‎src/ipc/collisions/tangential/tangential_collisions.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -97,7 +97,7 @@ class TangentialCollisions {
9797
/// @param blend_mu Function to blend vertex-based coefficients of friction. Defaults to average.
9898
void build(
9999
const CollisionMesh& mesh,
100-
const Eigen::MatrixXd& vertices,
100+
Eigen::ConstRef<Eigen::MatrixXd> vertices,
101101
const SmoothCollisions& collisions,
102102
const SmoothContactParameters& params,
103103
const double normal_stiffness,

‎src/ipc/math/math.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ double opposite_direction_penalty(
4949
return Math<double>::smooth_heaviside(d.dot(t) / t.norm(), alpha, beta);
5050
}
5151

52-
GradType<6> opposite_direction_penalty_grad(
52+
GradientType<6> opposite_direction_penalty_grad(
5353
Eigen::ConstRef<Eigen::Vector3d> t,
5454
Eigen::ConstRef<Eigen::Vector3d> d,
5555
const double alpha,
@@ -109,7 +109,7 @@ double negative_orientation_penalty(
109109
return opposite_direction_penalty(n, d, alpha, beta);
110110
}
111111

112-
GradType<9> negative_orientation_penalty_grad(
112+
GradientType<9> negative_orientation_penalty_grad(
113113
Eigen::ConstRef<Eigen::Vector3d> t1,
114114
Eigen::ConstRef<Eigen::Vector3d> t2,
115115
Eigen::ConstRef<Eigen::Vector3d> d,

‎src/ipc/math/math.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ double opposite_direction_penalty(
119119
const double alpha,
120120
const double beta);
121121

122-
GradType<6> opposite_direction_penalty_grad(
122+
GradientType<6> opposite_direction_penalty_grad(
123123
Eigen::ConstRef<Eigen::Vector3d> t,
124124
Eigen::ConstRef<Eigen::Vector3d> d,
125125
const double alpha,
@@ -139,7 +139,7 @@ double negative_orientation_penalty(
139139
const double alpha,
140140
const double beta);
141141

142-
GradType<9> negative_orientation_penalty_grad(
142+
GradientType<9> negative_orientation_penalty_grad(
143143
Eigen::ConstRef<Eigen::Vector3d> t1,
144144
Eigen::ConstRef<Eigen::Vector3d> t2,
145145
Eigen::ConstRef<Eigen::Vector3d> d,

‎src/ipc/potentials/potential.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ template <class TCollisions> class Potential {
1212
using TCollision = typename TCollisions::value_type;
1313
/// @brief Maximum degrees of freedom per collision
1414
static constexpr int STENCIL_NDOF = 3 * TCollision::STENCIL_SIZE;
15-
using VectorMaxNd = Vector<double, Eigen::Dynamic, STENCIL_NDOF>;
15+
using VectorMaxNd = VectorMax<double, STENCIL_NDOF>;
1616
using MatrixMaxNd = MatrixMax<double, STENCIL_NDOF, STENCIL_NDOF>;
1717

1818
public:

‎src/ipc/potentials/tangential_potential.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -683,7 +683,7 @@ TangentialPotential::VectorMaxNd TangentialPotential::smooth_contact_force(
683683
assert(rest_positions.size() == velocities.size());
684684

685685
// const VectorMaxNd x = dof(rest_positions, edges, faces);
686-
// const Vector<double, -1, STENCIL_NDOF> u =
686+
// const VectorMax<double, STENCIL_NDOF> u =
687687
// dof(lagged_displacements, edges, faces);
688688
// const VectorMaxNd v = dof(velocities, edges, faces);
689689
const VectorMaxNd lagged_positions =

‎src/ipc/smooth_contact/collisions/smooth_collision.cpp‎

Lines changed: 49 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -42,9 +42,9 @@ Eigen::VectorXd SmoothCollision::dof(Eigen::ConstRef<Eigen::MatrixXd> X) const
4242

4343
template <typename PrimitiveA, typename PrimitiveB>
4444
auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::get_core_indices() const
45-
-> Vector<int, N_CORE_DOFS>
45+
-> Eigen::Vector<int, N_CORE_DOFS>
4646
{
47-
Vector<int, N_CORE_DOFS> core_indices;
47+
Eigen::Vector<int, N_CORE_DOFS> core_indices;
4848
core_indices << Eigen::VectorXi::LinSpaced(
4949
N_CORE_DOFS_A, 0, N_CORE_DOFS_A - 1),
5050
Eigen::VectorXi::LinSpaced(
@@ -61,7 +61,7 @@ SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::SmoothCollisionTemplate(
6161
const CollisionMesh& mesh,
6262
const SmoothContactParameters& params,
6363
const double _dhat,
64-
const Eigen::MatrixXd& V)
64+
Eigen::ConstRef<Eigen::MatrixXd> V)
6565
: SmoothCollision(_primitive0, _primitive1, _dhat, mesh)
6666
{
6767
VectorMax3d d =
@@ -104,16 +104,16 @@ SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::SmoothCollisionTemplate(
104104

105105
template <typename PrimitiveA, typename PrimitiveB>
106106
double SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::operator()(
107-
Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
107+
Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
108108
const SmoothContactParameters& params) const
109109
{
110-
Vector<double, N_CORE_POINTS * DIM> x;
110+
Eigen::Vector<double, N_CORE_POINTS * DIM> x;
111111
x << positions.head(PrimitiveA::N_CORE_POINTS * DIM),
112112
positions.segment(
113113
primitive_a->n_dofs(), PrimitiveB::N_CORE_POINTS * DIM);
114114

115115
// grad of "d" wrt. points
116-
const Vector<double, DIM> closest_direction =
116+
const Eigen::Vector<double, DIM> closest_direction =
117117
PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, double>::
118118
compute_closest_direction(x, DTYPE::AUTO);
119119
const double dist = closest_direction.norm();
@@ -143,19 +143,19 @@ double SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::operator()(
143143

144144
template <typename PrimitiveA, typename PrimitiveB>
145145
auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
146-
Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
146+
Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
147147
const SmoothContactParameters& params) const
148-
-> Vector<double, -1, ELEMENT_SIZE>
148+
-> VectorMax<double, ELEMENT_SIZE>
149149
{
150150
const auto core_indices = get_core_indices();
151151

152-
Vector<double, N_CORE_DOFS> x;
152+
Eigen::Vector<double, N_CORE_DOFS> x;
153153
x = positions(core_indices);
154154

155155
const auto dtype =
156156
PrimitiveDistance<PrimitiveA, PrimitiveB>::compute_distance_type(x);
157157

158-
Vector<double, DIM> closest_direction;
158+
Eigen::Vector<double, DIM> closest_direction;
159159
Eigen::Matrix<double, DIM, N_CORE_DOFS> closest_direction_grad;
160160
std::tie(closest_direction, closest_direction_grad) = PrimitiveDistance<
161161
PrimitiveA, PrimitiveB>::compute_closest_direction_gradient(x, dtype);
@@ -171,42 +171,45 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
171171

172172
// gradient of barrier potential
173173
double barrier = 0;
174-
Vector<double, N_CORE_DOFS> gBarrier = Vector<double, N_CORE_DOFS>::Zero();
174+
Eigen::Vector<double, N_CORE_DOFS> gBarrier =
175+
Eigen::Vector<double, N_CORE_DOFS>::Zero();
175176
{
176177
barrier = Math<double>::inv_barrier(dist / dhat(), params.r);
177178

178-
const Vector<double, DIM> closest_direction_normalized =
179+
const Eigen::Vector<double, DIM> closest_direction_normalized =
179180
closest_direction / dist;
180181
const double barrier_1st_deriv =
181182
Math<double>::inv_barrier_grad(dist / dhat(), params.r) / dhat();
182-
const Vector<double, DIM> gBarrier_wrt_d =
183+
const Eigen::Vector<double, DIM> gBarrier_wrt_d =
183184
barrier_1st_deriv * closest_direction_normalized;
184185
gBarrier = closest_direction_grad.transpose() * gBarrier_wrt_d;
185186
}
186187

187188
// gradient of mollifier
188189
{
189190
double mollifier = 0;
190-
Vector<double, N_CORE_DOFS> gMollifier =
191-
Vector<double, N_CORE_DOFS>::Zero();
191+
Eigen::Vector<double, N_CORE_DOFS> gMollifier =
192+
Eigen::Vector<double, N_CORE_DOFS>::Zero();
192193
#ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
193194
ScalarBase::setVariableCount(N_CORE_DOFS);
194195
using T = ADGrad<N_CORE_DOFS>;
195-
Vector<T, N_CORE_DOFS> xAD = slice_positions<T, N_CORE_DOFS, 1>(x);
196-
Vector<T, DIM> closest_direction_autodiff = PrimitiveDistanceTemplate<
197-
PrimitiveA, PrimitiveB, T>::compute_closest_direction(xAD, dtype);
196+
Eigen::Vector<T, N_CORE_DOFS> xAD =
197+
slice_positions<T, N_CORE_DOFS, 1>(x);
198+
Eigen::Vector<T, DIM> closest_direction_autodiff =
199+
PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::
200+
compute_closest_direction(xAD, dtype);
198201
const auto dist_sqr_AD = closest_direction_autodiff.squaredNorm();
199202
auto mollifier_autodiff =
200203
PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::mollifier(
201204
xAD, dist_sqr_AD);
202205
mollifier = mollifier_autodiff.val;
203206
gMollifier = mollifier_autodiff.grad;
204207
#else
205-
Vector<double, N_CORE_DOFS + 1> mollifier_grad;
208+
Eigen::Vector<double, N_CORE_DOFS + 1> mollifier_grad;
206209
std::tie(mollifier, mollifier_grad) = PrimitiveDistance<
207210
PrimitiveA, PrimitiveB>::compute_mollifier_gradient(x, dist * dist);
208211

209-
const Vector<double, N_CORE_DOFS> dist_sqr_grad =
212+
const Eigen::Vector<double, N_CORE_DOFS> dist_sqr_grad =
210213
2 * closest_direction_grad.transpose() * closest_direction;
211214
mollifier_grad.head(N_CORE_DOFS) +=
212215
mollifier_grad(N_CORE_DOFS) * dist_sqr_grad;
@@ -219,11 +222,11 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
219222

220223
// grad of tangent/normal terms
221224
double orient = 0;
222-
Vector<double, -1, ELEMENT_SIZE> gOrient;
225+
VectorMax<double, ELEMENT_SIZE> gOrient;
223226
{
224-
Vector<double, -1, ELEMENT_SIZE>
225-
gA = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs()),
226-
gB = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs());
227+
VectorMax<double, ELEMENT_SIZE>
228+
gA = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs()),
229+
gB = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs());
227230
{
228231
gA(core_indices) =
229232
closest_direction_grad.transpose() * gA_reduced.head(DIM);
@@ -254,19 +257,19 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
254257

255258
template <typename PrimitiveA, typename PrimitiveB>
256259
auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
257-
Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
260+
Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
258261
const SmoothContactParameters& params) const
259262
-> MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE>
260263
{
261264
const auto core_indices = get_core_indices();
262265

263-
Vector<double, N_CORE_DOFS> x;
266+
Eigen::Vector<double, N_CORE_DOFS> x;
264267
x = positions(core_indices);
265268

266269
const auto dtype =
267270
PrimitiveDistance<PrimitiveA, PrimitiveB>::compute_distance_type(x);
268271

269-
Vector<double, DIM> closest_direction;
272+
Eigen::Vector<double, DIM> closest_direction;
270273
Eigen::Matrix<double, DIM, N_CORE_DOFS> closest_direction_grad;
271274
std::array<Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>, DIM>
272275
closest_direction_hess;
@@ -289,17 +292,18 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
289292

290293
// hessian of barrier potential
291294
double barrier = 0;
292-
Vector<double, N_CORE_DOFS> gBarrier = Vector<double, N_CORE_DOFS>::Zero();
295+
Eigen::Vector<double, N_CORE_DOFS> gBarrier =
296+
Eigen::Vector<double, N_CORE_DOFS>::Zero();
293297
Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS> hBarrier =
294298
Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>::Zero();
295299
{
296300
barrier = Math<double>::inv_barrier(dist / dhat(), params.r);
297301

298-
const Vector<double, DIM> closest_direction_normalized =
302+
const Eigen::Vector<double, DIM> closest_direction_normalized =
299303
closest_direction / dist;
300304
const double barrier_1st_deriv =
301305
Math<double>::inv_barrier_grad(dist / dhat(), params.r) / dhat();
302-
const Vector<double, DIM> gBarrier_wrt_d =
306+
const Eigen::Vector<double, DIM> gBarrier_wrt_d =
303307
barrier_1st_deriv * closest_direction_normalized;
304308
gBarrier = closest_direction_grad.transpose() * gBarrier_wrt_d;
305309

@@ -322,16 +326,18 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
322326
// hessian of mollifier
323327
{
324328
double mollifier = 0;
325-
Vector<double, N_CORE_DOFS> gMollifier =
326-
Vector<double, N_CORE_DOFS>::Zero();
329+
Eigen::Vector<double, N_CORE_DOFS> gMollifier =
330+
Eigen::Vector<double, N_CORE_DOFS>::Zero();
327331
Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS> hMollifier =
328332
Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>::Zero();
329333
#ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
330334
ScalarBase::setVariableCount(N_CORE_DOFS);
331335
using T = ADHessian<N_CORE_DOFS>;
332-
Vector<T, N_CORE_DOFS> xAD = slice_positions<T, N_CORE_DOFS, 1>(x);
333-
Vector<T, DIM> closest_direction_autodiff = PrimitiveDistanceTemplate<
334-
PrimitiveA, PrimitiveB, T>::compute_closest_direction(xAD, dtype);
336+
Eigen::Vector<T, N_CORE_DOFS> xAD =
337+
slice_positions<T, N_CORE_DOFS, 1>(x);
338+
Eigen::Vector<T, DIM> closest_direction_autodiff =
339+
PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::
340+
compute_closest_direction(xAD, dtype);
335341
const auto dist_sqr_AD = closest_direction_autodiff.squaredNorm();
336342
auto mollifier_autodiff =
337343
PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::mollifier(
@@ -341,12 +347,12 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
341347
gMollifier = mollifier_autodiff.grad;
342348
hMollifier = mollifier_autodiff.Hess;
343349
#else
344-
Vector<double, N_CORE_DOFS + 1> mollifier_grad;
350+
Eigen::Vector<double, N_CORE_DOFS + 1> mollifier_grad;
345351
Eigen::Matrix<double, N_CORE_DOFS + 1, N_CORE_DOFS + 1> mollifier_hess;
346352
std::tie(mollifier, mollifier_grad, mollifier_hess) = PrimitiveDistance<
347353
PrimitiveA, PrimitiveB>::compute_mollifier_hessian(x, dist * dist);
348354

349-
const Vector<double, N_CORE_DOFS> dist_sqr_grad =
355+
const Eigen::Vector<double, N_CORE_DOFS> dist_sqr_grad =
350356
2 * closest_direction_grad.transpose() * closest_direction;
351357
mollifier_grad.head(N_CORE_DOFS) +=
352358
mollifier_grad(N_CORE_DOFS) * dist_sqr_grad;
@@ -378,12 +384,12 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
378384

379385
// grad of tangent/normal terms
380386
double orient = 0;
381-
Vector<double, -1, ELEMENT_SIZE> gOrient;
387+
VectorMax<double, ELEMENT_SIZE> gOrient;
382388
MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE> hOrient;
383389
{
384-
Vector<double, -1, ELEMENT_SIZE>
385-
gA = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs()),
386-
gB = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs());
390+
VectorMax<double, ELEMENT_SIZE>
391+
gA = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs()),
392+
gB = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs());
387393
MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE>
388394
hA = MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE>::Zero(
389395
n_dofs(), n_dofs()),
@@ -467,9 +473,9 @@ template <typename PrimitiveA, typename PrimitiveB>
467473
double SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::compute_distance(
468474
Eigen::ConstRef<Eigen::MatrixXd> vertices) const
469475
{
470-
Vector<double, -1, ELEMENT_SIZE> positions = dof(vertices);
476+
VectorMax<double, ELEMENT_SIZE> positions = dof(vertices);
471477

472-
Vector<double, N_CORE_POINTS * DIM> x;
478+
Eigen::Vector<double, N_CORE_POINTS * DIM> x;
473479
x << positions.head(PrimitiveA::N_CORE_POINTS * DIM),
474480
positions.segment(
475481
primitive_a->n_dofs(), PrimitiveB::N_CORE_POINTS * DIM);

0 commit comments

Comments
 (0)