mpm_updated_lagrangian.hpp

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:        BSD License
//                  Kratos default license: kratos/license.txt
//
//  Main authors:    Ilaria Iaconeta, Bodhinanda Chandra
//
 
 
#if !defined(KRATOS_UPDATED_LAGRANGIAN_H_INCLUDED )
#define  KRATOS_UPDATED_LAGRANGIAN_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/element.h"
#include "includes/serializer.h"
#include "includes/ublas_interface.h"
#include "includes/variables.h"
#include "includes/constitutive_law.h"
 
namespace Kratos
{
 
 
class MPMUpdatedLagrangian
    : public Element
{
public:
 
    typedef ConstitutiveLaw ConstitutiveLawType;
    typedef ConstitutiveLawType::Pointer ConstitutiveLawPointerType;
    typedef ConstitutiveLawType::StressMeasure StressMeasureType;
    typedef GeometryData::IntegrationMethod IntegrationMethod;
 
    typedef typename GeometryType::CoordinatesArrayType CoordinatesArrayType;
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( MPMUpdatedLagrangian );
 
protected:
 
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_RHS_VECTOR );
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_LHS_MATRIX );
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_RHS_VECTOR_WITH_COMPONENTS );
    KRATOS_DEFINE_LOCAL_FLAG( COMPUTE_LHS_MATRIX_WITH_COMPONENTS );
 
    struct MaterialPointVariables
    {
    public:
        // Particle Position
        CoordinatesArrayType xg;
        // MP_MASS
        double mass;
        // MP_DENSITY
        double density;
        // MP_VOLUME
        double volume;
 
        // MP_DISPLACEMENT
        array_1d<double, 3> displacement;
        // MP_VELOCITY
        array_1d<double, 3> velocity;
        // MP_ACCELERATION
        array_1d<double, 3> acceleration;
 
        // MP_VOLUME_ACCELERATION
        array_1d<double, 3> volume_acceleration;
 
        // MP_CAUCHY_STRESS_VECTOR
        Vector cauchy_stress_vector;
        // MP_ALMANSI_STRAIN_VECTOR
        Vector almansi_strain_vector;
 
        // MP_DELTA_PLASTIC_STRAIN
        double delta_plastic_strain;
        // MP_DELTA_PLASTIC_VOLUMETRIC_STRAIN
        double delta_plastic_volumetric_strain;
        // MP_DELTA_PLASTIC_DEVIATORIC_STRAIN
        double delta_plastic_deviatoric_strain;
        // MP_EQUIVALENT_PLASTIC_STRAIN
        double equivalent_plastic_strain;
        // MP_ACCUMULATED_PLASTIC_VOLUMETRIC_STRAIN
        double accumulated_plastic_volumetric_strain;
        // MP_ACCUMULATED_PLASTIC_DEVIATORIC_STRAIN
        double accumulated_plastic_deviatoric_strain;
 
        explicit MaterialPointVariables()
        {
            // MP_MASS
            mass = 1.0;
            // MP_DENSITY
            density = 1.0;
            // MP_VOLUME
            volume = 1.0;
 
            // MP_DELTA_PLASTIC_STRAIN
            delta_plastic_strain = 1.0;
            // MP_DELTA_PLASTIC_VOLUMETRIC_STRAIN
            delta_plastic_volumetric_strain = 1.0;
            // MP_DELTA_PLASTIC_DEVIATORIC_STRAIN
            delta_plastic_deviatoric_strain = 1.0;
            // MP_EQUIVALENT_PLASTIC_STRAIN
            equivalent_plastic_strain = 1.0;
            // MP_ACCUMULATED_PLASTIC_VOLUMETRIC_STRAIN
            accumulated_plastic_volumetric_strain = 1.0;
            // MP_ACCUMULATED_PLASTIC_DEVIATORIC_STRAIN
            accumulated_plastic_deviatoric_strain = 1.0;
        }
 
    private:
 
        friend class Serializer;
 
        void save( Serializer& rSerializer ) const
        {
            rSerializer.save("xg",xg);
            rSerializer.save("mass",mass);
            rSerializer.save("density",density);
            rSerializer.save("volume",volume);
            rSerializer.save("displacement",displacement);
            rSerializer.save("velocity",velocity);
            rSerializer.save("acceleration",acceleration);
            rSerializer.save("volume_acceleration",volume_acceleration);
            rSerializer.save("cauchy_stress_vector",cauchy_stress_vector);
            rSerializer.save("almansi_strain_vector",almansi_strain_vector);
            rSerializer.save("delta_plastic_strain",delta_plastic_strain);
            rSerializer.save("delta_plastic_volumetric_strain",delta_plastic_volumetric_strain);
            rSerializer.save("delta_plastic_deviatoric_strain",delta_plastic_deviatoric_strain);
            rSerializer.save("equivalent_plastic_strain",equivalent_plastic_strain);
            rSerializer.save("accumulated_plastic_volumetric_strain",accumulated_plastic_volumetric_strain);
            rSerializer.save("accumulated_plastic_deviatoric_strain",accumulated_plastic_deviatoric_strain);
        }
 
        void load( Serializer& rSerializer )
        {
            rSerializer.load("xg",xg);
            rSerializer.load("mass",mass);
            rSerializer.load("density",density);
            rSerializer.load("volume",volume);
            rSerializer.load("displacement",displacement);
            rSerializer.load("velocity",velocity);
            rSerializer.load("acceleration",acceleration);
            rSerializer.load("volume_acceleration",volume_acceleration);
            rSerializer.load("cauchy_stress_vector",cauchy_stress_vector);
            rSerializer.load("almansi_strain_vector",almansi_strain_vector);
            rSerializer.load("delta_plastic_strain",delta_plastic_strain);
            rSerializer.load("delta_plastic_volumetric_strain",delta_plastic_volumetric_strain);
            rSerializer.load("delta_plastic_deviatoric_strain",delta_plastic_deviatoric_strain);
            rSerializer.load("equivalent_plastic_strain",equivalent_plastic_strain);
            rSerializer.load("accumulated_plastic_volumetric_strain",accumulated_plastic_volumetric_strain);
            rSerializer.load("accumulated_plastic_deviatoric_strain",accumulated_plastic_deviatoric_strain);
        }
    };
 
    struct GeneralVariables
    {
    public:
 
        StressMeasureType StressMeasure;
 
        // For axisymmetric use only
        double  CurrentRadius;
        double  ReferenceRadius;
 
        // General variables for large displacement use
        double  detF;
        double  detF0;
        double  detFT;
        Vector  StrainVector;
        Vector  StressVector;
        Matrix  B;
        Matrix  F;
        Matrix  FT;
        Matrix  F0;
        Matrix  DN_DX;
        Matrix  ConstitutiveMatrix;
 
        // Variables including all integration points
        Matrix CurrentDisp;
    };
 
public:
 
 
 
    MPMUpdatedLagrangian();
 
 
    MPMUpdatedLagrangian(IndexType NewId, GeometryType::Pointer pGeometry);
 
    MPMUpdatedLagrangian(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties);
 
    MPMUpdatedLagrangian(MPMUpdatedLagrangian const& rOther);
 
    ~MPMUpdatedLagrangian() override;
 
 
    MPMUpdatedLagrangian& operator=(MPMUpdatedLagrangian const& rOther);
 
 
    Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const override;
 
    Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override;
 
    Element::Pointer Clone(IndexType NewId, NodesArrayType const& ThisNodes) const override;
 
 
    //************* GETTING METHODS
 
    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;
 
    virtual void InitializeSolutionStep(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 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;
 
    //************************************************************************************
    //************************************************************************************
    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "MPM Element #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "MPM Element #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        GetGeometry().PrintData(rOStream);
    }
 
 
    void CalculateOnIntegrationPoints(const Variable<bool>& rVariable,
        std::vector<bool>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateOnIntegrationPoints(const Variable<int>& rVariable,
        std::vector<int>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateOnIntegrationPoints(const Variable<double>& rVariable,
        std::vector<double>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateOnIntegrationPoints(const Variable<array_1d<double, 3 > >& rVariable,
        std::vector<array_1d<double, 3 > >& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
        std::vector<Vector>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
 
    void SetValuesOnIntegrationPoints(const Variable<int>& rVariable,
        const std::vector<int>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void SetValuesOnIntegrationPoints(const Variable<double>& rVariable,
        const std::vector<double>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void SetValuesOnIntegrationPoints(const Variable<array_1d<double, 3 > >& rVariable,
        const std::vector<array_1d<double, 3 > >& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void SetValuesOnIntegrationPoints(const Variable<Vector>& rVariable,
        const std::vector<Vector>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
 
protected:
 
    MaterialPointVariables mMP;
 
    Matrix mDeformationGradientF0;
    double mDeterminantF0;
 
    ConstitutiveLaw::Pointer mConstitutiveLawVector;
 
 
    bool mFinalizedStep;
 
 
    virtual SizeType GetNumberOfDofs() {
        return GetGeometry().WorkingSpaceDimension();
    }
 
    virtual void CalculateElementalSystem(
        MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo,
        const bool CalculateStiffnessMatrixFlag,
        const bool CalculateResidualVectorFlag);
 
 
 
    virtual void CalculateAndAddLHS(
        MatrixType& rLeftHandSideMatrix,
        GeneralVariables& rVariables,
        const double& rIntegrationWeight,
        const ProcessInfo& rCurrentProcessInfo);
 
    virtual void CalculateAndAddRHS(
        VectorType& rRightHandSideVector,
        GeneralVariables& rVariables,
        Vector& rVolumeForce,
        const double& rIntegrationWeight,
        const ProcessInfo& rCurrentProcessInfo);
 
 
    virtual void CalculateAndAddKuum(MatrixType& rLeftHandSideMatrix,
                                     GeneralVariables& rVariables,
                                     const double& rIntegrationWeight);
 
    virtual void CalculateAndAddKuug(MatrixType& rLeftHandSideMatrix,
                                     GeneralVariables& rVariables,
                                     const double& rIntegrationWeight,
                                     const bool IsAxisymmetric = false);
 
 
    virtual void CalculateAndAddExternalForces(VectorType& rRightHandSideVector,
            GeneralVariables& rVariables,
            Vector& rVolumeForce,
            const double& rIntegrationWeight);
 
 
    virtual void CalculateAndAddInternalForces(VectorType& rRightHandSideVector,
            GeneralVariables & rVariables,
            const double& rIntegrationWeight);
 
    virtual void CalculateExplicitStresses(const ProcessInfo& rCurrentProcessInfo,
        GeneralVariables& rVariables);
 
 
    virtual void SetGeneralVariables(GeneralVariables& rVariables,
                                     ConstitutiveLaw::Parameters& rValues, const Vector& rN);
 
 
    virtual void InitializeMaterial (const ProcessInfo& rCurrentProcessInfo);
 
 
    void ResetConstitutiveLaw() override;
 
 
    void ClearNodalForces ();
 
    virtual void CalculateKinematics(GeneralVariables& rVariables, const ProcessInfo& rCurrentProcessInfo);
 
 
    Matrix& CalculateCurrentDisp(Matrix & rCurrentDisp, const ProcessInfo& rCurrentProcessInfo);
 
    void DecimalCorrection(Vector& rVector);
 
 
    virtual void InitializeGeneralVariables(GeneralVariables & rVariables, const ProcessInfo& rCurrentProcessInfo);
 
 
    virtual void FinalizeStepVariables(GeneralVariables & rVariables, const ProcessInfo& rCurrentProcessInfo);
 
    virtual void UpdateGaussPoint(GeneralVariables & rVariables, const ProcessInfo& rCurrentProcessInfo);
 
    virtual void GetHistoricalVariables( GeneralVariables& rVariables);
 
 
    virtual void CalculateGreenLagrangeStrain(const Matrix& rF,
            Vector& rStrainVector);
 
    virtual void CalculateAlmansiStrain(const Matrix& rF,
                                        Vector& rStrainVector);
 
 
    virtual void CalculateDeformationMatrix(Matrix& rB,
                                            const Matrix& rDN_DX,
                                            const Matrix& rN,
                                            const bool IsAxisymmetric = false);
 
    virtual double& CalculateIntegrationWeight(double& rIntegrationWeight);
 
 
    virtual double& CalculateVolumeChange(double& rVolumeChange, GeneralVariables& rVariables);
 
 
    void CalculateDeformationGradient(const Matrix& rDN_DX, Matrix& rF, Matrix& rDeltaPosition,
        const bool IsAxisymmetric = false);
 
 
private:
 
 
 
 
 
 
 
 
    friend class Serializer;
 
    // A private default constructor necessary for serialization
 
    void save(Serializer& rSerializer) const override;
 
    void load(Serializer& rSerializer) override;
 
 
 
}; // Class MPMUpdatedLagrangian
 
 
} // namespace Kratos.
#endif // KRATOS_UPDATED_LAGRANGIAN_H_INCLUDED  defined