1s_21_rcm.h

/* Copyright (C) 2004-2014 MBSim Development Team
 *
 * 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: thorsten.schindler@mytum.de
 */
 
#ifndef _FINITE_ELEMENT_1S_21_RCM_H_
#define _FINITE_ELEMENT_1S_21_RCM_H_
 
#include "mbsimFlexibleBody/discretization_interface.h"
#include "fmatvec/fmatvec.h"
 
namespace MBSimFlexibleBody {
 
  class FiniteElement1s21RCM : public DiscretizationInterface {
    public:
      FiniteElement1s21RCM(double l0_, double  Arho_, double EA_, double EI_, const fmatvec::Vec &g_);
 
       ~FiniteElement1s21RCM() override = default;
 
      /* INHERITED INTERFACE */
       const fmatvec::SymMat& getM() const override { return M; }
       const fmatvec::Vec& geth() const override { return h; }
       const fmatvec::SqrMat& getdhdq() const override { return Dhq; }    
       const fmatvec::SqrMat& getdhdu() const override { return Dhqp; }
       int getqSize() const override { return 8; }
       int getuSize() const override { return 8; }
       void computeM(const fmatvec::Vec& qElement) override;
       void computeh(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement) override;
       void computedhdz(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement) override;
       double computeKineticEnergy(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement) override;
       double computeGravitationalEnergy(const fmatvec::Vec& qElement) override;
      virtual double computePhysicalStrain(const fmatvec::Vec& qElement);
       double computeElasticEnergy(const fmatvec::Vec& qElement) override;
      virtual fmatvec::Vec3 getPosition(const fmatvec::Vec& qElement, double s);
      virtual fmatvec::SqrMat3 getOrientation(const fmatvec::Vec& qElement, double s);
      virtual fmatvec::Vec3 getVelocity (const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement, double s);
      virtual fmatvec::Vec3 getAngularVelocity(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement, double s);
      virtual fmatvec::Mat getJacobianOfMotion(const fmatvec::Vec& qElement, double s) { return JGeneralized(qElement,s); }
      /***************************************************/
      fmatvec::Vec computeAdditionalElementData(const fmatvec::Vec &qElement, const fmatvec::Vec &qpElement);
 
      /* GETTER / SETTER */
      void setCurlRadius(double R);
      void setMaterialDamping(double depsilons);
      void setLehrDamping(double D);
      /***************************************************/
 
      fmatvec::Vec3 getPositions(const fmatvec::Vec& qElement, double s);
 
      fmatvec::Vec3 getVelocities(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement, double s);
 
      fmatvec::Mat JGeneralizedInternal(const fmatvec::Vec& qElement, double s);
 
      fmatvec::Mat JGeneralized (const fmatvec::Vec& qElement, double s);
 
      fmatvec::Mat JpGeneralized(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement, double s, double sp);
 
      fmatvec::Vec ElementData(fmatvec::Vec qElement, fmatvec::Vec qpElement);
 
    protected:
      double l0, Arho, EA, EI;
 
      double wss0;
 
      double depsilon;
 
      fmatvec::Vec g;
 
      fmatvec::SymMat M;
 
      fmatvec::Vec h;
 
      fmatvec::SqrMat Dhq, Dhqp;
 
      fmatvec::SqrMat Damp;
 
      fmatvec::Vec qLocal, qpLocal, hIntermediate;
 
      fmatvec::SymMat MLocal;
 
      fmatvec::SqrMat Jeg, Jegp;
 
      fmatvec::Vec qElement_Old, qpElement_Old;
 
      double tol_comp;
 
    private:
      void BuildqLocal(const fmatvec::Vec& qGlobal, fmatvec::Vec& qLocal);
 
      void BuildJacobi(const fmatvec::Vec& qLocal, fmatvec::SqrMat& Jeg);
 
      void BuildJacobi(const fmatvec::Vec& qLocal, const fmatvec::Vec& qpIntern, fmatvec::SqrMat& Jeg, fmatvec::SqrMat& Jegp);
 
      fmatvec::Vec3 evalLocalPositions(const fmatvec::Vec& qLocal, double s);
 
      fmatvec::Vec3 evalLocalVelocities(const fmatvec::Vec& qLocal, const fmatvec::Vec& qpLocal, double s);
 
      fmatvec::Mat hFullJacobi(const fmatvec::Vec& qElement, const fmatvec::Vec& qpElement, const fmatvec::Vec& qLocal, const fmatvec::Vec& qpLocal, const fmatvec::SqrMat& Jeg, const fmatvec::SqrMat& Jegp, const fmatvec::SymMat& MLocal, const fmatvec::Vec& hIntermediate);
 
      double l0h2, l0h3, l0h4, l0h5, l0h7, l0h8;
      
      FiniteElement1s21RCM() = default;
  };
 
  inline double Sec(double alpha) { return 1.0/cos(alpha); }
  inline double Power(double base, int exponent) { return pow(base,exponent); }
}
 
#endif /* _FINITE_ELEMENT_1S_21_RCM_H_ */