Cable.h

#pragma once
#include <fhsim_environment/hydrodynamics/CylinderHydro.h>
 
#include <Eigen/Eigen>
#include <cmath>
 
#include "CableJointMoment.h"
#include "cable/GaussPointDrag.h"
#include "ElementInertia.h"
 
#include <fhsim/PrintDuringExec.h>
#include <fhsim/simobject/SimObject.h>
#include <fhsim_environment/EnvironmentProvider.h>
#include <fhsim_environment/SeafloorForcesInterface.h>
#include <fhsim_environment/ElementFlowInteraction.h>
#ifdef FH_VISUALIZATION
#   include "sfh/ogre/C3DLine.h"
#endif
 
namespace RbCable
{
inline bool IsBelowSurface(double z, double surfaceElevation)
{
   return z + surfaceElevation > 0.0;
}
 
class CableRM : public SimObject
{
public:
   CableRM(const std::string& simObjectName, ISimObjectCreator* const creator); 
   ~CableRM();
   void OdeFcn(const double T, const double* const X, double* const XDot) const;
 
   void InitialConditionSetup(const double T, const double* const currentIC, double* const updatedIC, ISimObjectCreator* const creator);
   void FinalSetup(const double T, const double* const X, ISimObjectCreator* const creator);
   void PreOdeFcn(const double T, const double* const X) override;
 
   const double* forceA(const double T, const double* const X);
   const double* forceB(const double T, const double* const X);
   void calculations(const double T, const double* const X);
 
   void BuildChain(wet_cylinder_water::ChainJoints& chain, const double* X, int firstElement, int numElements) const;
   void QueryJointWater(const double T, const double* const X);
   void calculationsCommon(const double T, const double* const X) const;
 
#ifdef FH_VISUALIZATION
   void RenderInit(Ogre::Root* const ogreRoot, ISimObjectCreator* const creator);
   void RenderUpdate(const double T, const double* const X);
#endif
 
protected:
   typedef Eigen::Matrix<double, 3, 3> mat3;
   typedef Eigen::Matrix<double, 3, 1> vec3;
 
   void DistributeCatenary(Eigen::Matrix<double, 3, 1> P1, Eigen::Matrix<double, 3, 1> P2, double L, double lengthA, double lengthB, const double* currentIC, double* states, int i1, int i2);
 
   void DistributeOnArcMidpoints(const vec3& P1, const Eigen::Matrix<double, 2, 1>& CatenaryU, double HorizontalDistance, double VerticalDifference, double L, const double* currentIC, double* states, int i1, int i2);
 
   void SetElementUnlessUserSet(int i, const vec3& centre, double angle, const Eigen::Matrix<double, 2, 1>& CatenaryU, const double* currentIC, double* states);
 
   [[noreturn]] void ReportError(const std::string& message) const;
 
   PrintDuringExec*                  m_print; 
   environment::SeafloorForcesInterface*     m_seafloor = nullptr; 
   environment::wake::ElementFlowInteraction m_flow; 
 
   int    m_numElements;
   double m_totalLength;
   double m_radius;
   double m_weight;
   double m_alphaN;
   double m_betaN;
   double m_epsilonN;
 
   JointLaw m_joint; 
 
   // double  m_rho;
   // double  m_surfacePosZ;
   // double    m_rhoWater;
   // double  m_Ct; // tangential damping coefficient
   // double  m_Cn; // normal damping coefficient
 
   double m_length;
   double m_mass;
   double m_Ixy;
   double m_Iz;
   // bool TEMP_BOOL_VAR__ISMAJORTIMESTEP;
 
   // Sugar kelp rope application - start
   bool   m_kelp;  // 1 if used as a kelp rope, 0 otherwise
   double m_wKelp; // Biomass of kelp per meter rope/cable in Newtons
   // Sugar kelp rope application - end
 
   bool                           m_gaussPointDrag = true;  
   wet_cylinder_water::ChainJoints m_chainJoints;            
   wet_cylinder_water::ChainJoints m_takeJoints;             
 
   struct element
   {
      int    p;
      int    q;
      int    v;
      int    w;
      mat3   Mi;   
      mat3   MiT;  
      vec3   k;
      vec3   eDot;
      double nuDot;
   };
   element* m_el;
 
   ISignalPort* m_posA;
   ISignalPort* m_posB;
   ISignalPort* m_velA;
   ISignalPort* m_velB;
   ISignalPort* m_retractedLengthA;
   ISignalPort* m_retractedLengthB;
   ISignalPort* m_retractedSpeedA;
   ISignalPort* m_retractedSpeedB;
 
   int m_retractedNodesA;
   int m_retractedNodesB;
   int m_numFreeNodes;
 
   struct Retraction
   {
      double lengthA; 
      double lengthB; 
      double speedA;  
      double speedB;  
      int    nodesA;  
      int    nodesB;  
   };
 
   Retraction ClampRetraction(double lengthA, double lengthB, double speedA, double speedB) const;
 
   vec3 m_ka;
   vec3 m_kb;
 
   ICommonComputation*                      m_calcDynamics;
   Eigen::Matrix<double, Eigen::Dynamic, 1> m_lambda;
   Eigen::Matrix<double, Eigen::Dynamic, 1> m_F_MDotV; //< external force F, minus Coriolis term dM/dt*V
   Eigen::Matrix<double, Eigen::Dynamic, 1> m_B;       //< right-hand side of the constraint system in OdeFcn; a member so that OdeFcn does not allocate it
 
   double m_forceA[3];
   double m_forceB[3];
 
#ifdef FH_VISUALIZATION
   double            m_scale;
   Ogre::SceneNode** m_ManualObjectNodes;
#endif
};
} // namespace RbCable