Loading [MathJax]/jax/output/HTML-CSS/config.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FrictionForm.cpp
Go to the documentation of this file.
1#include "FrictionForm.hpp"
4
7
8namespace polyfem::solver
9{
11 const ipc::CollisionMesh &collision_mesh,
12 const std::shared_ptr<time_integrator::ImplicitTimeIntegrator> time_integrator,
13 const double epsv,
14 const double mu,
15 const ipc::BroadPhaseMethod broad_phase_method,
16 const ContactForm &contact_form,
17 const int n_lagging_iters)
18 : collision_mesh_(collision_mesh),
19 time_integrator_(time_integrator),
20 epsv_(epsv),
21 mu_(mu),
22 broad_phase_method_(broad_phase_method),
23 n_lagging_iters_(n_lagging_iters < 0 ? std::numeric_limits<int>::max() : n_lagging_iters),
24 contact_form_(contact_form),
25 friction_potential_(epsv)
26 {
27 assert(epsv_ > 0);
28 }
29
31 const Eigen::MatrixXd &prev_solution,
32 const Eigen::MatrixXd &solution,
33 const Eigen::MatrixXd &adjoint,
34 const ipc::TangentialCollisions &friction_constraints_set,
35 Eigen::VectorXd &term)
36 {
37 Eigen::MatrixXd U = collision_mesh_.vertices(utils::unflatten(solution, collision_mesh_.dim()));
38 Eigen::MatrixXd U_prev = collision_mesh_.vertices(utils::unflatten(prev_solution, collision_mesh_.dim()));
39
40 // TODO: use the time integration to compute the velocity
41 const Eigen::MatrixXd velocities = (U - U_prev) / time_integrator_->dt();
42
43 StiffnessMatrix hess;
44 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm*>(&contact_form_))
45 {
46 hess = -friction_potential_.force_jacobian(
47 friction_constraints_set,
48 collision_mesh_, collision_mesh_.rest_positions(),
49 /*lagged_displacements=*/U_prev, velocities,
50 barrier_contact->barrier_potential(),
51 barrier_contact->barrier_stiffness(),
52 ipc::FrictionPotential::DiffWRT::REST_POSITIONS);
53 }
54
55 term = collision_mesh_.to_full_dof(hess).transpose() * adjoint;
56 }
57
58 Eigen::MatrixXd FrictionForm::compute_displaced_surface(const Eigen::VectorXd &x) const
59 {
61 }
62
63 Eigen::MatrixXd FrictionForm::compute_surface_velocities(const Eigen::VectorXd &x) const
64 {
65 // In the case of a static problem, the velocity is the displacement
66 const Eigen::VectorXd v = time_integrator_ != nullptr ? time_integrator_->compute_velocity(x) : x;
67 return collision_mesh_.map_displacements(utils::unflatten(v, collision_mesh_.dim()));
68 }
69
70 double FrictionForm::dv_dx() const
71 {
72 return time_integrator_ != nullptr ? time_integrator_->dv_dx() : 1;
73 }
74
79
80 void FrictionForm::first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
81 {
82 const Eigen::VectorXd grad_friction = friction_potential_.gradient(
84 gradv = collision_mesh_.to_full_dof(grad_friction);
85 }
86
87 void FrictionForm::second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const
88 {
89 POLYFEM_SCOPED_TIMER("friction hessian");
90
91 ipc::PSDProjectionMethod psd_projection_method;
92
94 {
95 psd_projection_method = ipc::PSDProjectionMethod::CLAMP;
96 }
97 else
98 {
99 psd_projection_method = ipc::PSDProjectionMethod::NONE;
100 }
101
102 hessian = dv_dx() * friction_potential_.hessian( //
104
105 hessian = collision_mesh_.to_full_dof(hessian);
106 }
107
108 void FrictionForm::update_lagging(const Eigen::VectorXd &x, const int iter_num)
109 {
110 const Eigen::MatrixXd displaced_surface = compute_displaced_surface(x);
111
112 auto broad_phase = ipc::build_broad_phase(broad_phase_method_);
113 if (const auto barrier_contact = dynamic_cast<const BarrierContactForm*>(&contact_form_))
114 {
115 ipc::NormalCollisions collision_set;
116 collision_set.set_use_area_weighting(barrier_contact->use_area_weighting());
117 collision_set.set_use_improved_max_approximator(barrier_contact->use_improved_max_operator());
118
119 collision_set.set_enable_shape_derivatives(barrier_contact->enable_shape_derivatives());
120 collision_set.build(
121 collision_mesh_, displaced_surface, barrier_contact->dhat(), /*dmin=*/0, broad_phase);
122
124 collision_mesh_, displaced_surface, collision_set,
125 barrier_contact->barrier_potential(), barrier_contact->barrier_stiffness(), Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_);
126 }
127 else if (const auto smooth_contact = dynamic_cast<const SmoothContactForm*>(&contact_form_))
128 {
129 ipc::SmoothCollisions collision_set;
130 if (smooth_contact->using_adaptive_dhat())
131 collision_set.compute_adaptive_dhat(collision_mesh_, collision_mesh_.rest_positions(), smooth_contact->get_params(), broad_phase);
132 collision_set.build(
133 collision_mesh_, displaced_surface, smooth_contact->get_params(),
134 smooth_contact->using_adaptive_dhat(), broad_phase);
135
136 friction_collision_set_.build_for_smooth_contact(
137 collision_mesh_, displaced_surface,
138 collision_set, smooth_contact->get_params(), contact_form_.barrier_stiffness(), Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_);
139 }
140 else
141 {
142 throw std::runtime_error("Unknown contact form");
143 }
144 }
145} // namespace polyfem::solver
int x
#define POLYFEM_SCOPED_TIMER(...)
Definition Timer.hpp:10
Form representing the contact potential and forces.
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.
bool project_to_psd_
If true, the form's second derivative is projected to be positive semidefinite.
Definition Form.hpp:147
const ipc::BroadPhaseMethod broad_phase_method_
Broad-phase method used for distance computation and collision detection.
FrictionForm(const ipc::CollisionMesh &collision_mesh, const std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator, const double epsv, const double mu, const ipc::BroadPhaseMethod broad_phase_method, const ContactForm &contact_form, const int n_lagging_iters)
Construct a new Friction Form object.
void update_lagging(const Eigen::VectorXd &x, const int iter_num) override
Update lagged fields.
void second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override
Compute the second derivative of the value wrt x.
const ipc::CollisionMesh & collision_mesh_
Reference to the collision mesh.
void force_shape_derivative(const Eigen::MatrixXd &prev_solution, const Eigen::MatrixXd &solution, const Eigen::MatrixXd &adjoint, const ipc::TangentialCollisions &friction_constraints_set, Eigen::VectorXd &term)
const std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator_
Pointer to the time integrator.
Eigen::MatrixXd compute_surface_velocities(const Eigen::VectorXd &x) const
Compute the surface velocities.
const ipc::FrictionPotential friction_potential_
Eigen::MatrixXd compute_displaced_surface(const Eigen::VectorXd &x) const
Compute the displaced positions of the surface nodes.
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form.
double dv_dx() const
Compute the derivative of the velocities wrt x.
void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
Compute the first derivative of the value wrt x.
const ContactForm & contact_form_
necessary to have the barrier stiffnes, maybe clean me
ipc::TangentialCollisions friction_collision_set_
Lagged friction constraint set.
const double epsv_
Smoothing factor between static and dynamic friction.
const double mu_
Global coefficient of friction.
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22