newmark_method.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:            JMCarbonell $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:            November 2017 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(KRATOS_NEWMARK_METHOD_H_INCLUDED)
#define  KRATOS_NEWMARK_METHOD_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_solvers/time_integration_methods/time_integration_method.hpp"
 
namespace Kratos
{
 
 
 
 
 
 
 
 
  template<class TVariableType, class TValueType>
  class KRATOS_API(SOLID_MECHANICS_APPLICATION) NewmarkMethod : public TimeIntegrationMethod<TVariableType,TValueType>
  {
  protected:
 
    struct NewmarkParameters
    {
      double beta;
      double gamma;
 
      double delta_time;
 
      //system constants
      double c0;
      double c1;
      double c2;
      double c3;
      double c4;
      double c5;
 
      void SetParameters(const double& rbeta,
             const double& rgamma,
             const double& rdelta_time)
      {
    beta  = rbeta;
    gamma = rgamma;
 
    delta_time = rdelta_time;
 
    c0 = ( 1.0 / (beta * delta_time * delta_time) );
        c1 = ( gamma / (beta * delta_time) );
        c2 = ( 1.0 / (beta * delta_time) );
        c3 = ( 0.5 / (beta) - 1.0 );
        c4 = ( (gamma / beta) - 1.0  );
        c5 = ( delta_time * 0.5 * ( ( gamma / beta ) - 2.0 ) );
      }
 
 
    private:
 
      friend class Serializer;
 
      void save(Serializer& rSerializer) const
      {
    rSerializer.save("beta", beta);
    rSerializer.save("gamma", gamma);
    rSerializer.save("delta_time", delta_time);
    rSerializer.save("c0", c0);
    rSerializer.save("c1", c1);
    rSerializer.save("c2", c2);
    rSerializer.save("c3", c3);
    rSerializer.save("c4", c4);
    rSerializer.save("c5", c5);
      };
 
      void load(Serializer& rSerializer)
      {
    rSerializer.load("beta", beta);
    rSerializer.load("gamma", gamma);
    rSerializer.load("delta_time", delta_time);
    rSerializer.load("c0", c0);
    rSerializer.load("c1", c1);
    rSerializer.load("c2", c2);
    rSerializer.load("c3", c3);
    rSerializer.load("c4", c4);
    rSerializer.load("c5", c5);
      };
 
    };
 
  public:
 
 
    typedef TimeIntegrationMethod<TVariableType,TValueType>  BaseType;
 
    typedef typename BaseType::Pointer                BasePointerType;
 
    typedef typename BaseType::NodeType                      NodeType;
 
    typedef typename BaseType::VariablePointer        VariablePointer;
 
    KRATOS_CLASS_POINTER_DEFINITION( NewmarkMethod );
 
 
 
    NewmarkMethod() : BaseType() {}
 
    NewmarkMethod(const TVariableType& rVariable) : BaseType(rVariable) {}
 
    NewmarkMethod(const TVariableType& rVariable, const TVariableType& rFirstDerivative, const TVariableType& rSecondDerivative) : BaseType(rVariable,rFirstDerivative,rSecondDerivative) {}
 
    NewmarkMethod(const TVariableType& rVariable, const TVariableType& rFirstDerivative, const TVariableType& rSecondDerivative, const TVariableType& rPrimaryVariable) : BaseType(rVariable,rFirstDerivative,rSecondDerivative,rPrimaryVariable) {}
 
    NewmarkMethod(NewmarkMethod& rOther)
      :BaseType(rOther)
      ,mNewmark(rOther.mNewmark)
    {
    }
 
    BasePointerType Clone() override
    {
      return BasePointerType( new NewmarkMethod(*this) );
    }
 
    ~NewmarkMethod() override{}
 
 
 
    //calculate parameters (to call it once with the original input parameters)
    void CalculateParameters(ProcessInfo& rCurrentProcessInfo) override
    {
     KRATOS_TRY
 
     double beta = 0.25;
     if (rCurrentProcessInfo.Has(NEWMARK_BETA))
       {
     beta = rCurrentProcessInfo[NEWMARK_BETA];
       }
 
     double gamma = 0.5;
     if (rCurrentProcessInfo.Has(NEWMARK_GAMMA))
       {
     gamma = rCurrentProcessInfo[NEWMARK_GAMMA];
       }
 
     rCurrentProcessInfo[NEWMARK_BETA]  = beta;
     rCurrentProcessInfo[NEWMARK_GAMMA] = gamma;
 
     this->SetParameters(rCurrentProcessInfo);
 
     KRATOS_CATCH( "" )
    }
 
    // set parameters (do not calculate parameters here, only read them)
    void SetParameters(const ProcessInfo& rCurrentProcessInfo) override
    {
     KRATOS_TRY
 
     double delta_time = rCurrentProcessInfo[DELTA_TIME];
 
     if (delta_time < 1.0e-24)
        {
      KRATOS_ERROR << " ERROR: detected delta_time = 0 in the Solution Method DELTA_TIME. PLEASE : check if the time step is created correctly for the current model part " << std::endl;
        }
 
     double beta = 0.25;
     if (rCurrentProcessInfo.Has(NEWMARK_BETA))
       {
     beta = rCurrentProcessInfo[NEWMARK_BETA];
       }
     double gamma = 0.5;
     if (rCurrentProcessInfo.Has(NEWMARK_GAMMA))
       {
     gamma = rCurrentProcessInfo[NEWMARK_GAMMA];
       }
 
     mNewmark.SetParameters(beta,gamma,delta_time);
 
     KRATOS_CATCH( "" )
    }
 
    // set parameters to process info
    void SetProcessInfoParameters(ProcessInfo& rCurrentProcessInfo) override
    {
     KRATOS_TRY
 
     rCurrentProcessInfo[NEWMARK_BETA]  = mNewmark.beta;
     rCurrentProcessInfo[NEWMARK_GAMMA] = mNewmark.gamma;
 
     KRATOS_CATCH( "" )
    }
 
 
    int Check( const ProcessInfo& rCurrentProcessInfo ) override
    {
      KRATOS_TRY
 
      // Perform base integration method checks
      int ErrorCode = 0;
      ErrorCode = BaseType::Check(rCurrentProcessInfo);
 
      // Check that all required variables have been registered
      if( this->mpFirstDerivative == nullptr ){
        KRATOS_ERROR << " time integration method FirstDerivative not set " <<std::endl;
      }
 
      if( this->mpSecondDerivative == nullptr ){
        KRATOS_ERROR << " time integration method SecondDerivative not set " <<std::endl;
      }
 
      return ErrorCode;
 
      KRATOS_CATCH("")
    }
 
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "NewmarkMethod";
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "NewmarkMethod";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
      rOStream << "NewmarkMethod Data";
    }
 
 
 
 
 
  protected:
 
 
 
    NewmarkParameters   mNewmark;
 
 
    void AssignFromVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      // predict variable from variable
      TValueType& CurrentFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      TValueType& CurrentSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentFirstDerivative  -= CurrentFirstDerivative;
      CurrentSecondDerivative -= CurrentSecondDerivative;
 
 
      KRATOS_CATCH( "" )
    }
 
    void AssignFromFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      // predict variable from first derivative
      TValueType& CurrentVariable                = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      TValueType& CurrentFirstDerivative         = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
 
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      // newmark consistent
      CurrentVariable = PreviousVariable + this->mNewmark.delta_time * ( ( 1.0 - (this->mNewmark.beta/this->mNewmark.gamma) ) * PreviousFirstDerivative + (this->mNewmark.beta/this->mNewmark.gamma) * CurrentFirstDerivative +  this->mNewmark.delta_time *  0.5 * ( 1.0 - 2.0 * (this->mNewmark.beta/this->mNewmark.gamma) ) * PreviousSecondDerivative );
      // uniform accelerated movement
      //CurrentVariable = PreviousVariable + 0.5 * mNewmark.delta_time * (PreviousFirstDerivative + CurrentFirstDerivative) + 0.5 * mNewmark.delta_time * mNewmark.delta_time * PreviousSecondDerivative;
 
      TValueType& CurrentSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      // variable prediction :: uniform accelerated movement **
      CurrentVariable -= this->mNewmark.delta_time * PreviousFirstDerivative;
 
      //CurrentFirstDerivative   = PreviousFirstDerivative;
      CurrentSecondDerivative -= CurrentSecondDerivative;
 
      //std::cout<<*this->mpVariable<<" Assign Node["<<rNode.Id()<<"]"<<CurrentVariable<<" "<<CurrentFirstDerivative<<" "<<CurrentSecondDerivative<<std::endl;
 
      KRATOS_CATCH( "" )
    }
 
    void AssignFromSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      // predict variable from second derivative
      TValueType& CurrentVariable                = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
 
      TValueType& CurrentSecondDerivative        = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      CurrentVariable = PreviousVariable + this->mNewmark.delta_time * PreviousFirstDerivative + this->mNewmark.delta_time * this->mNewmark.delta_time * ( 0.5 * ( 1.0 - 2.0 * this->mNewmark.beta ) * PreviousSecondDerivative + this->mNewmark.beta * CurrentSecondDerivative );
 
      // variable prediction :: uniform accelerated movement **
      CurrentVariable -= this->mNewmark.delta_time * PreviousFirstDerivative + 0.5 * this->mNewmark.delta_time * this->mNewmark.delta_time * PreviousSecondDerivative;
 
      //CurrentSecondDerivative = PreviousSecondDerivative;
 
      KRATOS_CATCH( "" )
    }
 
    void PredictFromVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      //if(this->Is(TimeIntegrationLocalFlags::NOT_PREDICT_PRIMARY_VARIABLE))
      this->PredictVariable(rNode);
      this->PredictFirstDerivative(rNode);
      this->PredictSecondDerivative(rNode);
 
      KRATOS_CATCH( "" )
    }
 
 
    void PredictFromFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      //if(this->Is(TimeIntegrationLocalFlags::NOT_PREDICT_PRIMARY_VARIABLE))
      //this->PredictFirstDerivative(rNode);
      this->PredictSecondDerivative(rNode);
      this->PredictVariable(rNode);
 
      KRATOS_CATCH( "" )
    }
 
    void PredictVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      TValueType& CurrentVariable               = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& CurrentFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      const TValueType& CurrentSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      // variable prediction :: uniform accelerated movement **
      CurrentVariable += this->mNewmark.delta_time * CurrentFirstDerivative + 0.5 * this->mNewmark.delta_time * this->mNewmark.delta_time * CurrentSecondDerivative;
 
      KRATOS_CATCH( "" )
    }
 
    void PredictFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      const TValueType& CurrentVariable          = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      TValueType& CurrentFirstDerivative         = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
 
      // consistent newmark
      CurrentFirstDerivative = this->mNewmark.c1 * (CurrentVariable-PreviousVariable) - this->mNewmark.c4 * PreviousFirstDerivative - this->mNewmark.c5 * PreviousSecondDerivative;
 
      // uniform accelerated movement
      // CurrentFirstDerivative = (2.0/this->mNewmark.delta_time) * (CurrentVariable-PreviousVariable) - this->mNewmark.delta_time * CurrentSecondDerivative - CurrentFirstDerivative;
 
 
      KRATOS_CATCH( "" )
    }
 
    void PredictSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      const TValueType& CurrentFirstDerivative   = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      TValueType& CurrentSecondDerivative        = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      CurrentSecondDerivative = ( 1.0 / (this->mNewmark.gamma * this->mNewmark.delta_time) ) * ( CurrentFirstDerivative - PreviousFirstDerivative - ( 1.0 - this->mNewmark.gamma ) * this->mNewmark.delta_time * PreviousSecondDerivative );
 
      KRATOS_CATCH( "" )
    }
 
 
    void UpdateFromVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      this->UpdateVariable(rNode);
      this->UpdateFirstDerivative(rNode);
      this->UpdateSecondDerivative(rNode);
 
      KRATOS_CATCH( "" )
    }
 
    void UpdateFromFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      TValueType& CurrentVariable                = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      TValueType& CurrentSecondDerivative        = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
      const TValueType& CurrentFirstDerivative   = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
 
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      CurrentSecondDerivative = (this->mNewmark.c0/this->mNewmark.c1) * (CurrentFirstDerivative + this->mNewmark.c4 * PreviousFirstDerivative + this->mNewmark.c5 * PreviousSecondDerivative) - this->mNewmark.c2 * PreviousFirstDerivative - this->mNewmark.c3 * PreviousSecondDerivative;
 
      CurrentVariable = PreviousVariable + (1.0/this->mNewmark.c1) * (CurrentFirstDerivative + this->mNewmark.c4 * PreviousFirstDerivative + this->mNewmark.c5 * PreviousSecondDerivative);
 
      KRATOS_CATCH( "" )
    }
 
 
    void UpdateFromSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      KRATOS_ERROR << " Calling UpdateFromSecondDerivative for Newmark time integration method : NOT IMPLEMENTED " <<std::endl;
 
      KRATOS_CATCH( "" )
    }
 
    void UpdateVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      KRATOS_CATCH( "" )
    }
 
 
    void UpdateFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      const TValueType& CurrentVariable          = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      TValueType& CurrentFirstDerivative         = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
 
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      CurrentFirstDerivative = (this->mNewmark.c1 * (CurrentVariable-PreviousVariable) - this->mNewmark.c4 * PreviousFirstDerivative - this->mNewmark.c5 * PreviousSecondDerivative);
 
      KRATOS_CATCH( "" )
    }
 
    void UpdateSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      const TValueType& CurrentVariable          = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      TValueType& CurrentSecondDerivative        = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
      const TValueType& PreviousSecondDerivative = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      CurrentSecondDerivative = (this->mNewmark.c0 * (CurrentVariable-PreviousVariable) - this->mNewmark.c2 * PreviousFirstDerivative - this->mNewmark.c3 * PreviousSecondDerivative);
 
      KRATOS_CATCH( "" )
    }
 
 
 
    // get parameters
    double& GetFirstDerivativeInertialParameter(double& rParameter) override
    {
      rParameter = mNewmark.c1;
      return rParameter;
    }
 
    double& GetSecondDerivativeInertialParameter(double& rParameter) override
    {
      rParameter = mNewmark.c0;
      return rParameter;
    }
 
 
 
 
  private:
 
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
      KRATOS_SERIALIZE_SAVE_BASE_CLASS( rSerializer, BaseType )
      rSerializer.save("NewmarkParameters", mNewmark);
    };
 
    void load(Serializer& rSerializer) override
    {
      KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, BaseType )
      rSerializer.load("NewmarkParameters", mNewmark);
    };
 
 
 
 
  }; // Class NewmarkMethod
 
 
 
 
 
  template<class TVariableType, class TValueType>
  inline std::istream & operator >> (std::istream & rIStream, NewmarkMethod<TVariableType,TValueType>& rThis)
  {
    return rIStream;
  }
 
  template<class TVariableType, class TValueType>
  inline std::ostream & operator << (std::ostream & rOStream, const NewmarkMethod<TVariableType,TValueType>& rThis)
  {
    return rOStream << rThis.Info();
  }
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_NEWMARK_METHOD_H_INCLUDED defined