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/system/CTimeLogger.h>
#include <mrpt/viz/CSetOfObjects.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;
 
  Eigen::Matrix<num_t, Eigen::Dynamic, 1> F_b_applied;
};
 
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;
 
  num_t rcond = std::numeric_limits<num_t>::quiet_NaN();
 
  std::vector<std::string> warnings;
};
 
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::Sparse_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::viz::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; }
 
  bool addNodeConstraint(node_index_t node, DoF_index dof, num_t value = 0);
 
  const std::map<std::pair<node_index_t, DoF_index>, num_t>& getNodeConstraintRequests() const
  {
    return m_node_constraint_requests;
  }
 
  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;
 
  std::map<std::pair<node_index_t, DoF_index>, num_t> m_node_constraint_requests;
 
  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();
 
  std::string describeSingularStiffness(
      const DynMatrix& Kff, const std::vector<size_t>& free_dof_indices) const;
 
  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::viz::CSetOfObjects& gl,
      const DrawStructureOptions& options,
      const StaticSolveProblemInfo& solver_info,
      const MeshOutputInfo* meshing_info,
      num_t DEFORMED_SCALE_FACTOR) const;
 
  void internal_getVisualization_constraints(
      mrpt::viz::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::viz::CSetOfObjects& gl,
      const DrawStructureOptions& options,
      const StaticSolveProblemInfo& solver_info,
      const MeshOutputInfo* meshing_info,
      num_t DEFORMED_SCALE_FACTOR) const;
 
  void internal_getVisualization_stressDiagrams(
      mrpt::viz::CSetOfObjects& gl,
      const DrawStructureOptions& options,
      const StaticSolveProblemInfo& solverInfo,
      const MeshOutputInfo* meshingInfo,
      num_t DEFORMED_SCALE_FACTOR,
      const StressInfo& stressInfo) const;
};
}  // namespace openbeam