AMIPSEnergy.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/assembler/GenericElastic.hpp>
#include <polyfem/assembler/MatParams.hpp>
#include <polyfem/utils/Types.hpp>
#include <polyfem/utils/MatrixUtils.hpp>
 
#include <polyfem/Common.hpp>
 
#include <Eigen/Dense>
 
#include <functional>
#include <iostream>
#include <vector>
 
namespace polyfem::assembler
{
    class AMIPSEnergy : public GenericElastic<AMIPSEnergy>
    {
    public:
        AMIPSEnergy();
 
        // sets material params
        void add_multimaterial(const int index, const json &params, const Units &units, const io::ResourceIO &resources) override;
 
        std::string name() const override { return "AMIPS"; }
        std::map<std::string, ParamFunc> parameters() const override;
 
        bool allow_inversion() const override { return false; }
 
        bool real_def_grad(const RowVectorNd &p, const double t, const int el_id) const override { return use_rest_pose(p, t, el_id); }
 
        template <typename T>
        T elastic_energy(
            const RowVectorNd &p,
            const double t,
            const int el_id,
            DefGradMatrix<T> &def_grad) const
        {
            typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
 
            const double power = 2. / size();
            const bool rest_pose = use_rest_pose(p, t, el_id);
 
            AutoDiffGradMat standard;
 
            if (size() == 2)
                standard = get_standard<2, T>(size(), rest_pose);
            else
                standard = get_standard<3, T>(size(), rest_pose);
 
            if (!rest_pose)
                def_grad = def_grad * standard;
 
            const T det = polyfem::utils::determinant(def_grad);
            if (det <= 0)
            {
                return T(std::nan(""));
            }
 
            const T powJ = pow(det, power);
            const double weight = weight_(p, t, el_id);
            return T(weight) * (def_grad.transpose() * def_grad).trace() / powJ; //+ barrier<T>::value(det);
        }
 
    private:
        bool use_rest_pose(const RowVectorNd &p, const double t, const int el_id) const;
        GenericMatParam use_rest_pose_;
        GenericMatParam weight_;
 
        template <int dimt, class T>
        static Eigen::Matrix<T, dimt, dimt> get_standard(const int dim, const bool use_rest_pose)
        {
            Eigen::Matrix<double, dimt, dimt> standard(dim, dim);
            if (use_rest_pose)
            {
                standard.setIdentity();
            }
            else
            {
                if (dim == 2)
                    standard << 1, 0,
                        0.5, std::sqrt(3) / 2;
                else
                    standard << 1, 0, 0,
                        0.5, std::sqrt(3) / 2., 0,
                        0.5, 0.5 / std::sqrt(3), std::sqrt(2. / 3.);
                standard = standard.inverse().transpose().eval();
            }
 
            Eigen::Matrix<T, dimt, dimt> res(dim, dim);
            for (int i = 0; i < dim; ++i)
            {
                for (int j = 0; j < dim; ++j)
                {
                    res(i, j) = T(standard(i, j));
                }
            }
 
            return res;
        }
 
    public:
        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> gradient(
            const RowVectorNd &p,
            const double t,
            const int el_id,
            const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &F) const override;
 
        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9> hessian(
            const RowVectorNd &p,
            const double t,
            const int el_id,
            const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &F) const override;
    };
 
} // namespace polyfem::assembler