phem56_integrator.h

/* Copyright (C) 2004-2018  Martin Förg
 
 * This library is free software; you can redistribute it and/or 
 * modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation; either 
 * version 2.1 of the License, or (at your option) any later version. 
 *  
 * This library is distributed in the hope that it will be useful, 
 * but WITHOUT ANY WARRANTY; without even the implied warranty of 
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
 * Lesser General Public License for more details. 
 *  
 * You should have received a copy of the GNU Lesser General Public 
 * License along with this library; if not, write to the Free Software 
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA
 
 *
 * Contact:
 *   martin.o.foerg@googlemail.com
 *
 */
 
#ifndef _PHEM56_INTEGRATOR_H_
#define _PHEM56_INTEGRATOR_H_
 
#include "root_finding_integrator.h"
 
namespace MBSim {
 
  class PHEM56Integrator : public RootFindingIntegrator {
 
    public:
 
      enum LinearAlgebra {
        DEC=0,
        DGETRF,
        unknown
      };
 
    private:
 
      static void fprob(int* ifcn, int* nq, int* nu, int* nx, int* nla, int* nzg, int* nzf, int* lrda, int* nblk, int* nmrc, int* npgp, int* npfl, int* indgr, int* indgc, int* indflr, int* indflc,  double* t, double* q, double* u, double* x, double* xl, double* g, double* WT_, double* f, double* wb, double* deltagd_, double* V_, double* qd_, double* xd_, double* M_);
      static void solout(int* nr, int* nq, int* nu, int* nx, int* nla, int* lrdo, double* q, double* u, double* x, double* a, double* rlam, double* dowk, int* irtrn);
 
      double tPlot{0};
      double dtOut{0};
      double s0; 
      double time{0};
 
      LinearAlgebra linearAlgebra{DGETRF};
      bool generalVMatrix{true};
      bool initialProjection{false};
      int numberOfStepsBetweenProjections{0};
      bool projectOntoIndex1ConstraintManifold{false};
      fmatvec::Vec aTol;
      fmatvec::Vec rTol;
      double dt0{1e-3};
      int maxSteps{std::numeric_limits<int>::max()};
      double dtMax{0};
 
      fmatvec::RangeV Rv;
 
      bool drift { false };
 
      std::exception_ptr exception;
 
    public:
 
      ~PHEM56Integrator() override = default;
 
      void setAbsoluteTolerance(const fmatvec::Vec &aTol_) { aTol <<= aTol_; }
      void setAbsoluteTolerance(double aTol_) { aTol.resize(1,fmatvec::INIT,aTol_); }
      void setRelativeTolerance(const fmatvec::Vec &rTol_) { rTol <<= rTol_; }
      void setRelativeTolerance(double rTol_) { rTol.resize(1,fmatvec::INIT,rTol_); }
      void setInitialStepSize(double dt0_) { dt0 = dt0_; }
      void setMaximumStepSize(double dtMax_) { dtMax = dtMax_; }
      void setStepLimit(int maxSteps_) { maxSteps = maxSteps_; }
      void setLinearAlgebra(LinearAlgebra linearAlgebra_) { linearAlgebra = linearAlgebra_; }
      void setGeneralVMatrix(bool generalVMatrix_) { generalVMatrix = generalVMatrix_; }
      void setInitialProjection(bool initialProjection_) { initialProjection = initialProjection_; }
      void setNumberOfStepsBetweenProjections(int numberOfStepsBetweenProjections_) { numberOfStepsBetweenProjections = numberOfStepsBetweenProjections_; }
      void setProjectOntoIndex1ConstraintManifold(bool projectOntoIndex1ConstraintManifold_) { projectOntoIndex1ConstraintManifold = projectOntoIndex1ConstraintManifold_; }
 
      const fmatvec::Vec& getAbsoluteTolerance() const { return aTol; }
      const fmatvec::Vec& getRelativeTolerance() const { return rTol; }
 
      using Integrator::integrate;
      void integrate() override;
 
      void initializeUsingXML(xercesc::DOMElement *element) override;
  };
 
}
 
#endif