CFiniteElementProblem.h

/* +---------------------------------------------------------------------------+
   |              OpenBeam - C++ Finite Element Analysis library               |
   |                                                                           |
   |   Copyright (C) 2010-2021  Jose Luis Blanco Claraco                       |
   |                              University of Malaga                         |
   |                                                                           |
   | OpenBeam is free software: you can redistribute it and/or modify          |
   |     it under the terms of the GNU General Public License as published by  |
   |     the Free Software Foundation, either version 3 of the License, or     |
   |     (at your option) any later version.                                   |
   |                                                                           |
   | OpenBeam is distributed in the hope that it will be useful,               |
   |     but WITHOUT ANY WARRANTY; without even the implied warranty of        |
   |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         |
   |     GNU General Public License for more details.                          |
   |                                                                           |
   |     You should have received a copy of the GNU General Public License     |
   |     along with OpenBeam.  If not, see <http://www.gnu.org/licenses/>.     |
   |                                                                           |
   +---------------------------------------------------------------------------+
 */
 
#pragma once
 
#include <mrpt/containers/yaml.h>
#include <mrpt/core/optional_ref.h>
#include <mrpt/opengl/CSetOfObjects.h>
#include <mrpt/system/CTimeLogger.h>
#include <openbeam/CElement.h>
#include <openbeam/types.h>
 
#include <Eigen/Sparse>
#include <cstdint>
#include <fstream>
#include <iostream>
 
namespace openbeam
{
struct RenderInitData;  // Fwd. decl. (defined in
                        // CFiniteElementProblem::saveAsImage)
 
enum class DoF_index : uint8_t
{
    DX = 0,
    DY,
    DZ,
    RX,
    RY,
    RZ
};
 
struct NodeDoF
{
    NodeDoF(node_index_t _node_id, DoF_index _dof) : nodeId(_node_id), dof(_dof)
    {
    }
    NodeDoF(node_index_t _node_id, uint8_t _dof)
        : nodeId(_node_id), dof(static_cast<DoF_index>(_dof))
    {
    }
 
    node_index_t nodeId;
    DoF_index    dof;  
 
    uint8_t dofAsInt() const { return static_cast<uint8_t>(dof); }
};
 
enum class StaticSolverAlgorithm : uint8_t
{
    Dense_LLT = 0,
    Sparse_LLT
};
 
struct BuildProblemInfo
{
    struct TDoFType
    {
        size_t bounded_index;  
        size_t free_index;  
    };
 
    Eigen::SparseMatrix<num_t> K_bb;
    Eigen::SparseMatrix<num_t> K_ff;
    Eigen::SparseMatrix<num_t> K_bf;
 
    std::vector<size_t>   free_dof_indices;
    std::vector<size_t>   bounded_dof_indices;
    std::vector<TDoFType> dof_types;
 
    Eigen::Matrix<num_t, Eigen::Dynamic, 1> U_b;
 
    Eigen::Matrix<num_t, Eigen::Dynamic, 1> F_f;
};
 
struct StaticSolveProblemInfo
{
    BuildProblemInfo build_info;
 
    Eigen::Matrix<num_t, Eigen::Dynamic, 1> F_b;
 
    Eigen::Matrix<num_t, Eigen::Dynamic, 1> U_f;
 
    Eigen::Matrix<num_t, Eigen::Dynamic, 1> F;
 
    Eigen::Matrix<num_t, Eigen::Dynamic, 1> U;
};
 
struct MeshParams
{
    MeshParams();
 
    double max_element_length;  
};
 
struct MeshOutputInfo
{
    MeshOutputInfo() = default;
 
    size_t num_original_nodes = 0;
 
    size_t num_original_elements = 0;
 
    std::deque<std::vector<node_index_t>> element2nodes;
 
    std::deque<std::vector<element_index_t>> element2elements;
};
 
struct ImageSaveOutputInfo
{
    ImageSaveOutputInfo() = default;
 
    unsigned int img_width = 0, img_height = 0;
};
 
struct StaticSolverOptions
{
    StaticSolverOptions() = default;
 
    StaticSolverAlgorithm algorithm          = StaticSolverAlgorithm::Dense_LLT;
    bool                  nonLinearIterative = false;
};
 
struct StressInfo
{
    StressInfo() {}
 
    std::vector<ElementStress> element_stress;
};
 
class CFiniteElementProblem
{
   public:
    using constraint_list_t = std::map<size_t, num_t>;
    using load_list_t       = std::map<size_t, num_t>;
 
    CFiniteElementProblem() = default;
    virtual ~CFiniteElementProblem();
 
    // ----------------------------------------------------------------------------
    virtual void clear();  
 
    bool loadFromStream(
        std::istream&                              is,
        const mrpt::optional_ref<vector_string_t>& errMsg  = std::nullopt,
        const mrpt::optional_ref<vector_string_t>& warnMsg = std::nullopt);
 
    bool loadFromFile(
        const std::string&                         file,
        const mrpt::optional_ref<vector_string_t>& errMsg  = std::nullopt,
        const mrpt::optional_ref<vector_string_t>& warnMsg = std::nullopt);
 
    bool saveAsImageSVG(
        const std::string& file, const DrawStructureOptions& options,
        const StaticSolveProblemInfo* solver_info  = nullptr,
        const MeshOutputInfo*         meshing_info = nullptr,
        ImageSaveOutputInfo*          out_img_info = nullptr) const;
    bool saveAsImagePNG(
        const std::string& file, const DrawStructureOptions& options,
        const StaticSolveProblemInfo* solver_info  = nullptr,
        const MeshOutputInfo*         meshing_info = nullptr,
        ImageSaveOutputInfo*          out_img_info = nullptr) const;
    bool saveAsImage(
        const std::string& file, const bool is_svg,
        const DrawStructureOptions&   options,
        const StaticSolveProblemInfo* solver_info  = nullptr,
        const MeshOutputInfo*         meshing_info = nullptr,
        ImageSaveOutputInfo*          out_img_info = nullptr) const;
 
    bool renderToCairoContext(
        void* _cairo_context, const RenderInitData& ri,
        const DrawStructureOptions&   options,
        const StaticSolveProblemInfo* solver_info,
        const MeshOutputInfo*         meshing_info) const;
 
    mrpt::opengl::CSetOfObjects::Ptr getVisualization(
        const DrawStructureOptions&   options,
        const StaticSolveProblemInfo& solver_info,
        const MeshOutputInfo*         meshing_info = nullptr,
        const StressInfo*             stressInfo   = nullptr) const;
 
    // ----------------------------------------------------------------------------
    size_t insertElement(CElement::Ptr el);
 
    template <typename ElementClass, typename... _Args>
    size_t createElement(_Args&&... __args)
    {
        return insertElement(
            std::make_shared<ElementClass>(std::forward<_Args>(__args)...));
    }
 
    CElement::ConstPtr getElement(size_t i) const;
 
    const CElement::Ptr& getElement(size_t i);
 
    size_t getNumberOfElements() const { return m_elements.size(); }
 
    // ----------------------------------------------------------------------------
    void insertConstraint(const size_t dof_index, const num_t value = 0);
 
    const constraint_list_t& getAllConstraints() const
    {
        return m_DoF_constraints;
    }
 
    void setLoadAtDOF(const size_t dof_index, const num_t f);
 
    void addLoadAtDOF(const size_t dof_index, const num_t f);
 
    const load_list_t& getOverallLoadsOnDOFs() const
    {
        return m_loads_at_each_dof;
    }
 
    // ----------------------------------------------------------------------------
    void setNumberOfNodes(size_t N);
 
    size_t getNumberOfNodes() const { return m_node_poses.size(); }
 
    node_index_t insertNode(const TRotationTrans3D& p);
 
    void setNodePose(size_t idx, const TRotationTrans3D& p);
 
    void setNodePose(size_t idx, const num_t x, const num_t y, const num_t z);
 
    TRotationTrans3D& getNodePose(size_t i)
    {
        ASSERT_(i < m_node_poses.size());
        return m_node_poses[i];
    }
    const TRotationTrans3D& getNodePose(size_t i) const
    {
        ASSERT_(i < m_node_poses.size());
        return m_node_poses[i];
    }
 
    void getNodeDeformedPosition(
        size_t i, Vector3& out_final_point,
        const StaticSolveProblemInfo& solver_info,
        const num_t                   exageration_factor = 1) const;
 
    num_t getMaximumDeformedDisplacement(
        const StaticSolveProblemInfo& solver_info) const;
 
    void getBoundingBox(
        num_t& min_x, num_t& max_x, num_t& min_y, num_t& max_y,
        bool                          deformed              = false,
        const StaticSolveProblemInfo* solver_info           = nullptr,
        num_t                         deformed_scale_factor = 1.0) const;
 
    // ----------------------------------------------------------------------------
    const std::vector<NodeDoF>& getProblemDoFs()
    {
        updateListDoFs();
        return m_problem_DoFs;
    }
 
    std::string getProblemDoFsDescription();
 
    static std::vector<size_t> complementaryDoFs(
        const std::vector<size_t>& ds, const size_t nTotalDOFs);
 
    size_t getDOFIndex(const size_t nNode, const DoF_index n) const;
 
    // ----------------------------------------------------------------------------
    virtual void updateAll();
 
    void assembleProblem(BuildProblemInfo& out_info);
 
    void solveStatic(
        StaticSolveProblemInfo&    out_info,
        const StaticSolverOptions& opts = StaticSolverOptions());
 
    void postProcCalcStress(
        StressInfo& out_stress, const StaticSolveProblemInfo& solver_info);
 
    std::string getNodeLabel(
        const size_t idx) const;  
 
    // ----------------------------------------------------------------------------
   protected:
    struct node_used_t
    {
        node_used_t() = default;
 
        bool used = false;
    };
 
    std::deque<node_used_t>      m_node_defined;
    std::deque<TRotationTrans3D> m_node_poses;
    std::vector<std::string>     m_node_labels;  
    std::deque<CElement::Ptr>    m_elements;
 
    constraint_list_t m_DoF_constraints;
 
    load_list_t m_loads_at_each_dof;
 
    load_list_t m_loads_at_each_dof_equivs;
 
    std::map<size_t, ElementStress> m_extra_stress_for_each_element;
 
    virtual void internalComputeStressAndEquivalentLoads() {}
 
    bool internal_loadFromYaml(
        const mrpt::containers::yaml&              f,
        const mrpt::optional_ref<vector_string_t>& errMsg,
        const mrpt::optional_ref<vector_string_t>& warnMsg);
 
    void internal_parser1_Parameters(
        const mrpt::containers::yaml& f, EvaluationContext& ctx) const;
 
    void internal_parser2_BeamSections(
        const mrpt::containers::yaml& f, EvaluationContext& ctx) const;
 
    void internal_parser3_nodes(
        const mrpt::containers::yaml& f, EvaluationContext& ctx);
    void internal_parser4_elements(
        const mrpt::containers::yaml& f, EvaluationContext& ctx);
    void internal_parser5_constraints(
        const mrpt::containers::yaml& f, EvaluationContext& ctx);
 
    void internal_parser6_node_loads(
        const mrpt::containers::yaml& f, EvaluationContext& ctx);
    void internal_parser7_element_loads(
        const mrpt::containers::yaml& f, EvaluationContext& ctx);
 
    void updateListDoFs();
 
    void updateNodeConnections();
 
    void updateElementsOrientation();
 
    void updateNodesMainOrientation();
 
    struct TProblemDOFIndicesForNode
    {
        TProblemDOFIndicesForNode() = default;
 
        std::array<int, 6> dof_index = {-1, -1, -1, -1, -1, -1};
    };
 
    std::vector<NodeDoF> m_problem_DoFs;
 
    std::vector<TProblemDOFIndicesForNode> m_problem_DoFs_inverse_list;
 
    struct TNodeElementConnection
    {
        unsigned char elementFaceId;  
        used_DoFs_t   dofs;  
    };
 
    using TNodeConnections = std::map<element_index_t, TNodeElementConnection>;
 
    std::vector<TNodeConnections> m_node_connections;
 
    std::vector<TRotation3D> m_nodeMainDirection;
 
    // Visualization subroutines:
    void internal_getVisualization_nodeLoads(
        mrpt::opengl::CSetOfObjects& gl, const DrawStructureOptions& options,
        const StaticSolveProblemInfo& solver_info,
        const MeshOutputInfo* meshing_info, num_t DEFORMED_SCALE_FACTOR) const;
 
    void internal_getVisualization_constraints(
        mrpt::opengl::CSetOfObjects& gl, const DrawStructureOptions& options,
        const StaticSolveProblemInfo& solver_info,
        const MeshOutputInfo* meshing_info, num_t DEFORMED_SCALE_FACTOR) const;
 
    void internal_getVisualization_distributedLoads(
        const CStructureProblem& str, mrpt::opengl::CSetOfObjects& gl,
        const DrawStructureOptions&   options,
        const StaticSolveProblemInfo& solver_info,
        const MeshOutputInfo* meshing_info, num_t DEFORMED_SCALE_FACTOR) const;
 
    void internal_getVisualization_stressDiagrams(
        mrpt::opengl::CSetOfObjects& gl, const DrawStructureOptions& options,
        const StaticSolveProblemInfo& solverInfo,
        const MeshOutputInfo* meshingInfo, num_t DEFORMED_SCALE_FACTOR,
        const StressInfo& stressInfo) const;
};
}  // namespace openbeam