beam_element.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:            JMCarbonell $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:              August 2017 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(KRATOS_BEAM_ELEMENT_H_INCLUDED)
#define  KRATOS_BEAM_ELEMENT_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "includes/checks.h"
#include "includes/element.h"
#include "utilities/beam_math_utilities.hpp"
#include "custom_utilities/element_utilities.hpp"
 
namespace Kratos
{
 
 
class KRATOS_API(SOLID_MECHANICS_APPLICATION) BeamElement
    :public Element
{
public:
 
    typedef ConstitutiveLaw                         ConstitutiveLawType;
    typedef ConstitutiveLawType::Pointer     ConstitutiveLawPointerType;
    typedef ConstitutiveLawType::StressMeasure        StressMeasureType;
    typedef GeometryData::IntegrationMethod           IntegrationMethod;
    typedef BeamMathUtils<double>                     BeamMathUtilsType;
    typedef Quaternion<double>                           QuaternionType;
    typedef GeometryData::SizeType                             SizeType;
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( BeamElement );
 
 
protected:
 
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_RHS_VECTOR );
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_LHS_MATRIX );
    KRATOS_DEFINE_LOCAL_FLAG( FINALIZED_STEP );
    struct SectionProperties
    {
      double Area;                            // Area or the beam section
      double Inertia_z;                       // Moment of Inertia about the local z axis, Iz local
      double Inertia_y;                       // Moment of Inertia about the local y axis, Iy local
      double Polar_Inertia;                   // Polar Moment of Inertia, measures an object's ability to resist twisting, when acted upon by differences of torque along its length.
      double Rotational_Inertia;              // Moment of Inertia about the local x axis, measures an object's resistance to changes in its rotational velocity when acted by a net resultant torque
    };
 
    struct DirectorsVariables
    {
      //Current
      std::vector<Vector> Current;
      std::vector<Vector> CurrentDerivatives;
 
      //Initial
      std::vector<Vector> Initial;
      std::vector<Vector> InitialDerivatives;
 
      //Previous
      std::vector<Vector> Previous;
      std::vector<Vector> PreviousDerivatives;
 
 
      //NODAL VARIABLES
      std::vector<Matrix> InitialNode;
      std::vector<Matrix> CurrentNode;
      std::vector<Matrix> PreviousNode;
 
      std::vector<Matrix> CurrentNodeVelocities;
      std::vector<Matrix> PreviousNodeVelocities;
    };
 
 
    struct ElementData
    {
      private:
 
        //variables including all integration points
        const GeometryType::ShapeFunctionsGradientsType* pDN_De;
        const Matrix* pNcontainer;
 
        //directors variables
        DirectorsVariables* Directors;
 
      public:
 
        //integration point
        unsigned int PointNumber;
 
        //delta time
        double DeltaTime;
 
        //element length
        double  Length;
        double  detJ;
 
        //equilibrium point
        double Alpha;
 
        //general variables
        Vector  StrainVector;
        Vector  StressVector;
        Vector  N;
 
        Matrix  B;
        Matrix  DN_DX;
        Matrix  ConstitutiveMatrix;
        Matrix  DeltaPosition;
 
        //large displacement
        Vector  CurrentCurvatureVector;
        Vector  PreviousCurvatureVector;
 
        Vector  CurrentStepRotationVector;
 
        Vector  CurrentStrainResultantsVector;
        Vector  PreviousStrainResultantsVector;
 
        Vector  InitialAxisPositionDerivatives;
        Vector  CurrentAxisPositionDerivatives;
        Vector  PreviousAxisPositionDerivatives;
 
 
        Matrix  AlphaRotationMatrix;
        Matrix  AlphaRotationMatrixAsterisk;
 
        Matrix  CurrentRotationMatrix;
        Matrix  PreviousRotationMatrix;
 
 
        //variables including all integration points
        GeometryType::JacobiansType J;
        GeometryType::JacobiansType j;
 
        //section properties
        SectionProperties Section;
 
        void SetShapeFunctionsGradients(const GeometryType::ShapeFunctionsGradientsType &rDN_De)
        {
            pDN_De=&rDN_De;
        };
 
        void SetShapeFunctions(const Matrix& rNcontainer)
        {
            pNcontainer=&rNcontainer;
        };
 
 
        void SetDirectors(DirectorsVariables& rDirectors)
        {
            Directors=&rDirectors;
        };
 
        const GeometryType::ShapeFunctionsGradientsType& GetShapeFunctionsGradients()
        {
            return *pDN_De;
        };
 
        const Matrix& GetShapeFunctions()
        {
            return *pNcontainer;
        };
 
        DirectorsVariables& GetDirectors()
        {
            return *Directors;
        };
 
        void Initialize( const unsigned int& voigt_size,
             const unsigned int& dimension,
             const unsigned int& number_of_nodes )
        {
      //scalars
      PointNumber = 0;
      Length = 0;
      detJ  = 1;
      Alpha = 1;
          DeltaTime = 0;
 
      //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-1) * 3 * number_of_nodes,false);
      DN_DX.resize(number_of_nodes, 1,false); //Local 1D
      ConstitutiveMatrix.resize(voigt_size, voigt_size,false);
      DeltaPosition.resize(number_of_nodes, dimension,false);
 
      noalias(DN_DX) = ZeroMatrix(number_of_nodes, 1);
      noalias(ConstitutiveMatrix) = ZeroMatrix(voigt_size, voigt_size);
      noalias(DeltaPosition) = ZeroMatrix(number_of_nodes, dimension);
 
      //large displacement
      CurrentRotationMatrix.resize(dimension, dimension,false);
      InitialAxisPositionDerivatives.resize(dimension,false);
      CurrentAxisPositionDerivatives.resize(dimension,false);
      PreviousAxisPositionDerivatives.resize(dimension,false);
 
      noalias(CurrentRotationMatrix) = ZeroMatrix(dimension, dimension);
      noalias(InitialAxisPositionDerivatives) = ZeroVector(dimension);
      noalias(CurrentAxisPositionDerivatives) = ZeroVector(dimension);
      noalias(PreviousAxisPositionDerivatives) = ZeroVector(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;
      Directors = 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;
 
 
    BeamElement(IndexType NewId, GeometryType::Pointer pGeometry);
 
    BeamElement(IndexType NewId, GeometryType::Pointer pGeometry,  PropertiesType::Pointer pProperties);
 
    BeamElement(BeamElement const& rOther);
 
    ~BeamElement() override;
 
 
 
 
    Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes,  PropertiesType::Pointer pProperties) 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;
 
 
    //************* 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 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 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< array_1d<double, 3 > >& rVariable,
                                       std::vector< array_1d<double, 3 > >& Output,
                                       const ProcessInfo& rCurrentProcessInfo) override;
 
 
    //************************************************************************************
    //************************************************************************************
    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Large Displacement Beam Element #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "Large Displacement Beam Element #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
      GetGeometry().PrintData(rOStream);
    }
 
protected:
 
 
    IntegrationMethod mThisIntegrationMethod;
 
    QuaternionType  mInitialLocalQuaternion;
 
    BeamElement() {};
 
    //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 );
 
    virtual unsigned int GetDofsSize() const;
 
    void InitializeSystemMatrices(MatrixType& rLeftHandSideMatrix,
                  VectorType& rRightHandSideVector,
                  Flags& rCalculationFlags);
 
 
    Vector& MapToInitialLocalFrame(Vector& rVariable, unsigned int PointNumber = 0);
 
    virtual void MapLocalToGlobal(ElementDataType& rVariables, Matrix& rVariable);
 
    virtual void MapLocalToGlobal(ElementDataType& rVariables, VectorType& rVector);
 
 
    Vector& GetLocalCurrentValue(const Variable<array_1d<double,3> >&rVariable,
                 Vector& rValue,
                 const Vector& rN);
 
 
    Vector& GetLocalPreviousValue(const Variable<array_1d<double,3> >&rVariable,
                  Vector& rValue,
                  const Vector& rN);
 
 
    Vector& GetNodalCurrentValue(const Variable<array_1d<double,3> >&rVariable,
                 Vector& rValue,
                 const unsigned int& rNode);
 
    Vector& GetNodalPreviousValue(const Variable<array_1d<double,3> >&rVariable,
                  Vector& rValue,
                  const unsigned int& rNode);
 
    virtual void InitializeElementData(ElementDataType & rVariables, const ProcessInfo& rCurrentProcessInfo);
 
 
    void CalculateSectionProperties(SectionProperties & rSection);
 
 
    virtual void CalculateKinematics(ElementDataType& rVariables,
                                     const unsigned int& rPointNumber);
 
    virtual void CalculateConstitutiveMatrix(ElementDataType& rVariables);
 
 
    void CalculateMaterialConstitutiveMatrix(Matrix& rConstitutiveMatrix, ElementDataType& rVariables);
 
 
    virtual void CalculateStressResultants(ElementDataType& rVariables, const unsigned int& rPointNumber);
 
    virtual void CalculateAndAddLHS(LocalSystemComponents& rLocalSystem, ElementDataType& rVariables, double& rIntegrationWeight);
 
    virtual void CalculateAndAddRHS(LocalSystemComponents& rLocalSystem, ElementDataType& rVariables, Vector& VolumeForce, double& rIntegrationWeight);
 
    virtual double& CalculateIntegrationWeight(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 CalculateAndAddInertiaLHS(MatrixType& rLeftHandSideMatrix,
                       ElementDataType& rVariables,
                       const ProcessInfo& rCurrentProcessInfo,
                       double& rIntegrationWeight);
 
    virtual void CalculateAndAddInertiaRHS(VectorType& rRightHandSideVector,
                       ElementDataType& rVariables,
                       const ProcessInfo& rCurrentProcessInfo,
                       double& rIntegrationWeight);
 
    virtual void CalculateAndAddInternalForces(VectorType& rRightHandSideVector,
                           ElementDataType & rVariables,
                           double& rIntegrationWeight);
 
 
    void CalculateLocalAxesMatrix(Matrix& rRotationMatrix);
 
 
    virtual Vector& CalculateVolumeForce(Vector& rVolumeForce, const Vector& rN);
 
    double& CalculateTotalMass( SectionProperties& Section, double& rTotalMass );
 
    void CalculateInertiaDyadic(SectionProperties& rSection, Matrix& rInertiaDyadic);
 
 
 
private:
 
    friend class Serializer;
 
 
    // A private default constructor necessary for serialization
 
 
    void save(Serializer& rSerializer) const override;
 
    void load(Serializer& rSerializer) override;
 
 
 
}; // Class BeamElement
 
} // namespace Kratos.
#endif //  KRATOS_BEAM_ELEMENT_H_INCLUDED defined