StateHomogenization.cpp

Go to the documentation of this file.

#include <polyfem/legacy/State.hpp>
#include <polyfem/Common.hpp>
#include <polyfem/assembler/Mass.hpp>
#include <polyfem/assembler/ViscousDamping.hpp>
#include <polyfem/utils/Logger.hpp>
#include <polyfem/utils/StringUtils.hpp>
#include <polyfem/utils/MaybeParallelFor.hpp>
#include <polyfem/utils/Timer.hpp>
#include <polyfem/utils/Types.hpp>
#include <polyfem/solver/NLHomoProblem.hpp>
#include <polyfem/solver/forms/PeriodicContactForm.hpp>
#include <polyfem/solver/forms/lagrangian/MacroStrainLagrangianForm.hpp>
 
#include <polyfem/io/MshWriter.hpp>
#include <polyfem/io/OBJWriter.hpp>
#include <polyfem/io/Evaluator.hpp>
 
#include <Eigen/Core>
#include <unsupported/Eigen/SparseExtra>
#include <polysolve/linear/FEMSolver.hpp>
#include <polysolve/nonlinear/Solver.hpp>
#include <ipc/ipc.hpp>
#include <spdlog/fmt/fmt.h>
 
#include <algorithm>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace polyfem::legacy
{
 
    using namespace assembler;
    using namespace mesh;
    using namespace solver;
    using namespace utils;
    using namespace quadrature;
 
    void State::init_homogenization_solve(const double t)
    {
        const int dim = mesh->dimension();
        const int ndof = n_bases * dim;
 
        const std::vector<std::shared_ptr<Form>> forms = solve_data.init_forms(
            // General
            units,
            mesh->dimension(), t, in_node_to_node,
            // Elastic form
            n_bases, bases, geom_bases(), *assembler, ass_vals_cache, mass_ass_vals_cache,
            args["solver"]["advanced"]["jacobian_threshold"], args["solver"]["advanced"]["check_inversion"],
            args["solver"]["advanced"]["conservative_max_iter"],
            // Body form
            n_pressure_bases, boundary_nodes, local_boundary, local_neumann_boundary,
            n_boundary_samples(), rhs, Eigen::VectorXd::Zero(ndof) /* only to set neumann BC, not used*/, mass_matrix_assembler->density(),
            // Pressure form
            local_pressure_boundary, local_pressure_cavity, elasticity_pressure_assembler,
            // Inertia form
            args.value("/time/quasistatic"_json_pointer, true), mass,
            nullptr,
            // Lagged regularization form
            args["solver"]["advanced"]["lagged_regularization_weight"],
            args["solver"]["advanced"]["lagged_regularization_iterations"],
            // Augmented lagrangian form
            obstacle.ndof(), args["constraints"]["hard"], args["constraints"]["soft"], args["constraints"]["zero_mean"],
            // Contact form
            args["contact"]["enabled"], args["contact"]["periodic"].get<bool>() ? periodic_collision_mesh : collision_mesh, args["contact"]["dhat"],
            avg_mass, args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_area_weighting"]) : false,
            args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_improved_max_operator"]) : false,
            args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_physical_barrier"]) : false,
            args["solver"]["contact"]["barrier_stiffness"],
            args["solver"]["contact"]["initial_barrier_stiffness"],
            args["solver"]["contact"]["CCD"]["broad_phase"],
            args["solver"]["contact"]["CCD"]["tolerance"],
            args["solver"]["contact"]["CCD"]["max_iterations"],
            optimization_enabled,
            // Smooth Contact Form
            args["contact"]["use_gcp_formulation"],
            args["contact"]["alpha_t"],
            args["contact"]["alpha_n"],
            args["contact"]["use_adaptive_dhat"],
            args["contact"]["min_distance_ratio"],
            // Normal Adhesion Form
            is_adhesion_enabled(),
            args["contact"]["adhesion"]["dhat_p"],
            args["contact"]["adhesion"]["dhat_a"],
            args["contact"]["adhesion"]["adhesion_strength"],
            // Tangential Adhesion Form
            args["contact"]["adhesion"]["tangential_adhesion_coefficient"],
            args["contact"]["adhesion"]["epsa"],
            args["solver"]["contact"]["tangential_adhesion_iterations"],
            // Homogenization
            macro_strain_constraint,
            // Periodic contact
            args["contact"]["periodic"], periodic_collision_mesh_to_basis,
            // Friction form
            args["contact"]["friction_coefficient"],
            args["contact"]["epsv"],
            args["solver"]["contact"]["friction_iterations"],
            // Rayleigh damping form
            args["solver"]["rayleigh_damping"],
            // BC AL lumping
            args["solver"]["augmented_lagrangian"]["lumping"],
            // Boundary-ID periodic constraints
            mesh.get(), &total_local_boundary,
            args["boundary_conditions"]["periodic"], /*fe_space_id=*/-1);
 
        for (const auto &[name, form] : solve_data.named_forms())
        {
            if (name == "augmented_lagrangian")
            {
                form->set_weight(0);
                form->disable();
            }
        }
 
        bool solve_symmetric_flag = false;
        {
            const auto &fixed_entry = macro_strain_constraint.get_fixed_entry();
            for (int i = 0; i < dim; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    if (std::find(fixed_entry.data(), fixed_entry.data() + fixed_entry.size(), i + j * dim) == fixed_entry.data() + fixed_entry.size() && std::find(fixed_entry.data(), fixed_entry.data() + fixed_entry.size(), j + i * dim) == fixed_entry.data() + fixed_entry.size())
                    {
                        logger().info("Strain entry [{},{}] and [{},{}] are not fixed, solve for symmetric strain...", i, j, j, i);
                        solve_symmetric_flag = true;
                        break;
                    }
                }
                if (solve_symmetric_flag)
                    break;
            }
        }
 
        std::shared_ptr<solver::NLHomoProblem> homo_problem = std::make_shared<solver::NLHomoProblem>(
            ndof,
            macro_strain_constraint,
            n_bases, mesh_nodes, t, forms, solve_data.al_form, solve_symmetric_flag,
            polysolve::linear::Solver::create(args["solver"]["linear"], logger()),
            characteristic_length, characteristic_force_density, pure_mass, mesh->dimension());
        if (solve_data.periodic_contact_form)
            homo_problem->add_form(solve_data.periodic_contact_form);
        if (solve_data.strain_al_lagr_form)
            homo_problem->add_form(solve_data.strain_al_lagr_form);
 
        solve_data.nl_problem = homo_problem;
        solve_data.nl_problem->init(Eigen::VectorXd::Zero(homo_problem->reduced_size() + homo_problem->macro_reduced_size()));
        solve_data.nl_problem->update_quantities(t, Eigen::VectorXd::Zero(homo_problem->reduced_size() + homo_problem->macro_reduced_size()));
    }
 
    void State::solve_homogenization_step(int step, Eigen::MatrixXd &sol, bool adaptive_initial_weight, UserPostStepCallback user_post_step)
    {
        const int dim = mesh->dimension();
        const int ndof = n_bases * dim;
 
        auto homo_problem = std::dynamic_pointer_cast<solver::NLHomoProblem>(solve_data.nl_problem);
 
        Eigen::VectorXd extended_sol;
        extended_sol.setZero(ndof + dim * dim);
 
        if (sol.size() == extended_sol.size())
            extended_sol = sol;
 
        const auto &fixed_entry = macro_strain_constraint.get_fixed_entry();
        homo_problem->set_fixed_entry({});
        {
            std::shared_ptr<polysolve::nonlinear::Solver> nl_solver = make_nl_solver(true);
 
            Eigen::VectorXi al_indices = fixed_entry.array() + homo_problem->full_size();
            Eigen::VectorXd al_values = utils::flatten(macro_strain_constraint.eval(step))(fixed_entry);
 
            std::shared_ptr<MacroStrainLagrangianForm> lagr_form = solve_data.strain_al_lagr_form;
            lagr_form->enable();
 
            const double initial_weight = args["solver"]["augmented_lagrangian"]["initial_weight"];
            const double max_weight = args["solver"]["augmented_lagrangian"]["max_weight"];
            const double eta_tol = args["solver"]["augmented_lagrangian"]["eta"];
            const double scaling = args["solver"]["augmented_lagrangian"]["scaling"];
            double al_weight = initial_weight;
 
            Eigen::VectorXd tmp_sol = homo_problem->extended_to_reduced(extended_sol);
            const Eigen::VectorXd initial_sol = tmp_sol;
            const double initial_error = lagr_form->compute_error(extended_sol);
            double current_error = initial_error;
 
            // try to enforce fixed values on macro strain
            extended_sol(al_indices) = al_values;
            Eigen::VectorXd reduced_sol = homo_problem->extended_to_reduced(extended_sol);
 
            homo_problem->line_search_begin(tmp_sol, reduced_sol);
            int al_steps = 0;
            bool force_al = true;
 
            lagr_form->set_initial_weight(al_weight);
 
            while (force_al
                   || !std::isfinite(homo_problem->value(reduced_sol))
                   || !homo_problem->is_step_valid(tmp_sol, reduced_sol)
                   || !homo_problem->is_step_collision_free(tmp_sol, reduced_sol))
            {
                force_al = false;
                homo_problem->line_search_end();
 
                logger().info("Solving AL Problem with weight {}", al_weight);
 
                homo_problem->init(tmp_sol);
                try
                {
                    homo_problem->normalize_forms();
                    nl_solver->minimize(*homo_problem, tmp_sol);
                }
                catch (const std::runtime_error &e)
                {
                    logger().error("AL solve failed!");
                }
 
                extended_sol = homo_problem->reduced_to_extended(tmp_sol);
                logger().debug("Current macro strain: {}", extended_sol.tail(dim * dim));
 
                current_error = lagr_form->compute_error(extended_sol);
                const double eta = 1 - sqrt(current_error / initial_error);
 
                logger().info("Current eta = {}, current error = {}, initial error = {}", eta, current_error, initial_error);
 
                if (eta < eta_tol && al_weight < max_weight)
                    al_weight *= scaling;
                else
                    lagr_form->update_lagrangian(extended_sol, al_weight);
 
                if (eta <= 0)
                {
                    if (adaptive_initial_weight)
                    {
                        args["solver"]["augmented_lagrangian"]["initial_weight"] = args["solver"]["augmented_lagrangian"]["initial_weight"].get<double>() * scaling;
                        {
                            json tmp = json::object();
                            tmp["/solver/augmented_lagrangian/initial_weight"_json_pointer] = args["solver"]["augmented_lagrangian"]["initial_weight"];
                        }
                        logger().warn("AL weight too small, increase weight and revert solution, new initial weight is {}", args["solver"]["augmented_lagrangian"]["initial_weight"].get<double>());
                    }
                    tmp_sol = initial_sol;
                }
 
                // try to enforce fixed values on macro strain
                extended_sol(al_indices) = al_values;
                reduced_sol = homo_problem->extended_to_reduced(extended_sol);
 
                homo_problem->line_search_begin(tmp_sol, reduced_sol);
            }
            homo_problem->line_search_end();
            lagr_form->disable();
        }
 
        homo_problem->set_fixed_entry(fixed_entry);
 
        Eigen::VectorXd reduced_sol = homo_problem->extended_to_reduced(extended_sol);
 
        homo_problem->init(reduced_sol);
        std::shared_ptr<polysolve::nonlinear::Solver> nl_solver = make_nl_solver(false);
        homo_problem->normalize_forms();
        nl_solver->minimize(*homo_problem, reduced_sol);
 
        logger().info("Macro Strain: {}", extended_sol.tail(dim * dim).transpose());
 
        // check saddle point
        {
            json linear_args = args["solver"]["linear"];
            std::string solver_name = linear_args["solver"];
            if (solver_name.find("Pardiso") != std::string::npos)
            {
                linear_args["solver"] = "Eigen::PardisoLLT";
                std::unique_ptr<polysolve::linear::Solver> solver =
                    polysolve::linear::Solver::create(linear_args, logger());
 
                StiffnessMatrix A;
                homo_problem->hessian(reduced_sol, A);
                Eigen::VectorXd x, b = Eigen::VectorXd::Zero(A.rows());
                try
                {
                    dirichlet_solve(
                        *solver, A, b, {}, x, A.rows(), args["output"]["data"]["stiffness_mat"], false, false, false);
                }
                catch (const std::runtime_error &error)
                {
                    logger().error("The solution is a saddle point!");
                }
            }
        }
 
        sol = homo_problem->reduced_to_extended(reduced_sol);
 
        if (user_post_step)
        {
            Eigen::MatrixXd disp_grad = utils::unflatten(sol.bottomRows(dim * dim), dim);
            user_post_step(step, *this, homo_problem->reduced_to_full(reduced_sol), &disp_grad, nullptr);
        }
    }
 
    void State::solve_homogenization(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, UserPostStepCallback user_post_step)
    {
        bool is_static = !is_param_valid(args, "time");
        if (!is_static && !args["time"]["quasistatic"])
            log_and_throw_error("Transient homogenization can only do quasi-static!");
 
        init_homogenization_solve(t0);
 
        const int t_offset = args["output"]["data"]["file_index_offset"].get<int>();
        const int dim = mesh->dimension();
        Eigen::MatrixXd extended_sol;
        for (int t = 0; t <= time_steps; ++t)
        {
            double forward_solve_time = 0, remeshing_time = 0, global_relaxation_time = 0;
 
            {
                POLYFEM_SCOPED_TIMER(forward_solve_time);
                solve_homogenization_step(t, extended_sol, false, user_post_step);
            }
            sol = extended_sol.topRows(extended_sol.size() - dim * dim) + polyfem::io::Evaluator::generate_linear_field(n_bases, mesh_nodes, utils::unflatten(extended_sol.bottomRows(dim * dim), dim));
 
            if (is_static)
                return;
 
            // Always save the solution for consistency
            save_timestep(t0 + dt * t, t + t_offset, t0, dt, sol, Eigen::MatrixXd()); // no pressure
 
            {
                POLYFEM_SCOPED_TIMER("Update quantities");
 
                //     solve_data.time_integrator->update_quantities(sol);
 
                solve_data.nl_problem->update_quantities(t0 + (t + 1) * dt, sol);
 
                solve_data.update_dt();
                solve_data.update_barrier_stiffness(sol);
            }
 
            logger().info("{}/{}  t={}", t, time_steps, t0 + dt * t);
 
            // const std::string rest_mesh_path = args["output"]["data"]["rest_mesh"].get<std::string>();
            // if (!rest_mesh_path.empty())
            // {
            //     Eigen::MatrixXd V;
            //     Eigen::MatrixXi F;
            //     build_mesh_matrices(V, F);
            //     polyfem::io::MshWriter::write(
            //         resolve_output_path(fmt::format(args["output"]["data"]["rest_mesh"], t)),
            //         V, F, mesh->get_body_ids(), mesh->is_volume(), /*binary=*/true);
            // }
 
            // const std::string &state_path = resolve_output_path(fmt::format(args["output"]["data"]["state"], t));
            // if (!state_path.empty())
            //     solve_data.time_integrator->save_state(state_path);
 
            // save restart file
            save_restart_json(t0, dt, t);
            // stats_csv.write(t, forward_solve_time, remeshing_time, global_relaxation_time, sol);
        }
    }
 
} // namespace polyfem::legacy