RSVScalc.hpp

Go to the documentation of this file.

 
#ifndef RSVSCALC_H_INCLUDED
#define RSVSCALC_H_INCLUDED
 
//=================================
// forward declared dependencies
//      class foo; //when you only need a pointer not the actual object
//      and to avoid circular dependencies
 
namespace param
{
namespace dev
{
class devparam;
}
} // namespace param
 
//=================================
// included dependencies
#include <Eigen>
#include <fstream>
#include <iostream>
#include <vector>
 
#include "SnakeVelocityCalculator.hpp"
#include "arraystructures.hpp"
#include "mesh.hpp"
#include "snake.hpp"
#include "triangulate.hpp"
#include "warning.hpp"
 
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> MatrixXd_sparse;
typedef HashedVectorPair<int, int, int, double> TripletMap;
 
template <typename T> class sparsetripletelement
{
    int rowVal = 0;
    int colVal = 0;
    T valueVal = 0;
 
  public:
    sparsetripletelement<T>(int r, int c, T v)
    {
#ifdef SAFE_ACCESS
        if (r < 0 || c < 0)
        {
            RSVS3D_ERROR_RANGE("Indices out of range.");
        }
#endif
        this->rowVal = r;
        this->colVal = c;
        this->valueVal = v;
    }
    int row() const
    {
        return this->rowVal;
    }
    int col() const
    {
        return this->colVal;
    }
    T value() const
    {
        return this->valueVal;
    }
};
 
class SparseMatrixTriplet : protected TripletMap
{
    int nRow = 0;
    int nCol = 0;
 
  public:
    double &operator()(int a, int b)
    {
#ifdef SAFE_ACCESS
        if (a < 0 || a >= this->nRow || b < 0 || b >= this->nCol)
        {
            RSVS3D_ERROR_RANGE("Indices out of range.");
        }
#endif
        return TripletMap::operator()(a + this->nRow * b);
    }
    const double &operator()(int a, int b) const
    {
#ifdef SAFE_ACCESS
        if (a < 0 || a >= this->nRow || b < 0 || b >= this->nCol)
        {
            RSVS3D_ERROR_RANGE("Indices out of range.");
        }
#endif
        return TripletMap::operator()(a + this->nRow * b);
    }
    sparsetripletelement<double> operator[](int a);
    double &coeffRef(int a, int b)
    {
        return this->operator()(a, b);
    }
    void resize(int row, int col)
    {
        this->nRow = row;
        this->nCol = col;
    }
    void setZero()
    {
        this->TripletMap::clear();
    }
    int nonZeros() const
    {
        return this->TripletMap::size();
    }
    void reserve(size_t a);
    void SetEqual(MatrixXd_sparse &targ);
    void SetEqual(Eigen::MatrixXd &targ);
};
 
class RSVScalc : public rsvs3d::SnakeVelocityCalculator
{
  protected:
    int nDv = 0;
    int nConstr = 0;
    int falseaccess = 0;
    int sparseDVcutoff = 200;
    int nonZeroPerDV = 10;
    bool returnDeriv = true;
    bool useSurfCentreDeriv = true;
    bool useSurfCentreHessian = false;
 
    int timer1 = 0;
    int timer2 = 0;
    int timer3 = 0;
    void ZeroTimers()
    {
        this->timer1 = 0;
        this->timer2 = 0;
        this->timer3 = 0;
    }
    void PrintTimers() const;
 
  public:
    Eigen::MatrixXd dConstr;
    SparseMatrixTriplet dConstr_sparse;
    Eigen::MatrixXd HConstr;
    SparseMatrixTriplet HConstr_sparse;
    Eigen::MatrixXd HObj;
    SparseMatrixTriplet HObj_sparse;
    Eigen::MatrixXd HLag;
    MatrixXd_sparse HLag_sparse;
    Eigen::RowVectorXd dLag;
    Eigen::RowVectorXd dObj;
    Eigen::VectorXd constr;
    Eigen::VectorXd lagMult;
    Eigen::VectorXd deltaDV;
    Eigen::VectorXd constrTarg;
    Eigen::MatrixXd dvCallConstr;
    Eigen::MatrixXd sensDv;
    double obj = 0.0;
    double limLag = INFINITY;
    std::vector<bool> isConstrAct;
    std::vector<bool> isDvAct;
    HashedVector<int, int> subConstrAct;
    HashedVector<int, int> subDvAct;
    HashedVector<int, int> dvMap;
    HashedMap<int, int, int> constrMap;
    std::vector<std::pair<int, int>> constrList;
 
    bool UseFullMath() const
    {
        return (this->nDv < this->sparseDVcutoff);
    }
    void BuildMathArrays(int nDv, int nConstr);
 
    void BuildConstrMap(const triangulation &triangleRSVS);
 
    void BuildConstrMap(const mesh &meshin);
 
    int BuildDVMap(const std::vector<int> &vecin);
 
    bool SnakDVcond(const triangulation &triRSVS, int ii);
 
    void PrepTriangulationCalc(const triangulation &triRSVS);
 
    void CalculateMesh(mesh &meshin);
 
    void CalculateTriangulation(const triangulation &triRSVS, int derivMethod = 0);
 
    void CalcTriangle(const triangle &triIn, const triangulation &triRSVS, bool isObj = true, bool isConstr = true,
                      bool isDeriv = true);
 
    void CalcTriangleFD(const triangle &triIn, const triangulation &triRSVS, bool isObj = true, bool isConstr = true,
                        bool isDeriv = true);
 
    void CalcTriangleDirectVolume(const triangle &triIn, const triangulation &triRSVS, bool isObj = true,
                                  bool isConstr = true, bool isDeriv = true);
 
    void CalcTriangleEdgeLength(const triangle &triIn, const triangulation &triRSVS, bool isObj = true,
                                bool isConstr = true, bool isDeriv = true);
 
    void ReturnConstrToMesh(triangulation &triRSVS) const;
 
    void ReturnConstrToMesh(mesh &meshin, double volu::*mp = &volu::volume) const;
 
    void CheckAndCompute(int calcMethod = 0, bool sensCalc = false);
 
    void ComputeSQPstep(int calcMethod, Eigen::MatrixXd &dConstrAct, Eigen::RowVectorXd &dObjAct,
                        Eigen::VectorXd &constrAct, Eigen::VectorXd &lagMultAct);
    void ComputeSQPstep(int calcMethod, MatrixXd_sparse &dConstrAct, Eigen::RowVectorXd &dObjAct,
                        Eigen::VectorXd &constrAct, Eigen::VectorXd &lagMultAct);
 
    void ComputeSQPsens(int calcMethod, const Eigen::MatrixXd &sensMult, const Eigen::MatrixXd &sensInv,
                        Eigen::MatrixXd &sensRes);
    void ComputeSQPsens(int calcMethod, Eigen::MatrixXd &sensMult, MatrixXd_sparse &sensInv, Eigen::MatrixXd &sensRes);
 
    bool PrepareMatricesForSQP(Eigen::MatrixXd &dConstrAct, Eigen::MatrixXd &HConstrAct, Eigen::MatrixXd &HObjAct,
                               MatrixXd_sparse &dConstrAct_sparse, MatrixXd_sparse &HConstrAct_sparse,
                               MatrixXd_sparse &HObjAct_sparse, Eigen::RowVectorXd &dObjAct, Eigen::VectorXd &constrAct,
                               Eigen::VectorXd &lagMultAct);
 
    void PrepareMatricesForSQPFull(Eigen::MatrixXd &dConstrAct, Eigen::MatrixXd &HConstrAct, Eigen::MatrixXd &HObjAct);
    void PrepareMatricesForSQPSparse(MatrixXd_sparse &dConstrAct_sparse, MatrixXd_sparse &HConstrAct_sparse,
                                     MatrixXd_sparse &HObjAct_sparse);
    bool PrepareMatricesForSQPSensitivity(const Eigen::MatrixXd &dConstrAct, const Eigen::MatrixXd &HConstrAct,
                                          const Eigen::MatrixXd &HObjAct, Eigen::MatrixXd &sensMult,
                                          Eigen::MatrixXd &sensInv, Eigen::MatrixXd &sensRes) const;
    bool PrepareMatricesForSQPSensitivity(const MatrixXd_sparse &dConstrAct, const MatrixXd_sparse &HConstrAct,
                                          MatrixXd_sparse &HObjAct, Eigen::MatrixXd &sensMult, MatrixXd_sparse &sensInv,
                                          Eigen::MatrixXd &sensRes) const;
 
    void ReturnVelocities(triangulation &triRSVS);
    void ReturnSensitivities(const triangulation &triRSVS, std::vector<double> &sensVec, int constrNum) const;
    void ReturnGradient(const triangulation &triRSVS, std::vector<double> &sensVec, int constrNum) const;
    int numConstr() const
    {
        return (this->nConstr);
    }
    // Output functions
 
    void Print2Screen(int outType = 0) const;
 
    void ConvergenceLog(std::ofstream &out, int loglvl = 3) const;
    void CalculateVelocities(triangulation &triRSVS, int calculationMethod = 0, bool calculateDerivatives = true,
                             int derivativeMethod = 1) override;
    void setDevParameters(const param::dev::devparam &devset) override;
    void SetUseSurfCentreDeriv(int in)
    {
        this->useSurfCentreDeriv = in;
    }
    bool GetUseSurfCentreDeriv() const
    {
        return (this->useSurfCentreDeriv);
    }
    void SetUseSurfCentreHessian(int in)
    {
        this->useSurfCentreHessian = in;
    }
    bool SetUseSurfCentreHessian() const
    {
        return (this->useSurfCentreHessian);
    }
    void SetSparseDVcutoff(int in)
    {
        this->sparseDVcutoff = in;
    }
    int GetSparseDVcutoff() const
    {
        return (this->sparseDVcutoff);
    }
};
 
#endif