Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SmoothContactForm.cpp
Go to the documentation of this file.
8
9#include <ipc/utils/eigen_ext.hpp>
10#include <ipc/barrier/adaptive_stiffness.hpp>
11#include <ipc/utils/world_bbox_diagonal_length.hpp>
12
13namespace polyfem::solver
14{
15 SmoothContactForm::SmoothContactForm(const ipc::CollisionMesh &collision_mesh,
16 const double dhat,
17 const double avg_mass,
18 const double alpha_t,
19 const double alpha_n,
20 const bool use_adaptive_dhat,
21 const double min_distance_ratio,
22 const bool use_adaptive_barrier_stiffness,
23 const bool is_time_dependent,
24 const bool enable_shape_derivatives,
25 const ipc::BroadPhaseMethod broad_phase_method,
26 const double ccd_tolerance,
27 const int ccd_max_iterations) : ContactForm(collision_mesh, dhat, avg_mass, use_adaptive_barrier_stiffness, is_time_dependent, enable_shape_derivatives, broad_phase_method, ccd_tolerance, ccd_max_iterations), params(dhat_, alpha_t, 0, alpha_n, 0, collision_mesh.dim() - 1), use_adaptive_dhat(use_adaptive_dhat), barrier_potential_(params)
28 {
29 params.set_adaptive_dhat_ratio(min_distance_ratio);
31 {
32 collision_set_.compute_adaptive_dhat(collision_mesh, collision_mesh.rest_positions(), params, broad_phase_);
34 logger().error("Adaptive dhat is not compatible with adaptive barrier stiffness");
35 }
36 }
37
38 void SmoothContactForm::update_barrier_stiffness(const Eigen::VectorXd &x, const Eigen::MatrixXd &grad_energy)
39 {
41 return;
42
43 log_and_throw_error("Adaptive barrier stiffness not implemented for SmoothContactForm!");
44 }
45
46 void SmoothContactForm::force_shape_derivative(const ipc::SmoothCollisions &collision_set, const Eigen::MatrixXd &solution, const Eigen::VectorXd &adjoint_sol, Eigen::VectorXd &term) const
47 {
48 StiffnessMatrix hessian = barrier_potential_.hessian(collision_set, collision_mesh_, compute_displaced_surface(solution), ipc::PSDProjectionMethod::NONE);
49 term = barrier_stiffness() * collision_mesh_.to_full_dof(hessian) * adjoint_sol;
50 }
51
52 void SmoothContactForm::update_collision_set(const Eigen::MatrixXd &displaced_surface)
53 {
54 // Store the previous value used to compute the constraint set to avoid duplicate computation.
55 static Eigen::MatrixXd cached_displaced_surface;
56 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
57 return;
58
60 collision_set_.build(
62 else
63 collision_set_.build(
65 cached_displaced_surface = displaced_surface;
66 }
67
72
73 Eigen::VectorXd SmoothContactForm::value_per_element_unweighted(const Eigen::VectorXd &x) const
74 {
75 log_and_throw_error("value_per_element_unweighted not implemented!");
76 }
77
78 void SmoothContactForm::first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
79 {
81 gradv = collision_mesh_.to_full_dof(gradv);
82 }
83
84 void SmoothContactForm::second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const
85 {
86 // {
87 // io::OBJWriter::write(
88 // "collision_mesh.obj",
89 // compute_displaced_surface(x),
90 // collision_mesh_.edges(), collision_mesh_.faces());
91 // }
92 POLYFEM_SCOPED_TIMER("barrier hessian");
93 hessian = barrier_potential_.hessian(collision_set_, collision_mesh_, compute_displaced_surface(x), project_to_psd_ ? ipc::PSDProjectionMethod::CLAMP : ipc::PSDProjectionMethod::NONE);
94 hessian = collision_mesh_.to_full_dof(hessian);
95 }
96
97 void SmoothContactForm::post_step(const polysolve::nonlinear::PostStepData &data)
98 {
99 const Eigen::MatrixXd displaced_surface = compute_displaced_surface(data.x);
100
101 const double curr_distance = collision_set_.compute_minimum_distance(collision_mesh_, displaced_surface);
102 const double curr_active_distance = collision_set_.compute_active_minimum_distance(collision_mesh_, displaced_surface);
103 if (!std::isinf(curr_distance))
104 {
105 const double ratio = sqrt(curr_distance) / dhat();
106 const auto log_level = (ratio < 1e-6) ? spdlog::level::err : ((ratio < 1e-4) ? spdlog::level::warn : spdlog::level::debug);
107 polyfem::logger().log(log_level, "Minimum distance during solve: {}, active distance: {}, dhat: {}", sqrt(curr_distance), sqrt(curr_active_distance), dhat());
108 }
109
110 if (data.iter_num == 0)
111 return;
112
114 {
116 {
117 const double prev_barrier_stiffness = barrier_stiffness();
118
119 barrier_stiffness_ = ipc::update_barrier_stiffness(
121 barrier_stiffness(), ipc::world_bbox_diagonal_length(displaced_surface), 1e-7);
122
123 if (barrier_stiffness() != prev_barrier_stiffness)
124 {
125 polyfem::logger().debug(
126 "updated barrier stiffness from {:g} to {:g}",
127 prev_barrier_stiffness, barrier_stiffness());
128 }
129 }
130 else
131 {
132 // TODO: missing feature
133 // update_barrier_stiffness(data.x);
134 }
135 }
136
137 prev_distance_ = curr_distance;
138 }
139} // namespace polyfem::solver
int x
#define POLYFEM_SCOPED_TIMER(...)
Definition Timer.hpp:10
Form representing the contact potential and forces.
const bool use_adaptive_barrier_stiffness_
If true, use an adaptive barrier stiffness.
bool use_cached_candidates_
If true, use the cached candidate set for the current solution.
Eigen::MatrixXd compute_displaced_surface(const Eigen::VectorXd &x) const
Compute the displaced positions of the surface nodes.
double barrier_stiffness() const
Get the current barrier stiffness.
double barrier_stiffness_
Barrier stiffness.
double max_barrier_stiffness_
Maximum barrier stiffness to use when using adaptive barrier stiffness.
double prev_distance_
Previous minimum distance between all elements.
const ipc::CollisionMesh & collision_mesh_
Collision mesh.
bool use_adaptive_barrier_stiffness() const
Get use_adaptive_barrier_stiffness.
ipc::Candidates candidates_
Cached candidate set for the current solution.
const std::shared_ptr< ipc::BroadPhase > broad_phase_
const bool is_time_dependent_
Is the simulation time dependent?
bool project_to_psd_
If true, the form's second derivative is projected to be positive semidefinite.
Definition Form.hpp:147
void update_collision_set(const Eigen::MatrixXd &displaced_surface) override
Update the cached candidate set for the current solution.
const ipc::SmoothCollisions & collision_set() const
ipc::SmoothContactPotential barrier_potential_
Contact potential.
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 second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override
Compute the second derivative of the value wrt x.
Eigen::VectorXd value_per_element_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form multiplied per element.
void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
Compute the first derivative of the value wrt x.
SmoothContactForm(const ipc::CollisionMesh &collision_mesh, const double dhat, const double avg_mass, const double alpha_t, const double alpha_n, const bool use_adaptive_dhat, const double min_distance_ratio, 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)
void post_step(const polysolve::nonlinear::PostStepData &data) override
Update fields after a step in the optimization.
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the contact barrier potential value.
void force_shape_derivative(const ipc::SmoothCollisions &collision_set, const Eigen::MatrixXd &solution, const Eigen::VectorXd &adjoint_sol, Eigen::VectorXd &term) const
ipc::SmoothCollisions collision_set_
Cached constraint set for the current solution.
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22