OutData.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/Common.hpp>
 
#include <polyfem/assembler/Problem.hpp>
 
#include <polyfem/basis/ElementBases.hpp>
 
#include <polyfem/mesh/Mesh.hpp>
 
#include <polyfem/solver/SolveData.hpp>
 
#include <paraviewo/ParaviewWriter.hpp>
#include <paraviewo/VTUWriter.hpp>
#include <paraviewo/HDF5VTUWriter.hpp>
 
#include <polyfem/utils/RefElementSampler.hpp>
 
#include <Eigen/Dense>
 
namespace polyfem
{
    class State;
}
 
namespace polyfem::io
{
    class SolutionFrame
    {
    public:
        std::string name;
        Eigen::MatrixXd points;
        Eigen::MatrixXi connectivity;
        Eigen::MatrixXd solution;
        Eigen::MatrixXd pressure;
        Eigen::MatrixXd exact;
        Eigen::MatrixXd error;
        Eigen::MatrixXd scalar_value;
        Eigen::MatrixXd scalar_value_avg;
    };
 
    class OutGeometryData
    {
    public:
        struct ExportOptions
        {
            bool volume;
            bool surface;
            bool wire;
            bool points;
            bool contact_forces;
            bool friction_forces;
            bool forces;
            bool jacobian_validity;
 
            bool use_sampler;
            bool boundary_only;
            bool material_params;
            bool body_ids;
            bool sol_on_grid;
            bool velocity;
            bool acceleration;
            bool tensor_values;
            bool scalar_values;
            bool discretization_order;
            bool nodes;
 
            bool use_spline;
            bool reorder_output;
 
            bool solve_export_to_file;
 
            bool use_hdf5;
 
            ExportOptions(const json &args,
                          const bool is_mesh_linear,
                          const bool is_problem_scalar,
                          const bool solve_export_to_file);
 
            inline std::string file_extension() const { return use_hdf5 ? ".hdf" : ".vtu"; }
        };
 
        static void extract_boundary_mesh(
            const mesh::Mesh &mesh,
            const int n_bases,
            const std::vector<basis::ElementBases> &bases,
            const std::vector<mesh::LocalBoundary> &total_local_boundary,
            Eigen::MatrixXd &node_positions,
            Eigen::MatrixXi &boundary_edges,
            Eigen::MatrixXi &boundary_triangles,
            std::vector<Eigen::Triplet<double>> &displacement_map_entries);
 
        void init_sampler(const polyfem::mesh::Mesh &mesh, const double vismesh_rel_area);
 
        void build_grid(const polyfem::mesh::Mesh &mesh, const double spacing);
 
        void export_data(
            const State &state,
            const Eigen::MatrixXd &sol,
            const Eigen::MatrixXd &pressure,
            const bool is_time_dependent,
            const double tend_in,
            const double dt,
            const ExportOptions &opts,
            const std::string &vis_mesh_path,
            const std::string &nodes_path,
            const std::string &solution_path,
            const std::string &stress_path,
            const std::string &mises_path,
            const bool is_contact_enabled,
            std::vector<SolutionFrame> &solution_frames) const;
 
        void save_vtu(const std::string &path,
                      const State &state,
                      const Eigen::MatrixXd &sol,
                      const Eigen::MatrixXd &pressure,
                      const double t,
                      const double dt,
                      const ExportOptions &opts,
                      const bool is_contact_enabled,
                      std::vector<SolutionFrame> &solution_frames) const;
 
        void save_volume(const std::string &path,
                         const State &state,
                         const Eigen::MatrixXd &sol,
                         const Eigen::MatrixXd &pressure,
                         const double t,
                         const double dt,
                         const ExportOptions &opts,
                         std::vector<SolutionFrame> &solution_frames) const;
 
        void save_surface(const std::string &export_surface,
                          const State &state,
                          const Eigen::MatrixXd &sol,
                          const Eigen::MatrixXd &pressure,
                          const double t,
                          const double dt_in,
                          const ExportOptions &opts,
                          const bool is_contact_enabled,
                          std::vector<SolutionFrame> &solution_frames) const;
 
        void save_contact_surface(
            const std::string &export_surface,
            const State &state,
            const Eigen::MatrixXd &sol,
            const Eigen::MatrixXd &pressure,
            const double t,
            const double dt_in,
            const ExportOptions &opts,
            const bool is_contact_enabled,
            std::vector<SolutionFrame> &solution_frames) const;
 
        void save_wire(const std::string &name,
                       const State &state,
                       const Eigen::MatrixXd &sol,
                       const double t,
                       const ExportOptions &opts,
                       std::vector<SolutionFrame> &solution_frames) const;
 
        void save_points(
            const std::string &path,
            const State &state,
            const Eigen::MatrixXd &sol,
            const ExportOptions &opts,
            std::vector<SolutionFrame> &solution_frames) const;
 
        void save_pvd(const std::string &name, const std::function<std::string(int)> &vtu_names,
                      int time_steps, double t0, double dt, int skip_frame = 1) const;
 
    private:
        utils::RefElementSampler ref_element_sampler;
 
        Eigen::MatrixXd grid_points;
        Eigen::MatrixXi grid_points_to_elements;
        Eigen::MatrixXd grid_points_bc;
 
        void build_vis_boundary_mesh(
            const mesh::Mesh &mesh,
            const std::vector<basis::ElementBases> &bases,
            const std::vector<basis::ElementBases> &gbases,
            const std::vector<mesh::LocalBoundary> &total_local_boundary,
            const Eigen::MatrixXd &solution,
            const int problem_dim,
            Eigen::MatrixXd &boundary_vis_vertices,
            Eigen::MatrixXd &boundary_vis_local_vertices,
            Eigen::MatrixXi &boundary_vis_elements,
            Eigen::MatrixXi &boundary_vis_elements_ids,
            Eigen::MatrixXi &boundary_vis_primitive_ids,
            Eigen::MatrixXd &boundary_vis_normals,
            Eigen::MatrixXd &displaced_boundary_vis_normals) const;
 
        void build_vis_mesh(
            const mesh::Mesh &mesh,
            const Eigen::VectorXi &disc_orders,
            const std::vector<basis::ElementBases> &gbases,
            const std::map<int, Eigen::MatrixXd> &polys,
            const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
            const bool boundary_only,
            Eigen::MatrixXd &points,
            Eigen::MatrixXi &tets,
            Eigen::MatrixXi &el_id,
            Eigen::MatrixXd &discr) const;
 
        void build_high_order_vis_mesh(
            const mesh::Mesh &mesh,
            const Eigen::VectorXi &disc_orders,
            const std::vector<basis::ElementBases> &bases,
            Eigen::MatrixXd &points,
            std::vector<std::vector<int>> &elements,
            Eigen::MatrixXi &el_id,
            Eigen::MatrixXd &discr) const;
 
        void save_volume_vector_field(
            const State &state,
            const Eigen::MatrixXd &points,
            const ExportOptions &opts,
            const std::string &name,
            const Eigen::VectorXd &field,
            paraviewo::ParaviewWriter &writer) const;
    };
 
    class OutRuntimeData
    {
    public:
        double building_basis_time;
        double loading_mesh_time;
        double computing_poly_basis_time;
        double assembling_stiffness_mat_time;
        double assembling_mass_mat_time;
        double assigning_rhs_time;
        double solving_time;
 
        double total_time()
        {
            return building_basis_time + assembling_mass_mat_time + assembling_stiffness_mat_time + solving_time;
        }
    };
 
    class OutStatsData
    {
    public:
        Eigen::Vector4d spectrum;
 
        json solver_info;
 
        double mesh_size;
        double min_edge_length;
        double average_edge_length;
 
        double l2_err, linf_err, lp_err, h1_err, h1_semi_err, grad_max_err;
 
        long long nn_zero, mat_size, num_dofs;
 
        double max_angle;
        double sigma_max, sigma_min, sigma_avg;
 
        int n_flipped;
 
        int simplex_count;
        int regular_count;
        int regular_boundary_count;
        int simple_singular_count;
        int multi_singular_count;
        int boundary_count;
        int non_regular_boundary_count;
        int non_regular_count;
        int undefined_count;
        int multi_singular_boundary_count;
 
        void compute_errors(const int n_bases,
                            const std::vector<polyfem::basis::ElementBases> &bases,
                            const std::vector<polyfem::basis::ElementBases> &gbases,
                            const polyfem::mesh::Mesh &mesh,
                            const assembler::Problem &problem,
                            const double tend,
                            const Eigen::MatrixXd &sol);
 
        void compute_mesh_stats(const polyfem::mesh::Mesh &mesh);
 
        void compute_mesh_size(const polyfem::mesh::Mesh &mesh_in, const std::vector<polyfem::basis::ElementBases> &bases_in, const int n_samples, const bool use_curved_mesh_size);
 
        void reset();
 
        void count_flipped_elements(const polyfem::mesh::Mesh &mesh, const std::vector<polyfem::basis::ElementBases> &gbases);
 
 
        void save_json(const nlohmann::json &args,
                       const int n_bases, const int n_pressure_bases,
                       const Eigen::MatrixXd &sol,
                       const mesh::Mesh &mesh,
                       const Eigen::VectorXi &disc_orders,
                       const assembler::Problem &problem,
                       const OutRuntimeData &runtime,
                       const std::string &formulation,
                       const bool isoparametric,
                       const int sol_at_node_id,
                       nlohmann::json &j);
    };
 
    class EnergyCSVWriter
    {
    public:
        EnergyCSVWriter(const std::string &path, const solver::SolveData &solve_data);
        ~EnergyCSVWriter();
 
        void write(const int i, const Eigen::MatrixXd &sol);
 
    protected:
        std::ofstream file;
        const solver::SolveData &solve_data;
    };
 
    class RuntimeStatsCSVWriter
    {
    public:
        RuntimeStatsCSVWriter(const std::string &path, const State &state, const double t0, const double dt);
        ~RuntimeStatsCSVWriter();
 
        void write(const int t, const double forward, const double remeshing, const double global_relaxation, const Eigen::MatrixXd &sol);
 
    protected:
        std::ofstream file;
        const State &state;
        const double t0;
        const double dt;
        double total_forward_solve_time = 0;
        double total_remeshing_time = 0;
        double total_global_relaxation_time = 0;
    };
} // namespace polyfem::io