MoorDyn2.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 "MoorDynAPI.h"
#include "IO.hpp"
#include "Misc.hpp"
 
#include "Time.hpp"
#include "Waves.hpp"
#include "MoorDyn.h"
#include "Line.hpp"
#include "Point.hpp"
#include "Rod.hpp"
#include "Body.hpp"
#include "Seafloor.hpp"
#include <limits>
 
namespace moordyn {
 
class MoorDyn final : public io::IO
{
  public:
    DECLDIR MoorDyn(const char* infilename = NULL,
                    int log_level = MOORDYN_MSG_LEVEL);
 
    DECLDIR ~MoorDyn();
 
    moordyn::error_id DECLDIR Init(const double* x,
                                   const double* xd,
                                   bool skip_ic = false);
 
    moordyn::error_id DECLDIR
    Step(const double* x, const double* xd, double* f, double& t, double& dt);
 
    inline vector<Body*> GetBodies() const { return BodyList; }
 
    inline vector<Rod*> GetRods() const { return RodList; }
 
    inline vector<Point*> GetPoints() const { return PointList; }
 
    inline vector<Line*> GetLines() const { return LineList; }
 
    inline void SetDisableOutput(bool disable) { disableOutput = disable; }
 
    inline unsigned int NCoupledDOF() const
    {
        std::size_t n = 6 * CpldBodyIs.size() + 3 * CpldPointIs.size();
        for (auto rodi : CpldRodIs) {
            if (RodList[rodi]->type == Rod::COUPLED)
                n += 6; // cantilevered rods
            else
                n += 3; // pinned rods
        }
        return static_cast<unsigned int>(n);
    }
 
    inline moordyn::WavesRef DECLDIR GetWaves() const { return waves; }
 
    inline moordyn::SeafloorRef GetSeafloor() const { return seafloor; }
 
    inline unsigned int ExternalWaveKinInit()
    {
        const auto& points = waves->getWaveKinematicsPoints();
        npW = static_cast<unsigned int>(points.size());
 
        return npW;
    }
 
    inline unsigned int ExternalWaveKinGetN() const { return npW; }
 
    inline moordyn::error_id DEPRECATED GetWaveKinCoordinates(double* r) const
    {
        const auto& points = ExternalWaveKinGetPoints();
        unsigned int i = 0;
        for (const auto& vec : points) {
            moordyn::vec2array(vec, r + i);
            i += 3;
        }
 
        return MOORDYN_SUCCESS;
    }
 
    inline std::vector<vec> ExternalWaveKinGetPoints() const
    {
        return waves->getWaveKinematicsPoints();
    }
 
    inline void DEPRECATED SetWaveKin(std::vector<vec> const& U,
                                      std::vector<vec> const& Ud,
                                      double t)
    {
        ExternalWaveKinSet(U, Ud, t);
    }
 
    inline void ExternalWaveKinSet(std::vector<vec> const& U,
                                   std::vector<vec> const& Ud,
                                   double t)
 
    {
        waves->setWaveKinematics(U, Ud);
    }
 
    std::vector<uint64_t> Serialize(void);
 
    uint64_t* Deserialize(const uint64_t* data);
 
    void saveVTK(const char* filename) const;
 
    inline real GetDt() const { return dtM0; }
 
    inline void SetDt(real dt)
    {
        this->dtM0 = dt;
        this->cfl = 0.0;
        for (auto obj : LineList)
            cfl = (std::max)(cfl, obj->dt2cfl(dtM0));
        for (auto obj : PointList)
            cfl = (std::max)(cfl, obj->dt2cfl(dtM0));
        for (auto obj : RodList)
            cfl = (std::max)(cfl, obj->dt2cfl(dtM0));
        for (auto obj : BodyList)
            cfl = (std::max)(cfl, obj->dt2cfl(dtM0));
    }
 
    inline real GetCFL() const { return cfl; }
 
    inline void SetCFL(real cfl)
    {
        this->cfl = cfl;
        this->dtM0 = (std::numeric_limits<real>::max)();
        for (auto obj : LineList)
            dtM0 = (std::min)(dtM0, obj->cfl2dt(cfl));
        for (auto obj : PointList)
            dtM0 = (std::min)(dtM0, obj->cfl2dt(cfl));
        for (auto obj : RodList)
            dtM0 = (std::min)(dtM0, obj->cfl2dt(cfl));
        for (auto obj : BodyList)
            dtM0 = (std::min)(dtM0, obj->cfl2dt(cfl));
    }
 
    inline time::Scheme* GetTimeScheme() const { return _t_integrator; }
 
    inline void SetTimeScheme(time::Scheme* tscheme)
    {
        if (_t_integrator)
            delete _t_integrator;
        _t_integrator = tscheme;
        _t_integrator->SetGround(GroundBody);
        for (auto obj : BodyList)
            _t_integrator->AddBody(obj);
        for (auto obj : RodList)
            _t_integrator->AddRod(obj);
        for (auto obj : PointList)
            _t_integrator->AddPoint(obj);
        for (auto obj : LineList)
            _t_integrator->AddLine(obj);
        _t_integrator->SetCFL(cfl);
        _t_integrator->Init();
    }
 
  protected:
    moordyn::error_id ReadInFile();
 
    moordyn::error_id readFileIntoBuffers(vector<string>& in_txt);
 
    int findStartOfSection(vector<string>& in_txt, vector<string> sectionName);
 
    LineProps* readLineProps(string inputText, int lineNum);
 
    RodProps* readRodProps(string inputText, int lineNum);
 
    Rod* readRod(string inputText, int lineNum);
 
    Body* readBody(string inputText, int lineNum);
 
    void readOptionsLine(vector<string>& in_txt, int index);
 
    bool checkNumberOfEntriesInLine(vector<string> entries,
                                    int supposedNumberOfEntries,
                                    int lineNum);
 
    moordyn::error_id icStationary();
 
    moordyn::error_id icLegacy();
 
    inline moordyn::error_id GetForces(double* f) const
    {
        if (!NCoupledDOF()) {
            if (f)
                _log->Cout(MOORDYN_WRN_LEVEL)
                    << "Warning: Forces have been asked on "
                    << "the coupled entities, but there are no such entities"
                    << std::endl;
            return MOORDYN_SUCCESS;
        }
        if (NCoupledDOF() && !f) {
            _log->Cout(MOORDYN_ERR_LEVEL)
                << "Error: " << __PRETTY_FUNC_NAME__
                << " called with a NULL forces pointer, but there are "
                << NCoupledDOF() << " coupled Degrees Of Freedom" << std::endl;
            return MOORDYN_INVALID_VALUE;
        }
        unsigned int ix = 0;
        for (auto l : CpldBodyIs) {
            // BUG: These conversions will not be needed in the future
            const vec6 f_i = BodyList[l]->getFnet();
            if (BodyList[l]->type == Body::COUPLED) {
                moordyn::vec62array(f_i, f + ix);
                ix += 6; // for coupled bodies 6 entries will be taken
            } else {
                moordyn::vec2array(f_i(Eigen::seqN(0, 3)), f + ix);
                ix += 3; // for pinned bodies 3 entries will be taken
            }
        }
        for (auto l : CpldRodIs) {
            const vec6 f_i = RodList[l]->getFnet();
            if (RodList[l]->type == Rod::COUPLED) {
                moordyn::vec62array(f_i, f + ix);
                ix += 6; // for cantilevered rods 6 entries will be taken
            } else {
                moordyn::vec2array(f_i(Eigen::seqN(0, 3)), f + ix);
                ix += 3; // for pinned rods 3 entries will be taken
            }
        }
        for (auto l : CpldPointIs) {
            vec fnet;
            PointList[l]->getFnet(fnet);
            moordyn::vec2array(fnet, f + ix);
            ix += 3;
        }
        return MOORDYN_SUCCESS;
    }
 
  private:
    string _filepath;
    string _basename;
    string _basepath;
 
    // factor by which to boost drag coefficients during dynamic relaxation IC
    // generation
    real ICDfac;
    // convergence analysis time step for IC generation
    real ICdt;
    // max time for IC generation
    real ICTmax;
    // threshold for relative change in tensions to call it converged
    real ICthresh;
    // use dynamic (true) or stationary (false) initial condition solver
    bool ICgenDynamic;
    // Load the initial state (before the initial condition solver) from a
    // file. Empty to do not load it.
    string ICfile;
    // temporary wave kinematics flag used to store input value while keeping
    // env.WaveKin=0 for IC gen
    moordyn::waves::waves_settings WaveKinTemp;
    real dtM0;
    real cfl;
    real dtOut;
 
    time::Scheme* _t_integrator;
 
    EnvCondRef env;
    Body* GroundBody;
    WavesRef waves{};
    moordyn::SeafloorRef seafloor;
 
    vector<LineProps*> LinePropList;
    vector<RodProps*> RodPropList;
    vector<FailProps*> FailList;
    vector<Body*> BodyList;
    vector<Rod*> RodList;
    vector<moordyn::Point*> PointList;
    vector<moordyn::Line*> LineList;
 
    vector<unsigned int> FreeBodyIs;
    vector<unsigned int> FixedBodyIs;
    vector<unsigned int> CpldBodyIs;
 
    vector<unsigned int> FreeRodIs;
    vector<unsigned int> CpldRodIs;
 
    vector<unsigned int> FreePointIs;
    vector<unsigned int> CpldPointIs;
 
    unsigned int npW;
 
    bool disableOutput = false;
 
    bool disableOutTime = false;
 
    ofstream outfileMain;
 
    vector<shared_ptr<ofstream>> outfiles;
 
    vector<OutChanProps> outChans;
 
    inline moordyn::error_id SetupLog()
    {
        // env.writeLog = 0 -> MOORDYN_NO_OUTPUT
        // env.writeLog = 1 -> MOORDYN_WRN_LEVEL
        // env.writeLog = 2 -> MOORDYN_MSG_LEVEL
        // env.writeLog >= 3 -> MOORDYN_DBG_LEVEL
        int log_level = MOORDYN_ERR_LEVEL - env->writeLog;
        if (log_level >= MOORDYN_ERR_LEVEL)
            log_level = MOORDYN_NO_OUTPUT;
        if (log_level < MOORDYN_DBG_LEVEL)
            log_level = MOORDYN_DBG_LEVEL;
        GetLogger()->SetLogLevel(log_level);
 
        if (env->writeLog > 0) {
            moordyn::error_id err = MOORDYN_SUCCESS;
            string err_msg;
            stringstream filepath;
            filepath << _basepath << _basename << ".log";
            try {
                GetLogger()->SetFile(filepath.str().c_str());
            }
            MOORDYN_CATCHER(err, err_msg);
            if (err != MOORDYN_SUCCESS)
                LOGERR << "Unable to create the log at '" << filepath.str()
                       << "': " << std::endl
                       << err_msg << std::endl;
            else
                LOGMSG << "MoorDyn v2 log file with output level "
                       << log_level_name(GetLogger()->GetLogLevel()) << " at '"
                       << filepath.str() << "'" << std::endl;
            return err;
        }
 
        return MOORDYN_SUCCESS;
    }
 
    moordyn::error_id read_curve(const char* entry,
                                 vector<double>& x,
                                 vector<double>& y)
    {
        try {
            y.push_back((real)std::stold(entry));
            // If we reach this point, it is a valid number
            x.push_back(0.0);
            return MOORDYN_SUCCESS;
        } catch (std::out_of_range) {
            LOGERR << "" << std::endl;
            return MOORDYN_INVALID_INPUT;
        } catch (std::invalid_argument) {
            // Do nothing, just proceed to read the curve file
        }
 
        string fpath = _basepath + entry;
        LOGMSG << "Loading a curve from '" << fpath << "'..." << std::endl;
        ifstream f(fpath);
        if (!f.is_open()) {
            LOGERR << "Cannot read the file '" << fpath << "'" << std::endl;
            return MOORDYN_INVALID_INPUT_FILE;
        }
 
        vector<string> flines;
        int i = 0;
        while (f.good()) {
            string fline;
            getline(f, fline);
            if (i > 2) { // skip first three lines as headers
                moordyn::str::rtrim(fline);
                flines.push_back(fline);
            }
            i++;
        }
        f.close();
 
        if (i < 5) {
            LOGERR << "Error: Not enough curve data in curve file" << endl;
            return MOORDYN_INVALID_INPUT;
        }
 
        for (auto fline : flines) {
            vector<string> entries = moordyn::str::split(fline, ' ');
            if (entries.size() < 2) {
                LOGERR << "Error: Bad curve point" << std::endl
                       << "\t'" << fline << "'" << std::endl
                       << "\t2 fields required, but just " << entries.size()
                       << " are provided" << std::endl;
                return MOORDYN_INVALID_INPUT;
            }
            x.push_back(atof(entries[0].c_str()));
            y.push_back(atof(entries[1].c_str()));
            LOGDBG << "(" << x.back() << ", " << y.back() << ")" << std::endl;
        }
 
        LOGMSG << (i - 3) << " lines of curve successfully loaded" << std::endl;
        return MOORDYN_SUCCESS;
    }
 
    moordyn::error_id read_curve(const char* entry,
                                 double* c,
                                 int* n,
                                 double* x,
                                 double* y)
    {
        vector<double> xv, yv;
        const moordyn::error_id error = read_curve(entry, xv, yv);
        if (error != MOORDYN_SUCCESS)
            return error;
 
        if (xv.size() == 1) {
            *c = yv.back();
            n = 0;
            return MOORDYN_SUCCESS;
        }
 
        if (xv.size() > nCoef) {
            _log->Cout(MOORDYN_ERR_LEVEL)
                << "Error: Too much points in the curve" << std::endl
                << "\t" << xv.size() << " points given, but just " << nCoef
                << " are accepted" << std::endl;
            return MOORDYN_INVALID_INPUT;
        }
 
        *c = 0.0;
        *n = static_cast<unsigned int>(xv.size());
        memcpy(x, xv.data(), xv.size() * sizeof(double));
        memcpy(y, yv.data(), yv.size() * sizeof(double));
 
        return MOORDYN_SUCCESS;
    }
 
    moordyn::error_id read_syrope_working_curves(const char* entry,
                                                 string& owc_path,
                                                 string& wc_formula,
                                                 double& k1,
                                                 double& k2);
 
    inline double GetOutput(const OutChanProps channel) const
    {
        if (channel.OType == 1)
            return LineList[channel.ObjID - 1]->GetLineOutput(channel);
        else if (channel.OType == 2)
            return PointList[channel.ObjID - 1]->GetPointOutput(channel);
        else if (channel.OType == 3)
            return RodList[channel.ObjID - 1]->GetRodOutput(channel);
        else if (channel.OType == 4)
            return BodyList[channel.ObjID - 1]->GetBodyOutput(channel);
        stringstream s;
        s << "Error: output type of " << channel.Name
          << " does not match a supported object type";
        MOORDYN_THROW(MOORDYN_INVALID_VALUE, s.str().c_str());
        return 0.0;
    }
 
    void detachLines(FailProps* failure);
 
    moordyn::error_id WriteOutputs(double t, double dt);
};
 
} // ::moordyn