mortar_contact_condition.h

Go to the documentation of this file.

// KRATOS    ______            __             __  _____ __                  __                   __
//          / ____/___  ____  / /_____ ______/ /_/ ___// /________  _______/ /___  ___________ _/ /
//         / /   / __ \/ __ \/ __/ __ `/ ___/ __/\__ \/ __/ ___/ / / / ___/ __/ / / / ___/ __ `/ /
//        / /___/ /_/ / / / / /_/ /_/ / /__/ /_ ___/ / /_/ /  / /_/ / /__/ /_/ /_/ / /  / /_/ / /
//        \____/\____/_/ /_/\__/\__,_/\___/\__//____/\__/_/   \__,_/\___/\__/\__,_/_/   \__,_/_/  MECHANICS
//
//  License:         BSD License
//                   license: ContactStructuralMechanicsApplication/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
#pragma once
 
// System includes
 
// External includes
 
// Project includes
#include "contact_structural_mechanics_application_variables.h"
#include "custom_conditions/paired_condition.h"
#include "utilities/math_utils.h"
#include "includes/kratos_flags.h"
#include "includes/checks.h"
#include "includes/mortar_classes.h"
 
/* Utilities */
#include "utilities/exact_mortar_segmentation_utility.h"
#include "custom_utilities/derivatives_utilities.h"
// #include "custom_utilities/logging_settings.hpp"
 
/* Geometries */
#include "geometries/line_2d_2.h"
#include "geometries/triangle_3d_3.h"
 
namespace Kratos
{
 
 
 
    using SizeType = std::size_t;
 
 
 
 
template< const SizeType TDim, const SizeType TNumNodes, const FrictionalCase TFrictional, const bool TNormalVariation, const SizeType TNumNodesMaster = TNumNodes>
class KRATOS_API(CONTACT_STRUCTURAL_MECHANICS_APPLICATION) MortarContactCondition
    : public PairedCondition
{
public:
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( MortarContactCondition );
 
    using BaseType = PairedCondition;
 
    using VectorType = typename BaseType::VectorType;
 
    using MatrixType = typename BaseType::MatrixType;
 
    using IndexType = typename BaseType::IndexType;
 
    using GeometryPointerType = typename BaseType::GeometryType::Pointer;
 
    using NodesArrayType = typename BaseType::NodesArrayType;
 
    using PropertiesPointerType = typename BaseType::PropertiesType::Pointer;
 
    using PointType = Point;
 
    using GeometryType = Geometry<Node>;
 
    // Type definition for integration methods
    using IntegrationPointsType = typename GeometryType::IntegrationPointsArrayType;
 
    using BelongType = typename std::conditional<TNumNodes == 2, PointBelongsLine2D2N, typename std::conditional<TNumNodes == 3, typename std::conditional<TNumNodesMaster == 3, PointBelongsTriangle3D3N, PointBelongsTriangle3D3NQuadrilateral3D4N>::type, typename std::conditional<TNumNodesMaster == 3, PointBelongsQuadrilateral3D4NTriangle3D3N, PointBelongsQuadrilateral3D4N>::type>::type>::type;
 
    using PointBelongType = PointBelong<TNumNodes, TNumNodesMaster>;
 
    using GeometryPointBelongType = Geometry<PointBelongType>;
 
    using ConditionArrayType = array_1d<PointBelongType, TDim>;
 
    using ConditionArrayListType = typename std::vector<ConditionArrayType>;
 
    using LineType = Line2D2<PointType>;
 
    using TriangleType = Triangle3D3<PointType>;
 
    using DecompositionType = typename std::conditional<TDim == 2, LineType, TriangleType >::type;
 
    using DerivativeDataType = typename std::conditional<TFrictional == FrictionalCase::FRICTIONAL || TFrictional == FrictionalCase::FRICTIONAL_PENALTY, DerivativeDataFrictional<TDim, TNumNodes, TNumNodesMaster>, DerivativeData<TDim, TNumNodes, TNumNodesMaster> >::type;
 
    static constexpr IndexType MatrixSize = (TFrictional == FrictionalCase::FRICTIONLESS) ? TDim * (TNumNodesMaster + TNumNodes) + TNumNodes : (TFrictional == FrictionalCase::FRICTIONLESS_COMPONENTS || TFrictional == FrictionalCase::FRICTIONAL) ? TDim * (TNumNodesMaster + TNumNodes + TNumNodes) :  TDim * (TNumNodesMaster + TNumNodes);
 
    static constexpr bool IsFrictional  = (TFrictional == FrictionalCase::FRICTIONAL || TFrictional == FrictionalCase::FRICTIONAL_PENALTY) ? true: false;
 
    using GeneralVariables = MortarKinematicVariablesWithDerivatives<TDim, TNumNodes, TNumNodesMaster>;
 
    using AeData = DualLagrangeMultiplierOperatorsWithDerivatives<TDim, TNumNodes, IsFrictional, TNumNodesMaster>;
 
    using MortarConditionMatrices = MortarOperatorWithDerivatives<TDim, TNumNodes, IsFrictional, TNumNodesMaster>;
 
    using IntegrationUtility = ExactMortarIntegrationUtility<TDim, TNumNodes, true, TNumNodesMaster>;
 
    using DerivativesUtilitiesType = DerivativesUtilities<TDim, TNumNodes, IsFrictional, TNormalVariation, TNumNodesMaster>;
 
    // The threshold coefficient considered for checking
    static constexpr double CheckThresholdCoefficient = 1.0e-12;
 
 
    MortarContactCondition()
        : PairedCondition()
    {}
 
    // Constructor 1
    MortarContactCondition(
        IndexType NewId,
        GeometryType::Pointer pGeometry
        ) :PairedCondition(NewId, pGeometry)
    {}
 
    // Constructor 2
    MortarContactCondition(
        IndexType NewId,
        GeometryType::Pointer pGeometry,
        PropertiesType::Pointer pProperties
        ) :PairedCondition( NewId, pGeometry, pProperties )
    {}
 
    // Constructor 3
    MortarContactCondition(
        IndexType NewId,
        GeometryType::Pointer pGeometry,
        PropertiesType::Pointer pProperties,
        GeometryType::Pointer pMasterGeometry
        )
        :PairedCondition( NewId, pGeometry, pProperties, pMasterGeometry)
    {}
 
    MortarContactCondition( MortarContactCondition const& rOther){}
 
    ~MortarContactCondition() override;
 
 
 
 
    void Initialize(const ProcessInfo& rCurrentProcessInfo) override;
 
    void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override;
 
    void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateMassMatrix(
        MatrixType& rMassMatrix,
        const ProcessInfo& rCurrentProcessInfo
        ) override;
 
    void CalculateDampingMatrix(
        MatrixType& rDampingMatrix,
        const ProcessInfo& rCurrentProcessInfo
        ) override;
 
    Condition::Pointer Create(
        IndexType NewId,
        NodesArrayType const& rThisNodes,
        PropertiesType::Pointer pProperties
        ) const override;
 
    Condition::Pointer Create(
        IndexType NewId,
        GeometryType::Pointer pGeom,
        PropertiesType::Pointer pProperties
        ) const override;
 
    Condition::Pointer Create(
        IndexType NewId,
        GeometryType::Pointer pGeom,
        PropertiesType::Pointer pProperties,
        GeometryType::Pointer pMasterGeom
        ) const override;
 
    void AddExplicitContribution(const ProcessInfo& rCurrentProcessInfo) override;
 
    void AddExplicitContribution(
        const VectorType& rRHSVector,
        const Variable<VectorType>& rRHSVariable,
        const Variable<double >& rDestinationVariable,
        const ProcessInfo& rCurrentProcessInfo
        ) override;
 
    void AddExplicitContribution(const VectorType& rRHSVector,
        const Variable<VectorType>& rRHSVariable,
        const Variable<array_1d<double, 3> >& rDestinationVariable,
        const ProcessInfo& rCurrentProcessInfo
        ) override;
 
    /******************************************************************/
    /********** AUXILLIARY METHODS FOR GENERAL CALCULATIONS ***********/
    /******************************************************************/
 
    void EquationIdVector(
        EquationIdVectorType& rResult,
        const ProcessInfo& rCurrentProcessInfo
        ) const override;
 
    void GetDofList(
        DofsVectorType& rConditionalDofList,
        const ProcessInfo& rCurrentProcessInfo
        ) 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;
 
    void CalculateOnIntegrationPoints(
        const Variable<Vector>& rVariable,
        std::vector<Vector>& rOutput,
        const ProcessInfo& rCurrentProcessInfo
        ) override;
 
    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "MortarContactCondition #" << this->Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "MortarContactCondition #" << this->Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        PrintInfo(rOStream);
        this->GetParentGeometry().PrintData(rOStream);
        this->GetPairedGeometry().PrintData(rOStream);
    }
 
 
 
protected:
 
 
 
 
    /******************************************************************/
    /*********************** 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 CalculateConditionSystem(
        MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        const ProcessInfo& CurrentProcessInfo,
        const bool ComputeLHS = true,
        const bool ComputeRHS = true
        );
 
    /********************************************************************************/
    /**************** METHODS TO CALCULATE MORTAR CONDITION MATRICES ****************/
    /********************************************************************************/
 
    virtual void CalculateLocalLHS(
        Matrix& rLocalLHS,
        const MortarConditionMatrices& rMortarConditionMatrices,
        const DerivativeDataType& rDerivativeData,
        const IndexType rActiveInactive,
        const ProcessInfo& rCurrentProcessInfo
        );
 
    virtual void CalculateLocalRHS(
        Vector& rLocalRHS,
        const MortarConditionMatrices& rMortarConditionMatrices,
        const DerivativeDataType& rDerivativeData,
        const IndexType rActiveInactive,
        const ProcessInfo& rCurrentProcessInfo
        );
 
    /******************************************************************/
    /********** AUXILLIARY METHODS FOR GENERAL CALCULATIONS ***********/
    /******************************************************************/
 
    virtual IndexType GetActiveInactiveValue(const GeometryType& CurrentGeometry) const
    {
        KRATOS_ERROR << "You are calling to the base class method GetActiveInactiveValue, you are evil, and your seed must be eradicated from the face of the earth" << std::endl;
 
        return 0;
    }
 
    bool CheckIsolatedElement(
        const double DeltaTime,
        const bool HalfJump = true
        );
 
    IntegrationMethod GetIntegrationMethod() const override
    {
        // Setting the auxiliary integration points
        const IndexType integration_order = GetProperties().Has(INTEGRATION_ORDER_CONTACT) ? GetProperties().GetValue(INTEGRATION_ORDER_CONTACT) : 2;
        switch (integration_order) {
        case 1:
            return GeometryData::IntegrationMethod::GI_GAUSS_1;
        case 2:
            return GeometryData::IntegrationMethod::GI_GAUSS_2;
        case 3:
            return GeometryData::IntegrationMethod::GI_GAUSS_3;
        case 4:
            return GeometryData::IntegrationMethod::GI_GAUSS_4;
        case 5:
            return GeometryData::IntegrationMethod::GI_GAUSS_5;
        default:
            return GeometryData::IntegrationMethod::GI_GAUSS_2;
        }
    }
 
    virtual bool IsAxisymmetric() const;
 
    virtual double GetAxisymmetricCoefficient(const GeneralVariables& rVariables) const;
 
    virtual void ResizeLHS(MatrixType& rLeftHandSideMatrix);
 
    virtual void ResizeRHS(VectorType& rRightHandSideVector);
 
    virtual void ZeroLHS(MatrixType& rLeftHandSideMatrix);
 
    virtual void ZeroRHS(VectorType& rRightHandSideVector);
 
 
 
 
private:
 
 
 
 
 
 
 
    // Serialization
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS( rSerializer, PairedCondition );
    }
 
    void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, PairedCondition );
    }
 
 
}; // Class MortarContactCondition
 
 
 
 
 
}// namespace Kratos.