AMIPSForm.cpp

Go to the documentation of this file.

#include <polyfem/optimization/forms/AMIPSForm.hpp>
 
#include <polyfem/varforms/diff/DifferentiableVarForm.hpp>
 
#include <polyfem/mesh/mesh2D/Mesh2D.hpp>
#include <polyfem/mesh/mesh3D/Mesh3D.hpp>
 
#include <polyfem/assembler/AssemblyValsCache.hpp>
#include <polyfem/basis/ElementBases.hpp>
#include <polyfem/assembler/Assembler.hpp>
#include <polyfem/assembler/AssemblerUtils.hpp>
#include <polyfem/assembler/GenericElastic.hpp>
#include <polyfem/assembler/AMIPSEnergy.hpp>
#include <polyfem/utils/GeometryUtils.hpp>
#include <polyfem/optimization/AdjointTools.hpp>
 
#include <cmath>
 
namespace polyfem::solver
{
    namespace
    {
        void scaled_jacobian(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, Eigen::VectorXd &quality)
        {
            const int dim = F.cols() - 1;
 
            quality.setZero(F.rows());
            if (dim == 2)
            {
                for (int i = 0; i < F.rows(); i++)
                {
                    Eigen::RowVector3d e0;
                    e0(2) = 0;
                    e0.head(2) = V.row(F(i, 2)) - V.row(F(i, 1));
                    Eigen::RowVector3d e1;
                    e1(2) = 0;
                    e1.head(2) = V.row(F(i, 0)) - V.row(F(i, 2));
                    Eigen::RowVector3d e2;
                    e2(2) = 0;
                    e2.head(2) = V.row(F(i, 1)) - V.row(F(i, 0));
 
                    double l0 = e0.norm();
                    double l1 = e1.norm();
                    double l2 = e2.norm();
 
                    double A = 0.5 * (e0.cross(e1)).norm();
                    double Lmax = std::max(l0 * l1, std::max(l1 * l2, l0 * l2));
 
                    quality(i) = 2 * A * (2 / sqrt(3)) / Lmax;
                }
            }
            else
            {
                for (int i = 0; i < F.rows(); i++)
                {
                    Eigen::RowVector3d e0 = V.row(F(i, 1)) - V.row(F(i, 0));
                    Eigen::RowVector3d e1 = V.row(F(i, 2)) - V.row(F(i, 1));
                    Eigen::RowVector3d e2 = V.row(F(i, 0)) - V.row(F(i, 2));
                    Eigen::RowVector3d e3 = V.row(F(i, 3)) - V.row(F(i, 0));
                    Eigen::RowVector3d e4 = V.row(F(i, 3)) - V.row(F(i, 1));
                    Eigen::RowVector3d e5 = V.row(F(i, 3)) - V.row(F(i, 2));
 
                    double l0 = e0.norm();
                    double l1 = e1.norm();
                    double l2 = e2.norm();
                    double l3 = e3.norm();
                    double l4 = e4.norm();
                    double l5 = e5.norm();
 
                    double J = std::abs((e0.cross(e3)).dot(e2));
 
                    double a1 = l0 * l2 * l3;
                    double a2 = l0 * l1 * l4;
                    double a3 = l1 * l2 * l5;
                    double a4 = l3 * l4 * l5;
 
                    double a = std::max({a1, a2, a3, a4, J});
                    quality(i) = J * sqrt(2) / a;
                }
            }
        }
    } // namespace
 
    double MinJacobianForm::value_unweighted(const Eigen::VectorXd &x) const
    {
        const bool is_volume = varform_->get_mesh().is_volume();
        double min_jacs = std::numeric_limits<double>::max();
        for (size_t e = 0; e < varform_->primary_space().geometry_basis_list().size(); ++e)
        {
            if (varform_->get_mesh().is_polytope(e))
                continue;
 
            const auto &gbasis = varform_->primary_space().geometry_basis_list()[e];
            const int n_local_bases = int(gbasis.bases.size());
 
            quadrature::Quadrature quad;
            gbasis.compute_quadrature(quad);
 
            std::vector<assembler::AssemblyValues> tmp;
 
            Eigen::MatrixXd dx = Eigen::MatrixXd::Zero(quad.points.rows(), quad.points.cols());
            Eigen::MatrixXd dy = Eigen::MatrixXd::Zero(quad.points.rows(), quad.points.cols());
            Eigen::MatrixXd dz;
            if (is_volume)
                dz = Eigen::MatrixXd::Zero(quad.points.rows(), quad.points.cols());
 
            gbasis.evaluate_grads(quad.points, tmp);
 
            for (int j = 0; j < n_local_bases; ++j)
            {
                const basis::Basis &b = gbasis.bases[j];
 
                for (std::size_t ii = 0; ii < b.global().size(); ++ii)
                {
                    dx += tmp[j].grad.col(0) * b.global()[ii].node * b.global()[ii].val;
                    dy += tmp[j].grad.col(1) * b.global()[ii].node * b.global()[ii].val;
                    if (is_volume)
                        dz += tmp[j].grad.col(2) * b.global()[ii].node * b.global()[ii].val;
                }
            }
 
            for (long i = 0; i < dx.rows(); ++i)
            {
                if (is_volume)
                {
                    Eigen::Matrix3d tmp;
                    tmp << dx.row(i), dy.row(i), dz.row(i);
                    min_jacs = std::min(min_jacs, tmp.determinant());
                }
                else
                {
                    Eigen::Matrix2d tmp;
                    tmp << dx.row(i), dy.row(i);
                    min_jacs = std::min(min_jacs, tmp.determinant());
                }
            }
        }
 
        return min_jacs;
    }
 
    void MinJacobianForm::compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
    {
        log_and_throw_adjoint_error("{} is not differentiable!", name());
    }
 
    AMIPSForm::AMIPSForm(const VariableToSimulationGroup &variable_to_simulation, std::shared_ptr<const varform::DifferentiableVarForm> varform)
        : AdjointForm(variable_to_simulation),
          varform_(std::move(varform))
    {
        amips_energy_ = assembler::AssemblerUtils::make_assembler("AMIPS");
        amips_energy_->set_size(varform_->get_mesh().dimension());
 
        json use_rest = {};
        use_rest["use_rest_pose"] = true;
        amips_energy_->add_multimaterial(0, use_rest, varform_->get_units(), varform_->get_root_path());
 
        Eigen::MatrixXd V;
        varform_->get_vertices(V);
        varform_->get_elements(F);
        X_rest = utils::flatten(V);
        init_geom_bases_ = varform_->primary_space().geometry_basis_list();
    }
 
    double AMIPSForm::value_unweighted(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd X = get_updated_mesh_nodes(x);
 
        return amips_energy_->assemble_energy(varform_->get_mesh().is_volume(), init_geom_bases_, init_geom_bases_, init_ass_vals_cache_, 0, 0, AdjointTools::map_primitive_to_node_order(*varform_, X - X_rest), Eigen::VectorXd());
    }
 
    void AMIPSForm::compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
    {
        gradv = weight() * variable_to_simulations_.apply_parametrization_jacobian(ParameterType::Shape, *varform_, x, [this, &x]() {
            const Eigen::VectorXd X = get_updated_mesh_nodes(x);
            Eigen::MatrixXd grad;
            amips_energy_->assemble_gradient(varform_->get_mesh().is_volume(), varform_->primary_space().geometry->n_bases, init_geom_bases_, init_geom_bases_, init_ass_vals_cache_, 0, 0, AdjointTools::map_primitive_to_node_order(*varform_, X - X_rest), Eigen::VectorXd(), grad); // grad wrt. gbases
            return AdjointTools::map_node_to_primitive_order(*varform_, grad);
        });
    }
 
    bool AMIPSForm::is_step_valid(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
    {
        Eigen::VectorXd X = get_updated_mesh_nodes(x1);
        Eigen::MatrixXd V1 = utils::unflatten(X, varform_->get_mesh().dimension());
        bool flipped = utils::is_flipped(V1, F);
 
        if (flipped)
            adjoint_logger().trace("[{}] Step flips elements.", name());
 
        return !flipped;
    }
} // namespace polyfem::solver