shell_3p_element.h

Go to the documentation of this file.

#if !defined(KRATOS_SHELL_3P_ELEMENT_H_INCLUDED )
#define  KRATOS_SHELL_3P_ELEMENT_H_INCLUDED
 
 
// System includes
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/element.h"
#include "utilities/math_utils.h"
 
// Application includes
#include "iga_application_variables.h"
 
namespace Kratos
{
 
class Shell3pElement
    : public Element
{
protected:
 
    struct KinematicVariables
    {
        // covariant metric
        array_1d<double, 3> a_ab_covariant;
        array_1d<double, 3> b_ab_covariant;
 
        //base vector 1
        array_1d<double, 3> a1;
        //base vector 2
        array_1d<double, 3> a2;
        //base vector 3 normalized
        array_1d<double, 3> a3;
        //not-normalized base vector 3
        array_1d<double, 3> a3_tilde;
 
        //differential area
        double dA;
 
        KinematicVariables(SizeType Dimension)
        {
            noalias(a_ab_covariant) = ZeroVector(Dimension);
            noalias(b_ab_covariant) = ZeroVector(Dimension);
 
            noalias(a1) = ZeroVector(Dimension);
            noalias(a2) = ZeroVector(Dimension);
            noalias(a3) = ZeroVector(Dimension);
 
            noalias(a3_tilde) = ZeroVector(Dimension);
 
            dA = 1.0;
        }
    };
 
    struct ConstitutiveVariables
    {
        Vector StrainVector;
        Vector StressVector;
        Matrix ConstitutiveMatrix;
 
        ConstitutiveVariables(SizeType StrainSize)
        {
            StrainVector = ZeroVector(StrainSize);
            StressVector = ZeroVector(StrainSize);
            ConstitutiveMatrix = ZeroMatrix(StrainSize, StrainSize);
        }
    };
 
    struct SecondVariations
    {
        Matrix B11;
        Matrix B22;
        Matrix B12;
 
        SecondVariations(const int& mat_size)
        {
            B11 = ZeroMatrix(mat_size, mat_size);
            B22 = ZeroMatrix(mat_size, mat_size);
            B12 = ZeroMatrix(mat_size, mat_size);
        }
    };
 
public:
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(Shell3pElement);
 
    typedef std::size_t SizeType;
    typedef std::size_t IndexType;
 
    // GometryType
    typedef Geometry<Node> GeometryType;
 
 
    Shell3pElement(
        IndexType NewId,
        GeometryType::Pointer pGeometry)
        : Element(NewId, pGeometry)
    {};
 
    Shell3pElement(
        IndexType NewId,
        GeometryType::Pointer pGeometry,
        PropertiesType::Pointer pProperties)
        : Element(NewId, pGeometry, pProperties)
    {};
 
    Shell3pElement()
        : Element()
    {};
 
    virtual ~Shell3pElement() = default;
 
 
    Element::Pointer Create(
        IndexType NewId,
        GeometryType::Pointer pGeom,
        PropertiesType::Pointer pProperties
    ) const override
    {
        return Kratos::make_intrusive<Shell3pElement>(
            NewId, pGeom, pProperties);
    };
 
    Element::Pointer Create(
        IndexType NewId,
        NodesArrayType const& ThisNodes,
        PropertiesType::Pointer pProperties
    ) const override
    {
        return Kratos::make_intrusive< Shell3pElement >(
            NewId, GetGeometry().Create(ThisNodes), pProperties);
    };
 
 
    void CalculateRightHandSide(
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo) override
    {
        const SizeType number_of_nodes = GetGeometry().size();
        const SizeType mat_size = number_of_nodes * 3;
 
        if (rRightHandSideVector.size() != mat_size)
            rRightHandSideVector.resize(mat_size);
        noalias(rRightHandSideVector) = ZeroVector(mat_size);
 
        MatrixType left_hand_side_matrix;
 
        CalculateAll(left_hand_side_matrix, rRightHandSideVector,
            rCurrentProcessInfo, false, true);
    }
 
    void CalculateLeftHandSide(
        MatrixType& rLeftHandSideMatrix,
        const ProcessInfo& rCurrentProcessInfo) override
    {
        const SizeType number_of_nodes = GetGeometry().size();
        const SizeType mat_size = number_of_nodes * 3;
 
        VectorType right_hand_side_vector;
 
        if (rLeftHandSideMatrix.size1() != mat_size)
            rLeftHandSideMatrix.resize(mat_size, mat_size);
        noalias(rLeftHandSideMatrix) = ZeroMatrix(mat_size, mat_size);
 
        CalculateAll(rLeftHandSideMatrix, right_hand_side_vector,
            rCurrentProcessInfo, true, false);
    }
 
    void CalculateLocalSystem(
        MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo) override
    {
        const SizeType number_of_nodes = GetGeometry().size();
        const SizeType mat_size = number_of_nodes * 3;
 
        if (rRightHandSideVector.size() != mat_size)
            rRightHandSideVector.resize(mat_size);
        noalias(rRightHandSideVector) = ZeroVector(mat_size);
 
        if (rLeftHandSideMatrix.size1() != mat_size)
            rLeftHandSideMatrix.resize(mat_size, mat_size);
        noalias(rLeftHandSideMatrix) = ZeroMatrix(mat_size, mat_size);
 
        CalculateAll(rLeftHandSideMatrix, rRightHandSideVector,
            rCurrentProcessInfo, true, true);
    }
 
    void CalculateMassMatrix(
        MatrixType& rMassMatrix,
        const ProcessInfo& rCurrentProcessInfo
    ) override;
 
    void CalculateDampingMatrix(
        MatrixType& rDampingMatrix,
        const ProcessInfo& rCurrentProcessInfo
    ) override;
 
    void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void EquationIdVector(
        EquationIdVectorType& rResult,
        const ProcessInfo& rCurrentProcessInfo
    ) const override;
 
    void GetDofList(
        DofsVectorType& rElementalDofList,
        const ProcessInfo& rCurrentProcessInfo
    ) const override;
 
 
    void Initialize(const ProcessInfo& rCurrentProcessInfo) override;
 
    void GetValuesVector(
        Vector& rValues,
        int Step) const override;
 
    void GetFirstDerivativesVector(
        Vector& rValues,
        int Step) const override;
 
    void GetSecondDerivativesVector(
        Vector& rValues,
        int Step) const override;
 
    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>>& rOutput,
        const ProcessInfo& rCurrentProcessInfo
    ) override;
 
 
    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Kirchhoff-Love Shell3pElement #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "Kirchhoff-Love Shell3pElement #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        pGetGeometry()->PrintData(rOStream);
    }
 
 
private:
 
    // Components of the metric coefficient tensor on the contravariant basis
    std::vector<array_1d<double, 3>> m_A_ab_covariant_vector;
    // Components of the curvature coefficient tensor on the contravariant basis
    std::vector<array_1d<double, 3>> m_B_ab_covariant_vector;
 
    // Determinant of the geometrical Jacobian.
    Vector m_dA_vector;
 
    /* Transformation the strain tensor from the curvilinear system
    *  to the local cartesian in voigt notation including a 2 in the
    *  shear part. */
    std::vector<Matrix> m_T_vector;
 
    std::vector<ConstitutiveLaw::Pointer> mConstitutiveLawVector;
 
 
    void CalculateAll(
        MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo,
        const bool CalculateStiffnessMatrixFlag,
        const bool CalculateResidualVectorFlag
    ) const;
 
    void InitializeMaterial();
 
    void CalculateKinematics(
        const IndexType IntegrationPointIndex,
        KinematicVariables& rKinematicVariables) const;
 
    // Computes transformation
    void CalculateTransformation(
        const KinematicVariables& rKinematicVariables,
        Matrix& rT) const;
 
    // Computes transformation for the stress tensor 
    void CalculateTransformationFromCovariantToCartesian(
        const KinematicVariables& rKinematicVariables,
        Matrix& rTCovToCar) const;
 
    void CalculateBMembrane(
        const IndexType IntegrationPointIndex,
        Matrix& rB,
        const KinematicVariables& rActualKinematic) const;
 
    void CalculateBCurvature(
        const IndexType IntegrationPointIndex,
        Matrix& rB,
        const KinematicVariables& rActualKinematic) const;
 
    void CalculateSecondVariationStrainCurvature(
        const IndexType IntegrationPointIndex,
        SecondVariations& rSecondVariationsStrain,
        SecondVariations& rSecondVariationsCurvature,
        const KinematicVariables& rActualKinematic) const;
 
    void CalculateConstitutiveVariables(
        const IndexType IntegrationPointIndex,
        KinematicVariables& rActualMetric,
        ConstitutiveVariables& rThisConstitutiveVariablesMembrane,
        ConstitutiveVariables& rThisConstitutiveVariablesCurvature,
        ConstitutiveLaw::Parameters& rValues,
        const ConstitutiveLaw::StressMeasure ThisStressMeasure
    ) const;
 
    inline void CalculateAndAddKm(
        MatrixType& rLeftHandSideMatrix,
        const Matrix& B,
        const Matrix& D,
        const double IntegrationWeight) const;
 
    inline void CalculateAndAddNonlinearKm(
        Matrix& rLeftHandSideMatrix,
        const SecondVariations& rSecondVariationsStrain,
        const Vector& rSD,
        const double IntegrationWeight) const;
 
    // Calculation of the PK2 stress
    void CalculatePK2Stress(
        const IndexType IntegrationPointIndex,
        array_1d<double, 3>& rPK2MembraneStressCartesian,
        array_1d<double, 3>& rPK2BendingStressCartesian,
        const ProcessInfo& rCurrentProcessInfo) const;
 
    // Calculation of the Cauchy stress by transforming the PK2 stress
    void CalculateCauchyStress(
        const IndexType IntegrationPointIndex,
        array_1d<double, 3>& rCauchyMembraneStressesCartesian, 
        array_1d<double, 3>& rCauchyBendingStressesCartesian, 
        const ProcessInfo& rCurrentProcessInfo) const;
 
    // Calculation of the shear force, shear force = derivative of moment
    void CalculateShearForce(
        const IndexType IntegrationPointIndex,
        array_1d<double, 2>& rq, 
        const ProcessInfo& rCurrentProcessInfo) const;
 
    void CalculateDerivativeOfCurvatureInitial(
        const IndexType IntegrationPointIndex,
        array_1d<double, 3>& rDCurvature_D1,
        array_1d<double, 3>& rDCurvature_D2,
        const Matrix& rHessian) const;
 
    void CalculateDerivativeOfCurvatureActual(
        const IndexType IntegrationPointIndex,
        array_1d<double, 3>& rDCurvature_D1,
        array_1d<double, 3>& rDCurvature_D2,
        const Matrix& rHessian,
        const KinematicVariables& rKinematicVariables) const;
 
    void CalculateDerivativeTransformationMatrices(
        const IndexType IntegrationPointIndex,
        std::vector<Matrix>& rDQ_Dalpha_init,
        std::vector<Matrix>& rDTransCartToCov_Dalpha_init,
        const Matrix& rHessian) const;
 
    template<class TType>
    void GetValueOnConstitutiveLaw(
        const Variable<TType>& rVariable,
        std::vector<TType>& rOutput
        )
    {
        const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints(this->GetIntegrationMethod());
 
        for (IndexType point_number = 0; point_number < integration_points.size(); ++point_number) {
            mConstitutiveLawVector[point_number]->GetValue(rVariable, rOutput[point_number]);
        }
    }
 
 
    void CalculateHessian(
        Matrix& Hessian,
        const Matrix& rDDN_DDe) const;
 
    void CalculateSecondDerivativesOfBaseVectors(
        const Matrix& rDDDN_DDDe,
        array_1d<double, 3>& rDDa1_DD11,
        array_1d<double, 3>& rDDa1_DD12,
        array_1d<double, 3>& rDDa2_DD21,
        array_1d<double, 3>& rDDa2_DD22) const;
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Element);
        rSerializer.save("A_ab_covariant_vector", m_A_ab_covariant_vector);
        rSerializer.save("B_ab_covariant_vector", m_B_ab_covariant_vector);
        rSerializer.save("dA_vector", m_dA_vector);
        rSerializer.save("T_vector", m_T_vector);
        rSerializer.save("constitutive_law_vector", mConstitutiveLawVector);
    }
 
    void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element);
        rSerializer.load("A_ab_covariant_vector", m_A_ab_covariant_vector);
        rSerializer.load("B_ab_covariant_vector", m_B_ab_covariant_vector);
        rSerializer.load("dA_vector", m_dA_vector);
        rSerializer.load("T_vector", m_T_vector);
        rSerializer.load("constitutive_law_vector", mConstitutiveLawVector);
    }
 
 
};     // Class Shell3pElement
 
}  // namespace Kratos.
 
#endif // KRATOS_MESHLESS_SHELL_3P_ELEMENT_H_INCLUDED  defined