Problem.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/Units.hpp>
 
#include <polyfem/assembler/Assembler.hpp>
#include <polyfem/basis/ElementBases.hpp>
#include <polyfem/mesh/LocalBoundary.hpp>
#include <polyfem/mesh/Mesh.hpp>
 
namespace polyfem
{
    namespace assembler
    {
        enum class BoundaryKind
        {
            Dirichlet,
            Neumann,
            NormalAlignedNeumann, // Follower traction from pressure_boundary.
            PressureCavity
        };
 
        class Problem
        {
        public:
            Problem(const std::string &name);
            virtual void set_units(const assembler::Assembler &assembler, const Units &units) {}
 
            virtual ~Problem() {}
 
            virtual void init(const mesh::Mesh &mesh) {}
            inline const std::string &name() const { return name_; }
 
            virtual bool is_scalar() const = 0;
 
            virtual void rhs(const assembler::Assembler &assembler, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val, const int fe_space_id = -1) const = 0;
            virtual void rhs(const assembler::Assembler &assembler, const mesh::Mesh &mesh, const int element_id, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val, const int fe_space_id = -1) const
            {
                (void)mesh;
                (void)element_id;
                rhs(assembler, pts, t, val, fe_space_id);
            }
            virtual bool is_rhs_zero(const int fe_space_id = -1) const = 0;
 
            virtual void dirichlet_bc(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val, const int fe_space_id = -1) const = 0;
            virtual void neumann_bc(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, const double t, Eigen::MatrixXd &val, const int fe_space_id = -1) const {}
            virtual void pressure_bc(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, const double t, Eigen::MatrixXd &val) const {}
            virtual double pressure_cavity_bc(const int boundary_id, const double t) const { return 0; }
 
            virtual bool is_boundary_pressure(const int boundary_id) const { return std::find(pressure_boundary_ids_.begin(), pressure_boundary_ids_.end(), boundary_id) != pressure_boundary_ids_.end(); }
 
            virtual void dirichlet_nodal_value(const mesh::Mesh &mesh, const int node_id, const RowVectorNd &pt, const double t, Eigen::MatrixXd &val, const int fe_space_id = -1) const {}
            virtual void neumann_nodal_value(const mesh::Mesh &mesh, const int node_id, const RowVectorNd &pt, const Eigen::MatrixXd &normal, const double t, Eigen::MatrixXd &val, const int fe_space_id = -1) const {}
            virtual bool is_nodal_dirichlet_boundary(const int n_id, const int tag, const int fe_space_id = -1) { return false; }
            virtual bool is_nodal_neumann_boundary(const int n_id, const int tag, const int fe_space_id = -1) { return false; }
            virtual bool has_nodal_dirichlet(const int fe_space_id = -1) { return false; }
            virtual bool has_nodal_neumann(const int fe_space_id = -1) { return false; }
 
            virtual bool has_exact_sol() const = 0;
            virtual void exact(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const {};
            virtual void exact_grad(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const {};
 
            virtual void clear() {}
 
            virtual bool is_time_dependent() const { return false; }
            virtual bool is_constant_in_time() const { return true; }
 
            virtual void initial_solution(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val, const int fe_space_id = -1) const {}
            virtual void initial_velocity(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val, const int fe_space_id = -1) const {}
            virtual void initial_acceleration(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val, const int fe_space_id = -1) const {}
            virtual void initial_density(const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const {}
 
            virtual void set_parameters(const json &params, const io::ResourceIO &resources) {}
 
            virtual bool might_have_no_dirichlet() { return false; }
            virtual bool is_dimension_dirichet(const int tag, const int dim, const int fe_space_id = -1) const { return true; }
            virtual bool is_nodal_dimension_dirichlet(const int n_id, const int tag, const int dim, const int fe_space_id = -1) const { return true; }
 
            virtual bool all_dimensions_dirichlet(const int fe_space_id) const
            {
                (void)fe_space_id;
                return true;
            } // here for efficiency reasons
 
            // Selects local boundary primitives for a BC kind and FE space.
            void setup_bc(
                const mesh::Mesh &mesh,
                const BoundaryKind kind,
                const int fe_space_id,
                const std::vector<mesh::LocalBoundary> &local_boundary,
                std::vector<mesh::LocalBoundary> &selected_local_boundary);
 
            // Also collects constrained DOF ids for Dirichlet conditions.
            void setup_bc(
                const mesh::Mesh &mesh,
                const BoundaryKind kind,
                const int fe_space_id,
                const std::vector<basis::ElementBases> &bases,
                const std::vector<mesh::LocalBoundary> &local_boundary,
                std::vector<mesh::LocalBoundary> &selected_local_boundary,
                std::vector<int> &boundary_nodes,
                const int value_dim = -1);
 
            // Groups selected primitives by boundary tag, e.g. for cavity pressure.
            void setup_bc(
                const mesh::Mesh &mesh,
                const BoundaryKind kind,
                const int fe_space_id,
                const std::vector<mesh::LocalBoundary> &local_boundary,
                std::unordered_map<int, std::vector<mesh::LocalBoundary>> &selected_by_tag);
 
            // Returns nodal BC node ids, not constrained DOF ids.
            void setup_nodal_bc(
                const mesh::Mesh &mesh,
                const BoundaryKind kind,
                const int fe_space_id,
                const int n_bases,
                std::vector<int> &nodes);
 
            virtual void update_nodes(const Eigen::VectorXi &in_node_to_node) {}
 
        protected:
            virtual bool has_boundary(const BoundaryKind kind, const int tag, const int fe_space_id);
 
            std::vector<int> boundary_ids_;
            std::vector<int> neumann_boundary_ids_;
            std::vector<int> normal_aligned_neumann_boundary_ids_;
            std::vector<int> pressure_boundary_ids_;
            std::vector<int> pressure_cavity_ids_;
 
            bool updated_dirichlet_node_ordering_ = false;
 
            std::string name_;
        };
    } // namespace assembler
} // namespace polyfem