fluid_element.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosPfemApplication    $
//   Created by:          $Author:     JMCarbonell $
//   Last modified by:    $Co-Author:              $
//   Date:                $Date:          May 2018 $
//   Revision:            $Revision:           0.0 $
//
//
 
#if !defined(KRATOS_FLUID_ELEMENT_H_INCLUDED)
#define  KRATOS_FLUID_ELEMENT_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "includes/checks.h"
#include "includes/element.h"
 
namespace Kratos
{
 
 
class KRATOS_API(PFEM_APPLICATION) FluidElement
    : 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( FluidElement );
 
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  Alpha;
        double  Tau;
        double  IntegrationWeight;
 
        //general variables for large displacement use
        double  detF;
        double  detJ;
 
        Vector  StrainVector;
        Vector  StressVector;
        Vector  N;
 
        Matrix  B;
        Matrix  L;    //Velocity Gradient
        Matrix  F;    //Incremental Deformation Gradient (n to n+1)
 
        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;
 
          //integration
          Alpha = 1.0;
 
          //stabilization
          Tau = 0;
 
          //time step
          IntegrationWeight = 1;
 
      //jacobians
      detF  = 1;
      detJ  = 1;
 
      //vectors
      StrainVector.resize(voigt_size,false);
          StressVector.resize(voigt_size,false);
      N.resize(number_of_nodes,false);
 
      //matrices
      B.resize(voigt_size,dimension*number_of_nodes,false);
      L.resize(dimension,dimension,false);
      F.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);
 
          //others
      J.resize(1,false);
      j.resize(1,false);
      J[0].resize(dimension,dimension,false);
      j[0].resize(dimension,dimension,false);
 
          //pointers
      pDN_De = NULL;
      pNcontainer = NULL;
          pProcessInfo = NULL;
 
          //fill containters
          // noalias(StrainVector) = ZeroVector(voigt_size);
      // noalias(StressVector) = ZeroVector(voigt_size);
      // noalias(N) = ZeroVector(number_of_nodes);
 
      // noalias(B)  = ZeroMatrix(voigt_size, dimension*number_of_nodes);
      // noalias(L)  = ZeroMatrix(dimension,dimension);
      // noalias(F)  = IdentityMatrix(dimension,dimension);
 
      // noalias(DN_DX) = ZeroMatrix(number_of_nodes, dimension);
      // noalias(ConstitutiveMatrix) = ZeroMatrix(voigt_size, voigt_size);
      // noalias(DeltaPosition) = ZeroMatrix(number_of_nodes, dimension);
 
          // noalias(J[0]) = ZeroMatrix(dimension,dimension);
      // noalias(j[0]) = ZeroMatrix(dimension,dimension);
 
    }
 
    };
 
 
    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;
 
 
    FluidElement();
 
    FluidElement(IndexType NewId, GeometryType::Pointer pGeometry);
 
    FluidElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties);
 
    FluidElement(FluidElement const& rOther);
 
    virtual ~FluidElement();
 
 
    FluidElement& operator=(FluidElement 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
    virtual 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
 
    virtual 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;
 
 
    virtual 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 << "Fluid Element #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "Fluid Element #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
      GetGeometry().PrintData(rOStream);
    }
 
protected:
 
    IntegrationMethod mThisIntegrationMethod;
 
    std::vector<ConstitutiveLaw::Pointer> mConstitutiveLawVector;
 
 
 
    virtual void SetProcessInformation(const ProcessInfo& rCurrentProcessInfo);
 
    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);
 
    virtual void CalculateAndAddRHS(LocalSystemComponents& rLocalSystem,
                                    ElementDataType& rVariables);
 
 
    virtual void CalculateAndAddDynamicLHS(MatrixType& rLeftHandSideMatrix,
                       ElementDataType& rVariables);
 
    virtual void CalculateAndAddDynamicRHS(VectorType& rRightHandSideVector,
                       ElementDataType& rVariables);
 
    virtual void CalculateAndAddKvvm(MatrixType& rLeftHandSideMatrix,
                                     ElementDataType& rVariables);
 
    virtual void CalculateAndAddKvvg(MatrixType& rLeftHandSideMatrix,
                                     ElementDataType& rVariables);
 
 
    virtual void CalculateAndAddExternalForces(VectorType& rRightHandSideVector,
                           ElementDataType& rVariables);
 
 
    virtual void CalculateAndAddInternalForces(VectorType& rRightHandSideVector,
                           ElementDataType & rVariables);
 
 
    virtual void SetElementData(ElementDataType& rVariables,
                                ConstitutiveLaw::Parameters& rValues,
                                const int & rPointNumber);
 
 
    virtual void CalculateMaterialResponse(ElementDataType& rVariables,
                                           ConstitutiveLaw::Parameters& rValues,
                                           const int & rPointNumber);
 
    virtual unsigned int GetDofsSize() const;
 
 
    bool IsSliver();
 
    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);
 
    void CalculateVelocityGradient(Matrix& rL,
                                   const Matrix& rDN_DX,
                                   unsigned int step = 0);
 
    void CalculateVelocityGradientVector(Vector& rVector,
                                         const Matrix& rL,
                                         const Matrix& rDN_DX,
                                         unsigned int step = 0);
    void CalculateVelocityGradientVector(Vector& rVector,
                                         const Matrix& rDN_DX,
                                         unsigned int step = 0);
 
    void CalculateSymmetricVelocityGradient(const Matrix& rL,
                                            Vector& rStrainVector);
 
    void CalculateSkewSymmetricVelocityGradient(const Matrix& rL,
                                                Vector& rStrainVector);
 
    virtual double& CalculateIntegrationWeight(double& rIntegrationWeight);
 
 
    virtual double& CalculateTotalMass(double& rTotalMass, const ProcessInfo& rCurrentProcessInfo);
 
 
    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
 
    virtual void save(Serializer& rSerializer) const override;
 
    virtual void load(Serializer& rSerializer) override;
 
 
 
}; // Class FluidElement
 
 
} // namespace Kratos.
#endif // KRATOS_FLUID_ELEMENT_H_INCLUDED  defined