solid_element.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:            JMCarbonell $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:                July 2013 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(KRATOS_SOLID_ELEMENT_H_INCLUDED)
#define  KRATOS_SOLID_ELEMENT_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "includes/kratos_flags.h"
#include "includes/checks.h"
#include "includes/element.h"
#include "custom_utilities/element_utilities.hpp"
 
 
namespace Kratos
{
 
 
class KRATOS_API(SOLID_MECHANICS_APPLICATION) SolidElement
    : public Element
{
public:
 
    typedef ConstitutiveLaw ConstitutiveLawType;
    typedef ConstitutiveLawType::Pointer ConstitutiveLawPointerType;
    typedef ConstitutiveLawType::StressMeasure StressMeasureType;
    typedef GeometryData::IntegrationMethod IntegrationMethod;
    typedef GeometryData::SizeType SizeType;
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( SolidElement );
 
protected:
 
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_RHS_VECTOR );
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_LHS_MATRIX );
    KRATOS_DEFINE_LOCAL_FLAG( FINALIZED_STEP );
 
    struct ElementData
    {
      private:
 
        //variables including all integration points
        const GeometryType::ShapeFunctionsGradientsType* pDN_De;
        const Matrix* pNcontainer;
        const ProcessInfo* pProcessInfo;
 
      public:
 
        StressMeasureType StressMeasure;
 
        double  Tau;
        double  IntegrationWeight;
 
        //for axisymmetric use only
        double  CurrentRadius;
        double  ReferenceRadius;
 
        //general variables for large displacement use
        double  detF;
        double  detF0;
        double  detH; //Wildcard ( detF(0 to n+1) )
        double  detJ;
        Vector  StrainVector;
        Vector  StressVector;
        Vector  N;
        Matrix  B;
        Matrix  H;    //Wildcard ( Displacement Gradient, F(0 to n+1), B-bar, Velocity Gradient...)
        Matrix  F;    //Incremental Deformation Gradient (n to n+1)
        Matrix  F0;   //Historical Deformation Gradient  (0 to n)
        Matrix  DN_DX;
        Matrix  ConstitutiveMatrix;
 
        //variables including all integration points
        GeometryType::JacobiansType J;
        GeometryType::JacobiansType j;
        Matrix  DeltaPosition;
 
 
        void SetShapeFunctionsGradients(const GeometryType::ShapeFunctionsGradientsType &rDN_De)
        {
            pDN_De=&rDN_De;
        };
 
        void SetShapeFunctions(const Matrix& rNcontainer)
        {
            pNcontainer=&rNcontainer;
        };
 
        void SetProcessInfo(const ProcessInfo& rProcessInfo)
        {
            pProcessInfo=&rProcessInfo;
        };
 
 
 
        const GeometryType::ShapeFunctionsGradientsType& GetShapeFunctionsGradients()
        {
            return *pDN_De;
        };
 
        const Matrix& GetShapeFunctions()
        {
            return *pNcontainer;
        };
 
        const ProcessInfo& GetProcessInfo()
        {
            return *pProcessInfo;
        };
 
        void Initialize( const unsigned int& voigt_size,
             const unsigned int& dimension,
             const unsigned int& number_of_nodes )
        {
      StressMeasure = ConstitutiveLaw::StressMeasure_PK2;
 
          //stabilization
          Tau = 0;
 
          //time step
          IntegrationWeight = 1;
 
          //radius
      CurrentRadius = 0;
      ReferenceRadius = 0;
 
          //jacobians
      detF  = 1;
      detF0 = 1;
      detH  = 1;
      detJ  = 1;
 
          //vectors
      StrainVector.resize(voigt_size,false);
          StressVector.resize(voigt_size,false);
      N.resize(number_of_nodes,false);
      noalias(StrainVector) = ZeroVector(voigt_size);
      noalias(StressVector) = ZeroVector(voigt_size);
      noalias(N) = ZeroVector(number_of_nodes);
 
          //matrices
      B.resize(voigt_size, dimension*number_of_nodes,false);
      H.resize(dimension,dimension,false);
      F.resize(dimension,dimension,false);
      F0.resize(dimension,dimension,false);
      DN_DX.resize(number_of_nodes, dimension,false);
      ConstitutiveMatrix.resize(voigt_size, voigt_size,false);
      DeltaPosition.resize(number_of_nodes, dimension,false);
 
      noalias(B)  = ZeroMatrix(voigt_size, dimension*number_of_nodes);
      noalias(H)  = IdentityMatrix(dimension);
      noalias(F)  = IdentityMatrix(dimension);
      noalias(F0) = IdentityMatrix(dimension);
      noalias(DN_DX) = ZeroMatrix(number_of_nodes, dimension);
      noalias(ConstitutiveMatrix) = ZeroMatrix(voigt_size, voigt_size);
      noalias(DeltaPosition) = ZeroMatrix(number_of_nodes, dimension);
 
          //others
      J.resize(1,false);
      j.resize(1,false);
      J[0].resize(dimension,dimension,false);
      j[0].resize(dimension,dimension,false);
      noalias(J[0]) = ZeroMatrix(dimension,dimension);
      noalias(j[0]) = ZeroMatrix(dimension,dimension);
 
          //pointers
      pDN_De = NULL;
      pNcontainer = NULL;
          pProcessInfo = NULL;
    }
 
    };
 
 
    struct LocalSystemComponents
    {
    private:
 
      //for calculation local system with compacted LHS and RHS
      MatrixType *mpLeftHandSideMatrix;
      VectorType *mpRightHandSideVector;
 
    public:
 
      //calculation flags
      Flags  CalculationFlags;
 
      void SetLeftHandSideMatrix( MatrixType& rLeftHandSideMatrix ) { mpLeftHandSideMatrix = &rLeftHandSideMatrix; };
 
      void SetRightHandSideVector( VectorType& rRightHandSideVector ) { mpRightHandSideVector = &rRightHandSideVector; };
 
 
      MatrixType& GetLeftHandSideMatrix() { return *mpLeftHandSideMatrix; };
 
      VectorType& GetRightHandSideVector() { return *mpRightHandSideVector; };
 
    };
 
 
public:
 
 
    typedef ElementData ElementDataType;
 
 
    SolidElement();
 
    SolidElement(IndexType NewId, GeometryType::Pointer pGeometry);
 
    SolidElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties);
 
    SolidElement(SolidElement const& rOther);
 
    ~SolidElement() override;
 
 
    SolidElement& operator=(SolidElement const& rOther);
 
 
    Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const override;
 
    Element::Pointer Clone(IndexType NewId, NodesArrayType const& ThisNodes) const override;
 
 
    //************* GETTING METHODS
 
    IntegrationMethod GetIntegrationMethod() const override;
 
    void GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo& rCurrentProcessInfo) const override;
 
    void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo) const override;
 
    void GetValuesVector(Vector& rValues, int Step = 0) const override;
 
    void GetFirstDerivativesVector(Vector& rValues, int Step = 0) const override;
 
    void GetSecondDerivativesVector(Vector& rValues, int Step = 0) const override;
 
 
 
    //on integration points:
    //SET
    void SetValuesOnIntegrationPoints(const Variable<double>& rVariable, const std::vector<double>& rValues, const ProcessInfo& rCurrentProcessInfo) override;
 
    void SetValuesOnIntegrationPoints(const Variable<Vector>& rVariable, const std::vector<Vector>& rValues, const ProcessInfo& rCurrentProcessInfo) override;
 
    void SetValuesOnIntegrationPoints(const Variable<Matrix>& rVariable, const std::vector<Matrix>& rValues, const ProcessInfo& rCurrentProcessInfo) override;
 
    void SetValuesOnIntegrationPoints( const Variable<ConstitutiveLaw::Pointer>& rVariable,
                                      const std::vector<ConstitutiveLaw::Pointer>& rValues,
                                      const ProcessInfo& rCurrentProcessInfo ) override;
 
 
    //GET:
 
    void CalculateOnIntegrationPoints( const Variable<ConstitutiveLaw::Pointer>& rVariable,
                                      std::vector<ConstitutiveLaw::Pointer>& rValues,
                                      const ProcessInfo& rCurrentProcessInfo ) override;
 
 
 
    //************* STARTING - ENDING  METHODS
 
    void Initialize(const ProcessInfo& rCurrentProcessInfo) override;
 
    void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override;
 
    void FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override;
 
    void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
 
    //************* COMPUTING  METHODS
 
 
    void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
                  VectorType& rRightHandSideVector,
                  const ProcessInfo& rCurrentProcessInfo) override;
 
 
    void CalculateRightHandSide(VectorType& rRightHandSideVector,
                const ProcessInfo& rCurrentProcessInfo) override;
 
 
    void CalculateLeftHandSide (MatrixType& rLeftHandSideMatrix,
                const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateFirstDerivativesContributions(MatrixType& rLeftHandSideMatrix,
                        VectorType& rRightHandSideVector,
                        const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateSecondDerivativesContributions(MatrixType& rLeftHandSideMatrix,
                        VectorType& rRightHandSideVector,
                        const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateSecondDerivativesLHS(MatrixType& rLeftHandSideMatrix,
                       const ProcessInfo& rCurrentProcessInfo) override;
 
 
    void CalculateSecondDerivativesRHS(VectorType& rRightHandSideVector,
                       const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateMassMatrix(MatrixType& rMassMatrix,
                             const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateDampingMatrix(MatrixType& rDampingMatrix,
                                const ProcessInfo& rCurrentProcessInfo) override;
 
 
    void AddExplicitContribution(const VectorType& rRHSVector,
                     const Variable<VectorType>& rRHSVariable,
                     const Variable<array_1d<double,3> >& rDestinationVariable,
                     const ProcessInfo& rCurrentProcessInfo) override;
 
    //on integration points:
    void CalculateOnIntegrationPoints(const Variable<double>& rVariable, std::vector<double>& rOutput, const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateOnIntegrationPoints(const Variable<Vector>& rVariable, std::vector<Vector>& rOutput, const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateOnIntegrationPoints(const Variable<Matrix >& rVariable, std::vector< Matrix >& rOutput, const ProcessInfo& rCurrentProcessInfo) override;
 
 
    //************************************************************************************
    //************************************************************************************
    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Large Displacement Element #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "Large Displacement Element #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
      GetGeometry().PrintData(rOStream);
    }
 
protected:
 
    IntegrationMethod mThisIntegrationMethod;
 
    std::vector<ConstitutiveLaw::Pointer> mConstitutiveLawVector;
 
 
 
    //constexpr const std::size_t& Dimension() const {return GetGeometry().WorkingSpaceDimension();}
 
 
    void IncreaseIntegrationMethod(IntegrationMethod& rThisIntegrationMethod,
                   unsigned int increment) const;
 
    virtual void CalculateElementalSystem(LocalSystemComponents& rLocalSystem,
                                          const ProcessInfo& rCurrentProcessInfo);
 
    virtual void CalculateDynamicSystem(LocalSystemComponents& rLocalSystem,
                    const ProcessInfo& rCurrentProcessInfo);
 
    void PrintElementCalculation(LocalSystemComponents& rLocalSystem, ElementDataType& rVariables);
 
 
    void CalculatePerturbedLeftHandSide (MatrixType& rLeftHandSideMatrix,
                     const ProcessInfo& rCurrentProcessInfo);
 
    virtual void CalculateAndAddLHS(LocalSystemComponents& rLocalSystem,
                                    ElementDataType& rVariables,
                                    double& rIntegrationWeight);
 
    virtual void CalculateAndAddRHS(LocalSystemComponents& rLocalSystem,
                                    ElementDataType& rVariables,
                                    Vector& rVolumeForce,
                                    double& rIntegrationWeight);
 
 
 
    virtual void CalculateAndAddDynamicLHS(MatrixType& rLeftHandSideMatrix,
                       ElementDataType& rVariables,
                       const ProcessInfo& rCurrentProcessInfo,
                       double& rIntegrationWeight);
 
    virtual void CalculateAndAddDynamicRHS(VectorType& rRightHandSideVector,
                       ElementDataType& rVariables,
                       const ProcessInfo& rCurrentProcessInfo,
                       double& rIntegrationWeight);
 
 
 
    virtual void CalculateAndAddKuum(MatrixType& rLeftHandSideMatrix,
                                     ElementDataType& rVariables,
                                     double& rIntegrationWeight);
 
    virtual void CalculateAndAddKuug(MatrixType& rLeftHandSideMatrix,
                                     ElementDataType& rVariables,
                                     double& rIntegrationWeight);
 
 
    virtual void CalculateAndAddExternalForces(VectorType& rRightHandSideVector,
                           ElementDataType& rVariables,
                           Vector& rVolumeForce,
                           double& rIntegrationWeight);
 
 
    virtual void CalculateAndAddInternalForces(VectorType& rRightHandSideVector,
                           ElementDataType & rVariables,
                           double& rIntegrationWeight);
 
 
    virtual void SetElementData(ElementDataType& rVariables,
                                ConstitutiveLaw::Parameters& rValues,
                                const int & rPointNumber);
 
    virtual void CalculateMaterialResponse(ElementDataType& rVariables,
                                           ConstitutiveLaw::Parameters& rValues,
                                           const int & rPointNumber);
    virtual SizeType GetDofsSize() const;
 
 
    bool IsSliver()
    {
      if( this->IsDefined(SELECTED) )
        return this->Is(SELECTED);
      else
        return false;
    }
 
    virtual void InitializeSystemMatrices(MatrixType& rLeftHandSideMatrix,
                                          VectorType& rRightHandSideVector,
                                          Flags& rCalculationFlags);
 
 
 
    void InitializeConstitutiveLaw();
 
 
    void ResetConstitutiveLaw() override;
 
 
    void InitializeExplicitContributions();
 
 
    virtual void CalculateKinematics(ElementDataType& rVariables,
                                     const double& rPointNumber);
 
    virtual void CalculateKinetics(ElementDataType& rVariables,
                   const double& rPointNumber);
 
    virtual void InitializeElementData(ElementDataType & rVariables,
                        const ProcessInfo& rCurrentProcessInfo);
 
    virtual void TransformElementData(ElementDataType & rVariables,
                       const double& rPointNumber);
 
    virtual void FinalizeStepVariables(ElementDataType & rVariables,
                       const double& rPointNumber);
 
 
    virtual double& CalculateIntegrationWeight(double& rIntegrationWeight);
 
 
    virtual double& CalculateTotalMass(double& rTotalMass, const ProcessInfo& rCurrentProcessInfo);
 
    virtual Matrix& CalculateDeltaPosition(Matrix & rDeltaPosition);
 
    virtual Matrix& CalculateTotalDeltaPosition(Matrix & rDeltaPosition);
 
 
    virtual double& CalculateVolumeChange(double& rVolumeChange, ElementDataType& rVariables);
 
    virtual Vector& CalculateVolumeForce(Vector& rVolumeForce, ElementDataType& rVariables);
 
 
 
private:
 
 
 
 
 
 
 
 
    friend class Serializer;
 
    // A private default constructor necessary for serialization
 
    void save(Serializer& rSerializer) const override;
 
    void load(Serializer& rSerializer) override;
 
 
 
}; // Class SolidElement
 
 
} // namespace Kratos.
#endif // KRATOS_SOLID_ELEMENT_H_INCLUDED  defined