VarForm.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/assembler/Assembler.hpp>
#include <polyfem/assembler/Problem.hpp>
#include <polyfem/assembler/RhsAssembler.hpp>
#include <polyfem/basis/ElementBases.hpp>
#include <polyfem/basis/InterfaceData.hpp>
#include <polyfem/mesh/Mesh.hpp>
#include <polyfem/mesh/MeshNodes.hpp>
#include <polyfem/mesh/GeometryLoader.hpp>
#include <polyfem/io/OutputData.hpp>
#include <polyfem/io/OutData.hpp>
#include <polyfem/io/OutStatsData.hpp>
#include <polyfem/io/ResourceIO.hpp>
#include <polyfem/io/Checkpoint.hpp>
#include <polyfem/utils/Types.hpp>
#include <polyfem/varforms/FESpace.hpp>
 
#include <Eigen/Dense>
 
#include <functional>
#include <cassert>
#include <iosfwd>
#include <memory>
#include <optional>
#include <map>
#include <vector>
#include <string>
 
namespace polyfem
{
    namespace time_integrator
    {
        class BDF;
        class ImplicitTimeIntegrator;
    } // namespace time_integrator
 
    namespace test
    {
        class VarFormTestAccess;
    }
 
    namespace varform
    {
        struct InitialConditionOverride
        {
            Eigen::MatrixXd solution;
            Eigen::MatrixXd velocity;
            Eigen::MatrixXd acceleration;
            Eigen::MatrixXd displacement_gradient;
        };
 
        // Called after each completed forward-simulation step.
        using ForwardStepCallback = std::function<void(int step, const Eigen::MatrixXd &solution)>;
 
        class VarForm
        {
            friend class polyfem::test::VarFormTestAccess;
 
        public:
            explicit VarForm(const io::ResourceIO &resources)
                : resources_(resources) {}
            virtual ~VarForm() = default;
 
            virtual std::string name() const = 0;
 
            void set_args(const json &args)
            {
                this->args = args;
                prepared_ = false;
                output_sampler_initialized_ = false;
            }
 
            virtual void init(
                const std::string &formulation,
                const Units &units,
                const json &args,
                const std::string &out_path);
 
            void set_checkpoint_reader(const io::CheckpointReader &reader) { checkpoint_reader_ = std::cref(reader); }
 
            void set_mesh(std::unique_ptr<mesh::Mesh> mesh, const double loading_mesh_time = 0);
            void set_geometry(mesh::LoadedGeometry geometry, double loading_mesh_time = 0);
 
            void solve(
                Eigen::MatrixXd &sol,
                const InitialConditionOverride *initial_condition_override = nullptr,
                const ForwardStepCallback &post_step = {});
 
            void set_time_callback(const std::function<void(int, int, double, double)> &callback) { time_callback = callback; }
 
            int problem_dimension() const;
            virtual bool is_contact_enabled() const { return false; }
 
            virtual io::OutputSpace output_space() const = 0;
            virtual std::vector<io::OutputField> output_fields(
                const io::OutputSample &sample,
                const Eigen::MatrixXd &solution,
                const io::OutputFieldOptions &options) const = 0;
 
            const io::OutRuntimeData &output_timings() const { return timings; }
            virtual io::OutStatsData compute_errors(const Eigen::MatrixXd &solution) = 0;
 
            virtual void save_json(const Eigen::MatrixXd &solution, std::ostream &out) const = 0;
            void save_json(const Eigen::MatrixXd &solution) const;
            virtual void export_data(const Eigen::MatrixXd &solution) const = 0;
 
            virtual void serialize_checkpoint(
                io::CheckpointWriter &writer,
                const Eigen::MatrixXd &solution,
                const io::CheckpointMetadata &metadata) const;
            virtual void deserialize_checkpoint(
                const io::CheckpointReader &reader,
                Eigen::MatrixXd &solution);
 
        protected:
            void prepare();
 
            std::string resolve_output_path(const std::string &path) const;
 
            void set_materials(assembler::Assembler &assembler, const int size) const;
            void restore_checkpoint_integrator(
                const std::shared_ptr<time_integrator::ImplicitTimeIntegrator> &integrator,
                const std::string &group,
                double dt,
                const std::string &ordering,
                const FESpace &space,
                int block_size) const;
            void write_checkpoint_ordering(
                io::CheckpointWriter &writer,
                const std::string &name,
                const FESpace &space) const;
            void reorder_checkpoint_block(
                const io::CheckpointReader &reader,
                const std::string &name,
                const FESpace &space,
                int block_size,
                int row_offset,
                Eigen::MatrixXd &value) const;
            Eigen::VectorXi checkpoint_ordering(const FESpace &space) const;
            void validate_checkpoint_solution(const Eigen::MatrixXd &solution, int expected_rows) const;
            bool checkpoint_reorder_enabled() const;
            virtual void reset() = 0;
 
            virtual void load_mesh(const mesh::Mesh &mesh, const json &args) = 0;
            virtual void build_basis(mesh::Mesh &mesh, const bool iso_parametric, const json &args) = 0;
            void assign_discr_orders(const json &space_args, const mesh::Mesh &mesh, Eigen::VectorXi &disc_orders, Eigen::VectorXi &disc_ordersq);
            void assign_discr_orders(const json &space_args, const int fe_space_id, const mesh::Mesh &mesh, Eigen::VectorXi &disc_orders, Eigen::VectorXi &disc_ordersq);
            virtual void assemble_rhs(const mesh::Mesh &mesh) = 0;
            virtual void assemble_mass_mat(const mesh::Mesh &mesh, const json &args) = 0;
            virtual void solve_problem(
                Eigen::MatrixXd &sol,
                const InitialConditionOverride *initial_condition_override,
                const ForwardStepCallback &post_step) = 0;
            QuadratureOrders n_boundary_samples(const int discr_order, const int discr_orderq, const int gdiscr_order) const;
 
            void build_fe_space(
                mesh::Mesh &mesh,
                const bool iso_parametric,
                const Eigen::VectorXi &disc_orders,
                const Eigen::VectorXi &disc_ordersq,
                const std::string &basis_type,
                const std::string &poly_basis_type,
                const assembler::Assembler &space_assembler,
                const int value_dim,
                const int quadrature_order,
                const int mass_quadrature_order,
                const bool use_corner_quadrature,
                const int n_harmonic_samples,
                const int integral_constraints,
                FESpace &space,
                VarFormBoundaryState &boundary,
                std::shared_ptr<GeometryMapping> geometry = nullptr);
 
        private:
            void build_polygonal_basis(
                const mesh::Mesh &mesh,
                const std::string &poly_basis_type,
                const assembler::Assembler &space_assembler,
                bool iso_parametric,
                const int quadrature_order,
                const int mass_quadrature_order,
                const int n_harmonic_samples,
                const int integral_constraints,
                FESpace &space,
                VarFormBoundaryState &boundary);
 
            void build_node_mapping(
                const mesh::Mesh &mesh,
                const std::string &basis_type,
                const FESpace &space,
                Eigen::VectorXi &space_in_node_to_node,
                Eigen::VectorXi &space_in_primitive_to_primitive) const;
 
        protected:
            virtual void build_rhs_assembler() = 0;
            virtual void set_obstacle(mesh::Obstacle &&obstacle);
 
            void save_step_state(
                const double t0,
                const double dt,
                const int t,
                const Eigen::MatrixXd &solution,
                const time_integrator::ImplicitTimeIntegrator *time_integrator) const;
 
            void ensure_output_sampler() const;
            void save_timestep(const double time, const int t, const double t0, const double dt, const Eigen::MatrixXd &solution) const;
            bool save_timestep_to_vtm(
                const double time, const int t, const double dt,
                const Eigen::MatrixXd &solution, paraviewo::VTMWriter &vtm,
                const std::string &block_prefix) const;
            void save_subsolve(const int i, const int t, const Eigen::MatrixXd &solution) const;
            int output_file_index(const int t) const;
            void notify_time_step(const int t, const int time_steps, const double t0, const double dt) const;
            void init_child_varform(
                VarForm &child,
                const std::string &formulation,
                const Units &units,
                const json &args,
                const std::string &out_path) const;
 
            io::OutGeometryData::ExportOptions export_options(const io::OutputSpace &space) const;
            io::OutputFieldFunction output_field_function(const Eigen::MatrixXd &solution, const io::OutGeometryData::ExportOptions &opts) const;
 
            std::shared_ptr<assembler::Problem> problem;
            Units units;
            json args;
 
            io::OutStatsData stats;
 
            io::OutRuntimeData timings;
 
            std::string output_path;
            std::optional<std::reference_wrapper<const io::CheckpointReader>> checkpoint_reader_;
 
            std::unique_ptr<mesh::Mesh> mesh_;
 
            std::function<void(int, int, double, double)> time_callback;
 
            mutable io::OutGeometryData output_geometry_;
            mutable bool output_sampler_initialized_ = false;
            bool prepared_ = false;
            int output_index_offset_ = 0;
 
        protected:
            const io::ResourceIO &resources_;
 
            static void rebuild_node_positions(
                const std::vector<basis::ElementBases> &bases,
                const std::vector<int> &node_ids,
                std::vector<RowVectorNd> &positions);
        };
    } // namespace varform
} // namespace polyfem