FrictionForm.cpp

Go to the documentation of this file.

#include "FrictionForm.hpp"
#include "BarrierContactForm.hpp"
#include "GCPContactForm.hpp"
#include "ESPContactForm.hpp"
 
#include <polyfem/utils/Timer.hpp>
#include <polyfem/utils/MatrixUtils.hpp>
 
#include <ipc/broad_phase/create_broad_phase.hpp>
#include <ipc/collisions/normal/normal_collisions.hpp>
#include <ipc/gcp/gcp_collisions.hpp>
 
#include <Eigen/Core>
 
#include <cassert>
#include <limits>
#include <memory>
#include <stdexcept>
 
namespace polyfem::solver
{
    FrictionForm::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)
        : collision_mesh_(collision_mesh),
          time_integrator_(time_integrator),
          epsv_(epsv),
          mu_(mu),
          broad_phase_method_(broad_phase_method),
          n_lagging_iters_(n_lagging_iters < 0 ? std::numeric_limits<int>::max() : n_lagging_iters),
          contact_form_(contact_form),
          friction_potential_(epsv)
    {
        assert(epsv_ > 0);
    }
 
    Eigen::MatrixXd FrictionForm::compute_displaced_surface(const Eigen::VectorXd &x) const
    {
        return contact_form_.compute_displaced_surface(x);
    }
 
    Eigen::MatrixXd FrictionForm::compute_surface_velocities(const Eigen::VectorXd &x) const
    {
        // In the case of a static problem, the velocity is the displacement
        const Eigen::VectorXd v = time_integrator_ != nullptr ? time_integrator_->compute_velocity(x) : x;
        return collision_mesh_.map_displacements(utils::unflatten(v, collision_mesh_.dim()));
    }
 
    double FrictionForm::dv_dx() const
    {
        return time_integrator_ != nullptr ? time_integrator_->dv_dx() : 1;
    }
 
    double FrictionForm::value_unweighted(const Eigen::VectorXd &x) const
    {
        return friction_potential_(friction_collision_set_, collision_mesh_, compute_surface_velocities(x)) / dv_dx();
    }
 
    void FrictionForm::first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
    {
        const Eigen::VectorXd grad_friction = friction_potential_.gradient(
            friction_collision_set_, collision_mesh_, compute_surface_velocities(x));
        gradv = collision_mesh_.to_full_dof(grad_friction);
    }
 
    void FrictionForm::second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const
    {
        POLYFEM_SCOPED_TIMER("friction hessian");
 
        ipc::PSDProjectionMethod psd_projection_method;
 
        if (project_to_psd_)
        {
            psd_projection_method = ipc::PSDProjectionMethod::CLAMP;
        }
        else
        {
            psd_projection_method = ipc::PSDProjectionMethod::NONE;
        }
 
        hessian = dv_dx() * friction_potential_.hessian( //
                      friction_collision_set_, collision_mesh_, compute_surface_velocities(x), psd_projection_method);
 
        hessian = collision_mesh_.to_full_dof(hessian);
    }
 
    void FrictionForm::update_lagging(const Eigen::VectorXd &x, const int iter_num)
    {
        const Eigen::MatrixXd displaced_surface = compute_displaced_surface(x);
 
        auto broad_phase = ipc::create_broad_phase(broad_phase_method_);
        if (const auto barrier_contact = dynamic_cast<const BarrierContactForm *>(&contact_form_))
        {
            ipc::NormalCollisions collision_set;
            collision_set.set_use_area_weighting(barrier_contact->use_area_weighting());
            // Copy the collision set type directly rather than going through
            // use_improved_max_operator(), which collapses it to a bool and so
            // would silently downgrade OGC to IPC here.
            collision_set.set_collision_set_type(barrier_contact->collision_set().collision_set_type());
 
            collision_set.set_enable_shape_derivatives(barrier_contact->enable_shape_derivatives());
            collision_set.build(
                collision_mesh_, displaced_surface, barrier_contact->dhat(), /*dmin=*/0, broad_phase.get());
 
            ipc::BarrierPotential bp = barrier_contact->barrier_potential();
            bp.set_stiffness(barrier_contact->barrier_stiffness());
            friction_collision_set_.build(
                collision_mesh_, displaced_surface, collision_set,
                bp, Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_, Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_);
        }
        else if (const auto gcp = dynamic_cast<const GCPContactForm *>(&contact_form_))
        {
            ipc::GCPCollisions collision_set;
            if (gcp->using_adaptive_dhat())
                collision_set.compute_adaptive_dhat(collision_mesh_, collision_mesh_.rest_positions(), gcp->get_params(), broad_phase.get());
            collision_set.build(
                collision_mesh_, displaced_surface, gcp->get_params(),
                gcp->using_adaptive_dhat(), broad_phase.get());
 
            friction_collision_set_.build(
                collision_mesh_, displaced_surface,
                collision_set, gcp->get_params(), contact_form_.barrier_stiffness(), Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_, Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_);
        }
        else if (const auto esp_contact = dynamic_cast<const ESPContactForm *>(&contact_form_))
        {
            ipc::ESPCollisions collision_set;
            collision_set.build(
                collision_mesh_, displaced_surface, esp_contact->get_params(),
                broad_phase.get());
 
            friction_collision_set_.build(
                collision_mesh_, displaced_surface, collision_set,
                esp_contact->get_params(), contact_form_.barrier_stiffness(),
                Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_,
                Eigen::VectorXd::Ones(collision_mesh_.num_vertices()) * mu_);
        }
        else
        {
            throw std::runtime_error("Unknown contact form");
        }
    }
} // namespace polyfem::solver