mortar_classes.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
#pragma once
 
// System includes
 
// External includes
 
// Project includes
#include "utilities/mortar_utilities.h"
#include "includes/variables.h"
#include "includes/mapping_variables.h"
#include "includes/serializer.h"
 
namespace Kratos
{
 
 
 
    // Point and nodes defines
    using PointType = Point;
    using GeometryType = Geometry<Node>;
 
    // Type definition for integration methods
    using IntegrationPointsType = typename GeometryType::IntegrationPointsArrayType;
 
    using SizeType = std::size_t;
 
    using IndexType = std::size_t;
 
 
 
    enum class PointBelongs
    {
        Master = 0,      
        Slave = 1,       
        Intersection = 2 
    };
 
    enum class PointBelongsLine2D2N
    {
        SlaveLine2D2N0 = 0,       
        SlaveLine2D2N1 = 1,       
        MasterLine2D2N0 = 2,      
        MasterLine2D2N1 = 3,      
        IntersectionLine2D2N = 4  
    };
 
    enum class PointBelongsTriangle3D3N
    {
        SlaveTriangle3D3N0           = 0,
        SlaveTriangle3D3N1           = 1,
        SlaveTriangle3D3N2           = 2,
        MasterTriangle3D3N0          = 3,
        MasterTriangle3D3N1          = 4,
        MasterTriangle3D3N2          = 5,
        IntersectionTriangle3D3N     = 6,
        IntersectionTriangle3D3N0101 = 10106,
        IntersectionTriangle3D3N1001 = 10016,
        IntersectionTriangle3D3N1201 = 10216,
        IntersectionTriangle3D3N2101 = 10126,
        IntersectionTriangle3D3N0201 = 10206,
        IntersectionTriangle3D3N2001 = 10026,
        IntersectionTriangle3D3N0110 = 1106,
        IntersectionTriangle3D3N1010 = 1016,
        IntersectionTriangle3D3N1210 = 1216,
        IntersectionTriangle3D3N2110 = 1126,
        IntersectionTriangle3D3N0210 = 1206,
        IntersectionTriangle3D3N2010 = 1026,
        IntersectionTriangle3D3N0112 = 21106,
        IntersectionTriangle3D3N1012 = 21016,
        IntersectionTriangle3D3N1212 = 21216,
        IntersectionTriangle3D3N2112 = 21126,
        IntersectionTriangle3D3N0212 = 21206,
        IntersectionTriangle3D3N2012 = 21026,
        IntersectionTriangle3D3N0121 = 12106,
        IntersectionTriangle3D3N1021 = 12016,
        IntersectionTriangle3D3N1221 = 12216,
        IntersectionTriangle3D3N2121 = 12126,
        IntersectionTriangle3D3N0221 = 12206,
        IntersectionTriangle3D3N2021 = 12026,
        IntersectionTriangle3D3N0102 = 20106,
        IntersectionTriangle3D3N1002 = 20016,
        IntersectionTriangle3D3N1202 = 20216,
        IntersectionTriangle3D3N2102 = 20126,
        IntersectionTriangle3D3N0202 = 20206,
        IntersectionTriangle3D3N2002 = 20026,
        IntersectionTriangle3D3N0120 = 2106,
        IntersectionTriangle3D3N1020 = 2016,
        IntersectionTriangle3D3N1220 = 2216,
        IntersectionTriangle3D3N2120 = 2126,
        IntersectionTriangle3D3N0220 = 2206,
        IntersectionTriangle3D3N2020 = 2026
    };
 
    enum class PointBelongsQuadrilateral3D4N
    {
        SlaveQuadrilateral3D4N0           = 0,
        SlaveQuadrilateral3D4N1           = 1,
        SlaveQuadrilateral3D4N2           = 2,
        SlaveQuadrilateral3D4N3           = 3,
        MasterQuadrilateral3D4N0          = 4,
        MasterQuadrilateral3D4N1          = 5,
        MasterQuadrilateral3D4N2          = 6,
        MasterQuadrilateral3D4N3          = 7,
        IntersectionQuadrilateral3D4N     = 8,
        IntersectionQuadrilateral3D4N0101 = 10108,
        IntersectionQuadrilateral3D4N1001 = 10018,
        IntersectionQuadrilateral3D4N1201 = 10218,
        IntersectionQuadrilateral3D4N2101 = 10128,
        IntersectionQuadrilateral3D4N2301 = 10328,
        IntersectionQuadrilateral3D4N3201 = 10238,
        IntersectionQuadrilateral3D4N3001 = 10038,
        IntersectionQuadrilateral3D4N0301 = 10308,
        IntersectionQuadrilateral3D4N0110 = 1108,
        IntersectionQuadrilateral3D4N1010 = 1018,
        IntersectionQuadrilateral3D4N1210 = 1218,
        IntersectionQuadrilateral3D4N2110 = 1128,
        IntersectionQuadrilateral3D4N2310 = 1328,
        IntersectionQuadrilateral3D4N3210 = 1238,
        IntersectionQuadrilateral3D4N3010 = 1038,
        IntersectionQuadrilateral3D4N0310 = 1308,
        IntersectionQuadrilateral3D4N0112 = 21108,
        IntersectionQuadrilateral3D4N1012 = 21018,
        IntersectionQuadrilateral3D4N1212 = 21218,
        IntersectionQuadrilateral3D4N2112 = 21128,
        IntersectionQuadrilateral3D4N2312 = 21328,
        IntersectionQuadrilateral3D4N3212 = 21238,
        IntersectionQuadrilateral3D4N3012 = 21038,
        IntersectionQuadrilateral3D4N0312 = 21308,
        IntersectionQuadrilateral3D4N0121 = 12108,
        IntersectionQuadrilateral3D4N1021 = 12018,
        IntersectionQuadrilateral3D4N1221 = 12218,
        IntersectionQuadrilateral3D4N2121 = 12128,
        IntersectionQuadrilateral3D4N2321 = 12328,
        IntersectionQuadrilateral3D4N3221 = 12238,
        IntersectionQuadrilateral3D4N3021 = 12038,
        IntersectionQuadrilateral3D4N0321 = 12308,
        IntersectionQuadrilateral3D4N0123 = 32108,
        IntersectionQuadrilateral3D4N1023 = 32018,
        IntersectionQuadrilateral3D4N1223 = 32218,
        IntersectionQuadrilateral3D4N2123 = 32128,
        IntersectionQuadrilateral3D4N2323 = 32328,
        IntersectionQuadrilateral3D4N3223 = 32238,
        IntersectionQuadrilateral3D4N3023 = 32038,
        IntersectionQuadrilateral3D4N0323 = 32308,
        IntersectionQuadrilateral3D4N0132 = 23108,
        IntersectionQuadrilateral3D4N1032 = 23018,
        IntersectionQuadrilateral3D4N1232 = 23218,
        IntersectionQuadrilateral3D4N2132 = 23128,
        IntersectionQuadrilateral3D4N2332 = 23328,
        IntersectionQuadrilateral3D4N3232 = 23238,
        IntersectionQuadrilateral3D4N3032 = 23038,
        IntersectionQuadrilateral3D4N0332 = 23308,
        IntersectionQuadrilateral3D4N0130 = 3108,
        IntersectionQuadrilateral3D4N1030 = 3018,
        IntersectionQuadrilateral3D4N1230 = 3218,
        IntersectionQuadrilateral3D4N2130 = 3128,
        IntersectionQuadrilateral3D4N2330 = 3328,
        IntersectionQuadrilateral3D4N3230 = 3238,
        IntersectionQuadrilateral3D4N3030 = 3038,
        IntersectionQuadrilateral3D4N0330 = 3308,
        IntersectionQuadrilateral3D4N0103 = 30108,
        IntersectionQuadrilateral3D4N1003 = 30018,
        IntersectionQuadrilateral3D4N1203 = 30218,
        IntersectionQuadrilateral3D4N2103 = 30128,
        IntersectionQuadrilateral3D4N2303 = 30328,
        IntersectionQuadrilateral3D4N3203 = 30238,
        IntersectionQuadrilateral3D4N3003 = 30038,
        IntersectionQuadrilateral3D4N0303 = 30308
    };
 
    enum class PointBelongsQuadrilateral3D4NTriangle3D3N
    {
        SlaveQuadrilateral3D4N0                       = 0,
        SlaveQuadrilateral3D4N1                       = 1,
        SlaveQuadrilateral3D4N2                       = 2,
        SlaveQuadrilateral3D4N3                       = 3,
        MasterTriangle3D3N0                           = 4,
        MasterTriangle3D3N1                           = 5,
        MasterTriangle3D3N2                           = 6,
        IntersectionQuadrilateral3D4NTriangle3D3N     = 7,
        IntersectionQuadrilateral3D4NTriangle3D3N0101 = 10107,
        IntersectionQuadrilateral3D4NTriangle3D3N1001 = 1107,
        IntersectionQuadrilateral3D4NTriangle3D3N1201 = 21107,
        IntersectionQuadrilateral3D4NTriangle3D3N2101 = 12107,
        IntersectionQuadrilateral3D4NTriangle3D3N2301 = 32107,
        IntersectionQuadrilateral3D4NTriangle3D3N3201 = 23107,
        IntersectionQuadrilateral3D4NTriangle3D3N3001 = 3107,
        IntersectionQuadrilateral3D4NTriangle3D3N0301 = 30107,
        IntersectionQuadrilateral3D4NTriangle3D3N0110 = 10017,
        IntersectionQuadrilateral3D4NTriangle3D3N1010 = 1017,
        IntersectionQuadrilateral3D4NTriangle3D3N1210 = 21017,
        IntersectionQuadrilateral3D4NTriangle3D3N2110 = 12017,
        IntersectionQuadrilateral3D4NTriangle3D3N2310 = 32017,
        IntersectionQuadrilateral3D4NTriangle3D3N3210 = 23017,
        IntersectionQuadrilateral3D4NTriangle3D3N3010 = 3017,
        IntersectionQuadrilateral3D4NTriangle3D3N0310 = 30017,
        IntersectionQuadrilateral3D4NTriangle3D3N0112 = 10217,
        IntersectionQuadrilateral3D4NTriangle3D3N1012 = 1217,
        IntersectionQuadrilateral3D4NTriangle3D3N1212 = 21217,
        IntersectionQuadrilateral3D4NTriangle3D3N2112 = 12217,
        IntersectionQuadrilateral3D4NTriangle3D3N2312 = 32217,
        IntersectionQuadrilateral3D4NTriangle3D3N3212 = 23217,
        IntersectionQuadrilateral3D4NTriangle3D3N3012 = 3217,
        IntersectionQuadrilateral3D4NTriangle3D3N0312 = 30217,
        IntersectionQuadrilateral3D4NTriangle3D3N0121 = 10127,
        IntersectionQuadrilateral3D4NTriangle3D3N1021 = 1127,
        IntersectionQuadrilateral3D4NTriangle3D3N1221 = 21127,
        IntersectionQuadrilateral3D4NTriangle3D3N2121 = 12127,
        IntersectionQuadrilateral3D4NTriangle3D3N2321 = 32127,
        IntersectionQuadrilateral3D4NTriangle3D3N3221 = 23127,
        IntersectionQuadrilateral3D4NTriangle3D3N3021 = 3127,
        IntersectionQuadrilateral3D4NTriangle3D3N0321 = 30127,
        IntersectionQuadrilateral3D4NTriangle3D3N0120 = 10027,
        IntersectionQuadrilateral3D4NTriangle3D3N1020 = 1027,
        IntersectionQuadrilateral3D4NTriangle3D3N1220 = 21027,
        IntersectionQuadrilateral3D4NTriangle3D3N2120 = 12027,
        IntersectionQuadrilateral3D4NTriangle3D3N2320 = 32027,
        IntersectionQuadrilateral3D4NTriangle3D3N3220 = 23027,
        IntersectionQuadrilateral3D4NTriangle3D3N3020 = 3027,
        IntersectionQuadrilateral3D4NTriangle3D3N0320 = 30027,
        IntersectionQuadrilateral3D4NTriangle3D3N0102 = 10207,
        IntersectionQuadrilateral3D4NTriangle3D3N1002 = 1207,
        IntersectionQuadrilateral3D4NTriangle3D3N1202 = 21207,
        IntersectionQuadrilateral3D4NTriangle3D3N2102 = 12207,
        IntersectionQuadrilateral3D4NTriangle3D3N2302 = 32207,
        IntersectionQuadrilateral3D4NTriangle3D3N3202 = 23207,
        IntersectionQuadrilateral3D4NTriangle3D3N3002 = 3207,
        IntersectionQuadrilateral3D4NTriangle3D3N0302 = 30207
    };
 
    enum class PointBelongsTriangle3D3NQuadrilateral3D4N
    {
        SlaveTriangle3D3N0                            = 0,
        SlaveTriangle3D3N1                            = 1,
        SlaveTriangle3D3N2                            = 2,
        MasterQuadrilateral3D4N0                      = 3,
        MasterQuadrilateral3D4N1                      = 4,
        MasterQuadrilateral3D4N2                      = 5,
        MasterQuadrilateral3D4N3                      = 6,
        IntersectionTriangle3D3NQuadrilateral3D4N     = 7,
        IntersectionTriangle3D3NQuadrilateral3D4N0101 = 10107,
        IntersectionTriangle3D3NQuadrilateral3D4N0110 = 1107,
        IntersectionTriangle3D3NQuadrilateral3D4N0112 = 21107,
        IntersectionTriangle3D3NQuadrilateral3D4N0121 = 12107,
        IntersectionTriangle3D3NQuadrilateral3D4N0123 = 32107,
        IntersectionTriangle3D3NQuadrilateral3D4N0132 = 23107,
        IntersectionTriangle3D3NQuadrilateral3D4N0130 = 3107,
        IntersectionTriangle3D3NQuadrilateral3D4N0103 = 30107,
        IntersectionTriangle3D3NQuadrilateral3D4N1001 = 10017,
        IntersectionTriangle3D3NQuadrilateral3D4N1010 = 1017,
        IntersectionTriangle3D3NQuadrilateral3D4N1012 = 21017,
        IntersectionTriangle3D3NQuadrilateral3D4N1021 = 12017,
        IntersectionTriangle3D3NQuadrilateral3D4N1023 = 32017,
        IntersectionTriangle3D3NQuadrilateral3D4N1032 = 23017,
        IntersectionTriangle3D3NQuadrilateral3D4N1030 = 3017,
        IntersectionTriangle3D3NQuadrilateral3D4N1003 = 30017,
        IntersectionTriangle3D3NQuadrilateral3D4N1201 = 10217,
        IntersectionTriangle3D3NQuadrilateral3D4N1210 = 1217,
        IntersectionTriangle3D3NQuadrilateral3D4N1212 = 21217,
        IntersectionTriangle3D3NQuadrilateral3D4N1221 = 12217,
        IntersectionTriangle3D3NQuadrilateral3D4N1223 = 32217,
        IntersectionTriangle3D3NQuadrilateral3D4N1232 = 23217,
        IntersectionTriangle3D3NQuadrilateral3D4N1230 = 3217,
        IntersectionTriangle3D3NQuadrilateral3D4N1203 = 30217,
        IntersectionTriangle3D3NQuadrilateral3D4N2101 = 10127,
        IntersectionTriangle3D3NQuadrilateral3D4N2110 = 1127,
        IntersectionTriangle3D3NQuadrilateral3D4N2112 = 21127,
        IntersectionTriangle3D3NQuadrilateral3D4N2121 = 12127,
        IntersectionTriangle3D3NQuadrilateral3D4N2123 = 32127,
        IntersectionTriangle3D3NQuadrilateral3D4N2132 = 23127,
        IntersectionTriangle3D3NQuadrilateral3D4N2130 = 3127,
        IntersectionTriangle3D3NQuadrilateral3D4N2103 = 30127,
        IntersectionTriangle3D3NQuadrilateral3D4N2001 = 10027,
        IntersectionTriangle3D3NQuadrilateral3D4N2010 = 1027,
        IntersectionTriangle3D3NQuadrilateral3D4N2012 = 21027,
        IntersectionTriangle3D3NQuadrilateral3D4N2021 = 12027,
        IntersectionTriangle3D3NQuadrilateral3D4N2023 = 32027,
        IntersectionTriangle3D3NQuadrilateral3D4N2032 = 23027,
        IntersectionTriangle3D3NQuadrilateral3D4N2030 = 3027,
        IntersectionTriangle3D3NQuadrilateral3D4N2003 = 30027,
        IntersectionTriangle3D3NQuadrilateral3D4N0201 = 10207,
        IntersectionTriangle3D3NQuadrilateral3D4N0210 = 1207,
        IntersectionTriangle3D3NQuadrilateral3D4N0212 = 21207,
        IntersectionTriangle3D3NQuadrilateral3D4N0221 = 12207,
        IntersectionTriangle3D3NQuadrilateral3D4N0223 = 32207,
        IntersectionTriangle3D3NQuadrilateral3D4N0232 = 23207,
        IntersectionTriangle3D3NQuadrilateral3D4N0230 = 3207,
        IntersectionTriangle3D3NQuadrilateral3D4N0203 = 30207
    };
 
 
 
template< const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class MortarKinematicVariables
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( MortarKinematicVariables );
 
 
    MortarKinematicVariables(){}
 
    virtual ~MortarKinematicVariables(){}
 
    // Shape functions for contact pair
    Vector NMaster = Vector(TNumNodesMaster, 0.0);
    Vector NSlave = Vector(TNumNodes, 0.0);
    Vector PhiLagrangeMultipliers = Vector(TNumNodes, 0.0);
 
    // Determinant of slave cell's jacobian
    double DetjSlave = 0.0;
 
 
 
    void Initialize()
    {
        // Shape functions
        noalias(NMaster)                = ZeroVector(TNumNodesMaster);
        noalias(NSlave)                 = ZeroVector(TNumNodes);
        noalias(PhiLagrangeMultipliers) = ZeroVector(TNumNodes);
 
        // Jacobian of slave
        DetjSlave = 0.0;
    }
 
    void PrintInfo(std::ostream& rOStream) const
    {
        rOStream << "NSlave:" << NSlave << std::endl;
        rOStream << "NMaster: " <<NMaster << std::endl;
        rOStream << "PhiLagrangeMultipliers: "<< PhiLagrangeMultipliers << std::endl;
        rOStream << "DetjSlave: " << DetjSlave << std::endl;
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    virtual void save(Serializer& rSerializer) const
    {
        rSerializer.save("NMaster", NMaster);
        rSerializer.save("NSlave", NSlave);
        rSerializer.save("PhiLagrangeMultipliers", PhiLagrangeMultipliers);
        rSerializer.save("DetjSlave", DetjSlave);
    }
 
    virtual void load(Serializer& rSerializer)
    {
        rSerializer.load("NMaster", NMaster);
        rSerializer.load("NSlave", NSlave);
        rSerializer.load("PhiLagrangeMultipliers", PhiLagrangeMultipliers);
        rSerializer.load("DetjSlave", DetjSlave);
    }
 
 
 
 
 
}; // Class MortarKinematicVariables
 
template< const SizeType TDim, const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class MortarKinematicVariablesWithDerivatives
    : public MortarKinematicVariables<TNumNodes, TNumNodesMaster>
{
public:
 
    using BaseClassType = MortarKinematicVariables<TNumNodes, TNumNodesMaster>;
 
    KRATOS_CLASS_POINTER_DEFINITION( MortarKinematicVariablesWithDerivatives );
 
 
    MortarKinematicVariablesWithDerivatives()= default;
 
    ~MortarKinematicVariablesWithDerivatives() override= default;
 
    // Shape functions local derivatives for contact pair
    Matrix DNDeMaster = ScalarMatrix(TNumNodesMaster, TDim - 1, 0.0);
    Matrix DNDeSlave = ScalarMatrix(TNumNodes, TDim - 1, 0.0);
 
    /*
    * Jacobians in current configuration on all integration points of slave segment
    * Only those two variables contain info on all GP
    * other variables contain info only on the currently-calculated GP
    */
    Matrix jSlave = ScalarMatrix(TDim, TDim - 1, 0.0);
    Matrix jMaster = ScalarMatrix(TDim, TDim - 1, 0.0);
 
 
 
    void Initialize()
    {
        BaseClassType::Initialize();
 
        // Shape functions local derivatives
        noalias(DNDeMaster) = ZeroMatrix(TNumNodesMaster, TDim - 1);
        noalias(DNDeSlave)  = ZeroMatrix(TNumNodes, TDim - 1);
 
        // Jacobians on all integration points
        noalias(jSlave)  = ZeroMatrix(TDim, TDim - 1);
        noalias(jMaster) = ZeroMatrix(TDim, TDim - 1);
    }
 
    void PrintInfo(std::ostream& rOStream) const
    {
        BaseClassType::PrintInfo(rOStream);
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.save("DNDeMaster", DNDeMaster);
        rSerializer.save("DNDeSlave", DNDeSlave);
        rSerializer.save("jSlave", jSlave);
        rSerializer.save("jMaster", jMaster);
    }
 
    void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.load("DNDeMaster", DNDeMaster);
        rSerializer.load("DNDeSlave", DNDeSlave);
        rSerializer.load("jSlave", jSlave);
        rSerializer.load("jMaster", jMaster);
    }
 
 
 
 
 
}; // Class MortarKinematicVariablesWithDerivatives
 
template< const SizeType TDim, const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class DerivativeData
{
public:
 
    // Auxiliary types
    using DummyBoundedMatrixType = BoundedMatrix<int, 1, 1>;
 
    using GeometryArraySlaveType = array_1d<double, TNumNodes>;
 
    using GeometryArrayMasterType = array_1d<double, TNumNodesMaster>;
 
    using GeometryDoFMatrixSlaveType = BoundedMatrix<double, TNumNodes, TDim>;
 
    using GeometryDoFMatrixMasterType = BoundedMatrix<double, TNumNodesMaster, TDim>;
 
    using GeometryMatrixType = BoundedMatrix<double, TNumNodes, TNumNodes>;
 
    using VertexDerivativesMatrixType = typename std::conditional<TNumNodes == 2, DummyBoundedMatrixType, BoundedMatrix<double, 3, 3>>::type;
 
    // Auxiliary sizes
    static const SizeType DummySize = 1;
 
    static const SizeType DoFSizeSlaveGeometry = (TNumNodes * TDim);
 
    static const SizeType DoFSizeMasterGeometry = (TNumNodesMaster * TDim);
 
    static const SizeType DoFSizePairedGeometry = DoFSizeSlaveGeometry + DoFSizeMasterGeometry;
 
    static const SizeType DoFSizeDerivativesDependence = (TDim == 2) ? DoFSizeSlaveGeometry : DoFSizePairedGeometry;
 
    static const SizeType DoFSizeDerivativesVertex = (TDim == 2) ? DummySize : DoFSizePairedGeometry;
 
 
    DerivativeData()= default;
 
    virtual ~DerivativeData()= default;
 
    array_1d<double, TNumNodes> PenaltyParameter;
    double ScaleFactor;
 
    GeometryDoFMatrixSlaveType NormalSlave;
    GeometryDoFMatrixMasterType NormalMaster;
 
    GeometryDoFMatrixSlaveType X1, u1;
    GeometryDoFMatrixMasterType X2, u2;
 
    array_1d<double, DoFSizeDerivativesDependence> DeltaDetjSlave;
    array_1d<GeometryArraySlaveType, DoFSizeDerivativesDependence> DeltaPhi;
    array_1d<GeometryArraySlaveType, DoFSizeDerivativesDependence> DeltaN1;
    array_1d<GeometryArrayMasterType, DoFSizeDerivativesDependence> DeltaN2;
    array_1d<GeometryDoFMatrixSlaveType, DoFSizeSlaveGeometry> DeltaNormalSlave;
    array_1d<GeometryDoFMatrixMasterType, DoFSizeMasterGeometry> DeltaNormalMaster;
    array_1d<VertexDerivativesMatrixType, DoFSizeDerivativesVertex> DeltaCellVertex;
 
    GeometryMatrixType Ae;
 
    array_1d<GeometryMatrixType, DoFSizeDerivativesDependence> DeltaAe;
 
 
 
    virtual void Initialize(
        const GeometryType& SlaveGeometry,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        // The normals of the nodes
        noalias(NormalSlave) = MortarUtilities::GetVariableMatrix<TDim,TNumNodes>(SlaveGeometry,  NORMAL, 0);
 
        // Displacements and velocities of the slave
        const IndexType step = (rCurrentProcessInfo[STEP] == 1) ? 0 : 1;
        noalias(u1) = step == 0 ?
                      MortarUtilities::GetVariableMatrix<TDim,TNumNodes>(SlaveGeometry, DISPLACEMENT, 0) :
                      MortarUtilities::GetVariableMatrix<TDim,TNumNodes>(SlaveGeometry, DISPLACEMENT, 0)
                    - MortarUtilities::GetVariableMatrix<TDim,TNumNodes>(SlaveGeometry, DISPLACEMENT, 1);
        noalias(X1) = MortarUtilities::GetCoordinates<TDim,TNumNodes>(SlaveGeometry, false, step);
 
        // We get the ALM variables
        for (IndexType i = 0; i < TNumNodes; ++i)
            PenaltyParameter[i] = SlaveGeometry[i].GetValue(INITIAL_PENALTY);
        ScaleFactor = rCurrentProcessInfo[SCALE_FACTOR];
 
        // We initialize the derivatives
        const array_1d<double, TNumNodes> zero_vector_slave(TNumNodes, 0.0);
        const array_1d<double, TNumNodesMaster> zero_vector_master(TNumNodesMaster, 0.0);
        for (IndexType i = 0; i < DoFSizeSlaveGeometry; ++i) {
            DeltaDetjSlave[i] = 0.0;
            noalias(DeltaPhi[i]) = zero_vector_slave;
            noalias(DeltaN1[i]) = zero_vector_slave;
            noalias(DeltaN2[i]) = zero_vector_master;
        }
        if constexpr (TDim == 3) {
            for (IndexType i = 0; i < DoFSizeMasterGeometry; ++i) {
                DeltaDetjSlave[i + DoFSizeSlaveGeometry] = 0.0;
                noalias(DeltaPhi[i + DoFSizeSlaveGeometry]) = zero_vector_slave;
                noalias(DeltaN1[i + DoFSizeSlaveGeometry]) = zero_vector_slave;
                noalias(DeltaN2[i + DoFSizeSlaveGeometry]) = zero_vector_master;
            }
        }
    }
 
    virtual void ResetDerivatives()
    {
        // Derivatives
        if constexpr (TDim == 3) { // Derivative of the cell vertex
            // Auxiliary zero matrix
            const BoundedMatrix<double, 3, 3> aux_zero = ZeroMatrix(3, 3);
            for (IndexType i = 0; i < TNumNodes * TDim; ++i) {
                noalias(DeltaCellVertex[i]) = aux_zero;
            }
            for (IndexType i = 0; i < TNumNodesMaster * TDim; ++i) {
                noalias(DeltaCellVertex[i + TNumNodes * TDim]) = aux_zero;
            }
        }
    }
 
    void InitializeDeltaAeComponents()
    {
        // Auxiliary zero matrix
        const GeometryMatrixType aux_zero = ZeroMatrix(TNumNodes, TNumNodes);
 
        // Ae
        noalias(Ae) = aux_zero;
 
        // Derivatives Ae
        for (IndexType i = 0; i < DoFSizeDerivativesDependence; ++i)
            noalias(DeltaAe[i]) = aux_zero;
    }
 
    virtual void UpdateMasterPair(
        const GeometryType& MasterGeometry,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        NormalMaster = MortarUtilities::GetVariableMatrix<TDim,TNumNodesMaster>(MasterGeometry,  NORMAL, 0);
 
        // Displacements, coordinates and normals of the master
        const IndexType step = (rCurrentProcessInfo[STEP] == 1) ? 0 : 1;
        noalias(u2) = step == 0 ?
                      MortarUtilities::GetVariableMatrix<TDim,TNumNodesMaster>(MasterGeometry, DISPLACEMENT, 0) :
                      MortarUtilities::GetVariableMatrix<TDim,TNumNodesMaster>(MasterGeometry, DISPLACEMENT, 0)
                    - MortarUtilities::GetVariableMatrix<TDim,TNumNodesMaster>(MasterGeometry, DISPLACEMENT, 1);
        noalias(X2) = MortarUtilities::GetCoordinates<TDim,TNumNodesMaster>(MasterGeometry, false, step);
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    virtual void save(Serializer& rSerializer) const
    {
        rSerializer.save("PenaltyParameter", PenaltyParameter);
        rSerializer.save("ScaleFactor", ScaleFactor);
        rSerializer.save("NormalSlave", NormalSlave);
        rSerializer.save("NormalMaster", NormalMaster);
        rSerializer.save("X1", X1);
        rSerializer.save("u1", u1);
        rSerializer.save("X2", X2);
        rSerializer.save("u2", u2);
        rSerializer.save("DeltaDetjSlave", DeltaDetjSlave);
        rSerializer.save("DeltaPhi", DeltaPhi);
        rSerializer.save("DeltaN1", DeltaN1);
        rSerializer.save("DeltaN2", DeltaN2);
        rSerializer.save("DeltaNormalSlave", DeltaNormalSlave);
        rSerializer.save("DeltaNormalMaster", DeltaNormalMaster);
        rSerializer.save("DeltaCellVertex", DeltaCellVertex);
        rSerializer.save("Ae", Ae);
        rSerializer.save("DeltaAe", DeltaAe);
    }
 
    virtual void load(Serializer& rSerializer)
    {
        rSerializer.load("PenaltyParameter", PenaltyParameter);
        rSerializer.load("ScaleFactor", ScaleFactor);
        rSerializer.load("NormalSlave", NormalSlave);
        rSerializer.load("NormalMaster", NormalMaster);
        rSerializer.load("X1", X1);
        rSerializer.load("u1", u1);
        rSerializer.load("X2", X2);
        rSerializer.load("u2", u2);
        rSerializer.load("DeltaDetjSlave", DeltaDetjSlave);
        rSerializer.load("DeltaPhi", DeltaPhi);
        rSerializer.load("DeltaN1", DeltaN1);
        rSerializer.load("DeltaN2", DeltaN2);
        rSerializer.load("DeltaNormalSlave", DeltaNormalSlave);
        rSerializer.load("DeltaNormalMaster", DeltaNormalMaster);
        rSerializer.load("DeltaCellVertex", DeltaCellVertex);
        rSerializer.load("Ae", Ae);
        rSerializer.load("DeltaAe", DeltaAe);
    }
 
 
 
 
 
};  // Class DerivativeData
 
template< const SizeType TDim, const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class DerivativeDataFrictional
    : public DerivativeData<TDim, TNumNodes, TNumNodesMaster>
{
public:
 
    using BaseClassType = DerivativeData<TDim, TNumNodes, TNumNodesMaster>;
 
    using GeometryDoFMatrixSlaveType = BoundedMatrix<double, TNumNodes, TDim>;
    using GeometryDoFMatrixMasterType = BoundedMatrix<double, TNumNodesMaster, TDim>;
 
    // Size of DoFs of a not paired dependency
    static const SizeType DoFSizeSlaveGeometry = (TNumNodes * TDim);
    static const SizeType DoFSizeMasterGeometry = (TNumNodes * TDim);
 
    static const SizeType DoFSizePairedGeometry = DoFSizeSlaveGeometry + DoFSizeMasterGeometry;
 
 
    DerivativeDataFrictional()= default;
 
    virtual ~DerivativeDataFrictional()= default;
 
    double TangentFactor = 0.0;
 
    GeometryDoFMatrixSlaveType u1old;
    GeometryDoFMatrixMasterType u2old;
 
 
 
    void Initialize(
        const GeometryType& SlaveGeometry,
        const ProcessInfo& rCurrentProcessInfo
        ) override
    {
        BaseClassType::Initialize(SlaveGeometry, rCurrentProcessInfo);
 
        TangentFactor = rCurrentProcessInfo[TANGENT_FACTOR];
 
        noalias(u1old) = MortarUtilities::GetVariableMatrix<TDim,TNumNodes>(SlaveGeometry, DISPLACEMENT, 1) - MortarUtilities::GetVariableMatrix<TDim,TNumNodes>(SlaveGeometry, DISPLACEMENT, 2);
    }
 
    void UpdateMasterPair(
        const GeometryType& MasterGeometry,
        const ProcessInfo& rCurrentProcessInfo
        ) override
    {
        BaseClassType::UpdateMasterPair(MasterGeometry, rCurrentProcessInfo);
 
        noalias(u2old) = MortarUtilities::GetVariableMatrix<TDim,TNumNodesMaster>(MasterGeometry, DISPLACEMENT, 1) - MortarUtilities::GetVariableMatrix<TDim,TNumNodesMaster>(MasterGeometry, DISPLACEMENT, 2);
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.save("TangentFactor", TangentFactor);
        rSerializer.save("u1old", u1old);
        rSerializer.save("u2old", u2old);
    }
 
    void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.load("TangentFactor", TangentFactor);
        rSerializer.load("u1old", u1old);
        rSerializer.load("u2old", u2old);
    }
 
 
 
 
 
};  // Class DerivativeDataFrictional
 
template< const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class MortarOperator
{
public:
 
    using KinematicVariables = MortarKinematicVariables<TNumNodes, TNumNodesMaster>;
 
    using GeometryMatrixSlaveType = BoundedMatrix<double, TNumNodes, TNumNodes>;
    using GeometryMatrixMasterType = BoundedMatrix<double, TNumNodes, TNumNodesMaster>;
 
    KRATOS_CLASS_POINTER_DEFINITION( MortarOperator );
 
 
    MortarOperator(){}
 
    virtual ~MortarOperator(){}
 
    GeometryMatrixSlaveType DOperator;
    GeometryMatrixMasterType MOperator;
 
 
 
    void Initialize()
    {
        // We initialize the D and M operators
        noalias(DOperator) = ZeroMatrix(TNumNodes, TNumNodes);
        noalias(MOperator) = ZeroMatrix(TNumNodes, TNumNodesMaster);
    }
 
    void CalculateMortarOperators(
        KinematicVariables& rKinematicVariables,
        const double rIntegrationWeight
        )
    {
        /* DEFINITIONS */
        const double det_j_slave = rKinematicVariables.DetjSlave;
        const Vector& phi_vector  = rKinematicVariables.PhiLagrangeMultipliers;
        const Vector& n1_vector   = rKinematicVariables.NSlave;
        const Vector& n2_vector   = rKinematicVariables.NMaster;
 
        for (IndexType i_slave = 0; i_slave < TNumNodes; ++i_slave) {
            const double phi = phi_vector[i_slave];
            for (IndexType j_slave = 0; j_slave < TNumNodes; ++j_slave) {
                DOperator(i_slave, j_slave) += det_j_slave * rIntegrationWeight * phi * n1_vector[j_slave];
            }
            for (IndexType j_slave = 0; j_slave < TNumNodesMaster; ++j_slave) {
                MOperator(i_slave, j_slave) += det_j_slave * rIntegrationWeight * phi * n2_vector[j_slave];
            }
        }
    }
 
    GeometryMatrixMasterType ComputePOperator()
    {
        // We calculate the inverse of D operator
        double auxdet;
        GeometryMatrixSlaveType inv_D_operator;
        MathUtils<double>::InvertMatrix(DOperator, inv_D_operator, auxdet);
 
        // We calculate the P operator
        const GeometryMatrixMasterType POperator = prod(inv_D_operator, MOperator);
 
        return POperator;
    }
 
    void PrintInfo(std::ostream& rOStream) const
    {
        rOStream << "DOperator: " << DOperator << std::endl;
        rOStream << "MOperator: " << MOperator << std::endl;
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    virtual void save(Serializer& rSerializer) const
    {
        rSerializer.save("DOperator", DOperator);
        rSerializer.save("MOperator", MOperator);
    }
 
    virtual void load(Serializer& rSerializer)
    {
        rSerializer.load("DOperator", DOperator);
        rSerializer.load("MOperator", MOperator);
    }
 
 
 
 
 
}; // Class MortarOperator
 
template< const SizeType TDim, const SizeType TNumNodes, bool TFrictional, const SizeType TNumNodesMaster = TNumNodes>
class MortarOperatorWithDerivatives
    : public MortarOperator<TNumNodes, TNumNodesMaster>
{
public:
 
    using BaseClassType = MortarOperator<TNumNodes, TNumNodesMaster>;
 
    using KinematicVariables = MortarKinematicVariables<TNumNodes, TNumNodesMaster>;
 
    using DerivativeDataFrictionalType = DerivativeDataFrictional<TDim, TNumNodes, TNumNodesMaster>;
 
    using DerivativeFrictionalessDataType = DerivativeData<TDim, TNumNodes, TNumNodesMaster>;
 
    using DerivativeDataType = typename std::conditional<TFrictional, DerivativeDataFrictionalType, DerivativeFrictionalessDataType>::type;
 
    using GeometryMatrixSlaveType = BoundedMatrix<double, TNumNodes, TNumNodes>;
 
    using GeometryMatrixMasterType = BoundedMatrix<double, TNumNodes, TNumNodesMaster>;
 
    // Auxiliary sizes
    static const SizeType DoFSizeSlaveGeometry = (TNumNodes * TDim);
    static const SizeType DoFSizeMasterGeometry = (TNumNodesMaster * TDim);
 
    static const SizeType DoFSizePairedGeometry = DoFSizeSlaveGeometry + DoFSizeMasterGeometry;
 
    static const SizeType DoFSizeDerivativesDependence = (TDim == 2) ? DoFSizeSlaveGeometry : DoFSizePairedGeometry;
 
    KRATOS_CLASS_POINTER_DEFINITION( MortarOperatorWithDerivatives );
 
 
    MortarOperatorWithDerivatives(){}
 
    ~MortarOperatorWithDerivatives() override{}
 
    // D and M directional derivatives
    array_1d<GeometryMatrixSlaveType, DoFSizePairedGeometry> DeltaDOperator;
    array_1d<GeometryMatrixMasterType, DoFSizePairedGeometry> DeltaMOperator;
 
 
 
    void Initialize()
    {
        BaseClassType::Initialize();
 
        // Auxiliary zero matrix
        const GeometryMatrixSlaveType aux_zero_slave = ZeroMatrix(TNumNodes, TNumNodes);
        const GeometryMatrixMasterType aux_zero_master = ZeroMatrix(TNumNodes, TNumNodesMaster);
 
        // We initialize the D and M derivatives operators
        for (IndexType i = 0; i < DoFSizeSlaveGeometry; ++i) {
            noalias(DeltaDOperator[i]) = aux_zero_slave;
            noalias(DeltaMOperator[i]) = aux_zero_master;
        }
        for (IndexType i = 0; i < DoFSizeMasterGeometry; ++i) {
            noalias(DeltaDOperator[i + DoFSizeSlaveGeometry]) = aux_zero_slave;
            noalias(DeltaMOperator[i + DoFSizeSlaveGeometry]) = aux_zero_master;
        }
    }
 
    void CalculateDeltaMortarOperators(
        KinematicVariables& rKinematicVariables,
        DerivativeDataType& rDerivativeData,
        const double rIntegrationWeight
        )
    {
        /* DEFINITIONS */
        const double det_j_slave = rKinematicVariables.DetjSlave;
        const Vector& vector_phi = rKinematicVariables.PhiLagrangeMultipliers;
        const Vector& vector_n1  = rKinematicVariables.NSlave;
        const Vector& vector_n2  = rKinematicVariables.NMaster;
 
        // Derivatives
        const array_1d<double, DoFSizeDerivativesDependence>& delta_det_j_slave = rDerivativeData.DeltaDetjSlave;
        const array_1d<array_1d<double, TNumNodes >, DoFSizeDerivativesDependence>& delta_phi = rDerivativeData.DeltaPhi;
        const array_1d<array_1d<double, TNumNodes >, DoFSizeDerivativesDependence>& delta_n1  = rDerivativeData.DeltaN1;
        const array_1d<array_1d<double, TNumNodesMaster >, DoFSizeDerivativesDependence>& delta_n2  = rDerivativeData.DeltaN2;
 
        for (IndexType i_node = 0; i_node < TNumNodes; ++i_node) {
            const double phi = vector_phi[i_node];
 
            for (IndexType j_node = 0; j_node < TNumNodes; ++j_node) {
                const double n1 = vector_n1[j_node];
 
                BaseClassType::DOperator(i_node, j_node) += det_j_slave * rIntegrationWeight * phi * n1;
 
                for (IndexType i = 0; i < DoFSizeSlaveGeometry; ++i) {
                    DeltaDOperator[i](i_node, j_node) += delta_det_j_slave[i] * rIntegrationWeight * phi* n1
                                                       + det_j_slave * rIntegrationWeight * delta_phi[i][i_node] * n1
                                                       + det_j_slave * rIntegrationWeight * phi* delta_n1[i][j_node];
                }
                if constexpr (TDim == 3) {
                    for (IndexType i = DoFSizeSlaveGeometry; i < DoFSizePairedGeometry; ++i) {
                        DeltaDOperator[i](i_node, j_node) += det_j_slave * rIntegrationWeight * phi * delta_n1[i][j_node];
                        DeltaDOperator[i](i_node, j_node) += delta_det_j_slave[i] * rIntegrationWeight * phi * n1;
                        DeltaDOperator[i](i_node, j_node) += det_j_slave * rIntegrationWeight * delta_phi[i][i_node] * n1;
                    }
                }
            }
            for (IndexType j_node = 0; j_node < TNumNodesMaster; ++j_node) {
                const double n2 = vector_n2[j_node];
 
                BaseClassType::MOperator(i_node, j_node) += det_j_slave * rIntegrationWeight * phi * n2;
 
                for (IndexType i = 0; i < DoFSizeSlaveGeometry; ++i) {
                    DeltaMOperator[i](i_node, j_node) += delta_det_j_slave[i] * rIntegrationWeight * phi* n2
                                                       + det_j_slave * rIntegrationWeight * delta_phi[i][i_node] * n2
                                                       + det_j_slave * rIntegrationWeight * phi* delta_n2[i][j_node];
                }
                if constexpr (TDim == 3) {
                    for (IndexType i = DoFSizeSlaveGeometry; i < DoFSizePairedGeometry; ++i) {
                        DeltaMOperator[i](i_node, j_node) += det_j_slave * rIntegrationWeight * phi * delta_n2[i][j_node];
                        DeltaMOperator[i](i_node, j_node) += delta_det_j_slave[i] * rIntegrationWeight * phi * n2;
                        DeltaMOperator[i](i_node, j_node) += det_j_slave * rIntegrationWeight * delta_phi[i][i_node] * n2;
                    }
                }
            }
        }
    }
 
    void PrintInfo(std::ostream& rOStream) const
    {
        BaseClassType::PrintInfo(rOStream);
 
        for (IndexType i = 0; i < TNumNodes * TDim; ++i) {
            rOStream << "DeltaDOperator_" << i << ": " << DeltaDOperator[i] << std::endl;
            rOStream << "DeltaMOperator_" << i << ": " << DeltaMOperator[i] << std::endl;
        }
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.save("DeltaDOperator", DeltaDOperator);
        rSerializer.save("DeltaMOperator", DeltaMOperator);
    }
 
    void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.load("DeltaDOperator", DeltaDOperator);
        rSerializer.load("DeltaMOperator", DeltaMOperator);
    }
 
 
 
 
 
}; // Class MortarOperatorWithDerivatives
 
template< const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class DualLagrangeMultiplierOperators
{
public:
 
    using KinematicVariables = MortarKinematicVariables<TNumNodes, TNumNodesMaster>;
 
    using GeometryMatrixType = BoundedMatrix<double, TNumNodes, TNumNodes>;
 
    KRATOS_CLASS_POINTER_DEFINITION( DualLagrangeMultiplierOperators );
 
 
    DualLagrangeMultiplierOperators(){}
 
    virtual ~DualLagrangeMultiplierOperators(){}
 
    GeometryMatrixType Me, De;
 
 
 
    void Initialize()
    {
        // We initialize the De and Me operators
        noalias(Me) = ZeroMatrix(TNumNodes, TNumNodes);
        noalias(De) = ZeroMatrix(TNumNodes, TNumNodes);
    }
 
    void CalculateAeComponents(
        KinematicVariables& rKinematicVariables,
        const double rIntegrationWeight
        )
    {
        /* DEFINITIONS */
        const Vector& n1 = rKinematicVariables.NSlave;
        const double det_j = rKinematicVariables.DetjSlave;
 
        noalias(De) += rIntegrationWeight * (ComputeDe(n1, det_j));
        noalias(Me) += rIntegrationWeight * det_j * outer_prod(n1, n1);
    }
 
    bool CalculateAe(GeometryMatrixType& Ae)
    {
        const double tolerance = std::numeric_limits<double>::epsilon();
 
        // We compute the norm
        const double norm_me = norm_frobenius(Me);
 
        // Now we normalize the matrix
        if (norm_me >= tolerance) {
            const GeometryMatrixType normalized_Me = Me/norm_me;
 
            // We compute the normalized inverse
            double aux_det;
            GeometryMatrixType normalized_inv_Me;
            MathUtils<double>::InvertMatrix(normalized_Me, normalized_inv_Me, aux_det, -1.0);
            const bool good_condition_number = MathUtils<double>::CheckConditionNumber(normalized_Me, normalized_inv_Me, tolerance, false);
            if (good_condition_number) {
                noalias(Ae) = (1.0/norm_me) * prod(De, normalized_inv_Me);
                return true;
            }
        }
 
        noalias(Ae) = IdentityMatrix(TNumNodes);
        return false;
    }
 
    template<class TArray>
    GeometryMatrixType ComputeDe(
        const TArray& N1,
        const double detJ
        ) const
    {
        GeometryMatrixType De;
 
        for (IndexType i = 0; i < TNumNodes; ++i) {
            for (IndexType j = 0; j < TNumNodes; ++j) {
                if (i == j) De(i,i) = detJ * N1[i];
                else De(i,j) = 0.0;
            }
        }
 
        return De;
    }
 
    void PrintInfo(std::ostream& rOStream) const {
        rOStream << "Me: " << Me << std::endl;
        rOStream << "De: " << De << std::endl;
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    virtual void save(Serializer& rSerializer) const
    {
        rSerializer.save("Me", Me);
        rSerializer.save("De", De);
    }
 
    virtual void load(Serializer& rSerializer)
    {
        rSerializer.load("Me", Me);
        rSerializer.load("De", De);
    }
 
 
 
 
 
}; // Class DualLagrangeMultiplierOperators
 
template< const SizeType TDim, const SizeType TNumNodes, bool TFrictional, const SizeType TNumNodesMaster = TNumNodes>
class DualLagrangeMultiplierOperatorsWithDerivatives
    : public DualLagrangeMultiplierOperators<TNumNodes, TNumNodesMaster>
{
public:
 
    using BaseClassType = DualLagrangeMultiplierOperators<TNumNodes, TNumNodesMaster>;
 
    using KinematicVariablesType = MortarKinematicVariablesWithDerivatives<TDim, TNumNodes, TNumNodesMaster>;
 
    using DerivativeDataFrictionalType = DerivativeDataFrictional<TDim, TNumNodes, TNumNodesMaster>;
 
    using DerivativeFrictionalessDataType = DerivativeData<TDim, TNumNodes, TNumNodesMaster>;
 
    using DerivativeDataType = typename std::conditional<TFrictional, DerivativeDataFrictionalType, DerivativeFrictionalessDataType>::type;
 
    using GeometryMatrixType = BoundedMatrix<double, TNumNodes, TNumNodes>;
 
    // Auxiliary sizes
    static const SizeType DoFSizeSlaveGeometry = (TNumNodes * TDim);
 
    static const SizeType DoFSizeMasterGeometry = (TNumNodesMaster * TDim);
 
    static const SizeType DoFSizePairedGeometry = DoFSizeSlaveGeometry + DoFSizeMasterGeometry;
 
    static const SizeType DoFSizeDerivativesDependence = (TDim == 2) ? DoFSizeSlaveGeometry : DoFSizePairedGeometry;
 
    KRATOS_CLASS_POINTER_DEFINITION( DualLagrangeMultiplierOperatorsWithDerivatives );
 
 
    DualLagrangeMultiplierOperatorsWithDerivatives(){}
 
    ~DualLagrangeMultiplierOperatorsWithDerivatives() override{}
 
    // Derivatives matrices
    array_1d<GeometryMatrixType, DoFSizeDerivativesDependence> DeltaMe, DeltaDe;
 
 
 
    void Initialize()
    {
        BaseClassType::Initialize();
 
        // Derivatives matrices
        const BoundedMatrix<double, TNumNodes, TNumNodes> zeromatrix = ZeroMatrix(TNumNodes, TNumNodes);
        for (IndexType i = 0; i < DoFSizeDerivativesDependence; ++i) {
            noalias(DeltaMe[i]) = zeromatrix;
            noalias(DeltaDe[i]) = zeromatrix;
        }
    }
 
    void CalculateDeltaAeComponents(
        KinematicVariablesType& rKinematicVariables,
        DerivativeDataType& rDerivativeData,
        const double rIntegrationWeight
        )
    {
        /* DEFINITIONS */
        const double det_j_slave = rKinematicVariables.DetjSlave;
        const Vector& n1 = rKinematicVariables.NSlave;
 
        BaseClassType::CalculateAeComponents(rKinematicVariables, rIntegrationWeight);
 
        for (IndexType i = 0; i < DoFSizeDerivativesDependence; ++i) {
            const double delta_det_j = rDerivativeData.DeltaDetjSlave[i];
            const array_1d<double, TNumNodes>& delta_n1 = rDerivativeData.DeltaN1[i];
 
            noalias(DeltaDe[i]) += rIntegrationWeight * this->ComputeDe( n1, delta_det_j )
                                +  rIntegrationWeight * this->ComputeDe( delta_n1, det_j_slave );
 
            noalias(DeltaMe[i]) += rIntegrationWeight * delta_det_j * outer_prod(n1, n1)
                                +  rIntegrationWeight * det_j_slave * (outer_prod(delta_n1, n1) + outer_prod(n1, delta_n1));
        }
    }
 
    bool CalculateDeltaAe(DerivativeDataType& rDerivativeData)
    {
        const double tolerance = std::numeric_limits<double>::epsilon();
 
        // We compute the norm
        const double norm_Me = norm_frobenius(BaseClassType::Me);
 
        // Now we normalize the matrix
        const GeometryMatrixType normalized_Me = BaseClassType::Me/norm_Me;
 
        // We compute the normalized inverse
        double aux_det;
        GeometryMatrixType normalized_inv_Me;
        MathUtils<double>::InvertMatrix(normalized_Me, normalized_inv_Me, aux_det, -1.0);
        const bool good_condition_number = MathUtils<double>::CheckConditionNumber(normalized_Me, normalized_inv_Me, tolerance, false);
        if (!good_condition_number) {
            return false;
        }
 
        // Now we compute the inverse
        const GeometryMatrixType inv_Me = normalized_inv_Me/norm_Me;
 
        noalias(rDerivativeData.Ae) = prod(BaseClassType::De, inv_Me);
 
        array_1d<GeometryMatrixType , DoFSizeDerivativesDependence>& delta_Ae = rDerivativeData.DeltaAe;
 
        for (IndexType i = 0; i < DoFSizeDerivativesDependence; ++i) {
            const GeometryMatrixType aux_matrix = DeltaDe[i] - prod(rDerivativeData.Ae, DeltaMe[i]);
            noalias(delta_Ae[i]) = prod(aux_matrix, inv_Me);
        }
 
        return true;
    }
 
    void PrintInfo(std::ostream& rOStream) const
    {
        BaseClassType::PrintInfo(rOStream);
 
        // Derivatives matrices
        for (IndexType i = 0; i < DoFSizeDerivativesDependence; ++i) {
            rOStream << "DeltaMe_" << i << ": " << DeltaMe[i] << std::endl;
            rOStream << "DeltaDe_" << i << ": " << DeltaDe[i] << std::endl;
        }
    }
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.save("DeltaMe", DeltaMe);
        rSerializer.save("DeltaDe", DeltaDe);
    }
 
    void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseClassType );
        rSerializer.load("DeltaMe", DeltaMe);
        rSerializer.load("DeltaDe", DeltaDe);
    }
 
 
 
 
 
}; // Class DualLagrangeMultiplierOperatorsWithDerivatives
 
template<const SizeType TNumNodes, const SizeType TNumNodesMaster = TNumNodes>
class PointBelong
    : public Point
{
public:
 
    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;
 
    KRATOS_CLASS_POINTER_DEFINITION( PointBelong );
 
 
    PointBelong():
        Point()
    {}
 
    PointBelong(const array_1d<double, 3> Coords):
        Point(Coords)
    {}
 
    PointBelong(const array_1d<double, 3> Coords, const BelongType& ThisBelongs):
        Point(Coords),
        mBelongs(ThisBelongs)
    {}
 
    ~PointBelong() override= default;
 
 
 
    void SetBelong(BelongType ThisBelongs)
    {
        mBelongs = ThisBelongs;
    }
 
    BelongType GetBelong() const
    {
        return mBelongs;
    }
private:
 
    BelongType mBelongs; 
 
 
 
 
 
 
}; // Class PointBelong
 
 
 
 
 
}// namespace Kratos.