DirichletEnergy.hpp

#pragma once
 
#include <polysolve/nonlinear/Problem.hpp>
#include <wmtk/Types.hpp>
 
namespace wmtk::optimization {
 
class DirichletEnergy2D : public polysolve::nonlinear::Problem
{
public:
    using typename polysolve::nonlinear::Problem::Scalar;
    using typename polysolve::nonlinear::Problem::THessian;
    using typename polysolve::nonlinear::Problem::TVector;
 
    DirichletEnergy2D(std::vector<std::array<double, 4>>& cells);
 
    TVector initial_position() const;
 
    double value(const TVector& x) override;
    void gradient(const TVector& x, TVector& gradv) override;
    void hessian(const TVector& x, THessian& hessian) override
    {
        log_and_throw_error("Sparse functions do not exist, use dense solver");
    }
    void hessian(const TVector& x, MatrixXd& hessian) override;
 
    void solution_changed(const TVector& new_x) override {}
 
private:
    std::vector<std::array<double, 4>> m_cells;
};
 
class BiharmonicEnergy2D : public polysolve::nonlinear::Problem
{
public:
    using typename polysolve::nonlinear::Problem::Scalar;
    using typename polysolve::nonlinear::Problem::THessian;
    using typename polysolve::nonlinear::Problem::TVector;
 
    BiharmonicEnergy2D(
        const std::array<Vector2d, 3>& pts,
        const double& M,
        const Vector3d& L_w,
        const double weight = 1);
 
    TVector initial_position() const;
 
    double value(const TVector& x) override;
    void gradient(const TVector& x, TVector& gradv) override;
    void hessian(const TVector& x, THessian& hessian) override
    {
        log_and_throw_error("Sparse functions do not exist, use dense solver");
    }
    void hessian(const TVector& x, MatrixXd& hessian) override;
 
    void solution_changed(const TVector& new_x) override {}
 
    static void
    local_mass_and_stiffness(const std::array<Vector2d, 3>& pts, double& M, Vector3d& L_w);
 
    static void
    uniform_mass_and_stiffness(const std::array<Vector2d, 3>& pts, double& M, Vector3d& L_w);
 
private:
    std::array<Vector2d, 3> m_pts;
    Vector3d m_L_w_row;
    double m_M; // the mass of the optimized vertex
    double m_M_inv; // the inverse mass, 1 / M
 
    Vector3d m_LTML_row;
 
    double m_weight;
};
 
class BiharmonicEnergy3D : public polysolve::nonlinear::Problem
{
public:
    using typename polysolve::nonlinear::Problem::Scalar;
    using typename polysolve::nonlinear::Problem::THessian;
    using typename polysolve::nonlinear::Problem::TVector;
 
    BiharmonicEnergy3D(
        const MatrixXd& pts,
        const double& M,
        const VectorXd& L_w,
        const double weight = 1);
 
    TVector initial_position() const;
 
    double value(const TVector& x) override;
    void gradient(const TVector& x, TVector& gradv) override;
    void hessian(const TVector& x, THessian& hessian) override
    {
        log_and_throw_error("Sparse functions do not exist, use dense solver");
    }
    void hessian(const TVector& x, MatrixXd& hessian) override;
 
    void solution_changed(const TVector& new_x) override {}
 
    static void global_mass_and_stiffness(
        const MatrixXd& pts,
        const MatrixXi& tris,
        Eigen::SparseMatrix<double>& M,
        Eigen::SparseMatrix<double>& L_w);
 
    static void extract_local_mass_and_stiffness(
        const size_t vid,
        const Eigen::SparseMatrix<double>& M_glob,
        const Eigen::SparseMatrix<double>& L_w_glob,
        double& M_loc,
        VectorXd& L_w_loc);
 
    static void adjacency_from_stiffness(
        const size_t vid,
        const Eigen::SparseMatrix<double>& L_w_glob,
        std::vector<size_t>& adj);
 
    static void uniform_mass_and_stiffness(const MatrixXd& pts, double& M, VectorXd& L_w);
 
private:
    MatrixXd m_pts;
    VectorXd m_L_w_row;
    double m_M; // the mass of the optimized vertex
    double m_M_inv; // the inverse mass, 1 / M
 
    VectorXd m_LTML_row;
 
    double m_weight;
};
 
} // namespace wmtk::optimization