shell_5p_element.h

Go to the documentation of this file.

#if !defined(KRATOS_SHELL_5P_ELEMENT_H_INCLUDED )
#define  KRATOS_SHELL_5P_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 Shell5pElement final
    : public Element
{
protected:
 
    struct KinematicVariables
    {
        // covariant metric
        array_1d<double, 3> metricChange;
        array_1d<double, 3> curvature;
        array_1d<double, 2> transShear; //shear
 
        //base vectors current and reference
        BoundedVector<double, 3> a1;
        BoundedVector<double, 3> a2;
        BoundedVector<double, 3> A1;
        BoundedVector<double, 3> A2;
 
        BoundedVector<double, 3> dud1;
        BoundedVector<double, 3> dud2;
        //director
        BoundedVector<double, 3> t;
        //array_1d<double, 3> t;
 
        //director derivatives wrt physical space
        array_1d<double, 3> dtd1;
        array_1d<double, 3> dtd2;
 
        //differential area
        double dA;
 
        void setZero()
        {
            metricChange = ZeroVector(3);
            curvature = ZeroVector(3);
            transShear = ZeroVector(2);
 
            a1 = ZeroVector(3);
            a2 = ZeroVector(3);
            A1 = ZeroVector(3);
            A2 = ZeroVector(3);
            t = ZeroVector(3);
            dtd1 = ZeroVector(3);
            dtd2 = ZeroVector(3);
            dud1 = ZeroVector(3);
            dud2 = ZeroVector(3);
 
            dA = 0;
        }
    };
 
    struct ConstitutiveVariables
    {
        Vector StrainVector;
        Vector StressVector;
        Matrix ConstitutiveMatrix;
 
        ConstitutiveVariables(SizeType StrainSize)
        {
            StrainVector = ZeroVector(StrainSize);
            StressVector = ZeroVector(StrainSize);
            ConstitutiveMatrix = ZeroMatrix(StrainSize, StrainSize);
        }
    };
 
    using Matrix3d = BoundedMatrix<double, 3, 3>;
    using Matrix2d = BoundedMatrix<double, 2, 2>;
    using Matrix32d = BoundedMatrix<double, 3, 2>;
    using Matrix23d = BoundedMatrix<double, 2, 3>;
 
    struct VariationVariables
    {
        Matrix3d P;
        Matrix3d Q1;
        Matrix3d Q2;
        Matrix3d S1;
        Matrix3d S2;
        Matrix3d Chi11;
        Matrix3d Chi12Chi21;
        Matrix3d Chi22;
    };
 
public:
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(Shell5pElement);
 
    typedef std::size_t SizeType;
    typedef std::size_t IndexType;
 
    // GometryType
    typedef Node NodeType;
    typedef Geometry<NodeType> GeometryType;
 
 
 
    Shell5pElement(
        IndexType NewId,
        GeometryType::Pointer pGeometry)
        : Element(NewId, pGeometry)
    {};
 
    Shell5pElement(
        IndexType NewId,
        GeometryType::Pointer pGeometry,
        PropertiesType::Pointer pProperties)
        : Element(NewId, pGeometry, pProperties)
    {};
 
    Shell5pElement()
        : Element()
    {};
 
    virtual ~Shell5pElement() final
    {};
 
 
    Element::Pointer Create(
        IndexType NewId,
        GeometryType::Pointer pGeom,
        PropertiesType::Pointer pProperties
    ) const final
    {
        return Kratos::make_intrusive<Shell5pElement>(
            NewId, pGeom, pProperties);
    };
 
    Element::Pointer Create(
        IndexType NewId,
        NodesArrayType const& ThisNodes,
        PropertiesType::Pointer pProperties
    ) const final
    {
        return Kratos::make_intrusive< Shell5pElement >(
            NewId, GetGeometry().Create(ThisNodes), pProperties);
    };
 
 
    void CalculateRightHandSide(
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo) final
    {
        const SizeType number_of_nodes = GetGeometry().size();
        const SizeType mat_size = number_of_nodes * 5;
 
        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) final
    {
        const SizeType number_of_nodes = GetGeometry().size();
        const SizeType mat_size = number_of_nodes * 5;
 
        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) final
    {
        const SizeType number_of_nodes = GetGeometry().size();
        const SizeType mat_size = number_of_nodes * 5;
 
        if (rRightHandSideVector.size() != mat_size)
            rRightHandSideVector.resize(mat_size);
        noalias(rRightHandSideVector) = ZeroVector(mat_size);
 
        if (rLeftHandSideMatrix.size1() != mat_size || rLeftHandSideMatrix.size2() != mat_size)
            rLeftHandSideMatrix.resize(mat_size, mat_size);
        noalias(rLeftHandSideMatrix) = ZeroMatrix(mat_size, mat_size);
 
        CalculateAll(rLeftHandSideMatrix, rRightHandSideVector,
            rCurrentProcessInfo, true, true);
    }
 
    void EquationIdVector(
        EquationIdVectorType& rResult,
        const ProcessInfo& rCurrentProcessInfo
    ) const final;
 
    void GetDofList(
        DofsVectorType& rElementalDofList,
        const ProcessInfo& rCurrentProcessInfo
    ) const final;
 
 
    void Initialize(const ProcessInfo& rCurrentProcessInfo) final;
 
    void InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) final;
 
    void FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) final;
 
    void GetValuesVector(
        Vector& rValues,
        int Step) const final;
 
    void GetFirstDerivativesVector(
        Vector& rValues,
        int Step) const final;
 
    void GetSecondDerivativesVector(
        Vector& rValues,
        int Step) const final;
 
 
    int Check(const ProcessInfo& rCurrentProcessInfo) const final {
        auto& r_geometry = GetGeometry();
        for (IndexType i = 0; i < r_geometry.size(); ++i) {
            KRATOS_ERROR_IF_NOT(r_geometry[i].Has(DIRECTOR))
                << "DIRECTOR not provided at node #" << r_geometry[i].Id() << std::endl;
        }
 
        return 0;
    }
 
 
    std::string Info() const final
    {
        std::stringstream buffer;
        buffer << "RMElement #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const final
    {
        rOStream << "RMElement #" << Id();
    }
 
    void PrintData(std::ostream& rOStream) const final {
        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>> reference_Curvature;
    // Components of the shear coefficient tensor on the contravariant basis
    std::vector<array_1d<double, 2>> reference_TransShear;
 
    // Determinant of the geometrical Jacobian.
    Vector m_dA_vector;
 
    Kratos::Variable<Vector>::Type const& (NodeType::*m_GetValueFunctor)(const Kratos::Variable<Vector >&) const = &NodeType::GetValue;
    Kratos::Variable<array_1d<double, 3>>::Type const& (NodeType::*m_FastGetSolutionStepValueFunctor)(const Kratos::Variable<array_1d<double,3> >&) const = &NodeType::FastGetSolutionStepValue;
    NodeType::PointType::CoordinatesArrayType const& (NodeType::* m_GetCoordinatesFunctor)() const = &NodeType::Coordinates;
    NodeType::PointType const& (NodeType::* m_GetInitialPositionFunctor)() const = &NodeType::GetInitialPosition;
 
    // Transformed curvilinear derivatives into cartesian derivatives/
    std::vector<Matrix> m_cart_deriv;
 
    //The St. Venant Kirchhoff 5P Material Tangent
    BoundedMatrix<double, 8, 8> mC;
 
 
    void CalculateAll(
        MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        const ProcessInfo& rCurrentProcessInfo,
        const bool CalculateStiffnessMatrixFlag,
        const bool CalculateResidualVectorFlag
    );
 
    std::pair< Shell5pElement::KinematicVariables, Shell5pElement::VariationVariables>
        CalculateKinematics(const IndexType IntegrationPointIndex) const;
 
    // Computes the cartesian derivatives from curvilinear ones
    Matrix CalculateCartesianDerivatives(
        const IndexType IntegrationPointIndex);
 
    Matrix CalculateStrainDisplacementOperator(
        const IndexType IntegrationPointIndex,
        const KinematicVariables& rActualKinematic,
        const VariationVariables& rVariations) const ;
 
    Matrix CalculateGeometricStiffness(
        const IndexType IntegrationPointIndex,
        const KinematicVariables& rActKin,
        const VariationVariables& ractVar,
        const ConstitutiveVariables& rThisConstitutiveVariables) const;
 
    void CalculateConstitutiveVariables(
        const IndexType IntegrationPointIndex,
        const KinematicVariables& rActualMetric,
        ConstitutiveVariables& rThisConstitutiveVariables,
        ConstitutiveLaw::Parameters& rValues,
        const ConstitutiveLaw::StressMeasure ThisStressMeasure
    );
 
    void CalculateSVKMaterialTangent();
 
    template< typename ContainerType, typename NodeFunctor, typename ...Args>
    BoundedVector<double, 3> InterpolateNodalVariable(const ContainerType& vec, const NodeFunctor& funct, const Args&... args) const;
 
    public:
    static BoundedMatrix<double, 3, 2> TangentSpaceFromStereographicProjection(const array_1d<double, 3 >& director);
    private:
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const final
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Element);
        rSerializer.save("reference_Curvature", reference_Curvature);
        rSerializer.save("reference_TransShear", reference_TransShear);
        rSerializer.save("dA_vector", m_dA_vector);
        rSerializer.save("cart_deriv", m_cart_deriv);
    }
 
    void load(Serializer& rSerializer) final
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element);
        rSerializer.load("curvature", reference_Curvature);
        rSerializer.load("reference_TransShear", reference_TransShear);
        rSerializer.load("dA_vector", m_dA_vector);
        rSerializer.save("cart_deriv", m_cart_deriv);
    }
 
 
};     // Class Shell5pElement
 
}  // namespace Kratos.
 
#endif // KRATOS_MESHLESS_SHELL_5P_ELEMENT_H_INCLUDED  defined