VariableToSimulation.cpp

Go to the documentation of this file.

#include <polyfem/optimization/forms/VariableToSimulation.hpp>
 
#include <polyfem/State.hpp>
#include <polyfem/Common.hpp>
#include <polyfem/io/MatrixIO.hpp>
#include <polyfem/assembler/ViscousDamping.hpp>
#include <polyfem/optimization/Optimizations.hpp>
#include <polyfem/optimization/StateDiff.hpp>
#include <polyfem/optimization/parametrization/NodeCompositeParametrizations.hpp>
#include <polyfem/utils/Logger.hpp>
#include <polyfem/utils/Types.hpp>
#include <polyfem/mesh/mesh2D/Mesh2D.hpp>
#include <polyfem/mesh/mesh3D/Mesh3D.hpp>
 
#include <Eigen/Core>
#include <spdlog/fmt/fmt.h>
 
#include <cassert>
#include <functional>
#include <memory>
#include <string>
#include <vector>
 
namespace polyfem::solver
{
 
    void VariableToSimulation::set_output_indexing(const json &args)
    {
        const std::string composite_map_type = args["composite_map_type"];
        const State &state = *(states[0]);
        if (composite_map_type == "none")
        {
            output_indexing_.resize(0);
        }
        else if (composite_map_type == "indices")
        {
            if (args["composite_map_indices"].is_string())
            {
                Eigen::MatrixXi tmp_mat;
                polyfem::io::read_matrix(state.resolve_input_path(args["composite_map_indices"].get<std::string>()), tmp_mat);
                output_indexing_ = tmp_mat;
            }
            else if (args["composite_map_indices"].is_array())
                output_indexing_ = args["composite_map_indices"];
            else
                log_and_throw_adjoint_error("Invalid composite map indices type!");
        }
        else
            log_and_throw_adjoint_error("Unknown composite_map_type!");
    }
 
    Eigen::VectorXi VariableToSimulation::get_output_indexing(const Eigen::VectorXd &x) const
    {
        const int out_size = parametrization.size(x.size());
        if (output_indexing_.size() == out_size || out_size == 0)
            return output_indexing_;
        else if (output_indexing_.size() == 0)
        {
            Eigen::VectorXi ind;
            ind.setLinSpaced(out_size, 0, out_size - 1);
            return ind;
        }
        else
            log_and_throw_adjoint_error(fmt::format("[{}] Indexing size and output size of the Parametrization do not match! {} vs {}", name(), output_indexing_.size(), out_size));
        return Eigen::VectorXi();
    }
 
    Eigen::VectorXd VariableToSimulation::apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const
    {
        return parametrization.apply_jacobian(term(get_output_indexing(x)), x);
    }
 
    Eigen::VectorXd VariableToSimulation::inverse_eval()
    {
        log_and_throw_adjoint_error("[{}] inverse_eval not implemented!", name());
        return Eigen::VectorXd();
    }
 
    void VariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        log_and_throw_adjoint_error("[{}] update_state not implemented!", name());
    }
 
    Eigen::VectorXd VariableToSimulationGroup::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd adjoint_term = Eigen::VectorXd::Zero(x.size());
        for (const auto &v2s : data)
            adjoint_term += v2s->compute_adjoint_term(x);
        return adjoint_term;
    }
 
    void VariableToSimulationGroup::compute_state_variable(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, Eigen::VectorXd &state_variable) const
    {
        for (const auto &v2s : data)
        {
            if (v2s->get_parameter_type() != type)
                continue;
 
            const Eigen::VectorXd var = v2s->parametrization.eval(x);
            for (const auto &state : v2s->states)
            {
                if (state.get() != state_ptr)
                    continue;
 
                state_variable(v2s->get_output_indexing(x)) = var;
            }
        }
    }
 
    Eigen::VectorXd VariableToSimulationGroup::apply_parametrization_jacobian(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, const std::function<Eigen::VectorXd()> &grad) const
    {
        Eigen::VectorXd gradv = Eigen::VectorXd::Zero(x.size());
        for (const auto &v2s : data)
        {
            if (v2s->get_parameter_type() != type)
                continue;
 
            for (const auto &state : v2s->states)
            {
                if (state.get() != state_ptr)
                    continue;
 
                gradv += v2s->apply_parametrization_jacobian(grad(), x);
            }
        }
        return gradv;
    }
 
    void ShapeVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        for (auto state : states)
        {
            const int dim = state->mesh->dimension();
 
            // If indices include one vertex entry, we assume it include all entries of this vertex.
            // for (int i = 0; i < indices.size(); i += dim)
            //  for (int j = 0; j < dim; j++)
            //      assert(indices(i + j) == indices(i) + j);
 
            for (int i = 0; i < indices.size(); ++i)
            {
                const int vid = indices(i) / dim;
                Eigen::VectorXd p = state->mesh->point(vid);
                p(indices(i) - vid * dim) = state_variable(i);
                state->mesh->set_point(vid, p);
            }
        }
    }
    Eigen::VectorXd ShapeVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
 
            if (state->problem->is_time_dependent())
            {
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                Eigen::MatrixXd adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                AdjointTools::dJ_shape_transient_adjoint_term(*state, *diff_cache, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                if (!state->is_homogenization())
                    AdjointTools::dJ_shape_static_adjoint_term(*state, *diff_cache, diff_cache->u(0), get_adjoint_mat(*state, *diff_cache, 0), cur_term);
                else
                    AdjointTools::dJ_shape_homogenization_adjoint_term(*state, *diff_cache, diff_cache->u(0), get_adjoint_mat(*state, *diff_cache, 0), cur_term);
            }
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    Eigen::VectorXd ShapeVariableToSimulation::inverse_eval()
    {
        const int dim = states[0]->mesh->dimension();
        const int npts = states[0]->mesh->n_vertices();
 
        Eigen::VectorXd x;
        Eigen::VectorXi indices = get_output_indexing(x);
 
        if (indices.size() == 0)
            indices.setLinSpaced(npts * dim, 0, npts * dim - 1);
 
        Eigen::MatrixXd V;
        states[0]->get_vertices(V);
        if (indices.maxCoeff() >= V.size())
            log_and_throw_adjoint_error("Output indices larger than DoFs of vertices!");
        x = utils::flatten(V)(indices);
 
        return parametrization.inverse_eval(x);
    }
    void ShapeVariableToSimulation::set_output_indexing(const json &args)
    {
        const std::string composite_map_type = args["composite_map_type"];
        const State &state = *(states[0]);
 
        if (composite_map_type == "interior" || composite_map_type == "boundary" || composite_map_type == "boundary_excluding_surface")
        {
            std::vector<int> active_dimensions = args["active_dimensions"];
            if (active_dimensions.size() == 0)
                for (int d = 0; d < state.mesh->dimension(); d++)
                    active_dimensions.push_back(d);
 
            if (composite_map_type == "interior")
            {
                VariableToInteriorNodes variable_to_node(state, active_dimensions, args["volume_selection"]);
                output_indexing_ = variable_to_node.get_output_indexing();
            }
            else if (composite_map_type == "boundary")
            {
                VariableToBoundaryNodes variable_to_node(state, active_dimensions, args["surface_selection"]);
                output_indexing_ = variable_to_node.get_output_indexing();
            }
            else if (composite_map_type == "boundary_excluding_surface")
            {
                const std::vector<int> excluded_surfaces = args["surface_selection"];
                VariableToBoundaryNodesExclusive variable_to_node(state, active_dimensions, excluded_surfaces);
                output_indexing_ = variable_to_node.get_output_indexing();
            }
        }
        else
            VariableToSimulation::set_output_indexing(args);
    }
 
    void ElasticVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        for (auto state : states)
        {
            const int n_elem = state->bases.size();
            assert(n_elem * 2 == state_variable.size());
            state->assembler->update_lame_params(state_variable.segment(0, n_elem), state_variable.segment(n_elem, n_elem));
        }
    }
    Eigen::VectorXd ElasticVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
 
            if (state->problem->is_time_dependent())
            {
 
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                Eigen::MatrixXd adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                AdjointTools::dJ_material_transient_adjoint_term(*state, *diff_cache, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                AdjointTools::dJ_material_static_adjoint_term(*state, diff_cache->u(0), adjoint_p, cur_term);
            }
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    Eigen::VectorXd ElasticVariableToSimulation::inverse_eval()
    {
        auto &state = *(states[0]);
        auto params_map = state.assembler->parameters();
 
        auto search_lambda = params_map.find("lambda");
        auto search_mu = params_map.find("mu");
        if (search_lambda == params_map.end() || search_mu == params_map.end())
        {
            log_and_throw_adjoint_error("[{}] Failed to find Lame parameters!", name());
            return Eigen::VectorXd();
        }
 
        Eigen::VectorXd lambdas(state.mesh->n_elements());
        Eigen::VectorXd mus(state.mesh->n_elements());
        for (int e = 0; e < state.mesh->n_elements(); e++)
        {
            RowVectorNd barycenter;
            if (!state.mesh->is_volume())
            {
                const auto &mesh2d = *dynamic_cast<mesh::Mesh2D *>(state.mesh.get());
                barycenter = mesh2d.face_barycenter(e);
            }
            else
            {
                const auto &mesh3d = *dynamic_cast<mesh::Mesh3D *>(state.mesh.get());
                barycenter = mesh3d.cell_barycenter(e);
            }
            lambdas(e) = search_lambda->second(RowVectorNd::Zero(state.mesh->dimension()), barycenter, 0., e);
            mus(e) = search_mu->second(RowVectorNd::Zero(state.mesh->dimension()), barycenter, 0., e);
        }
        state.assembler->update_lame_params(lambdas, mus);
 
        Eigen::VectorXd params(lambdas.size() + mus.size());
        params << lambdas, mus;
 
        return parametrization.inverse_eval(params);
    }
 
    void FrictionCoeffientVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        assert(state_variable.size() == 1);
        assert(state_variable(0) >= 0);
        for (auto state : states)
            state->args["contact"]["friction_coefficient"] = state_variable(0);
    }
    Eigen::VectorXd FrictionCoeffientVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
 
            if (state->problem->is_time_dependent())
            {
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                Eigen::MatrixXd adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                AdjointTools::dJ_friction_transient_adjoint_term(*state, *diff_cache, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                log_and_throw_adjoint_error("[{}] Gradient in static simulations not implemented!", name());
            }
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    Eigen::VectorXd FrictionCoeffientVariableToSimulation::inverse_eval()
    {
        log_and_throw_adjoint_error("[{}] inverse_eval not implemented!", name());
        return Eigen::VectorXd();
    }
 
    void DampingCoeffientVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        assert(state_variable.size() == 2);
        json damping_param = {
            {"psi", state_variable(0)},
            {"phi", state_variable(1)},
        };
        for (auto state : states)
        {
            if (!state->args["materials"].is_array())
            {
                state->args["materials"]["psi"] = damping_param["psi"];
                state->args["materials"]["phi"] = damping_param["phi"];
            }
            else
            {
                for (auto &arg : state->args["materials"])
                {
                    arg["psi"] = damping_param["psi"];
                    arg["phi"] = damping_param["phi"];
                }
            }
 
            if (state->damping_assembler)
                state->damping_assembler->add_multimaterial(0, damping_param, state->units);
        }
        logger().info("[{}] Current params: {}, {}", name(), state_variable(0), state_variable(1));
    }
    Eigen::VectorXd DampingCoeffientVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
            if (state->problem->is_time_dependent())
            {
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                Eigen::MatrixXd adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                AdjointTools::dJ_damping_transient_adjoint_term(*state, *diff_cache, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                log_and_throw_adjoint_error("[{}] Static simulation not supported!", name());
            }
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    Eigen::VectorXd DampingCoeffientVariableToSimulation::inverse_eval()
    {
        log_and_throw_adjoint_error("[{}] inverse_eval not implemented!", name());
        return Eigen::VectorXd();
    }
 
    void InitialConditionVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = *states[i];
            auto &diff_cache = *diff_caches[i];
 
            // For initial condition var2sim, the state variable is the initial solution and velocity
            // of the forward simulation. So DOF should be 2x the state.
            assert(state_variable.size() == 2 * state.ndof());
 
            diff_cache.initial_condition_override.solution = state_variable.head(state.ndof());
            diff_cache.initial_condition_override.velocity = state_variable.tail(state.ndof());
            diff_cache.initial_condition_override.acceleration.resize(0, 0);
        }
    }
    Eigen::VectorXd InitialConditionVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
            if (state->problem->is_time_dependent())
            {
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                Eigen::MatrixXd adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                AdjointTools::dJ_initial_condition_adjoint_term(*state, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                log_and_throw_adjoint_error("[{}] Static simulation not supported!", name());
            }
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    Eigen::VectorXd InitialConditionVariableToSimulation::inverse_eval()
    {
        auto &state = *states[0];
        Eigen::MatrixXd sol, vel;
        state.initial_solution(sol);
        state.initial_velocity(vel);
 
        Eigen::VectorXd x(sol.size() + vel.size());
        x << sol, vel;
        return x;
    }
 
    void DirichletVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(states[0]->problem);
        assert(dirichlet_boundaries_.size() > 0);
        int dim = states[0]->mesh->dimension();
        int num_steps = indices.size() / dim;
        for (int i = 0; i < num_steps; ++i)
            for (const int &b : dirichlet_boundaries_)
                tensor_problem->update_dirichlet_boundary(b, indices(i * dim) + 1, state_variable.segment(i * dim, dim));
 
        logger().info("Current dirichlet boundary {} is {}.", dirichlet_boundaries_[0], state_variable.transpose());
    }
 
    Eigen::VectorXd DirichletVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
            if (state->problem->is_time_dependent())
            {
                Eigen::MatrixXd adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                Eigen::MatrixXd adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                AdjointTools::dJ_dirichlet_transient_adjoint_term(*state, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                log_and_throw_adjoint_error("[{}] Static dirichlet boundary optimization not supported!", name());
            }
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    std::string DirichletVariableToSimulation::variable_to_string(const Eigen::VectorXd &variable)
    {
        return "";
    }
    Eigen::VectorXd DirichletVariableToSimulation::inverse_eval()
    {
        assert(dirichlet_boundaries_.size() > 0);
        assert(states.size() > 0);
 
        int dim = states[0]->mesh->dimension();
        Eigen::VectorXd x;
        for (const auto &b : states[0]->args["boundary_conditions"]["dirichlet_boundary"])
            if (b["id"].get<int>() == dirichlet_boundaries_[0])
            {
                auto value = b["value"];
                if (value[0].is_array())
                {
                    if (!states[0]->problem->is_time_dependent())
                        log_and_throw_adjoint_error("Simulation must be time dependent for timestep wise dirichlet.");
                    Eigen::MatrixXd dirichlet = value;
                    x.setZero(dirichlet.rows() * (dirichlet.cols() - 1));
                    for (int j = 1; j < dirichlet.cols(); ++j)
                        x.segment((j - 1) * dim, dim) = dirichlet.col(j);
                }
                else if (value[0].is_number())
                {
                    if (states[0]->problem->is_time_dependent())
                        log_and_throw_adjoint_error("Simulation must be quasistatic for single value dirichlet.");
                    x.resize(dim);
                    x = value;
                }
                else if (value.is_string())
                    assert(false);
                break;
            }
 
        return parametrization.inverse_eval(x);
    }
    void DirichletVariableToSimulation::set_output_indexing(const json &args)
    {
        const std::string composite_map_type = args["composite_map_type"];
        const State &state = *(states[0]);
        if (composite_map_type == "time_step_indexing")
        {
            const int time_steps = state.args["time"]["time_steps"];
            const int dim = state.mesh->dimension();
 
            output_indexing_.setZero(time_steps * dim);
            for (int i = 0; i < time_steps; ++i)
                for (int k = 0; k < dim; ++k)
                    output_indexing_(i * dim + k) = i;
        }
        else
            VariableToSimulation::set_output_indexing(args);
 
        set_dirichlet_boundaries(args["surface_selection"]);
    }
 
    void DirichletNodesVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        for (auto state : states)
        {
            assert(state_variable.size() == (state->mesh->dimension() * dirichlet_nodes_.size()));
            auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(state->problem);
            assert(!state->problem->is_time_dependent());
            tensor_problem->update_dirichlet_nodes(state->in_node_to_node, dirichlet_nodes_, utils::unflatten(state_variable, state->mesh->dimension()));
 
            logger().info("Updated dirichlet nodes");
        }
    }
 
    Eigen::VectorXd DirichletNodesVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
            if (state->problem->is_time_dependent())
                log_and_throw_adjoint_error("[{}] Transient dirichlet node optimization not supported!", name());
            else
                AdjointTools::dJ_dirichlet_static_adjoint_term(*state, *diff_cache, get_adjoint_mat(*state, *diff_cache, 0), cur_term);
 
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
    std::string DirichletNodesVariableToSimulation::variable_to_string(const Eigen::VectorXd &variable)
    {
        return "";
    }
    Eigen::VectorXd DirichletNodesVariableToSimulation::inverse_eval()
    {
        log_and_throw_adjoint_error("Inverse eval not implemented!");
    }
    void DirichletNodesVariableToSimulation::set_output_indexing(const json &args)
    {
        json args_ = args;
        const std::string composite_map_type = args_["composite_map_type"];
        if (composite_map_type != "indices")
        {
            log_and_throw_adjoint_error("Must specify indices on which the nodal dirichlet is applied!");
        }
        else
        {
            set_dirichlet_nodes(args_["composite_map_indices"]);
            int dim = 3;
            std::vector<int> composite_map_indices = {};
            for (int i = 0; i < dirichlet_nodes_.size(); ++i)
                for (int k = 0; k < dim; ++k)
                    composite_map_indices.push_back(dirichlet_nodes_[i] * dim + k);
            args_["composite_map_indices"] = composite_map_indices;
        }
        VariableToSimulation::set_output_indexing(args_);
    }
 
    void PressureVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
    {
        auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(states[0]->problem);
        assert(pressure_boundaries_.size() > 0);
        for (int i = 0; i < indices.size(); ++i)
            for (const int &b : pressure_boundaries_)
                tensor_problem->update_pressure_boundary(b, indices(i) + 1, state_variable(i));
 
        logger().info("Current pressure boundary {} is {}.", pressure_boundaries_[0], state_variable.transpose());
    }
 
    Eigen::VectorXd PressureVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
 
            if (state->problem->is_time_dependent())
            {
                Eigen::MatrixXd adjoint_nu, adjoint_p;
                adjoint_nu = get_adjoint_mat(*state, *diff_cache, 1);
                adjoint_p = get_adjoint_mat(*state, *diff_cache, 0);
                AdjointTools::dJ_pressure_transient_adjoint_term(*state, *diff_cache, pressure_boundaries_, adjoint_nu, adjoint_p, cur_term);
            }
            else
            {
                AdjointTools::dJ_pressure_static_adjoint_term(*state, pressure_boundaries_, diff_cache->u(0), get_adjoint_mat(*state, *diff_cache, 0), cur_term);
            }
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return apply_parametrization_jacobian(term, x);
    }
 
    std::string PressureVariableToSimulation::variable_to_string(const Eigen::VectorXd &variable)
    {
        return "";
    }
 
    Eigen::VectorXd PressureVariableToSimulation::inverse_eval()
    {
        assert(pressure_boundaries_.size() > 0);
        assert(states.size() > 0);
 
        Eigen::VectorXd x;
        for (const auto &b : states[0]->args["boundary_conditions"]["pressure_boundary"])
            if (b["id"].get<int>() == pressure_boundaries_[0])
            {
                auto value = b["value"];
                if (value.is_array())
                {
                    if (!states[0]->problem->is_time_dependent())
                        log_and_throw_adjoint_error("Simulation must be time dependent for timestep wise pressure.");
                    Eigen::VectorXd pressures = value;
                    x = pressures.segment(1, pressures.size() - 1);
                }
                else if (value.is_number())
                {
                    if (states[0]->problem->is_time_dependent())
                        log_and_throw_adjoint_error("Simulation must be quasistatic for single value pressure.");
                    x.resize(1);
                    x(0) = value;
                }
                else if (value.is_string())
                    assert(false);
                break;
            }
 
        return parametrization.inverse_eval(x);
    }
 
    void PressureVariableToSimulation::set_output_indexing(const json &args)
    {
        const std::string composite_map_type = args["composite_map_type"];
        const State &state = *(states[0]);
        if (composite_map_type == "time_step_indexing")
        {
            const int time_steps = state.args["time"]["time_steps"];
            output_indexing_.setZero(time_steps);
            for (int i = 0; i < time_steps; ++i)
                output_indexing_(i) = i;
        }
        else
            VariableToSimulation::set_output_indexing(args);
 
        set_pressure_boundaries(args["surface_selection"]);
    }
 
    Eigen::VectorXd PeriodicShapeVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
    {
        Eigen::VectorXd term, cur_term;
 
        for (int i = 0; i < states.size(); ++i)
        {
            auto &state = states[i];
            auto &diff_cache = diff_caches[i];
 
            if (state->problem->is_time_dependent())
            {
                log_and_throw_error("Not implemented!");
            }
            else
            {
                AdjointTools::dJ_periodic_shape_adjoint_term(*state, *diff_cache, *periodic_mesh_map, periodic_mesh_representation, diff_cache->u(0), get_adjoint_mat(*state, *diff_cache, 0), cur_term);
            }
            if (term.size() != cur_term.size())
                term = cur_term;
            else
                term += cur_term;
        }
        return VariableToSimulation::apply_parametrization_jacobian(term, x);
    }
    void PeriodicShapeVariableToSimulation::update(const Eigen::VectorXd &x)
    {
        const int dim = states[0]->mesh->dimension();
        periodic_mesh_representation = parametrization.eval(x);
        const Eigen::MatrixXd V = utils::unflatten(periodic_mesh_map->eval(periodic_mesh_representation), dim);
 
        for (auto state : states)
        {
            const int n_verts = state->mesh->n_vertices();
 
            for (int i = 0; i < n_verts; i++)
                state->mesh->set_point(i, V.row(i));
        }
    }
    Eigen::VectorXd PeriodicShapeVariableToSimulation::inverse_eval()
    {
        const auto &state = *(states[0]);
 
        Eigen::MatrixXd V;
        state.get_vertices(V);
 
        if (!state.periodic_bc->all_direction_periodic())
            log_and_throw_error("Cannot inverse evaluate periodic shape!");
 
        periodic_mesh_map = std::make_unique<PeriodicMeshToMesh>(V);
        periodic_mesh_representation = periodic_mesh_map->inverse_eval(utils::flatten(V));
 
        return parametrization.inverse_eval(periodic_mesh_representation);
    }
    Eigen::VectorXd PeriodicShapeVariableToSimulation::apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const
    {
        const Eigen::VectorXd mid = periodic_mesh_map->apply_jacobian(term, periodic_mesh_representation);
        return parametrization.apply_jacobian(mid, x);
    }
} // namespace polyfem::solver