Line.hpp

Go to the documentation of this file.

/*
 * Copyright (c) 2022, Matt Hall
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its
 *    contributors may be used to endorse or promote products derived from
 *    this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */
 
#pragma once
 
#include "Misc.hpp"
#include "Instance.hpp"
#include "Seafloor.hpp"
#include "Util/CFL.hpp"
#include "leanvtk/leanvtk.hpp"
#include <utility>
 
using namespace std;
 
namespace moordyn {
 
class Waves;
typedef std::shared_ptr<Waves> WavesRef;
 
class DECLDIR Line final
  : public Instance
  , public NatFreqCFL
{
  public:
    Line(moordyn::Log* log, size_t lineId);
 
    ~Line();
 
  private:
    typedef enum
    {
        ELASTIC_CONSTANT = 1,
        ELASTIC_VISCO_CTE = 2,
        ELASTIC_VISCO_MEAN = 3,
        ELASTIC_SYROPE = 4,
    } elastic_model;
 
    typedef enum
    {
        SYROPE_LINEAR_WC = 1,
        SYROPE_QUADRATIC_WC = 2,
        SYROPE_EXP_WC = 3,
    } syrope_wc_formula;
 
    real getNonlinearEA(real l_stretched, real l_unstretched) const;
 
    real getNonlinearBA(real ld_stretched, real l_unstretched) const;
 
    real getNonlinearEI(real curv) const;
 
    inline real getWaterDepth(real x, real y)
    {
        return seafloor ? seafloor->getDepthAt(x, y) : -env->WtrDpth;
    }
 
    inline real avgWaterDepth()
    {
        return seafloor ? seafloor->getAverageDepth() : -env->WtrDpth;
    }
 
    inline void setWorkingCurve(real Tmax);
 
    // ENVIRONMENTAL STUFF
    EnvCondRef env;
    moordyn::WavesRef waves;
    moordyn::SeafloorRef seafloor;
 
    unsigned int N;
    moordyn::real UnstrLen;
    moordyn::real UnstrLen0;
    moordyn::real UnstrLend;
    moordyn::real d;
    moordyn::real rho;
    elastic_model ElasticMod;
    moordyn::real EA;
    moordyn::real EA_D;
    moordyn::real alphaMBL;
    moordyn::real vbeta;
    moordyn::real EA_2;
    std::vector<moordyn::real> dl_1;
    std::vector<moordyn::real> ld_1;
    moordyn::real EI;
    moordyn::real BA;
    moordyn::real BA_D;
    moordyn::real Can;
    moordyn::real Cat;
    moordyn::real Cdn;
    moordyn::real Cdt;
    moordyn::real Cl;
    moordyn::real dF;
    moordyn::real cF;
 
    moordyn::real BAin;
    moordyn::real A;
    unsigned int nEApoints;
    std::vector<moordyn::real> stiffXs;
    std::vector<moordyn::real> stiffYs;
    std::vector<moordyn::real> stiffZs;
    unsigned int nEIpoints;
    std::vector<moordyn::real> bstiffXs;
    std::vector<moordyn::real> bstiffYs;
    unsigned int nBApoints;
    std::vector<moordyn::real> dampXs;
    std::vector<moordyn::real> dampYs;
 
    std::vector<moordyn::real> stiffxs_;
    std::vector<moordyn::real> stiffys_;
    std::vector<moordyn::real> stiffzs_;
 
    // Externally provided data
    bool isPb;
    std::vector<moordyn::real> pin;
 
    // kinematics
    std::vector<vec> r;
    std::vector<vec> rd;
    std::vector<vec> q;
    std::vector<vec> pvec;
    std::vector<vec> qs;
    std::vector<moordyn::real> l;
    std::vector<moordyn::real> lstr;
    std::vector<moordyn::real> ldstr;
    std::vector<moordyn::real> Kurv;
 
    std::vector<mat> M;
    // line segment volumes
    std::vector<moordyn::real> V;
 
    // Syrope variables
    syrope_wc_formula syropeWCForm;
    moordyn::real k1;
    moordyn::real k2;
    std::vector<moordyn::real> Tmax;
    std::vector<moordyn::real> Tmean;
    const unsigned int nEApointsWC = 30;
 
    // forces
    std::vector<vec> T;
    std::vector<vec> Td;
    std::vector<vec> Bs;
    std::vector<vec> Pb;
    std::vector<vec> W;
    std::vector<vec> Dp;
    std::vector<vec> Dq;
    std::vector<vec> Ap;
    std::vector<vec> Aq;
    std::vector<vec> B;
    std::vector<vec> Lf;
    std::vector<vec> Fnet;
 
    // wave things
    std::vector<moordyn::real> F;
 
    // time
    moordyn::real t;
    moordyn::real dtm0;
 
    // VIV stuff
    // /// VIV amplitude updated every zero crossing of crossflow velcoity
    // std::vector<moordyn::real> Amp;
    const unsigned int n_m = 500;
    real phi_yd;
    std::vector<moordyn::real> phi;
    std::vector<moordyn::real> phi_dot;
 
    // back indexing one dtm for VIV
    moordyn::real t_old = 0.0;
    // /// running amplitude total, from previous zero crossing of yd
    // std::vector<moordyn::real> A_int_old;
    std::vector<moordyn::real> yd_rms_old;
    std::vector<moordyn::real> ydd_rms_old;
    std::vector<vec> rdd_old;
 
    // end conditions
    typedef enum
    {
        PINNED = 0,
        CANTILEVERED = 1,
        // Some aliases
        POINT = PINNED,
        ROD = CANTILEVERED,
    } endTypes;
 
    endTypes endTypeA;
    endTypes endTypeB;
    vec endMomentA;
    vec endMomentB;
 
    // file stuff
    ofstream* outfile;
    string channels;
 
    std::vector<std::vector<moordyn::real>> zetaTS;
    std::vector<std::vector<moordyn::real>> FTS;
    std::vector<std::vector<vec>> UTS;
    std::vector<std::vector<vec>> UdTS;
    unsigned int ntWater;
    moordyn::real dtWater;
 
  public:
    int number;
    size_t lineId;
 
    bool IC_gen = false;
 
    void setup(int number,
               LineProps* props,
               real l,
               unsigned int n,
               EnvCondRef env_in,
               shared_ptr<ofstream> outfile,
               string channels,
               real dtM0);
 
    inline void setWaves(moordyn::WavesRef waves_in,
                         moordyn::SeafloorRef seafloor_in)
    {
        waves = waves_in;
        seafloor = seafloor_in;
    }
 
    void initialize();
 
    inline void initialize(InstanceStateVarView state)
    {
        // ------ Initialize the line ------
        initialize();
 
        // ------ Assign the intialized values to the state (bascially
        // Line::setState but flipped) ------ Error check for number of columns
        // (if VIV and Visco need row.size() = 8, if VIV xor Visco need
        // row.size() = 7, if not VIV need row.size() = 6)
        if ((state.row(0).size() != 8 && Cl > 0 && ElasticMod != ELASTIC_CONSTANT) ||
            (state.row(0).size() != 7 && ((Cl > 0) ^ (ElasticMod != ELASTIC_CONSTANT))) ||
            (state.row(0).size() != 6 && Cl == 0 && ElasticMod == ELASTIC_CONSTANT)) {
            LOGERR << "Invalid state.row size for Line " << number << endl;
            throw moordyn::mem_error("Invalid state.row size");
        }
 
        // Error check for number of rows (if visco need N rows, if normal need
        // N-1 rows)
        if ((state.rows() != N && ElasticMod != ELASTIC_CONSTANT) ||
            (state.rows() != N - 1 && ElasticMod == ELASTIC_CONSTANT)) {
            LOGERR << "Invalid number of rows in state matrix for Line "
                   << number << endl;
            throw moordyn::mem_error("Invalid number of rows in state matrix");
        }
 
        // If using the viscoelastic model, iterate N rows, else iterate N-1
        // rows.
        for (unsigned int i = 0; i < (ElasticMod != ELASTIC_CONSTANT ? N : N - 1); i++) {
            // node number is i+1
            // segment number is i
            state.row(i).head<3>() = r[i + 1];
            state.row(i).segment<3>(3) = rd[i + 1];
 
            if (ElasticMod != ELASTIC_CONSTANT)
                state.row(i).tail<1>()[0] =
                    dl_1[i]; // [0] needed becasue tail<1> returns a one element
                             // vector. Viscoelastic state is always the last
                             // element in the row
 
            if (Cl > 0) {
                if (ElasticMod != ELASTIC_CONSTANT)
                    state.row(i).tail<2>()[0] =
                        phi[i + 1]; // if both VIV and viscoelastic second to
                                    // last element in the row
                else
                    state.row(i).tail<1>()[0] =
                        phi[i + 1]; // else last element in the row
            }
        }
    }
    inline void setInitialTmax(moordyn::real Tmax0) { Tmax.assign(N, Tmax0); }
 
    inline void setInitialTmean(moordyn::real Tmean0)
    {
        Tmean.assign(N, Tmean0);
    }
 
    inline elastic_model getElasticModel() const { return ElasticMod; }
 
    inline unsigned int getN() const { return N; }
 
    inline moordyn::real getUnstretchedLength() const { return UnstrLen; }
 
    inline void setUnstretchedLength(const moordyn::real len)
    {
        UnstrLen = len;
        for (unsigned int i = 0; i < N; i++) {
            l[i] = UnstrLen / double(N);
            V[i] = l[i] * A;
        }
    }
 
    inline void setUnstretchedLengthVel(const moordyn::real v)
    {
        UnstrLend = v;
    }
 
    inline void updateUnstretchedLength(const moordyn::real dt = 0.0)
    {
        if (!UnstrLend)
            return;
        if (!dt) {
            UnstrLen0 = UnstrLen;
            return;
        }
        setUnstretchedLength(UnstrLen0 + dt * UnstrLend);
    }
 
    inline bool isConstantEA() const { return nEApoints > 0; }
 
    inline moordyn::real getConstantEA() const { return EA; }
 
    inline void setConstantEA(moordyn::real EA_in)
    {
        if (ElasticMod != ELASTIC_CONSTANT) {
            LOGERR << "Cannot set constant EA for viscoelastic model" << endl;
            throw moordyn::invalid_value_error(
                "Cannot set constant EA for viscoelastic model");
        } else {
            EA = EA_in;
        }
    }
 
    inline const vec& getNodePos(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        if (isnan(r[i].sum())) {
            stringstream s;
            s << "NaN detected" << endl
              << "Line " << number << " node positions:" << endl;
            for (unsigned int j = 0; j <= N; j++)
                s << j << " : " << r[j] << ";" << endl;
            throw moordyn::nan_error(s.str().c_str());
        }
        return r[i];
    }
 
    inline const vec& getNodeVel(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        if (isnan(rd[i].sum())) {
            stringstream s;
            s << "NaN detected" << endl
              << "Line " << number << " node velocities:" << endl;
            for (unsigned int j = 0; j <= N; j++)
                s << j << " : " << rd[j] << ";" << endl;
            throw moordyn::nan_error(s.str().c_str());
        }
        return rd[i];
    }
 
    inline const vec& getNodeForce(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return Fnet[i];
    };
 
    inline const vec getNodeTen(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        if (i == 0) // bottom node, return ten and damping of bottom section
            return T[0] + Td[0];
        else if (i == N) // top node, return ten and damping of top section
            return T[N - 1] + Td[N - 1];
        // internal node, take average of tension in adjacent segments.
        return (0.5 * (T[i] + T[i - 1] + Td[i] + Td[i - 1]));
    };
 
    inline const vec& getNodeBendStiff(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return Bs[i];
    };
 
    inline const vec& getNodeWeight(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return W[i];
    };
 
    inline const vec getNodeDrag(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return Dp[i] + Dq[i];
    };
 
    inline const vec getNodeFroudeKrilov(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return Ap[i] + Aq[i];
    };
 
    inline const vec getNodeSeabedForce(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return B[i] + Bs[i];
    };
 
    inline const real& getNodeCurv(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return Kurv[i];
    }
 
    inline const mat& getNodeM(unsigned int i) const
    {
        if (i > N) {
            LOGERR << "Asking node " << i << " of line " << number
                   << ", which only has " << N + 1 << " nodes" << std::endl;
            throw moordyn::invalid_value_error("Invalid node index");
        }
        return M[i];
    }
 
    inline std::vector<vec> getNodeCoordinates() const { return r; }
 
    inline void getFASTtens(float* FairHTen,
                            float* FairVTen,
                            float* AnchHTen,
                            float* AnchVTen) const
    {
        *FairHTen = (float)(Fnet[N](Eigen::seqN(0, 2)).norm());
        *FairVTen = (float)(Fnet[N][2] + M[N](0, 0) * (-env->g));
        *AnchHTen = (float)(Fnet[0](Eigen::seqN(0, 2)).norm());
        *AnchVTen = (float)(Fnet[0][2] + M[0](0, 0) * (-env->g));
    }
 
    inline void getEndStuff(vec& Fnet_out,
                            vec& Moment_out,
                            mat& M_out,
                            EndPoints end_point) const
    {
        switch (end_point) {
            case ENDPOINT_TOP:
                Fnet_out = Fnet[N];
                Moment_out = endMomentB;
                M_out = M[N];
                break;
            case ENDPOINT_BOTTOM:
                Fnet_out = Fnet[0];
                Moment_out = endMomentA;
                M_out = M[0];
                break;
            default:
                LOGERR << "Invalid end point qualifier: " << end_point << endl;
                throw moordyn::invalid_value_error("Invalid end point");
        }
    }
 
    real GetLineOutput(OutChanProps outChan);
 
    void storeWaterKin(real dt,
                       std::vector<std::vector<moordyn::real>> zeta,
                       std::vector<std::vector<moordyn::real>> f,
                       std::vector<std::vector<vec>> u,
                       std::vector<std::vector<vec>> ud);
 
    real calcSubSeg(unsigned int firstNodeIdx,
                    unsigned int secondNodeIdx,
                    real surfaceHeight);
 
    inline std::pair<real, real> getDrag() const { return make_pair(Cdn, Cdt); }
 
    inline void setDrag(real cdn, real cdt)
    {
        Cdn = cdn;
        Cdt = cdt;
    }
 
    inline void scaleDrag(real scaler)
    {
        Cdn = Cdn * scaler;
        Cdt = Cdt * scaler;
    }
 
    inline void enablePb() { isPb = true; }
 
    inline void disablePb() { isPb = false; }
 
    inline bool enabledPb() const { return isPb; }
 
    void setPin(std::vector<real> p);
 
    inline void setTime(real time) { t = time; }
 
    void setState(const InstanceStateVarView r);
 
    void setEndKinematics(vec r, vec rd, EndPoints end_point);
 
    void setEndOrientation(vec q, EndPoints end_point, EndPoints rod_end_point);
 
    vec getEndSegmentMoment(EndPoints end_point, EndPoints rod_end_point) const;
 
    void getStateDeriv(InstanceStateVarView drdt);
 
    // void initiateStep(vector<double> &rFairIn, vector<double> &rdFairIn,
    // double time);
 
    void Output(real);
 
    inline const size_t stateN() const
    {
        if (ElasticMod != ELASTIC_CONSTANT)
            return getN(); // N rows for viscoelastic case
        else
            return getN() - 1; // N-1 rows for other cases
    }
 
    inline const size_t stateDims() const
    {
        if (Cl > 0 && ElasticMod != ELASTIC_CONSTANT)
            return 8; // 3 for position, 3 for velocity, 1 for VIV phase, 1 for
                      // viscoelasticity
        else if ((Cl > 0) ^ (ElasticMod != ELASTIC_CONSTANT))
            return 7; // 3 for position, 3 for velocity, 1 for VIV phase or
                      // viscoelasticity
        else
            return 6;
    }
 
    std::vector<uint64_t> Serialize(void);
 
    uint64_t* Deserialize(const uint64_t* data);
 
    void saveVTK(const char* filename);
 
    const leanvtk::VTPWriter* getVTK() const { return &vtk; }
private:
    leanvtk::VTPWriter vtk;
};
 
} // ::moordyn