generic_total_lagrangian_femdem_element.h

Go to the documentation of this file.

// KRATOS  ___|  |                   |                   |
//       \___ \  __|  __| |   |  __| __| |   |  __| _` | |
//             | |   |    |   | (    |   |   | |   (   | |
//       _____/ \__|_|   \__,_|\___|\__|\__,_|_|  \__,_|_| MECHANICS
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//                   Vicente Mataix Ferrandiz
//                   Alejandro Cornejo
//
 
 
#pragma once
 
 
// System includes
 
 
// External include
 
// Project includes
#include "includes/define.h"
#include "custom_elements/base_solid_element.h"
#include "includes/variables.h"
#include "custom_utilities/constitutive_law_utilities.h"
 
namespace Kratos
{
 
 
 
template<unsigned int TDim, unsigned int TyieldSurf>
class KRATOS_API(FEM_TO_DEM_APPLICATION) GenericTotalLagrangianFemDemElement
    : public BaseSolidElement
{
public:
 
    typedef ConstitutiveLaw ConstitutiveLawType;
 
    typedef ConstitutiveLawType::Pointer ConstitutiveLawPointerType;
 
    typedef GeometryData::IntegrationMethod IntegrationMethod;
 
    typedef BaseSolidElement BaseType;
 
    typedef std::size_t IndexType;
 
    typedef std::size_t SizeType;
 
    typedef Geometry<NodeType>::PointsArrayType NodesArrayType;
 
    typedef Properties PropertiesType;
 
    typedef Geometry<NodeType> GeometryType;
 
    static constexpr double tolerance = std::numeric_limits<double>::epsilon();
 
    static constexpr SizeType VoigtSize = (TDim == 3) ? 6 : 3;
 
    static constexpr SizeType NumberOfEdges = (TDim == 3) ? 6 : 3;
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(GenericTotalLagrangianFemDemElement);
 
 
    GenericTotalLagrangianFemDemElement(IndexType NewId, GeometryType::Pointer pGeometry);
    GenericTotalLagrangianFemDemElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties);
 
    // Copy constructor
    GenericTotalLagrangianFemDemElement(GenericTotalLagrangianFemDemElement const& rOther)
        : BaseSolidElement(rOther)
    {};
 
    ~GenericTotalLagrangianFemDemElement() override;
 
 
    Element::Pointer Create(
        IndexType NewId,
        GeometryType::Pointer pGeom,
        PropertiesType::Pointer pProperties
        ) const override;
 
    Element::Pointer Create(
        IndexType NewId,
        NodesArrayType const& ThisNodes,
        PropertiesType::Pointer pProperties
        ) const override;
 
    Element::Pointer Clone(
        IndexType NewId,
        NodesArrayType const& rThisNodes
        ) const override;
 
    // int Check(const ProcessInfo& rCurrentProcessInfo) override;
 
    //std::string Info() const;
 
    void CalculateSensitivityMatrix(const Variable<array_1d<double, 3>>& rDesignVariable,
                                    Matrix& rOutput,
                                    const ProcessInfo& rCurrentProcessInfo) override;
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Total Lagrangian FEMDEM Solid Element #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "Total Lagrangian FEMDEM Solid Element #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        pGetGeometry()->PrintData(rOStream);
    }
 
 
protected:
 
    Vector mThresholds;                 // Stress mThreshold on edge
    Vector mDamages;                    // Converged Damage on each edge
 
    double mThreshold = 0.0;            // Converged Threshold
    double mDamage = 0.0;               // Converged Damage
 
 
    GenericTotalLagrangianFemDemElement() : BaseSolidElement()
    {
    }
 
    void CalculateAll(
        MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo,
        const bool CalculateStiffnessMatrixFlag,
        const bool CalculateResidualVectorFlag) override;
 
    void CalculateKinematicVariables(
        KinematicVariables& rThisKinematicVariables,
        const IndexType PointNumber,
        const GeometryType::IntegrationMethod& rIntegrationMethod) override;
 
    // ************** Methods to compute the tangent constitutive tensor via numerical derivation ************** 
    void CalculateTangentTensor(Matrix& rTangentTensor, const Vector& rStrainVectorGP, const Vector& rStressVectorGP, const Matrix& rDeformationGradientGP, const Matrix& rElasticMatrix, ConstitutiveLaw::Parameters& rValues);
    void CalculateTangentTensorSecondOrder(Matrix& rTangentTensor, const Vector& rStrainVectorGP, const Vector& rStressVectorGP, const Matrix& rDeformationGradientGP, const Matrix& rElasticMatrix, ConstitutiveLaw::Parameters& rValues);
 
    void CalculatePerturbation(const Vector& rStrainVectorGP, double& rPerturbation, const int Component);
 
    void PerturbateStrainVector(Vector& rPerturbedStrainVector, const Vector& rStrainVectorGP, const double Perturbation, const int Component);
 
    virtual void IntegratePerturbedStrain(Vector& rPerturbedStressVector, const Vector& rPerturbedStrainVector, const Matrix& rElasticMatrix, ConstitutiveLaw::Parameters& rValues);
 
    void AssignComponentsToTangentTensor(Matrix& rTangentTensor, const Vector& rDeltaStress, const double Perturbation, const int Component);
 
    void PerturbateDeformationGradient(Matrix& rPerturbedDeformationGradient, const Matrix& rDeformationGradientGP, const double Perturbation, const int ComponentI, const int ComponentJ);
 
    int CalculateVoigtIndex(const SizeType VoigtSize, const int ComponentI, const int ComponentJ);
 
    double CalculateElementalDamage(const Vector& rEdgeDamages);
    double CalculateElementalDamage3D(const Vector& rEdgeDamages);
    double CalculateElementalDamage2D(const Vector& rEdgeDamages);
 
    void IntegrateStressDamageMechanics(double& rThreshold,double& rDamage, const Vector& rStrainVector,
        const Vector& rStressVector, const int Edge, const double CharacteristicLength,
        ConstitutiveLaw::Parameters& rValues, bool& rIsDamaging);
 
    virtual Vector IntegrateSmoothedConstitutiveLaw(const std::string &rYieldSurface, ConstitutiveLaw::Parameters &rValues,
                                             const ConstitutiveVariables &rThisConstVars, const KinematicVariables &rKinVariables, 
                                             Vector &rStrainVector, double& rDamageElement,  bool& rIsDamaging, const double CharacteristicLength,
                                             const bool SaveIntVars);
 
    void ComputeEdgeNeighbours(const ProcessInfo& rCurrentProcessInfo);
 
    void AuxComputeEdgeNeighbours(const ProcessInfo& rCurrentProcessInfo);
 
    std::vector<Element*> GetEdgeNeighbourElements(const int edge) {return mEdgeNeighboursContainer[edge];}
 
    void SaveEdgeNeighboursContainer(const std::vector<std::vector<Element*>>& rtoSave) {mEdgeNeighboursContainer = rtoSave;}
 
    void SetNodeIndexes(Matrix& rMatrix)
    {
        rMatrix.resize(6, 2);
        rMatrix(0, 0) = 0; rMatrix(0, 1) = 1; rMatrix(1, 0) = 0;
        rMatrix(1, 1) = 2; rMatrix(2, 0) = 0; rMatrix(2, 1) = 3;
        rMatrix(3, 0) = 1; rMatrix(3, 1) = 2; rMatrix(4, 0) = 1;
        rMatrix(4, 1) = 3; rMatrix(5, 0) = 2; rMatrix(5, 1) = 3;
    }
 
 
    void CalculateAverageVariableOnEdge(const Element* pCurrentElement,   const Variable<Vector>& rThisVariable, Vector& rAverageStress, const int edge);
    void CalculateAverageVariableOnEdge2D(const Element* pCurrentElement, const Variable<Vector>& rThisVariable, Vector& rAverageStress, const int edge);
    void CalculateAverageVariableOnEdge3D(const Element* pCurrentElement, const Variable<Vector>& rThisVariable, Vector& rAverageStress, const int edge);
 
    void CalculateEquivalentStress(const array_1d<double, VoigtSize>& rPredictiveStressVector, const Vector& rStrainVector, double& rEquivalentStress, ConstitutiveLaw::Parameters& rValues);
 
    void GetInitialUniaxialThreshold(ConstitutiveLaw::Parameters& rValues, double& rThreshold);
 
    void CalculateDamageParameter(ConstitutiveLaw::Parameters& rValues, double& rAParameter, const double CharacteristicLength);
 
    double CalculateCharacteristicLength(GenericTotalLagrangianFemDemElement *pCurrentElement);
 
    void CalculateExponentialDamage(double& rDamage, const double DamageParameter, const double UniaxialStress, const double InitialThrehsold);
 
    void CalculateGreenLagrangeStrainVector(Vector& rStrainVector, const Matrix& rF);
 
    std::size_t GetStrainSize() const;
 
    void SetValuesOnIntegrationPoints(
        const Variable<double> &rVariable,
        const std::vector<double> &rValues,
        const ProcessInfo &rCurrentProcessInfo) override;
 
    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;
 
 
    virtual void CheckIfEraseElement(
        const ProcessInfo &rCurrentProcessInfo, 
        const Properties& rProperties)
    {
        const double max_damage = 0.98;
        if (mDamage >= max_damage) {
            this->Set(ACTIVE, false);
            mDamage = max_damage;
            this->SetValue(GENERATE_DEM, true);
        }
    }
 
 
private:
 
 
 
 
    void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override;
 
 
    void CalculateB(Matrix& rB, Matrix const& rF, const Matrix& rDN_DX);
 
    void Calculate2DB(Matrix& rB, const Matrix& rF, const Matrix& rDN_DX);
 
    void Calculate3DB(Matrix& rB, const Matrix& rF, const Matrix& rDN_DX);
 
    void CalculateStress(Vector& rStrain,
                         std::size_t IntegrationPoint,
                         Vector& rStress,
                         ProcessInfo const& rCurrentProcessInfo);
 
    void CalculateStress(Matrix const& rF,
                         std::size_t IntegrationPoint,
                         Vector& rStress,
                         ProcessInfo const& rCurrentProcessInfo);
 
    void CalculateStrain(Matrix const& rF,
                         std::size_t IntegrationPoint,
                         Vector& rStrain,
                         ProcessInfo const& rCurrentProcessInfo);
 
    void CalculateShapeSensitivity(ShapeParameter Deriv,
                                   Matrix& rDN_DX0,
                                   Matrix& rDN_DX0_Deriv,
                                   Matrix& rF_Deriv,
                                   double& rDetJ0_Deriv,
                                   std::size_t IntegrationPointIndex);
 
    void CalculateBSensitivity(Matrix const& rDN_DX,
                               Matrix const& rF,
                               Matrix const& rDN_DX_Deriv,
                               Matrix const& rF_Deriv,
                               Matrix& rB_Deriv);
 
 
    // Vector to storage the neigh elements sharing a certain edge
    std::vector<std::vector<Element*>> mEdgeNeighboursContainer;
 
    friend class Serializer;
 
    // A private default constructor necessary for serialization
 
    void save(Serializer& rSerializer) const override;
 
    void load(Serializer& rSerializer) override;
 
}; // Class GenericTotalLagrangianFemDemElement
 
template<unsigned int TDim, unsigned int TyieldSurf> constexpr SizeType GenericTotalLagrangianFemDemElement<TDim, TyieldSurf>::VoigtSize;
template<unsigned int TDim, unsigned int TyieldSurf> constexpr SizeType GenericTotalLagrangianFemDemElement<TDim, TyieldSurf>::NumberOfEdges;
 
 
 
} // namespace Kratos.