dynamic_scheme.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_DYNAMIC_SCHEME_H_INCLUDED)
#define  KRATOS_DYNAMIC_SCHEME_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_solvers/solution_schemes/static_scheme.hpp"
 
namespace Kratos
{
 
  template<class TSparseSpace,  class TDenseSpace >
  class DynamicScheme: public StaticScheme<TSparseSpace,TDenseSpace>
  {
  protected:
 
    struct  GeneralMatrices
    {
      std::vector< Matrix > M;     // First derivative matrix  (usually mass matrix)
      std::vector< Matrix > D;     // Second derivative matrix (usually damping matrix)
    };
 
    struct GeneralVectors
    {
      std::vector< Vector > v;    // Velocity
      std::vector< Vector > a;    // Acceleration
      std::vector< Vector > c;    // Composed Variable
    };
 
 
  public:
 
    KRATOS_CLASS_POINTER_DEFINITION( DynamicScheme );
 
    typedef SolutionScheme<TSparseSpace,TDenseSpace>                             BaseType;
    typedef typename BaseType::SolutionSchemePointerType                  BasePointerType;
    typedef typename BaseType::LocalFlagType                                LocalFlagType;
 
    typedef StaticScheme<TSparseSpace,TDenseSpace>                           DerivedType;
 
    typedef typename BaseType::NodeType                                          NodeType;
    typedef typename BaseType::DofsArrayType                                DofsArrayType;
    typedef typename BaseType::SystemMatrixType                          SystemMatrixType;
    typedef typename BaseType::SystemVectorType                          SystemVectorType;
    typedef typename BaseType::LocalSystemVectorType                LocalSystemVectorType;
    typedef typename BaseType::LocalSystemMatrixType                LocalSystemMatrixType;
 
    typedef typename BaseType::NodesContainerType                      NodesContainerType;
    typedef typename BaseType::ElementsContainerType                ElementsContainerType;
    typedef typename BaseType::ConditionsContainerType            ConditionsContainerType;
 
    typedef typename BaseType::IntegrationMethodsVectorType  IntegrationMethodsVectorType;
    typedef typename BaseType::IntegrationMethodsScalarType  IntegrationMethodsScalarType;
 
 
    DynamicScheme(IntegrationMethodsVectorType& rTimeVectorIntegrationMethods, Flags& rOptions)
        :DerivedType(rTimeVectorIntegrationMethods, rOptions)
    {
    }
 
    DynamicScheme(IntegrationMethodsVectorType& rTimeVectorIntegrationMethods)
        :DerivedType(rTimeVectorIntegrationMethods)
    {
    }
 
    DynamicScheme(IntegrationMethodsScalarType& rTimeScalarIntegrationMethods, Flags& rOptions)
        :DerivedType(rTimeScalarIntegrationMethods, rOptions)
    {
    }
 
    DynamicScheme(IntegrationMethodsScalarType& rTimeScalarIntegrationMethods)
        :DerivedType(rTimeScalarIntegrationMethods)
    {
    }
 
 
    DynamicScheme(IntegrationMethodsVectorType& rTimeVectorIntegrationMethods,
                  IntegrationMethodsScalarType& rTimeScalarIntegrationMethods,
                  Flags& rOptions)
        :DerivedType(rTimeVectorIntegrationMethods, rTimeScalarIntegrationMethods, rOptions)
    {
    }
 
    DynamicScheme(IntegrationMethodsVectorType& rTimeVectorIntegrationMethods,
                  IntegrationMethodsScalarType& rTimeScalarIntegrationMethods)
        :DerivedType(rTimeVectorIntegrationMethods, rTimeScalarIntegrationMethods)
    {
    }
 
    DynamicScheme(Flags& rOptions)
      :DerivedType(rOptions)
    {
    }
 
    DynamicScheme(DynamicScheme& rOther)
      :DerivedType(rOther)
      ,mMatrix(rOther.mMatrix)
      ,mVector(rOther.mVector)
    {
    }
 
    BasePointerType Clone() override
    {
      return BasePointerType( new DynamicScheme(*this) );
    }
 
    ~DynamicScheme() override {}
 
 
 
    void Initialize(ModelPart& rModelPart) override
    {
        KRATOS_TRY
 
    DerivedType::Initialize(rModelPart);
 
    // Allocate auxiliary memory
    const unsigned int NumThreads = ParallelUtilities::GetNumThreads();
 
    mMatrix.M.resize(NumThreads);
    mMatrix.D.resize(NumThreads);
 
    mVector.v.resize(NumThreads);
    mVector.a.resize(NumThreads);
 
        double parameter = 0.0;
        if(this->mTimeVectorIntegrationMethods.size() != 0)
          parameter = this->mTimeVectorIntegrationMethods.front()->GetSecondDerivativeKineticFactor(parameter);
        else if(this->mTimeScalarIntegrationMethods.size() != 0)
          parameter = this->mTimeScalarIntegrationMethods.front()->GetSecondDerivativeKineticFactor(parameter);
 
        if( parameter != 0 )
          mVector.c.resize(NumThreads);
 
    KRATOS_CATCH("")
    }
 
 
    void CalculateSystemContributions(Element::Pointer pCurrentElement,
                      LocalSystemMatrixType& rLHS_Contribution,
                      LocalSystemVectorType& rRHS_Contribution,
                      Element::EquationIdVectorType& EquationId,
                      ProcessInfo& rCurrentProcessInfo) override
    {
      KRATOS_TRY;
 
      int thread = OpenMPUtils::ThisThread();
 
      (pCurrentElement) -> CalculateLocalSystem(rLHS_Contribution,rRHS_Contribution, rCurrentProcessInfo);
 
      (pCurrentElement) -> EquationIdVector(EquationId, rCurrentProcessInfo);
 
      if ( rCurrentProcessInfo[COMPUTE_DYNAMIC_TANGENT] == true ){
 
        (pCurrentElement) -> CalculateSecondDerivativesContributions(this->mMatrix.M[thread],this->mVector.a[thread],rCurrentProcessInfo);
 
        (pCurrentElement) -> CalculateFirstDerivativesContributions(this->mMatrix.D[thread],this->mVector.v[thread],rCurrentProcessInfo);
 
        AddDynamicTangentsToLHS(rLHS_Contribution,this->mMatrix.D[thread],this->mMatrix.M[thread],rCurrentProcessInfo);
 
        AddDynamicForcesToRHS(rRHS_Contribution,this->mVector.v[thread],this->mVector.a[thread],rCurrentProcessInfo);
 
      }
      else{
 
        (pCurrentElement) -> CalculateMassMatrix(mMatrix.M[thread], rCurrentProcessInfo);
 
        (pCurrentElement) -> CalculateDampingMatrix(mMatrix.D[thread], rCurrentProcessInfo);
 
        AddDynamicsToLHS (rLHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], rCurrentProcessInfo);
 
        AddDynamicsToRHS (pCurrentElement, rRHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], rCurrentProcessInfo);
 
      }
 
      KRATOS_CATCH("")
    }
 
    void Calculate_RHS_Contribution(Element::Pointer pCurrentElement,
                    LocalSystemVectorType& rRHS_Contribution,
                    Element::EquationIdVectorType& EquationId,
                    ProcessInfo& rCurrentProcessInfo) override
    {
 
      KRATOS_TRY;
 
      int thread = OpenMPUtils::ThisThread();
 
      // Basic operations for the element considered
      (pCurrentElement) -> CalculateRightHandSide(rRHS_Contribution, rCurrentProcessInfo);
 
      (pCurrentElement) -> EquationIdVector(EquationId, rCurrentProcessInfo);
 
      if ( rCurrentProcessInfo[COMPUTE_DYNAMIC_TANGENT] == true ){
 
        (pCurrentElement) -> CalculateSecondDerivativesRHS(this->mVector.a[thread],rCurrentProcessInfo);
 
        (pCurrentElement) -> CalculateFirstDerivativesRHS(this->mVector.v[thread],rCurrentProcessInfo);
 
        AddDynamicForcesToRHS(rRHS_Contribution,this->mVector.v[thread],this->mVector.a[thread],rCurrentProcessInfo);
 
      }
      else{
 
        (pCurrentElement) -> CalculateMassMatrix(mMatrix.M[thread], rCurrentProcessInfo);
 
        (pCurrentElement) -> CalculateDampingMatrix(mMatrix.D[thread], rCurrentProcessInfo);
 
        AddDynamicsToRHS (pCurrentElement, rRHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], rCurrentProcessInfo);
      }
 
      KRATOS_CATCH("")
    }
 
    void Condition_CalculateSystemContributions(Condition::Pointer pCurrentCondition,
                        LocalSystemMatrixType& rLHS_Contribution,
                        LocalSystemVectorType& rRHS_Contribution,
                        Element::EquationIdVectorType& EquationId,
                        ProcessInfo& rCurrentProcessInfo) override
    {
      KRATOS_TRY;
 
      int thread = OpenMPUtils::ThisThread();
 
      // Basic operations for the condition considered
      (pCurrentCondition) -> CalculateLocalSystem(rLHS_Contribution,rRHS_Contribution, rCurrentProcessInfo);
 
      (pCurrentCondition) -> EquationIdVector(EquationId, rCurrentProcessInfo);
 
      if ( rCurrentProcessInfo[COMPUTE_DYNAMIC_TANGENT] == true ){
 
        (pCurrentCondition) -> CalculateSecondDerivativesContributions(this->mMatrix.M[thread],this->mVector.a[thread],rCurrentProcessInfo);
 
        (pCurrentCondition) -> CalculateFirstDerivativesContributions(this->mMatrix.D[thread],this->mVector.v[thread],rCurrentProcessInfo);
 
        AddDynamicTangentsToLHS(rLHS_Contribution,this->mMatrix.D[thread],this->mMatrix.M[thread],rCurrentProcessInfo);
 
        AddDynamicForcesToRHS(rRHS_Contribution,this->mVector.v[thread],this->mVector.a[thread],rCurrentProcessInfo);
 
      }
      else{
 
        (pCurrentCondition) -> CalculateMassMatrix(mMatrix.M[thread], rCurrentProcessInfo);
 
        (pCurrentCondition) -> CalculateDampingMatrix(mMatrix.D[thread], rCurrentProcessInfo);
 
        AddDynamicsToLHS  (rLHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], rCurrentProcessInfo);
 
        AddDynamicsToRHS  (pCurrentCondition, rRHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], rCurrentProcessInfo);
      }
 
      KRATOS_CATCH("")
    }
 
    void Condition_Calculate_RHS_Contribution(Condition::Pointer pCurrentCondition,
                          LocalSystemVectorType& rRHS_Contribution,
                          Element::EquationIdVectorType& EquationId,
                          ProcessInfo& rCurrentProcessInfo) override
    {
      KRATOS_TRY;
 
      int thread = OpenMPUtils::ThisThread();
 
      // Basic operations for the condition considered
      (pCurrentCondition) -> CalculateRightHandSide(rRHS_Contribution, rCurrentProcessInfo);
 
      (pCurrentCondition) -> EquationIdVector(EquationId, rCurrentProcessInfo);
 
      if ( rCurrentProcessInfo[COMPUTE_DYNAMIC_TANGENT] == true ){
 
        (pCurrentCondition) -> CalculateSecondDerivativesRHS(this->mVector.a[thread],rCurrentProcessInfo);
 
        (pCurrentCondition) -> CalculateFirstDerivativesRHS(this->mVector.v[thread],rCurrentProcessInfo);
 
        AddDynamicForcesToRHS(rRHS_Contribution,this->mVector.v[thread],this->mVector.a[thread],rCurrentProcessInfo);
 
      }
      else{
 
        (pCurrentCondition) -> CalculateMassMatrix(mMatrix.M[thread], rCurrentProcessInfo);
 
        (pCurrentCondition) -> CalculateDampingMatrix(mMatrix.D[thread], rCurrentProcessInfo);
 
        AddDynamicsToRHS  (pCurrentCondition, rRHS_Contribution, mMatrix.D[thread], mMatrix.M[thread], rCurrentProcessInfo);
      }
 
      KRATOS_CATCH("")
    }
 
 
    int Check(ModelPart& rModelPart) override
    {
      KRATOS_TRY;
 
      // Perform base base checks
      int ErrorCode = 0;
      ErrorCode  = DerivedType::Check(rModelPart);
 
      if( !rModelPart.GetProcessInfo().Has(COMPUTE_DYNAMIC_TANGENT) )
        KRATOS_ERROR << "COMPUTE_DYNAMIC_TANGENT must be set to use a Dynamic Scheme" << std::endl;
 
 
      return ErrorCode;
 
      KRATOS_CATCH("")
    }
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "DynamicScheme";
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "DynamicScheme";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
      rOStream << "DynamicScheme Data";
    }
 
 
 
  protected:
 
 
 
    GeneralMatrices     mMatrix;
 
    GeneralVectors      mVector;
 
 
 
 
    virtual void AddDynamicsToLHS(LocalSystemMatrixType& rLHS_Contribution,
                  LocalSystemMatrixType& rD,
                  LocalSystemMatrixType& rM,
                  ProcessInfo& rCurrentProcessInfo)
    {
 
      double parameter = 0;
      // Adding mass contribution to the dynamic stiffness
      if (rM.size1() != 0) // if M matrix declared
        {
          if(this->mTimeVectorIntegrationMethods.size() != 0)
            parameter = this->mTimeVectorIntegrationMethods.front()->GetSecondDerivativeInertialFactor(parameter);
          else if(this->mTimeScalarIntegrationMethods.size() != 0)
            parameter = this->mTimeScalarIntegrationMethods.front()->GetSecondDerivativeInertialFactor(parameter);
 
      noalias(rLHS_Contribution) += rM * parameter;
        }
 
      // Adding  damping contribution
      if (rD.size1() != 0) // if D matrix declared
        {
          parameter = 0;
          if(this->mTimeVectorIntegrationMethods.size() != 0)
            parameter = this->mTimeVectorIntegrationMethods.front()->GetFirstDerivativeInertialFactor(parameter);
          else if(this->mTimeScalarIntegrationMethods.size() != 0)
            parameter = this->mTimeScalarIntegrationMethods.front()->GetFirstDerivativeInertialFactor(parameter);
 
      noalias(rLHS_Contribution) += rD * parameter;
        }
    }
 
 
    virtual void AddDynamicTangentsToLHS(LocalSystemMatrixType& rLHS_Contribution,
                                         LocalSystemMatrixType& rD,
                                         LocalSystemMatrixType& rM,
                                         ProcessInfo& rCurrentProcessInfo)
    {
 
      // Adding mass contribution to the dynamic stiffness
      if (rM.size1() != 0) // if M matrix declared
        {
      noalias(rLHS_Contribution) += rM ;
        }
 
      // Adding  damping contribution
      if (rD.size1() != 0) // if D matrix declared
        {
      noalias(rLHS_Contribution) += rD;
        }
    }
 
    virtual void AddDynamicsToRHS(Element::Pointer pCurrentElement,
                  LocalSystemVectorType& rRHS_Contribution,
                  LocalSystemMatrixType& rD,
                  LocalSystemMatrixType& rM,
                  ProcessInfo& rCurrentProcessInfo)
    {
      int thread = OpenMPUtils::ThisThread();
 
      // Adding inertia contribution
      if (rM.size1() != 0)
        {
      pCurrentElement->GetSecondDerivativesVector(mVector.a[thread], 0);
 
          double parameter = 0.0;
          if(this->mTimeVectorIntegrationMethods.size() != 0)
            parameter = this->mTimeVectorIntegrationMethods.front()->GetSecondDerivativeKineticFactor(parameter);
          else if(this->mTimeScalarIntegrationMethods.size() != 0)
            parameter = this->mTimeScalarIntegrationMethods.front()->GetSecondDerivativeKineticFactor(parameter);
 
          if( parameter != 0 ){
 
            (mVector.a[thread]) *= (1.00 - parameter);
 
            pCurrentElement->GetSecondDerivativesVector(mVector.c[thread], 1);
 
            noalias(mVector.a[thread]) += parameter * mVector.c[thread];
 
          }
 
      noalias(rRHS_Contribution) -= prod(rM, mVector.a[thread]);
        }
 
      // Adding damping contribution
      if (rD.size1() != 0)
        {
      pCurrentElement->GetFirstDerivativesVector(mVector.v[thread], 0);
 
      noalias(rRHS_Contribution) -= prod(rD, mVector.v[thread]);
        }
    }
 
    virtual void AddDynamicsToRHS(Condition::Pointer pCurrentCondition,
                  LocalSystemVectorType& rRHS_Contribution,
                  LocalSystemMatrixType& rD,
                  LocalSystemMatrixType& rM,
                  ProcessInfo& rCurrentProcessInfo)
    {
      int thread = OpenMPUtils::ThisThread();
 
      // Adding inertia contribution
      if (rM.size1() != 0)
        {
      pCurrentCondition->GetSecondDerivativesVector(mVector.a[thread], 0);
 
          double parameter = 0.0;
          if(this->mTimeVectorIntegrationMethods.size() != 0)
            parameter = this->mTimeVectorIntegrationMethods.front()->GetSecondDerivativeKineticFactor(parameter);
          else if(this->mTimeScalarIntegrationMethods.size() != 0)
            parameter = this->mTimeScalarIntegrationMethods.front()->GetSecondDerivativeKineticFactor(parameter);
 
          if( parameter != 0 ){
 
            (mVector.a[thread]) *= (1.00 - parameter);
 
            pCurrentCondition->GetSecondDerivativesVector(mVector.c[thread], 1);
 
            noalias(mVector.a[thread]) += parameter * mVector.c[thread];
 
          }
 
      noalias(rRHS_Contribution) -= prod(rM, mVector.a[thread]);
        }
 
      // Adding damping contribution
      if (rD.size1() != 0)
        {
      pCurrentCondition->GetFirstDerivativesVector(mVector.v[thread], 0);
 
      noalias(rRHS_Contribution) -= prod(rD, mVector.v[thread]);
        }
 
    }
 
 
    void AddDynamicForcesToRHS(LocalSystemVectorType& rRHS_Contribution,
                   LocalSystemVectorType& rfv ,
                   LocalSystemVectorType& rfa,
                   ProcessInfo& rCurrentProcessInfo)
    {
 
      // Adding inertia contribution
      if (rfa.size() != 0)
        {
      noalias(rRHS_Contribution) -=  rfa;
        }
 
      // Adding damping contribution
      if (rfv.size() != 0)
        {
      noalias(rRHS_Contribution) -=  rfv;
        }
    }
 
 
 
 
 
  private:
 
 
 
 
 
 
 
 
 
  }; // Class DynamicScheme
 
 
 
 
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_DYNAMIC_SCHEME_H_INCLUDED defined