PolyFEM
Loading...
Searching...
No Matches
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_convergent_formulation,
17 const bool use_adaptive_barrier_stiffness,
18 const bool is_time_dependent,
19 const bool enable_shape_derivatives,
20 const ipc::BroadPhaseMethod broad_phase_method,
21 const double ccd_tolerance,
22 const int ccd_max_iterations) : ContactForm(periodic_collision_mesh, dhat, avg_mass, use_convergent_formulation, 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)
23 {
24 assert(tiled_to_single_.size() == collision_mesh_.full_num_vertices());
25
27
28 // const Eigen::MatrixXd displaced = collision_mesh_.displace_vertices(
29 // Eigen::MatrixXd::Zero(collision_mesh_.full_num_vertices(), collision_mesh_.dim()));
30
31 // io::OBJWriter::write(
32 // "tiled-rest.obj", displaced,
33 // collision_mesh_.edges(), collision_mesh_.faces());
34 }
35
37 {
38 const int dim = collision_mesh_.dim();
39 const auto &boundary_vertices = collision_mesh_.rest_positions();
40
41 std::vector<Eigen::Triplet<double>> entries;
42 for (int i = 0; i < collision_mesh_.num_vertices(); i++)
43 {
44 const int i_full = collision_mesh_.to_full_vertex_id(i);
45 for (int d = 0; d < dim; d++)
46 entries.emplace_back(tiled_to_single_(i_full) * dim + d, i_full * dim + d, 1);
47
48 for (int p = 0; p < dim; p++)
49 for (int q = 0; q < dim; q++)
50 entries.emplace_back(n_single_dof_ * dim + p * dim + q, i_full * dim + p, boundary_vertices(i, q));
51 }
52
53 proj.resize(n_single_dof_ * dim + dim * dim, tiled_to_single_.size() * dim);
54 proj.setZero();
55 proj.setFromTriplets(entries.begin(), entries.end());
56 }
57
58 Eigen::VectorXd PeriodicContactForm::single_to_tiled(const Eigen::VectorXd &x) const
59 {
60 const int dim = collision_mesh_.dim();
61 const auto &boundary_vertices = collision_mesh_.rest_positions();
62 assert(x.size() == n_single_dof_ * dim + dim * dim);
63
64 Eigen::VectorXd tiled_x;
65 tiled_x.setZero(tiled_to_single_.size() * dim);
66 for (int i = 0; i < collision_mesh_.num_vertices(); i++)
67 {
68 const int i_full = collision_mesh_.to_full_vertex_id(i);
69 for (int d = 0; d < dim; d++)
70 tiled_x(i_full * dim + d) += x(tiled_to_single_(i_full) * dim + d);
71
72 for (int p = 0; p < dim; p++)
73 for (int q = 0; q < dim; q++)
74 tiled_x(i_full * dim + p) += x(n_single_dof_ * dim + p * dim + q) * boundary_vertices(i, q);
75 }
76
77 return tiled_x;
78 }
79 Eigen::VectorXd PeriodicContactForm::tiled_to_single_grad(const Eigen::VectorXd &grad) const
80 {
81 const int dim = collision_mesh_.dim();
82 const auto &boundary_vertices = collision_mesh_.rest_positions();
83 assert(grad.size() == tiled_to_single_.size() * dim);
84
85 Eigen::VectorXd reduced_grad;
86 reduced_grad.setZero(n_single_dof_ * dim + dim * dim);
87 for (int i = 0; i < collision_mesh_.num_vertices(); i++)
88 {
89 const int i_full = collision_mesh_.to_full_vertex_id(i);
90 for (int d = 0; d < dim; d++)
91 reduced_grad(tiled_to_single_(i_full) * dim + d) += grad(i_full * dim + d);
92
93 for (int p = 0; p < dim; p++)
94 for (int q = 0; q < dim; q++)
95 reduced_grad(n_single_dof_ * dim + p * dim + q) += boundary_vertices(i, q) * grad(i_full * dim + p);
96 }
97
98 return reduced_grad;
99 }
100
101 void PeriodicContactForm::init(const Eigen::VectorXd &x)
102 {
104 }
105
106 double PeriodicContactForm::value_unweighted(const Eigen::VectorXd &x) const
107 {
109 }
110
111 void PeriodicContactForm::first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
112 {
114 gradv = tiled_to_single_grad(gradv);
115 }
116
117 void PeriodicContactForm::second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const
118 {
119 StiffnessMatrix hessian_full;
121
123 hessian = proj * hessian_full * proj.transpose();
124
125 // const Eigen::MatrixXd displaced = collision_mesh_.displace_vertices(utils::unflatten(single_to_tiled(x), collision_mesh_.dim()));
126
127 // io::OBJWriter::write(
128 // "tiled.obj", displaced,
129 // collision_mesh_.edges(), collision_mesh_.faces());
130 }
131
132 void PeriodicContactForm::update_quantities(const double t, const Eigen::VectorXd &x)
133 {
135 }
136
137 double PeriodicContactForm::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
138 {
140 }
141
142 void PeriodicContactForm::line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
143 {
145 }
146
147 void PeriodicContactForm::solution_changed(const Eigen::VectorXd &new_x)
148 {
150 }
151
152 void PeriodicContactForm::post_step(const polysolve::nonlinear::PostStepData &data)
153 {
155 polysolve::nonlinear::PostStepData(
156 data.iter_num, data.solver_info, single_to_tiled(data.x), single_to_tiled(data.grad)));
157 }
158
159 bool PeriodicContactForm::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
160 {
162 }
163
164 void PeriodicContactForm::update_barrier_stiffness(const Eigen::VectorXd &x, const Eigen::MatrixXd &grad_energy)
165 {
167 }
168
169 void PeriodicContactForm::force_periodic_shape_derivative(const State& state, const PeriodicMeshToMesh &periodic_mesh_map, const Eigen::VectorXd &periodic_mesh_representation, const ipc::Collisions &contact_set, const Eigen::VectorXd &solution, const Eigen::VectorXd &adjoint_sol, Eigen::VectorXd &term)
170 {
171 const int dim = collision_mesh_.dim();
172 const Eigen::MatrixXd displaced_surface = compute_displaced_surface(single_to_tiled(solution));
173
174 Eigen::VectorXd tiled_term;
175
176 {
177 StiffnessMatrix dq_h = collision_mesh_.to_full_dof(barrier_potential_.shape_derivative(contact_set, collision_mesh_, displaced_surface));
178 tiled_term = dq_h.transpose() * single_to_tiled(adjoint_sol);
179 }
180
181 {
182 Eigen::VectorXd force;
183 force = barrier_potential_.gradient(contact_set, collision_mesh_, displaced_surface);
184 force = collision_mesh_.to_full_dof(force);
185 Eigen::MatrixXd adjoint_affine = utils::unflatten(adjoint_sol.tail(dim * dim), dim);
186 for (int k = 0; k < collision_mesh_.num_vertices(); k++)
187 {
188 const int k_full = collision_mesh_.to_full_vertex_id(k);
189 tiled_term.segment(k_full * dim, dim) += adjoint_affine.transpose() * force.segment(k_full * dim, dim);
190 }
191 }
192
193 {
194 StiffnessMatrix hessian_full;
196 Eigen::VectorXd tmp = single_to_tiled(adjoint_sol).transpose() * hessian_full;
197 Eigen::MatrixXd sol_affine = utils::unflatten(solution.tail(dim * dim), dim);
198 for (int k = 0; k < collision_mesh_.num_vertices(); k++)
199 {
200 const int k_full = collision_mesh_.to_full_vertex_id(k);
201 tiled_term.segment(k_full * dim, dim) += sol_affine.transpose() * tmp.segment(k_full * dim, dim);
202 }
203 }
204
205 // chain rule from tiled to periodic
206 Eigen::VectorXd unit_term;
207 Eigen::MatrixXd affine_term;
208 affine_term.setZero(dim, dim);
209 unit_term.setZero(tiled_term.size());
210 Eigen::MatrixXd affine = utils::unflatten(periodic_mesh_representation.tail(dim * dim), dim).transpose();
211 for (int k = 0; k < collision_mesh_.num_vertices(); k++)
212 {
213 const int k_full = collision_mesh_.to_full_vertex_id(k);
214 affine_term += tiled_term.segment(k_full * dim, dim) * collision_mesh_.rest_positions().row(k);
215 unit_term.segment(k_full * dim, dim) += affine.transpose() * tiled_term.segment(k_full * dim, dim);
216 }
217
218 Eigen::VectorXd single_term = state.basis_nodes_to_gbasis_nodes * proj.topRows(dim * n_single_dof_) * unit_term;
219 single_term = utils::flatten(utils::unflatten(single_term, dim)(state.primitive_to_node(), Eigen::all));
220
221 term.setZero(periodic_mesh_representation.size());
222 for (int i = 0; i < state.n_geom_bases; i++)
223 term.segment(periodic_mesh_map.full_to_periodic(i) * dim, dim).array() += single_term.segment(i * dim, dim).array();
224 term.tail(dim * dim) = Eigen::Map<Eigen::VectorXd>(affine_term.data(), dim * dim, 1);
225
226 term *= weight();
227 }
228}
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:228
int n_geom_bases
number of geometric bases
Definition State.hpp:182
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:699
Form representing the contact potential and forces.
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.
void post_step(const polysolve::nonlinear::PostStepData &data) override
Update fields after a step in the optimization.
virtual void second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override
Compute the second derivative of the value wrt x.
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 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::BarrierPotential barrier_potential_
virtual double value_unweighted(const Eigen::VectorXd &x) const override
Compute the contact barrier potential value.
const ipc::CollisionMesh & collision_mesh_
Collision mesh.
void line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
Initialize variables used during the line search.
virtual void update_barrier_stiffness(const Eigen::VectorXd &x, const Eigen::MatrixXd &grad_energy)
Update the barrier stiffness based on the current elasticity energy.
virtual void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
Compute the first derivative of the value wrt x.
void update_quantities(const double t, const Eigen::VectorXd &x) override
Update time-dependent fields.
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
void force_periodic_shape_derivative(const State &state, const PeriodicMeshToMesh &periodic_mesh_map, const Eigen::VectorXd &periodic_mesh_representation, const ipc::Collisions &contact_set, const Eigen::VectorXd &solution, const Eigen::VectorXd &adjoint_sol, Eigen::VectorXd &term)
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.
PeriodicContactForm(const ipc::CollisionMesh &periodic_collision_mesh, const Eigen::VectorXi &tiled_to_single, const double dhat, const double avg_mass, const bool use_convergent_formulation, 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.
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.
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.
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