Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PeriodicContactForm.cpp
Go to the documentation of this file.
2
6#include <polyfem/State.hpp>
7
8#include <iostream>
9
10namespace polyfem::solver
11{
12 PeriodicContactForm::PeriodicContactForm(const ipc::CollisionMesh &periodic_collision_mesh,
13 const Eigen::VectorXi &tiled_to_single,
14 const double dhat,
15 const double avg_mass,
16 const bool use_area_weighting,
17 const bool use_improved_max_operator,
18 const bool use_physical_barrier,
19 const bool use_adaptive_barrier_stiffness,
20 const bool is_time_dependent,
21 const bool enable_shape_derivatives,
22 const ipc::BroadPhaseMethod broad_phase_method,
23 const double ccd_tolerance,
24 const int ccd_max_iterations) : BarrierContactForm(periodic_collision_mesh, dhat, avg_mass, use_area_weighting, use_improved_max_operator, use_physical_barrier, use_adaptive_barrier_stiffness, is_time_dependent, enable_shape_derivatives, broad_phase_method, ccd_tolerance, ccd_max_iterations), tiled_to_single_(tiled_to_single), n_single_dof_(tiled_to_single_.maxCoeff() + 1)
25 {
26 assert(tiled_to_single_.size() == collision_mesh_.full_num_vertices());
27
29
30 // const Eigen::MatrixXd displaced = collision_mesh_.displace_vertices(
31 // Eigen::MatrixXd::Zero(collision_mesh_.full_num_vertices(), collision_mesh_.dim()));
32
33 // io::OBJWriter::write(
34 // "tiled-rest.obj", displaced,
35 // collision_mesh_.edges(), collision_mesh_.faces());
36 }
37
39 {
40 const int dim = collision_mesh_.dim();
41 const auto &boundary_vertices = collision_mesh_.rest_positions();
42
43 std::vector<Eigen::Triplet<double>> entries;
44 for (int i = 0; i < collision_mesh_.num_vertices(); i++)
45 {
46 const int i_full = collision_mesh_.to_full_vertex_id(i);
47 for (int d = 0; d < dim; d++)
48 entries.emplace_back(tiled_to_single_(i_full) * dim + d, i_full * dim + d, 1);
49
50 for (int p = 0; p < dim; p++)
51 for (int q = 0; q < dim; q++)
52 entries.emplace_back(n_single_dof_ * dim + p * dim + q, i_full * dim + p, boundary_vertices(i, q));
53 }
54
55 proj.resize(n_single_dof_ * dim + dim * dim, tiled_to_single_.size() * dim);
56 proj.setZero();
57 proj.setFromTriplets(entries.begin(), entries.end());
58 }
59
60 Eigen::VectorXd PeriodicContactForm::single_to_tiled(const Eigen::VectorXd &x) const
61 {
62 const int dim = collision_mesh_.dim();
63 const auto &boundary_vertices = collision_mesh_.rest_positions();
64 assert(x.size() == n_single_dof_ * dim + dim * dim);
65
66 Eigen::VectorXd tiled_x;
67 tiled_x.setZero(tiled_to_single_.size() * dim);
68 for (int i = 0; i < collision_mesh_.num_vertices(); i++)
69 {
70 const int i_full = collision_mesh_.to_full_vertex_id(i);
71 for (int d = 0; d < dim; d++)
72 tiled_x(i_full * dim + d) += x(tiled_to_single_(i_full) * dim + d);
73
74 for (int p = 0; p < dim; p++)
75 for (int q = 0; q < dim; q++)
76 tiled_x(i_full * dim + p) += x(n_single_dof_ * dim + p * dim + q) * boundary_vertices(i, q);
77 }
78
79 return tiled_x;
80 }
81 Eigen::VectorXd PeriodicContactForm::tiled_to_single_grad(const Eigen::VectorXd &grad) const
82 {
83 const int dim = collision_mesh_.dim();
84 const auto &boundary_vertices = collision_mesh_.rest_positions();
85 assert(grad.size() == tiled_to_single_.size() * dim);
86
87 Eigen::VectorXd reduced_grad;
88 reduced_grad.setZero(n_single_dof_ * dim + dim * dim);
89 for (int i = 0; i < collision_mesh_.num_vertices(); i++)
90 {
91 const int i_full = collision_mesh_.to_full_vertex_id(i);
92 for (int d = 0; d < dim; d++)
93 reduced_grad(tiled_to_single_(i_full) * dim + d) += grad(i_full * dim + d);
94
95 for (int p = 0; p < dim; p++)
96 for (int q = 0; q < dim; q++)
97 reduced_grad(n_single_dof_ * dim + p * dim + q) += boundary_vertices(i, q) * grad(i_full * dim + p);
98 }
99
100 return reduced_grad;
101 }
102
103 void PeriodicContactForm::init(const Eigen::VectorXd &x)
104 {
106 }
107
108 double PeriodicContactForm::value_unweighted(const Eigen::VectorXd &x) const
109 {
111 }
112
113 void PeriodicContactForm::first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
114 {
116 gradv = tiled_to_single_grad(gradv);
117 }
118
119 void PeriodicContactForm::second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const
120 {
121 StiffnessMatrix hessian_full;
123
125 hessian = proj * hessian_full * proj.transpose();
126
127 // const Eigen::MatrixXd displaced = collision_mesh_.displace_vertices(utils::unflatten(single_to_tiled(x), collision_mesh_.dim()));
128
129 // io::OBJWriter::write(
130 // "tiled.obj", displaced,
131 // collision_mesh_.edges(), collision_mesh_.faces());
132 }
133
134 void PeriodicContactForm::update_quantities(const double t, const Eigen::VectorXd &x)
135 {
137 }
138
139 double PeriodicContactForm::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
140 {
142 }
143
144 void PeriodicContactForm::line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
145 {
147 }
148
149 void PeriodicContactForm::solution_changed(const Eigen::VectorXd &new_x)
150 {
152 }
153
154 void PeriodicContactForm::post_step(const polysolve::nonlinear::PostStepData &data)
155 {
157 polysolve::nonlinear::PostStepData(
158 data.iter_num, data.solver_info, single_to_tiled(data.x), single_to_tiled(data.grad)));
159 }
160
161 bool PeriodicContactForm::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
162 {
164 }
165
166 void PeriodicContactForm::update_barrier_stiffness(const Eigen::VectorXd &x, const Eigen::MatrixXd &grad_energy)
167 {
169 }
170
171 void PeriodicContactForm::force_periodic_shape_derivative(const State& state, const PeriodicMeshToMesh &periodic_mesh_map, const Eigen::VectorXd &periodic_mesh_representation, const ipc::NormalCollisions &contact_set, const Eigen::VectorXd &solution, const Eigen::VectorXd &adjoint_sol, Eigen::VectorXd &term) const
172 {
173 const int dim = collision_mesh_.dim();
174 const Eigen::MatrixXd displaced_surface = compute_displaced_surface(single_to_tiled(solution));
175
176 Eigen::VectorXd tiled_term;
177
178 {
179 StiffnessMatrix dq_h = collision_mesh_.to_full_dof(barrier_potential().shape_derivative(contact_set, collision_mesh_, displaced_surface));
180 tiled_term = dq_h.transpose() * single_to_tiled(adjoint_sol);
181 }
182
183 {
184 Eigen::VectorXd force;
185 force = barrier_potential().gradient(contact_set, collision_mesh_, displaced_surface);
186 force = collision_mesh_.to_full_dof(force);
187 Eigen::MatrixXd adjoint_affine = utils::unflatten(adjoint_sol.tail(dim * dim), dim);
188 for (int k = 0; k < collision_mesh_.num_vertices(); k++)
189 {
190 const int k_full = collision_mesh_.to_full_vertex_id(k);
191 tiled_term.segment(k_full * dim, dim) += adjoint_affine.transpose() * force.segment(k_full * dim, dim);
192 }
193 }
194
195 {
196 StiffnessMatrix hessian_full;
198 Eigen::VectorXd tmp = single_to_tiled(adjoint_sol).transpose() * hessian_full;
199 Eigen::MatrixXd sol_affine = utils::unflatten(solution.tail(dim * dim), dim);
200 for (int k = 0; k < collision_mesh_.num_vertices(); k++)
201 {
202 const int k_full = collision_mesh_.to_full_vertex_id(k);
203 tiled_term.segment(k_full * dim, dim) += sol_affine.transpose() * tmp.segment(k_full * dim, dim);
204 }
205 }
206
207 // chain rule from tiled to periodic
208 Eigen::VectorXd unit_term;
209 Eigen::MatrixXd affine_term;
210 affine_term.setZero(dim, dim);
211 unit_term.setZero(tiled_term.size());
212 Eigen::MatrixXd affine = utils::unflatten(periodic_mesh_representation.tail(dim * dim), dim).transpose();
213 for (int k = 0; k < collision_mesh_.num_vertices(); k++)
214 {
215 const int k_full = collision_mesh_.to_full_vertex_id(k);
216 affine_term += tiled_term.segment(k_full * dim, dim) * collision_mesh_.rest_positions().row(k);
217 unit_term.segment(k_full * dim, dim) += affine.transpose() * tiled_term.segment(k_full * dim, dim);
218 }
219
220 Eigen::VectorXd single_term = state.basis_nodes_to_gbasis_nodes * proj.topRows(dim * n_single_dof_) * unit_term;
221 single_term = utils::flatten(utils::unflatten(single_term, dim)(state.primitive_to_node(), Eigen::all));
222
223 term.setZero(periodic_mesh_representation.size());
224 for (int i = 0; i < state.n_geom_bases; i++)
225 term.segment(periodic_mesh_map.full_to_periodic(i) * dim, dim).array() += single_term.segment(i * dim, dim).array();
226 term.tail(dim * dim) = Eigen::Map<Eigen::VectorXd>(affine_term.data(), dim * dim, 1);
227
228 term *= weight();
229 }
230}
std::vector< Eigen::Triplet< double > > entries
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
std::vector< int > primitive_to_node() const
Definition State.cpp:227
int n_geom_bases
number of geometric bases
Definition State.hpp:182
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
virtual void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
Compute the first derivative of the value wrt x.
virtual void second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override
Compute the second derivative of the value wrt x.
const ipc::BarrierPotential & barrier_potential() const
virtual double value_unweighted(const Eigen::VectorXd &x) const override
Compute the contact barrier potential value.
virtual void update_barrier_stiffness(const Eigen::VectorXd &x, const Eigen::MatrixXd &grad_energy) override
Update the barrier stiffness based on the current elasticity energy.
void post_step(const polysolve::nonlinear::PostStepData &data) override
Update fields after a step in the optimization.
virtual void init(const Eigen::VectorXd &x) override
Initialize the form.
Eigen::MatrixXd compute_displaced_surface(const Eigen::VectorXd &x) const
Compute the displaced positions of the surface nodes.
virtual double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Determine the maximum step size allowable between the current and next solution.
virtual void solution_changed(const Eigen::VectorXd &new_x) override
Update cached fields upon a change in the solution.
double weight() const override
Get the form's multiplicative constant weight.
const ipc::CollisionMesh & collision_mesh_
Collision mesh.
virtual void line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
Initialize variables used during the line search.
virtual void update_quantities(const double t, const Eigen::VectorXd &x) override
Update time-dependent fields.
virtual bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
void post_step(const polysolve::nonlinear::PostStepData &data) override
Update fields after a step in the optimization.
Eigen::VectorXd single_to_tiled(const Eigen::VectorXd &x) const
void solution_changed(const Eigen::VectorXd &new_x) override
Update cached fields upon a change in the solution.
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the contact barrier potential value.
void line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
Initialize variables used during the line search.
void second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override
Compute the second derivative of the value wrt x.
void force_periodic_shape_derivative(const State &state, const PeriodicMeshToMesh &periodic_mesh_map, const Eigen::VectorXd &periodic_mesh_representation, const ipc::NormalCollisions &contact_set, const Eigen::VectorXd &solution, const Eigen::VectorXd &adjoint_sol, Eigen::VectorXd &term) const
double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Determine the maximum step size allowable between the current and next solution.
void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
Compute the first derivative of the value wrt x.
void init(const Eigen::VectorXd &x) override
Initialize the form.
PeriodicContactForm(const ipc::CollisionMesh &periodic_collision_mesh, const Eigen::VectorXi &tiled_to_single, const double dhat, const double avg_mass, const bool use_area_weighting, const bool use_improved_max_operator, const bool use_physical_barrier, const bool use_adaptive_barrier_stiffness, const bool is_time_dependent, const bool enable_shape_derivatives, const ipc::BroadPhaseMethod broad_phase_method, const double ccd_tolerance, const int ccd_max_iterations)
Construct a new Contact Form object.
Eigen::VectorXd tiled_to_single_grad(const Eigen::VectorXd &grad) const
void update_quantities(const double t, const Eigen::VectorXd &x) override
Update time-dependent fields.
void update_barrier_stiffness(const Eigen::VectorXd &x, const Eigen::MatrixXd &grad_energy) override
Update the barrier stiffness based on the current elasticity energy.
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22