ShipSimObject.h

#ifndef _CSHIP_SIM_OBJECT_H_
#define _CSHIP_SIM_OBJECT_H_
 
#include <fhsim/simobject/ISimObjectCreator.h>
#include <fhsim/FhSimMgr.h>
#include <fhsim/PrintDuringExec.h>
#include <array>
#include <cmath>
#include <limits>
#include <memory>
#include <vector>
#include <string>
#include <queue>
#include <string>
#include <iostream>
#include <locale>
 
#ifdef FH_VISUALIZATION
#pragma warning(push)
#pragma warning( disable : 4251 )
    #include <Ogre.h>
#pragma warning(pop)
    //#include <fhsim/visual/renderer/FhVisualization.h>
    #include "sfh/ogre/C3DArrow.h"
    #include "sfh/ogre/DotSceneLoader.h"
#endif
 
#include <fhsim_environment/EnvironmentProvider.h>
 
#include <sfh/math/interpolation/SplineInterp.h>
 
// TODO: this will not work if someone inherits the simobject externally
#include "Ship/Resistance/ShipResistance.h"
#include "Ship/ViscousShip.h"
#include "Ship/HullData/HullData.h"
#include "Ship/HullData/VeresHullData.h"
#include "Ship/HullData/NapaCGALHullData.h"
#include "Ship/PotentialData/PotentialData.h"
#include "Ship/PotentialData/VeresData.h"
 
namespace Ship{
    struct ShipParameters
    {
        std::string forceInput = "NED";
        int numWaveComponents = 52;
        bool zeroSpeedRadiation = true;
        double stabilizeDOF[6] = {0,0,0,0,0,0};
        double alfaN = 0;
        double alfaE = 0;
        double alfaD = 0;
        double surgeVelocity = 0;
 
        std::string potentialFileLocation;
        std::string potentialFileName;
        std::string potentialType = "VERES";
        int flipPhaseForOSC = 0;
 
        int DOFFree[6] = {1,1,1,1,1,1};
        double forceScale[11] = {1,1,1,1,1,1,1,1,1,1,1};
        bool coriolisRollOnly = false;
        double waveForceTransient = 0.1;
        bool hydroStaticsType = false;
        int waveForceType = 0;
        bool waveForceHFrame = false;
        double hydroStaticsStationVerticalStep = 1;
        double hydroStaticsStationHorizontalStep = 1;
        double hydroStaticsStationTop = 15;
        double hydroStaticsStationBottom = -1;
 
        bool useCrossFlow = true;
        double damperLinear[3] = {0,0,0};
        double damperQuadratic[3] = {0,0,0};
        double damperRotLinear[3] = {0,0,0};
        double damperRotQuadratic[3] = {0,0,0};
        double wakeFraction = 0;
 
        double Yv_man = 0, Yr_man = 0, Nv_man = 0, Nr_man = 0;
        double Yvd_man = 0, Yrd_man = 0, Nvd_man = 0, Nrd_man = 0;
        double Xud_man = 0, Xvr_man = 0, Xrr_man = 0, Xvv_man = 0, Xvvvv_man = 0;
        double Xu_dp = 0, Xv_dp = 0, Xr_dp = 0;
        double Yu_dp = 0, Yv_dp = 0, Yr_dp = 0;
        double Nu_dp = 0, Nv_dp = 0, Nr_dp = 0;
        double dampingBlendSpeed = 8;
        double dampingBlendFraction = 0.9;
        double manoeuvreMinSpeed = 0.1;
 
        double bilgeKeelStartPosition = std::numeric_limits<double>::quiet_NaN();
        double bilgeKeelEndPosition = std::numeric_limits<double>::quiet_NaN();
        double bilgeKeelHeight = 0;
        int withSkeg = 0;
        bool skegLift = false;
        int skegStartSection = 5;
        int skegEndSection = 15;
        std::vector<int> skegStartPoints = {42, 27, 24, 24, 24, 24, 24, 26, 27, 29, 35};
        std::vector<int> skegEndPoints = std::vector<int>(11, 49);
        bool hullLift = true;
        int rollDampingTheory = 0;   
        double rollAmplitudeHorizon = 2.75*sqrt(10.0);   
 
        bool hasResistanceFriction = false;   
        std::string resistanceFriction;       
        std::string resistanceFrictionCurveX;
        std::string resistanceFrictionCurveY;
        bool hasResistanceResidual = false;   
        std::string resistanceResidual;       
        double hollenbachNumRudders = 0;
        double hollenbachNumBrackets = 0;
        double hollenbachNumBoss = 0;
        double hollenbachNumThrusters = 0;
        double hollenbachDiamPropeller = 0;
        double scaleResidualResistance = 1.0;
        std::string resistanceResidualCurveX;
        std::string resistanceResidualCurveY;
    };
 
    class ShipSimObject:public SimObject/*, public ICurrentField*/
    {
    public:
        EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 
    #ifdef FH_VISUALIZATION
        enum RenderType
        {
            BOX,
            STATIONS,
            SLICED,
            INTERPOLATED
        };
    #endif
 
        ShipSimObject( const std::string simObjectName, ISimObjectCreator* const creator);
 
        ~ShipSimObject();
 
        void FinalSetup( const double T, const double *const X, ISimObjectCreator* const creator );
 
        virtual void PreOdeFcn(const double T, const double *const X);
 
        virtual void OdeFcn(const double T, const double* const X, double* const XDot) const;
 
        void AcceptedStep(double T, const double* X) override;
 
    #ifdef FH_VISUALIZATION
 
        virtual void RenderInit(Ogre::Root* const ogreRoot, ISimObjectCreator* const creator);
 
        virtual void RenderUpdate( const double T, const double* const X );
    #endif
 
 
    public:
 
        const double* PositionNED(const double T, const double* const X);
 
        const double* PositionZeroCrossNED(const double T, const double* const X);
 
        const double* PositionAPNED(const double T, const double* const X);
 
        const double* VelocityNED(const double T, const double* const X);
 
        const double* QuaternionNED(const double T, const double* const X);
 
        const double* EulerAnglesNED(const double T, const double* const X);
 
        const double* RotationBody(const double T, const double* const X);
 
        const double* OmegaNED(const double T, const double* const X);
 
 
        void ComputeRollPitchYaw(const double T, const double* const X) const;
        void ComputeRelativeVelocity(const double T, const double* const X) const;
 
        const double* BodyForce(const double T, const double* const X);
 
        const double* VelocityBody(const double T, const double* const X);
 
        const double* FirstOrderWaveBody(const double T, const double* const X);
 
        const double* WaveDriftBody(const double T, const double* const X);
 
        const double* const FluidMemForceBody(const double T, const double* const X);
 
        const double* RestoringForceBody(const double T, const double* const X);
 
        const double* RelativeVelocityBody(const double T, const double* const X);
 
        const double* ViscousForceBody(const double T, const double* const X);
 
        const double* Encounter(const double T, const double* const X);
 
        const double* CoriolisRigidBody(const double T, const double* const X);
 
        const double* CoriolisFluidBody(const double T, const double* const X);
 
        const double* SpeedInfo( const double T, const double *const X);
 
        const double* RotationForceFrame( const double T, const double *const X);
 
        const double * KineticEnergy(const double T, const double* const X);
 
        const double * FluidEnergy(const double T, const double* const X);
        const double * FluidEnergy2(const double T, const double* const X);
        const double * UFluidMemSpeed(const double T, const double* const X);
 
        const double* FluidForceBySpeedBody(const double T, const double* const X, int index);
        const double* FluidMemoryForceBody( const double T, const double *const X);
 
        const double * CrossFlowBody(const double T, const double* const X);
 
        const double * ResistanceBody(const double T, const double* const X);
 
        const double * AdditionalDampeningBody(const double T, const double* const X);
 
        const double* ExternalForceNED(const double T, const double* const X);
 
        const double* ExternalForceBody(const double T, const double* const X);
        const double * ViscousRoll(const double T, const double* const X);
        double m_ViscousRoll[6];
 
        const double * DriftSpeed(const double T, const double* const X);
 
        ViscousShip& Viscous() { return m_viscous; }
        const ViscousShip& Viscous() const { return m_viscous; }
 
 
        double m_waveElevation;
        double m_waveSlope;
        /*const double * WaveElevation(const double T, const double* const X);
        const double * WaveSlope(const double T, const double* const X);*/
 
    private:
        static ShipParameters ReadParameters( ISimObjectCreator* const creator );
 
        void ComputeVelocityBody(const double T, const double* const X, double* const out) const;
        void ComputeFirstOrderWaveBody(const double T, const double* const X, double* const out) const;
        Eigen::Matrix<double,6,1> ComputeBodyForceTotal(const double T, const double* const X) const;
        void ComputeBodyForceTotalCommon(const double T, const double* const X) const { m_bodyForceTotal = ComputeBodyForceTotal(T,X); }
        void ComputeWaveDriftBody(const double T, const double* const X, double* const out) const;
        void ComputeRestoringForceBody(const double T, const double* const X, double* const out) const;
        void ComputeViscousForceBody(const double T, const double* const X, double* const out) const;
        void ComputeCrossFlowBody(const double T, const double* const X, double* const out) const;
        void ComputeAdditionalDampeningBody(const double T, const double* const X, double* const out) const;
        void ComputeResistanceBody(const double T, const double* const X, double* const out) const;
        void ComputeCoriolisRigidBody(const double T, const double* const X, double* const out) const;
        void ComputeCoriolisFluidBody(const double T, const double* const X, double* const out) const;
        void ComputeSpeedInfo(const double T, const double* const X, double* const out) const;
        void ComputeFluidForceBySpeedBody(const double T, const double* const X, int index, double* const out) const;
 
        void printShipInfo();
 
        PrintDuringExec* m_logg;
 
        Eigen::Vector3d m_forceExtNED;
        Eigen::Vector3d m_forceExtBody;
 
        Eigen::Vector3d m_PositionMidpoint;
        Eigen::Vector3d m_PositionAP;
        Eigen::Vector3d m_angles;
 
        Eigen::Matrix<double,6,6> m_mrb;
        Eigen::Matrix<double,6,6> m_invMrb;
 
        std::unique_ptr<HullData>       m_Hull;
        std::unique_ptr<ShipResistance> m_resitance;
                environment::EnvironmentProvider     *m_environment = nullptr;
 
 
        std::unique_ptr<sfh::math::SplineInterp> hoerner;
 
 
        ISignalPort *m_InForceNED;
 
        ISignalPort *m_InMomentNED;
 
        ISignalPort *m_InForceBODY;
 
        ISignalPort *m_InMomentBODY;
 
        ICommonComputation *m_computeRelativeVelocity;
        ICommonComputation *m_computeRotationAngles;
        ICommonComputation *m_computeBodyForceTotal = nullptr;
        mutable Eigen::Matrix<double,6,1> m_bodyForceTotal; 
 
        std::unique_ptr<PotentialData> m_TheoryData;
 
        int m_numWaveComponents = 52;
        bool RefreshWaveComponents(const double T);
 
        std::vector<double> m_waveHeight;
        std::vector<double> m_waveFreq;
        std::vector<double> m_wavePhase;
        std::vector<double> m_waveDirection;
        std::vector<double> m_waveNumber;
 
        std::vector<double> m_uHistory;
        std::vector<double> m_psiHistory;
        std::vector<double> m_tHistory;
 
        std::vector<std::array<double,6>> m_fluidForce;
        double m_fluidForceExperienced[6];
        double m_rotationForceFrame[4];
        double m_speedInformation[5];
        double m_kinetic;
        double m_fluid;
        double m_fluid2[36];
 
        double m_damping[6];
 
        mutable double m_rollPitchYaw[3];
        double m_speed[6];
        mutable double m_relativeSpeed[6];
 
        double m_firstOrderForce[6];
        double m_waveDriftForce[6];
        double m_CforcesBody[6];
        double m_CforcesFluid[6];
 
        double m_crossFlowBody[6];
        double m_resistance[6];
 
        double m_fluidMemForce[6];
        double m_restoringForce[6];
        double m_sumBodyForce[6];
        double m_encounterFreq[52];
        double m_viscousForce[6];
        ViscousShip m_viscous;
 
        double m_hydrostaticSpring[3];
 
        double m_driftSpeed;
 
        bool    m_zeroSpeedRadiation;
        bool    m_hydroStaticsType;
        bool    m_useCrossFlow;
        bool    m_waveforceHframe;
        int     m_waveForceType;
        int     m_DOFFree[6];
        double  m_alfaN;
        double  m_alfaE;
        double  m_alfaD;
        double  m_surgeVelocity;
        double  m_wakeFraction;
        double  m_hydroStaticsStationVerticalStep;
        double  m_hydroStaticsStationHorizontalStep;
        double  m_hydroStaticsStationTop;
        double  m_hydroStaticsStationBottom;
        double  m_bilgeKeelStartPosition;
        double  m_bilgeKeelEndPosition;
        double m_bilgeKeelHeight;
        double  m_damperLinear[3];
        double  m_damperQuadratic[3];
        double  m_damperRotLinear[3];
        double  m_damperRotQuadratic[3];
        double m_dampingBlendSpeed;
        double m_dampingBlendFraction;
 
        double m_waveForceTransient;
 
        double  m_ForceScale[11];
        std::string m_ForceInput;
        std::string m_PotentialFileLocation;
        std::string m_PotentialFileName;
        std::string m_PotentialType;
        double m_stablizeDOF[6];
 
        bool m_coriolisRollOnly;
        Eigen::Matrix<double, 6,1> m_coriolisDOFmangler;
 
        static inline Eigen::Matrix3d S(const Eigen::Vector3d &v){
            Eigen::Matrix3d retval = Eigen::Matrix3d::Zero();
            retval(0,1)=-v(2);
            retval(0,2)=v(1);
            retval(1,0)=v(2);
            retval(1,2)=-v(0);
            retval(2,0)=-v(1);
            retval(2,1)=v(0);
            return retval;
        }
 
#ifdef FH_VISUALIZATION
        bool m_ShowForceArrows = false;
        std::string m_Material;
        std::string m_Mesh;
        std::string m_Scene;
        RenderType m_RenderAs;
        std::string m_RenderType;
 
    #endif
 
        std::vector<std::array<std::array<int,6>,6>> m_IFluidState;
        std::vector<std::array<std::array<int,6>,6>> m_FluidStateSize;
 
        unsigned int m_ISR;
        unsigned int m_ISV;
        unsigned int m_ISQ;
        unsigned int m_ISW;
 
 
    #ifdef FH_VISUALIZATION
 
    public:
        void DrawShip(Ogre::SceneNode* renderNode, Ogre::SceneManager* sceneMgr);
        void DrawStations( Ogre::SceneNode* m_RenderNode, Ogre::SceneManager * sceneMgr );
    private:
        Ogre::ManualObject *manual ;
        Ogre::SceneNode*    m_RenderNodeHull;
        Ogre::SceneNode*    m_RenderNodeZeroCross;
        Ogre::SceneNode*    m_RenderNodeStations;
        Ogre::SceneNode**   m_RenderNodeSkegPoint;
        double m_HideUntil;
        double m_NodeOffset[3];
        C3DArrow *FKF_u,*FKF_v,*FKF_w;
        C3DArrow *FKT_u,*FKT_v,*FKT_w;
 
        C3DArrow *FMF_u,*FMF_v,*FMF_w;
        C3DArrow *FMT_u,*FMT_v,*FMT_w;
 
        C3DArrow *RSF_u,*RSF_v,*RSF_w;
        C3DArrow *RST_u,*RST_v,*RST_w;
 
        C3DArrow *FKArrow;
        C3DArrow *FMArrow;
        C3DArrow *RArrow;
        C3DArrow *Res;
 
        C3DArrow *FKArrow_M;
        C3DArrow *FMArrow_M;
        C3DArrow *RArrow_M;
        C3DArrow *Res_M;
 
        C3DArrow *VArrow;
    #endif
 
        double m_maxWavePeriod;
    };
}
 
#endif