Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BarrierForms.cpp
Go to the documentation of this file.
1#include "BarrierForms.hpp"
2#include <polyfem/State.hpp>
4
5namespace polyfem::solver
6{
7 namespace
8 {
9 class QuadraticBarrier : public ipc::Barrier
10 {
11 public:
12 QuadraticBarrier(const double weight = 1) : weight_(weight) {}
13
14 double operator()(const double d, const double dhat) const override
15 {
16 if (d > dhat)
17 return 0;
18 else
19 return weight_ * (d - dhat) * (d - dhat);
20 }
21 double first_derivative(const double d, const double dhat) const override
22 {
23 if (d > dhat)
24 return 0;
25 else
26 return 2 * weight_ * (d - dhat);
27 }
28 double second_derivative(const double d, const double dhat) const override
29 {
30 if (d > dhat)
31 return 0;
32 else
33 return 2 * weight_;
34 }
35 double units(const double dhat) const override
36 {
37 return dhat * dhat;
38 }
39
40 private:
41 const double weight_;
42 };
43
44 } // namespace
45
46 CollisionBarrierForm::CollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat, const double dmin)
47 : AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), dmin_(dmin), barrier_potential_(dhat)
48 {
52 [this](const std::string &p) { return this->state_.resolve_input_path(p); },
54
55 Eigen::MatrixXd V;
58
59 broad_phase_method_ = ipc::BroadPhaseMethod::HASH_GRID;
60 }
61
62 double CollisionBarrierForm::value_unweighted(const Eigen::VectorXd &x) const
63 {
64 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
65 return barrier_potential_(collision_set, collision_mesh_, displaced_surface);
66 }
67
68 void CollisionBarrierForm::compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
69 {
71 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
72 const Eigen::VectorXd grad = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
74 });
75 }
76
77 void CollisionBarrierForm::solution_changed(const Eigen::VectorXd &x)
78 {
80
81 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
82 build_collision_set(displaced_surface);
83 }
84
85 bool CollisionBarrierForm::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
86 {
87 const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
88 const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
89
90 // Skip CCD if the displacement is zero.
91 if ((V1 - V0).lpNorm<Eigen::Infinity>() == 0.0)
92 return true;
93
94 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
95 bool is_valid = ipc::is_step_collision_free(
97 collision_mesh_.vertices(V0),
98 collision_mesh_.vertices(V1),
99 dmin_,
100 ipc::build_broad_phase(broad_phase_method_),
101 tight_inclusion_ccd);
102
103 return is_valid;
104 }
105
106 double CollisionBarrierForm::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
107 {
108 const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
109 const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
110
111 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
112 double max_step = ipc::compute_collision_free_stepsize(
114 collision_mesh_.vertices(V0),
115 collision_mesh_.vertices(V1),
116 dmin_,
117 ipc::build_broad_phase(broad_phase_method_),
118 tight_inclusion_ccd);
119
120 adjoint_logger().info("Objective {}: max step size is {}.", name(), max_step);
121
122 return max_step;
123 }
124
125 void CollisionBarrierForm::build_collision_set(const Eigen::MatrixXd &displaced_surface)
126 {
127 static Eigen::MatrixXd cached_displaced_surface;
128 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
129 return;
130
131 collision_set.build(collision_mesh_, displaced_surface, dhat_, dmin_, ipc::build_broad_phase(broad_phase_method_));
132
133 cached_displaced_surface = displaced_surface;
134 }
135
136 Eigen::VectorXd CollisionBarrierForm::get_updated_mesh_nodes(const Eigen::VectorXd &x) const
137 {
138 Eigen::VectorXd X = X_init;
141 }
142
144 const State &state,
145 const std::vector<int> &boundary_ids,
146 const double dhat,
147 const bool use_log_barrier,
148 const double dmin)
149 : CollisionBarrierForm(variable_to_simulations, state, dhat, dmin),
150 boundary_ids_(boundary_ids)
151 {
152 for (const auto &id : boundary_ids_)
153 boundary_ids_to_dof_[id] = std::set<int>();
154
156
157 if (!use_log_barrier)
158 barrier_potential_.set_barrier(std::make_shared<QuadraticBarrier>());
159 }
160
162 {
163 Eigen::MatrixXd node_positions;
164 Eigen::MatrixXi boundary_edges, boundary_triangles;
165 std::vector<Eigen::Triplet<double>> displacement_map_entries;
167 node_positions, boundary_edges, boundary_triangles, displacement_map_entries);
168
169 std::vector<bool> is_on_surface;
170 is_on_surface.resize(node_positions.rows(), false);
171
173 Eigen::MatrixXd points, uv, normals;
174 Eigen::VectorXd weights;
175 Eigen::VectorXi global_primitive_ids;
176 for (const auto &lb : state_.total_local_boundary)
177 {
178 const int e = lb.element_id();
179 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, state_.n_boundary_samples(), *state_.mesh, false, uv, points, normals, weights, global_primitive_ids);
180
181 if (!has_samples)
182 continue;
183
184 const basis::ElementBases &gbs = state_.geom_bases()[e];
185
186 vals.compute(e, state_.mesh->is_volume(), points, gbs, gbs);
187
188 for (int i = 0; i < lb.size(); ++i)
189 {
190 const int primitive_global_id = lb.global_primitive_id(i);
191 const auto nodes = gbs.local_nodes_for_primitive(primitive_global_id, *state_.mesh);
192 const int boundary_id = state_.mesh->get_boundary_id(primitive_global_id);
193
194 if (!std::count(boundary_ids_.begin(), boundary_ids_.end(), boundary_id))
195 continue;
196
197 for (long n = 0; n < nodes.size(); ++n)
198 {
199 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
200 is_on_surface[v.global[0].index] = true;
201 assert(v.global[0].index < node_positions.rows());
202 boundary_ids_to_dof_[boundary_id].insert(v.global[0].index);
203 }
204 }
205 }
206
207 Eigen::SparseMatrix<double> displacement_map;
208 if (!displacement_map_entries.empty())
209 {
210 displacement_map.resize(node_positions.rows(), state_.n_geom_bases);
211 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
212 }
213
214 // Fix boundary edges and boundary triangles to exclude vertices not on triangles
215 Eigen::MatrixXi boundary_edges_alt(0, 2), boundary_triangles_alt(0, 3);
216 {
217 for (int i = 0; i < boundary_edges.rows(); ++i)
218 {
219 bool on_surface = true;
220 for (int j = 0; j < boundary_edges.cols(); ++j)
221 on_surface &= is_on_surface[boundary_edges(i, j)];
222 if (on_surface)
223 {
224 boundary_edges_alt.conservativeResize(boundary_edges_alt.rows() + 1, 2);
225 boundary_edges_alt.row(boundary_edges_alt.rows() - 1) = boundary_edges.row(i);
226 }
227 }
228
229 if (state_.mesh->is_volume())
230 {
231 for (int i = 0; i < boundary_triangles.rows(); ++i)
232 {
233 bool on_surface = true;
234 for (int j = 0; j < boundary_triangles.cols(); ++j)
235 on_surface &= is_on_surface[boundary_triangles(i, j)];
236 if (on_surface)
237 {
238 boundary_triangles_alt.conservativeResize(boundary_triangles_alt.rows() + 1, 3);
239 boundary_triangles_alt.row(boundary_triangles_alt.rows() - 1) = boundary_triangles.row(i);
240 }
241 }
242 }
243 else
244 boundary_triangles_alt.resize(0, 0);
245 }
246
247 collision_mesh_ = ipc::CollisionMesh(is_on_surface,
248 std::vector<bool>(is_on_surface.size(), false),
249 node_positions,
250 boundary_edges_alt,
251 boundary_triangles_alt,
252 displacement_map);
253
254 can_collide_cache_.resize(collision_mesh_.num_vertices(), collision_mesh_.num_vertices());
255 for (int i = 0; i < can_collide_cache_.rows(); ++i)
256 {
257 int dof_idx_i = collision_mesh_.to_full_vertex_id(i);
258 if (!is_on_surface[dof_idx_i])
259 continue;
260 for (int j = 0; j < can_collide_cache_.cols(); ++j)
261 {
262 int dof_idx_j = collision_mesh_.to_full_vertex_id(j);
263 if (!is_on_surface[dof_idx_j])
264 continue;
265
266 bool collision_allowed = true;
267 for (const auto &id : boundary_ids_)
268 if (boundary_ids_to_dof_[id].count(dof_idx_i) && boundary_ids_to_dof_[id].count(dof_idx_j))
269 collision_allowed = false;
270 can_collide_cache_(i, j) = collision_allowed;
271 }
272 }
273
274 collision_mesh_.can_collide = [&](size_t vi, size_t vj) {
275 return (bool)can_collide_cache_(vi, vj);
276 };
277
278 collision_mesh_.init_area_jacobians();
279 }
280
281 DeformedCollisionBarrierForm::DeformedCollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat)
282 : AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), barrier_potential_(dhat)
283 {
285 log_and_throw_adjoint_error("[{}] Should use linear FE basis!", name());
286
290 [this](const std::string &p) { return this->state_.resolve_input_path(p); },
292
293 Eigen::MatrixXd V;
296
297 broad_phase_method_ = ipc::BroadPhaseMethod::HASH_GRID;
298 }
299
300 double DeformedCollisionBarrierForm::value_unweighted(const Eigen::VectorXd &x) const
301 {
302 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
303
304 return barrier_potential_(collision_set, collision_mesh_, displaced_surface);
305 }
306
307 void DeformedCollisionBarrierForm::compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
308 {
310 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
311 const Eigen::VectorXd grad = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
313 });
314 }
315
317 {
319
320 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
321 build_collision_set(displaced_surface);
322 }
323
324 bool DeformedCollisionBarrierForm::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
325 {
326 // const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
327 // const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
328
329 // // Skip CCD if the displacement is zero.
330 // if ((V1 - V0).lpNorm<Eigen::Infinity>() == 0.0)
331 // return true;
332
333 // bool is_valid = ipc::is_step_collision_free(
334 // collision_mesh_,
335 // collision_mesh_.vertices(V0),
336 // collision_mesh_.vertices(V1),
337 // broad_phase_method_,
338 // 1e-6, 1e6);
339
340 return true; // is_valid;
341 }
342
343 double DeformedCollisionBarrierForm::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
344 {
345 // const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
346 // const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
347
348 // double max_step = ipc::compute_collision_free_stepsize(
349 // collision_mesh_,
350 // collision_mesh_.vertices(V0),
351 // collision_mesh_.vertices(V1),
352 // broad_phase_method_, 1e-6, 1e6);
353
354 return 1; // max_step;
355 }
356
357 void DeformedCollisionBarrierForm::build_collision_set(const Eigen::MatrixXd &displaced_surface)
358 {
359 static Eigen::MatrixXd cached_displaced_surface;
360 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
361 return;
362
363 collision_set.build(collision_mesh_, displaced_surface, dhat_, 0, ipc::build_broad_phase(broad_phase_method_));
364
365 cached_displaced_surface = displaced_surface;
366 }
367
374
376 const VariableToSimulationGroup &variable_to_simulations,
377 const State &state,
378 const json &args)
379 : StaticForm(variable_to_simulations),
380 state_(state),
381 params_(state.args["contact"]["dhat"], state.args["contact"]["alpha_t"], 0, state.args["contact"]["alpha_n"], 0, state.mesh->is_volume() ? 2 : 1),
382 potential_(params_)
383 {
384 auto tmp_ids = args["surface_selection"].get<std::vector<int>>();
385 boundary_ids_ = std::set(tmp_ids.begin(), tmp_ids.end());
386
388 }
389
391 {
392 // Deep copy and change the can_collide() function
394
395 // const int num_fe_nodes = state_.n_bases - state_.obstacle.n_vertices();
396
397 // collision_mesh_.can_collide = [this, num_fe_nodes](size_t vi, size_t vj) {
398 // return this->collision_mesh_.to_full_vertex_id(vi) >= num_fe_nodes || this->collision_mesh_.to_full_vertex_id(vj) >= num_fe_nodes;
399 // };
400
401 std::vector<int> is_obstacle(state_.n_bases);
402 for (int e = 0; e < state_.bases.size(); e++)
403 {
404 const auto &b = state_.bases[e];
405 if (state_.mesh->get_body_id(e) == 1)
406 for (const auto &bs : b.bases)
407 {
408 for (const auto &g : bs.global())
409 {
410 is_obstacle[g.index] = true;
411 }
412 }
413 }
414
415 collision_mesh_.can_collide = [this, is_obstacle](size_t vi, size_t vj) {
416 return is_obstacle[this->collision_mesh_.to_full_vertex_id(vi)] || is_obstacle[this->collision_mesh_.to_full_vertex_id(vj)];
417 };
418 }
419
420 ipc::SmoothCollisions SmoothContactForceForm::get_smooth_collision_set(const Eigen::MatrixXd &displaced_surface)
421 {
422 ipc::SmoothCollisions collisions;
423 const auto smooth_contact = dynamic_cast<const SmoothContactForm*>(state_.solve_data.contact_form.get());
424 collisions.build(collision_mesh_, displaced_surface, smooth_contact->get_params(), smooth_contact->using_adaptive_dhat(), smooth_contact->get_broad_phase());
425 return collisions;
426 }
427
428 double SmoothContactForceForm::value_unweighted_step(const int time_step, const Eigen::VectorXd &x) const
429 {
430 assert(state_.solve_data.contact_form != nullptr);
431
432 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
433
434 Eigen::VectorXd forces = collision_mesh_.to_full_dof(potential_.gradient(collisions_, collision_mesh_, displaced_surface));
435
436 // return forces.squaredNorm();
437
438 Eigen::VectorXd coeff(forces.size());
439 coeff.setZero();
440 coeff(Eigen::seq(1, coeff.size(), collision_mesh_.dim())).array() = 1;
441 return (coeff.array() * forces.array()).matrix().squaredNorm() / 2;
442 }
443
444 Eigen::VectorXd SmoothContactForceForm::compute_adjoint_rhs_step(const int time_step, const Eigen::VectorXd &x, const State &state) const
445 {
446 assert(state_.solve_data.contact_form != nullptr);
447
448 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
449
450 Eigen::VectorXd forces = potential_.gradient(collisions_, collision_mesh_, displaced_surface);
451 forces = collision_mesh_.to_full_dof(forces);
452
453 StiffnessMatrix hessian = potential_.hessian(collisions_, collision_mesh_, displaced_surface, ipc::PSDProjectionMethod::NONE);
454 hessian = collision_mesh_.to_full_dof(hessian);
455
456 Eigen::VectorXd coeff(forces.size());
457 coeff.setZero();
458 coeff(Eigen::seq(1, coeff.size(), collision_mesh_.dim())).array() = 1;
459 return weight() * (hessian * (coeff.array() * forces.array()).matrix());
460 }
461
462 void SmoothContactForceForm::compute_partial_gradient_step(const int time_step, const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
463 {
464 assert(state_.solve_data.contact_form != nullptr);
465
466 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
467
468 Eigen::VectorXd forces = potential_.gradient(collisions_, collision_mesh_, displaced_surface);
469 forces = collision_mesh_.to_full_dof(forces);
470
471 StiffnessMatrix hessian = potential_.hessian(collisions_, collision_mesh_, displaced_surface, ipc::PSDProjectionMethod::NONE);
472 hessian = collision_mesh_.to_full_dof(hessian);
473
474 gradv = weight() * variable_to_simulations_.apply_parametrization_jacobian(ParameterType::Shape, &state_, x, [this, &x, &forces, &hessian]() {
475 // Eigen::VectorXd grads = 2 * hessian.transpose() * forces;
476
477 Eigen::VectorXd coeff(forces.size());
478 coeff.setZero();
479 coeff(Eigen::seq(1, coeff.size(), collision_mesh_.dim())).array() = 1;
480 Eigen::VectorXd grads = (hessian * (coeff.array() * forces.array()).matrix());
481
482 grads = state_.basis_nodes_to_gbasis_nodes * grads;
484 });
485 }
486
487 void SmoothContactForceForm::solution_changed_step(const int time_step, const Eigen::VectorXd &x)
488 {
490 const Eigen::MatrixXd displaced_surface = collision_mesh_.displace_vertices(utils::unflatten(state_.diff_cached.u(time_step), collision_mesh_.dim()));
491 collisions_ = get_smooth_collision_set(displaced_surface);
492 }
493} // namespace polyfem::solver
int V
ElementAssemblyValues vals
Definition Assembler.cpp:22
const double weight_
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
void get_vertices(Eigen::MatrixXd &vertices) const
Definition StateDiff.cpp:72
int n_bases
number of bases
Definition State.hpp:178
int n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:263
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:455
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:473
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
json args
main input arguments containing all defaults
Definition State.hpp:101
solver::DiffCache diff_cached
Definition State.hpp:656
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
int n_geom_bases
number of geometric bases
Definition State.hpp:182
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
Definition State.cpp:1361
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
stores per local bases evaluations
std::vector< basis::Local2Global > global
stores per element basis values at given quadrature points and geometric mapping
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
Eigen::VectorXi local_nodes_for_primitive(const int local_index, const mesh::Mesh &mesh) const
static void extract_boundary_mesh(const mesh::Mesh &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &node_positions, Eigen::MatrixXi &boundary_edges, Eigen::MatrixXi &boundary_triangles, std::vector< Eigen::Triplet< double > > &displacement_map_entries)
extracts the boundary mesh
Definition OutData.cpp:145
virtual void solution_changed(const Eigen::VectorXd &new_x) override
Update cached fields upon a change in the solution.
const VariableToSimulationGroup variable_to_simulations_
ipc::BarrierPotential barrier_potential_
std::string name() const override
void build_collision_set(const Eigen::MatrixXd &displaced_surface)
ipc::BroadPhaseMethod broad_phase_method_
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
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.
CollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat, const double dmin=0)
void compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
void solution_changed(const Eigen::VectorXd &x) override
Update cached fields upon a change in the solution.
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form.
Eigen::VectorXd get_updated_mesh_nodes(const Eigen::VectorXd &x) const
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
void compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
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.
DeformedCollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat)
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form.
const ipc::BarrierPotential barrier_potential_
void build_collision_set(const Eigen::MatrixXd &displaced_surface)
Eigen::VectorXd get_updated_mesh_nodes(const Eigen::VectorXd &x) const
void solution_changed(const Eigen::VectorXd &x) override
Update cached fields upon a change in the solution.
Eigen::VectorXd u(int step) const
virtual double weight() const
Get the form's multiplicative constant weight.
Definition Form.hpp:128
LayerThicknessForm(const VariableToSimulationGroup &variable_to_simulations, const State &state, const std::vector< int > &boundary_ids, const double dhat, const bool use_log_barrier=false, const double dmin=0)
std::map< int, std::set< int > > boundary_ids_to_dof_
ipc::SmoothCollisions get_smooth_collision_set(const Eigen::MatrixXd &displaced_surface)
SmoothContactForceForm(const VariableToSimulationGroup &variable_to_simulations, const State &state, const json &args)
void solution_changed_step(const int time_step, const Eigen::VectorXd &x) override
Eigen::VectorXd compute_adjoint_rhs_step(const int time_step, const Eigen::VectorXd &x, const State &state) const override
double value_unweighted_step(const int time_step, const Eigen::VectorXd &x) const override
void compute_partial_gradient_step(const int time_step, const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
ipc::SmoothContactPotential potential_
std::shared_ptr< solver::ContactForm > contact_form
A collection of VariableToSimulation.
void compute_state_variable(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, Eigen::VectorXd &state_variable) const
Evaluate the variable to simulations and overwrite the state_variable based on x.
virtual Eigen::VectorXd apply_parametrization_jacobian(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, const std::function< Eigen::VectorXd()> &grad) const
Maps the partial gradient wrt.
static bool boundary_quadrature(const mesh::LocalBoundary &local_boundary, const int order, const mesh::Mesh &mesh, const bool skip_computation, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::MatrixXd &normals, Eigen::VectorXd &weights, Eigen::VectorXi &global_primitive_ids)
Eigen::VectorXd map_primitive_to_node_order(const State &state, const Eigen::VectorXd &primitives)
Eigen::VectorXd map_node_to_primitive_order(const State &state, const Eigen::VectorXd &nodes)
std::tuple< bool, int, Tree > is_valid(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u, const double threshold)
Definition Jacobian.cpp:297
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.
spdlog::logger & adjoint_logger()
Retrieves the current logger for adjoint.
Definition Logger.cpp:30
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22