derivatives_utilities.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"
 
/* Includes */
#include "includes/model_part.h"
#include "includes/mortar_classes.h"
 
/* 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, bool TFrictional, const bool TNormalVariation, const SizeType TNumNodesMaster = TNumNodes>
class KRATOS_API(CONTACT_STRUCTURAL_MECHANICS_APPLICATION) DerivativesUtilities
{
public:
 
    using IndexType = std::size_t;
 
    using GeometryType = Geometry<Node>;
 
    using NodesArrayType = Geometry<Node>::PointsArrayType;
 
    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, DerivativeDataFrictional<TDim, TNumNodes, TNumNodesMaster>, DerivativeData<TDim, TNumNodes, TNumNodesMaster> >::type;
 
    using GeneralVariables = MortarKinematicVariablesWithDerivatives<TDim, TNumNodes, TNumNodesMaster>;
 
    using AeData = DualLagrangeMultiplierOperatorsWithDerivatives<TDim, TNumNodes, TFrictional, TNumNodesMaster>;
 
    using MortarConditionMatrices = MortarOperatorWithDerivatives<TDim, TNumNodes, TFrictional, TNumNodesMaster>;
 
    // The threshold coefficient considered for checking
    static constexpr double CheckThresholdCoefficient = 1.0e-12;
 
    static constexpr double ZeroTolerance = std::numeric_limits<double>::epsilon();
 
    KRATOS_CLASS_POINTER_DEFINITION( DerivativesUtilities );
 
 
    static void CalculateDeltaDetjSlave(
        const DecompositionType& DecompGeom,
        const GeneralVariables& rVariables,
        DerivativeDataType& rDerivativeData
        );
 
    static inline array_1d<array_1d<double, 3>, TDim * TNumNodes> GPDeltaNormalSlave(
        const Matrix& rJacobian,
        const Matrix& rDNDe
        );
 
    static inline array_1d<array_1d<double, 3>, TDim * TNumNodesMaster> GPDeltaNormalMaster(
        const Matrix& rJacobian,
        const Matrix& rDNDe
        );
 
    static array_1d<array_1d<double, 3>, TDim * TNumNodes> DeltaNormalCenter(const GeometryType& rThisGeometry);
 
    static void CalculateDeltaNormalSlave(
        array_1d<BoundedMatrix<double, TNumNodes, TDim>, TNumNodes * TDim>& rDeltaNormal,
        GeometryType& rThisGeometry
        );
 
    static void CalculateDeltaNormalMaster(
        array_1d<BoundedMatrix<double, TNumNodesMaster, TDim>, TNumNodesMaster * TDim>& rDeltaNormal,
        const GeometryType& rThisGeometry
        );
 
    static void CalculateDeltaCellVertex(
        const GeneralVariables& rVariables,
        DerivativeDataType& rDerivativeData,
        const array_1d<BelongType, TDim>& rTheseBelongs,
        const NormalDerivativesComputation ConsiderNormalVariation,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const array_1d<double, 3>& rNormal
        );
 
    static inline void CalculateDeltaN1(
        const GeneralVariables& rVariables,
        DerivativeDataType& rDerivativeData,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const array_1d<double, 3>& rSlaveNormal,
        const DecompositionType& rDecompGeom,
        const PointType& rLocalPointDecomp,
        const PointType& rLocalPointParent,
        const NormalDerivativesComputation ConsiderNormalVariation = NormalDerivativesComputation::NO_DERIVATIVES_COMPUTATION
        );
 
    static void CalculateDeltaN(
        const GeneralVariables& rVariables,
        DerivativeDataType& rDerivativeData,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const array_1d<double, 3>& rSlaveNormal,
        const array_1d<double, 3>& rMasterNormal,
        const DecompositionType& rDecompGeom,
        const PointType& rLocalPointDecomp,
        const PointType& rLocalPointParent,
        const NormalDerivativesComputation ConsiderNormalVariation = NormalDerivativesComputation::NO_DERIVATIVES_COMPUTATION,
        const bool DualLM = false
        );
 
    Matrix& CalculateDeltaPosition(
        Matrix& rDeltaPosition,
        const GeometryType& rThisGeometry,
        const ConditionArrayType& rLocalCoordinates
        );
 
    static inline Matrix& CalculateDeltaPosition(
        Matrix& rDeltaPosition,
        const GeometryType& ThisGeometry
        );
 
    static inline void CalculateDeltaPosition(
        Vector& rDeltaPosition,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const IndexType IndexNode
        );
 
    static inline void CalculateDeltaPosition(
        Vector& rDeltaPosition,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const IndexType IndexNode,
        const IndexType iDoF
        );
 
    static inline void CalculateDeltaPosition(
        double& rDeltaPosition,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const IndexType IndexNode,
        const IndexType iDoF
        );
 
    static bool CalculateAeAndDeltaAe(
        const GeometryType& rSlaveGeometry,
        const array_1d<double, 3>& rSlaveNormal,
        const GeometryType& rMasterGeometry,
        DerivativeDataType& rDerivativeData,
        GeneralVariables& rVariables,
        const NormalDerivativesComputation ConsiderNormalVariation,
        ConditionArrayListType& rConditionsPointsSlave,
        GeometryData::IntegrationMethod ThisIntegrationMethod,
        const double AxiSymCoeff = 1.0
        );
 
private:
 
 
 
    static inline array_1d<double, 3> LocalDeltaVertex(
        const array_1d<double, 3>& rNormal,
        const array_1d<double, 3>& rDeltaNormal,
        const IndexType iDoF,
        const IndexType iBelong,
        const NormalDerivativesComputation ConsiderNormalVariation,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        double Coeff = 1.0
        );
 
    static inline BoundedMatrix<double, 3, 3> ComputeRenormalizerMatrix(
        const array_1d<double, 3>& rDiffVector,
        const array_1d<double, 3>& rDeltaNormal
        );
 
    static inline array_1d<double, 3> PreviousNormalGeometry(
        const GeometryType& rThisGeometry,
        const GeometryType::CoordinatesArrayType& rPointLocal
        );
 
    static inline void ConvertAuxHashIndex(
        const IndexType AuxIndex,
        IndexType& riBelongSlaveStart,
        IndexType& riBelongSlaveEnd,
        IndexType& riBelongMasterStart,
        IndexType& riBelongMasterEnd
        );
 
    static inline void DeltaPointLocalCoordinatesSlave(
        array_1d<double, 2>& rResult,
        const array_1d<double, 3>& rDeltaPoint,
        const Matrix& rDNDe,
        const GeometryType& rThisGeometry,
        const array_1d<double, 3>& rThisNormal
        );
    static inline void DeltaPointLocalCoordinatesMaster(
        array_1d<double, 2>& rResult,
        const array_1d<double, 3>& rDeltaPoint,
        const Matrix& rDNDe,
        const GeometryType& rThisGeometry,
        const array_1d<double, 3>& rThisNormal
        );
 
    static inline double LocalDeltaSegmentN1(
        const array_1d<array_1d<double, 3>, TDim * TNumNodes>& rDeltaNormal,
        const array_1d<double, 3>& rSlaveNormal,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const Vector& rN1,
        const Matrix& rDNDe1,
        const IndexType MortarNode,
        const IndexType iNode,
        const IndexType iDoF,
        const NormalDerivativesComputation ConsiderNormalVariation
        );
 
     static inline double LocalDeltaSegmentN2(
         const array_1d<array_1d<double, 3>, TDim * TNumNodes>& rDeltaNormal,
         const array_1d<double, 3>& rSlaveNormal,
         const GeometryType& rSlaveGeometry,
         const GeometryType& rMasterGeometry,
         const Vector& rN2,
         const Matrix& rDNDe2,
         const IndexType MortarNode,
         const IndexType iNode,
         const IndexType iDoF,
         const NormalDerivativesComputation ConsiderNormalVariation
         );
 
};// class DerivativesUtilities
 
class ImplementationDerivativesUtilities
{
public:
 
    static constexpr double ZeroTolerance = std::numeric_limits<double>::epsilon();
 
 
    template<const SizeType TNumNodes1, const SizeType TNumNodes2>
    static inline BoundedMatrix<double, 3, 3> ComputeRenormalizerMatrix(
        const BoundedMatrix<double, TNumNodes1, 3>& rDiffMatrix,
        const BoundedMatrix<double, TNumNodes2, 3>& rDeltaNormal,
        const IndexType iGeometry
        )
    {
        KRATOS_TRY
 
        BoundedMatrix<double, 3, 3> aux_matrix;
        for (IndexType itry = 0; itry < 3; ++itry) {
            if (rDeltaNormal(iGeometry, itry) > ZeroTolerance) {
 
                const IndexType aux_index_1 = itry == 2 ? 0 : itry + 1;
                const IndexType aux_index_2 = itry == 2 ? 1 : (itry == 1 ? 0 : 2);
 
                const double diff = rDeltaNormal(iGeometry, aux_index_1) + rDeltaNormal(iGeometry, aux_index_2);
                const double coeff = rDeltaNormal(iGeometry, itry);
 
                aux_matrix(0, aux_index_1) = 1.0;
                aux_matrix(0, aux_index_2) = 1.0;
                aux_matrix(1, aux_index_1) = 1.0;
                aux_matrix(1, aux_index_2) = 1.0;
                aux_matrix(2, aux_index_1) = 1.0;
                aux_matrix(2, aux_index_2) = 1.0;
 
                aux_matrix(0, itry) = (rDiffMatrix(iGeometry, 0) - diff)/coeff;
                aux_matrix(1, itry) = (rDiffMatrix(iGeometry, 1) - diff)/coeff;
                aux_matrix(2, itry) = (rDiffMatrix(iGeometry, 2) - diff)/coeff;
 
                return aux_matrix;
            }
        }
 
        return IdentityMatrix(3, 3);
 
        KRATOS_CATCH("")
    }
 
}; // class ImplementationDerivativesUtilities
 
} // namespace Kratos