backward_euler_method.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:            JMCarbonell $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:               April 2018 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(KRATOS_BACKWARD_EULER_METHOD_H_INCLUDED)
#define  KRATOS_BACKWARD_EULER_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 BackwardEulerMethod : public TimeIntegrationMethod<TVariableType,TValueType>
  {
  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( BackwardEulerMethod );
 
 
 
    BackwardEulerMethod() : BaseType() {}
 
    BackwardEulerMethod(const TVariableType& rVariable) : BaseType(rVariable) {}
 
    BackwardEulerMethod(const TVariableType& rVariable, const TVariableType& rFirstDerivative, const TVariableType& rSecondDerivative) : BaseType(rVariable,rFirstDerivative,rSecondDerivative) {}
 
    BackwardEulerMethod(const TVariableType& rVariable, const TVariableType& rFirstDerivative, const TVariableType& rSecondDerivative, const TVariableType& rPrimaryVariable) : BaseType(rVariable,rFirstDerivative,rSecondDerivative,rPrimaryVariable) {}
 
    BackwardEulerMethod(BackwardEulerMethod& rOther)
      :BaseType(rOther)
      ,mDeltaTime(rOther.mDeltaTime)
    {
    }
 
    BasePointerType Clone() override
    {
      return BasePointerType( new BackwardEulerMethod(*this) );
    }
 
    ~BackwardEulerMethod() override{}
 
 
 
    // set parameters (do not calculate parameters here, only read them)
    void SetParameters(const ProcessInfo& rCurrentProcessInfo) override
    {
     KRATOS_TRY
 
     const 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;
        }
 
     mDeltaTime = delta_time;
 
     KRATOS_CATCH( "" )
    }
 
    // set parameters to process info
    void SetProcessInfoParameters(ProcessInfo& rCurrentProcessInfo) override
    {
     KRATOS_TRY
 
 
     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 << "BackwardEulerMethod";
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "BackwardEulerMethod";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
      rOStream << "BackwardEulerMethod Data";
    }
 
 
 
 
 
  protected:
 
 
 
    double  mDeltaTime;
 
 
    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);
      const TValueType& PreviousVariable         = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& OldPreviousVariable      = rNode.FastGetSolutionStepValue(*this->mpVariable,         2);
 
      TValueType& CurrentFirstDerivative         = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      // backward euler consistent
      CurrentVariable = (1.0/3.0) * ( 4.0 * PreviousVariable - OldPreviousVariable + 2.0 * mDeltaTime * CurrentFirstDerivative);
 
      TValueType& CurrentSecondDerivative        = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentFirstDerivative   = PreviousFirstDerivative;
      CurrentSecondDerivative -= CurrentSecondDerivative;
 
      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);
 
      const TValueType& CurrentSecondDerivative        = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      const TValueType& PreviousFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      // backward euler consistent
      CurrentVariable = PreviousVariable + mDeltaTime * ( PreviousFirstDerivative + 0.5 * mDeltaTime * CurrentSecondDerivative );
 
      KRATOS_CATCH( "" )
    }
 
    void PredictFromVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      this->PredictVariable(rNode);
      this->PredictFirstDerivative(rNode);
      this->PredictSecondDerivative(rNode);
 
      KRATOS_CATCH( "" )
    }
 
    void PredictFromFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      this->PredictVariable(rNode);
      this->PredictFirstDerivative(rNode);
      this->PredictSecondDerivative(rNode);
 
      KRATOS_CATCH( "" )
    }
 
    void PredictVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      TValueType& CurrentVariable               = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable        = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
 
      const TValueType& PreviousFirstDerivative = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      // variable prediction
      CurrentVariable = PreviousVariable + mDeltaTime * PreviousFirstDerivative;
 
 
      KRATOS_CATCH( "" )
    }
 
    void PredictFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      // TValueType& CurrentFirstDerivative           = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      // const TValueType& PreviousFirstDerivative    = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      // CurrentFirstDerivative = PreviousFirstDerivative;
 
      const TValueType& CurrentVariable            = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable           = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& OldPreviousVariable        = rNode.FastGetSolutionStepValue(*this->mpVariable,         2);
 
      TValueType& CurrentFirstDerivative           = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
 
      CurrentFirstDerivative = (1.0/mDeltaTime) * ( 3.0 * 0.5 * CurrentVariable - 2.0 * PreviousVariable + 0.5 * OldPreviousVariable);
 
      KRATOS_CATCH( "" )
    }
 
    void PredictSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      // TValueType& CurrentSecondDerivative          = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
      // const TValueType& PreviousSecondDerivative   = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 1);
 
      // CurrentSecondDerivative = PreviousSecondDerivative;
 
      const TValueType& CurrentVariable            = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable           = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& PreviousFirstDerivative    = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      TValueType& CurrentSecondDerivative          = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentSecondDerivative = (2.0/ (mDeltaTime*mDeltaTime)) * (CurrentVariable - PreviousVariable - mDeltaTime * PreviousFirstDerivative);
 
      KRATOS_CATCH( "" )
    }
 
    void UpdateFromVariable(NodeType& rNode) override
    {
      KRATOS_TRY
 
      const TValueType& CurrentVariable            = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable           = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& OldPreviousVariable        = rNode.FastGetSolutionStepValue(*this->mpVariable,         2);
 
      TValueType& CurrentFirstDerivative           = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      const TValueType& PreviousFirstDerivative    = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      TValueType& CurrentSecondDerivative          = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentSecondDerivative = (2.0/ (mDeltaTime*mDeltaTime)) * (CurrentVariable - PreviousVariable - mDeltaTime * PreviousFirstDerivative);
 
      CurrentFirstDerivative = (1.0/mDeltaTime) * ( (3.0/2.0) * CurrentVariable - 2.0 * PreviousVariable + 0.5 * OldPreviousVariable);
 
      KRATOS_CATCH( "" )
    }
 
 
    void UpdateFromFirstDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
 
      TValueType& CurrentVariable               = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable        = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& OldPreviousVariable     = rNode.FastGetSolutionStepValue(*this->mpVariable,         2);
 
      const TValueType& CurrentFirstDerivative  = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      const TValueType& PreviousFirstDerivative = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      TValueType& CurrentSecondDerivative       = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentVariable = (2.0 / 3.0) * ( mDeltaTime * CurrentFirstDerivative + 2.0 * PreviousVariable - 0.5 * OldPreviousVariable);
      CurrentSecondDerivative = (2.0/ (mDeltaTime*mDeltaTime)) * (CurrentVariable - PreviousVariable - mDeltaTime * PreviousFirstDerivative);
 
 
      KRATOS_CATCH( "" )
    }
 
    void UpdateFromSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      TValueType& CurrentVariable                  = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable           = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& OldPreviousVariable        = rNode.FastGetSolutionStepValue(*this->mpVariable,         2);
 
      TValueType& CurrentFirstDerivative           = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
      const TValueType& PreviousFirstDerivative    = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      const TValueType& CurrentSecondDerivative    = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentVariable = PreviousVariable + mDeltaTime * PreviousFirstDerivative + 0.5 * mDeltaTime * mDeltaTime * CurrentSecondDerivative;
      CurrentFirstDerivative = (1.0/mDeltaTime) * ( 3.0 * 0.5 * CurrentVariable - 2.0 * PreviousVariable + 0.5 * OldPreviousVariable);
 
      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);
      const TValueType& PreviousVariable           = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& OldPreviousVariable        = rNode.FastGetSolutionStepValue(*this->mpVariable,         2);
 
      TValueType& CurrentFirstDerivative           = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  0);
 
      CurrentFirstDerivative = (1.0/mDeltaTime) * ( 3.0 * 0.5 * CurrentVariable - 2.0 * PreviousVariable + 0.5 * OldPreviousVariable);
 
      KRATOS_CATCH( "" )
    }
 
    void UpdateSecondDerivative(NodeType& rNode) override
    {
      KRATOS_TRY
 
      const TValueType& CurrentVariable            = rNode.FastGetSolutionStepValue(*this->mpVariable,         0);
      const TValueType& PreviousVariable           = rNode.FastGetSolutionStepValue(*this->mpVariable,         1);
      const TValueType& PreviousFirstDerivative    = rNode.FastGetSolutionStepValue(*this->mpFirstDerivative,  1);
 
      TValueType& CurrentSecondDerivative          = rNode.FastGetSolutionStepValue(*this->mpSecondDerivative, 0);
 
      CurrentSecondDerivative = (2.0/ (mDeltaTime*mDeltaTime)) * (CurrentVariable - PreviousVariable - mDeltaTime * PreviousFirstDerivative);
 
      KRATOS_CATCH( "" )
    }
 
 
 
    // get parameters
    double& GetFirstDerivativeInertialParameter(double& rParameter) override
    {
      rParameter = 3.0 / (2.0 * mDeltaTime);
      return rParameter;
    }
 
    double& GetSecondDerivativeInertialParameter(double& rParameter) override
    {
      rParameter = 2.0 / (mDeltaTime*mDeltaTime);
      return rParameter;
    }
 
 
 
 
  private:
 
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
      KRATOS_SERIALIZE_SAVE_BASE_CLASS( rSerializer, BaseType )
      rSerializer.save("DeltaTime", mDeltaTime);
    };
 
    void load(Serializer& rSerializer) override
    {
      KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, BaseType )
      rSerializer.load("DeltaTime", mDeltaTime);
    };
 
 
 
 
  }; // Class BackwardEulerMethod
 
 
 
 
 
  template<class TVariableType, class TValueType>
  inline std::istream & operator >> (std::istream & rIStream, BackwardEulerMethod<TVariableType,TValueType>& rThis)
  {
    return rIStream;
  }
 
  template<class TVariableType, class TValueType>
  inline std::ostream & operator << (std::ostream & rOStream, const BackwardEulerMethod<TVariableType,TValueType>& rThis)
  {
    return rOStream << rThis.Info();
  }
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_BACKWARD_EULER_METHOD_H_INCLUDED defined