exact_mortar_segmentation_utility.h

Go to the documentation of this file.

//
//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
#if !defined(KRATOS_EXACT_MORTAR_INTEGRATION_UTILITY_H_INCLUDED)
#define KRATOS_EXACT_MORTAR_INTEGRATION_UTILITY_H_INCLUDED
 
// System includes
#include <iostream>
#include <iomanip>
#include <cmath>
 
// External includes
 
// Project includes
#include "includes/mortar_classes.h"
 
// The geometry of the triangle for the "tessellation"
/* LINES */
#include "geometries/line_2d_2.h"
/* TRIANGLES */
#include "geometries/triangle_3d_3.h"
/* QUADRILATERALS */
#include "geometries/quadrilateral_3d_4.h"
 
// /* The integration points (we clip triangles in 3D, so with line and triangle is enough)*/
// #include "integration/line_gauss_legendre_integration_points.h"
#include "integration/triangle_gauss_legendre_integration_points.h"
 
/* Utilities */
#include "utilities/math_utils.h"
#include "utilities/mortar_utilities.h"
 
namespace Kratos {
 
 
    typedef Point PointType;
    typedef Node NodeType;
    typedef Geometry<NodeType> GeometryType;
    typedef Geometry<PointType> GeometryPointType;
 
    typedef GeometryData::IntegrationMethod IntegrationMethod;
    typedef IntegrationPoint<2> IntegrationPointType;
    typedef GeometryType::IntegrationPointsArrayType IntegrationPointsType;
 
    typedef std::size_t SizeType;
 
 
 
 
template <SizeType TDim, SizeType TNumNodes, bool TBelong = false, SizeType TNumNodesMaster = TNumNodes>
class KRATOS_API(KRATOS_CORE) ExactMortarIntegrationUtility
{
public:
 
    typedef 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 BelongType;
 
    typedef PointBelong<TNumNodes, TNumNodesMaster> PointBelongType;
 
    typedef std::vector<array_1d<PointBelongType, TDim>> VectorArrayPointsBelong;
 
    typedef std::vector<array_1d<PointType, TDim>> VectorArrayPoints;
 
    typedef typename std::conditional<TBelong, VectorArrayPointsBelong,VectorArrayPoints>::type ConditionArrayListType;
 
    typedef std::vector<PointBelongType> VectorPointsBelong;
 
    typedef std::vector<PointType> VectorPoints;
 
    typedef typename std::conditional<TBelong, VectorPointsBelong, VectorPoints>::type PointListType;
 
    typedef array_1d<PointBelongType, 3> ArrayPointsBelong;
 
    typedef array_1d<PointType, 3> ArrayPoints;
 
    typedef typename std::conditional<TBelong, ArrayPointsBelong, ArrayPoints>::type ArrayTriangleType;
 
    typedef Line2D2<Point> LineType;
 
    typedef Triangle3D3<Point> TriangleType;
 
    typedef typename std::conditional<TDim == 2, LineType, TriangleType>::type DecompositionType;
 
    typedef std::size_t IndexType;
 
    static constexpr double ZeroTolerance = std::numeric_limits<double>::epsilon();
 
    KRATOS_CLASS_POINTER_DEFINITION(ExactMortarIntegrationUtility);
 
 
    ExactMortarIntegrationUtility(
        const SizeType IntegrationOrder = 0,
        const double DistanceThreshold = std::numeric_limits<double>::max(),
        const SizeType EchoLevel = 0,
        const double ZeroToleranceFactor = 1.0,
        const bool ConsiderDelaunator = false
        ) :mIntegrationOrder(IntegrationOrder),
           mDistanceThreshold(DistanceThreshold),
           mEchoLevel(EchoLevel),
           mZeroToleranceFactor(ZeroToleranceFactor),
           mConsiderDelaunator(ConsiderDelaunator)
    {
        GetIntegrationMethod();
    }
 
    virtual ~ExactMortarIntegrationUtility() = default;
 
 
 
    bool GetExactIntegration(
        const GeometryType& rOriginalSlaveGeometry,
        const array_1d<double, 3>& rSlaveNormal,
        const GeometryType& rOriginalMasterGeometry,
        const array_1d<double, 3>& rMasterNormal,
        ConditionArrayListType& rConditionsPointsSlave
        );
 
    bool GetExactIntegration(
        const GeometryType& rOriginalSlaveGeometry,
        const array_1d<double, 3>& rSlaveNormal,
        const GeometryType& rOriginalMasterGeometry,
        const array_1d<double, 3>& rMasterNormal,
        IntegrationPointsType& rIntegrationPointsSlave
        );
 
    bool GetExactAreaIntegration(
        const GeometryType& rOriginalSlaveGeometry,
        const array_1d<double, 3>& rSlaveNormal,
        const GeometryType& rOriginalMasterGeometry,
        const array_1d<double, 3>& rMasterNormal,
        double& rArea
        );
 
    void GetTotalArea(
        const GeometryType& rOriginalSlaveGeometry,
        ConditionArrayListType& rConditionsPointsSlave,
        double& rArea
        );
 
    bool TestGetExactIntegration(
        Condition::Pointer pSlaveCond,
        Condition::Pointer pMasterCond,
        Matrix& rCustomSolution
        );
 
 
    double TestGetExactAreaIntegration(
        Condition::Pointer pSlaveCond,
        Condition::Pointer pMasterCond
        );
 
 
    double TestGetExactAreaIntegration(
        ModelPart& rMainModelPart,
        Condition::Pointer pSlaveCond
        );
 
    void TestIODebug(
        ModelPart& rMainModelPart,
        const std::string IOConsidered = "GiD"
        );
 
 
    SizeType& GetIntegrationOrder()
    {
        return mIntegrationOrder;
    }
 
    void SetIntegrationOrder(const SizeType IntegrationOrder)
    {
        mIntegrationOrder = IntegrationOrder;
        GetIntegrationMethod();
    }
 
    double& GetDistanceThreshold()
    {
        return mDistanceThreshold;
    }
 
    void SetDistanceThreshold(const double DistanceThreshold)
    {
        mDistanceThreshold = DistanceThreshold;
    }
 
    SizeType& GetEchoLevel()
    {
        return mEchoLevel;
    }
 
    void SetEchoLevel(const SizeType EchoLevel)
    {
        mEchoLevel = EchoLevel;
    }
 
    double& GetZeroToleranceFactor()
    {
        return mZeroToleranceFactor;
    }
 
    void SetZeroToleranceFactor(const double ZeroToleranceFactor)
    {
        mZeroToleranceFactor = ZeroToleranceFactor;
    }
 
    bool& GetConsiderDelaunator()
    {
        return mConsiderDelaunator;
    }
 
    void SetConsiderDelaunator(const bool ConsiderDelaunator)
    {
        mConsiderDelaunator = ConsiderDelaunator;
    }
 
protected:
 
 
 
 
    void GetIntegrationMethod();
 
    GeometryType::IntegrationPointsArrayType GetIntegrationTriangle();
 
    template<SizeType TSizeCheck = TNumNodes>
    static inline bool CheckAllInside(const array_1d<bool, TSizeCheck>& rAllInside)
    {
        for (IndexType i_node = 0; i_node < TSizeCheck; ++i_node)
            if (!rAllInside[i_node])
                return false;
 
        return true;
    }
 
    static inline bool Clipping2D(
        PointType& rPointIntersection,
        const PointType& rPointOrig1,
        const PointType& rPointOrig2,
        const PointType& rPointDest1,
        const PointType& rPointDest2
        )
    {
        const array_1d<double, 3>& r_coord_point_orig1 = rPointOrig1.Coordinates();
        const array_1d<double, 3>& r_coord_point_orig2 = rPointOrig2.Coordinates();
        const array_1d<double, 3>& r_coord_point_dest1 = rPointDest1.Coordinates();
        const array_1d<double, 3>& r_coord_point_dest2 = rPointDest2.Coordinates();
 
        const double s_orig1_orig2_x = r_coord_point_orig2[0] - r_coord_point_orig1[0];
        const double s_orig1_orig2_y = r_coord_point_orig2[1] - r_coord_point_orig1[1];
        const double s_dest1_dest2_x = r_coord_point_dest2[0] - r_coord_point_dest1[0];
        const double s_dest1_dest2_y = r_coord_point_dest2[1] - r_coord_point_dest1[1];
 
        const double denom = s_orig1_orig2_x * s_dest1_dest2_y -
                             s_dest1_dest2_x * s_orig1_orig2_y;
 
        const double tolerance = 1.0e-15;
//         const double tolerance = std::numeric_limits<double>::epsilon();
 
        if (std::abs(denom) < tolerance) // NOTE: Collinear
            return false;
 
        const double s_orig1_dest1_x = r_coord_point_orig1[0] - r_coord_point_dest1[0];
        const double s_orig1_dest1_y = r_coord_point_orig1[1] - r_coord_point_dest1[1];
 
        const double s = (s_orig1_orig2_x * s_orig1_dest1_y - s_orig1_orig2_y * s_orig1_dest1_x)/denom;
 
        const double t = (s_dest1_dest2_x * s_orig1_dest1_y - s_dest1_dest2_y * s_orig1_dest1_x)/denom;
 
        if (s >= -tolerance && s <= (1.0 + tolerance) && t >= -tolerance && t <= (1.0 + tolerance)) {
            rPointIntersection.Coordinates()[0] = r_coord_point_orig1[0] + t * s_orig1_orig2_x;
            rPointIntersection.Coordinates()[1] = r_coord_point_orig1[1] + t * s_orig1_orig2_y;
 
            return true;
        } else
            return false;
    }
 
    static inline array_1d<double, 3> GetNormalVector2D(const array_1d<double, 3>& rVector)
    {
        array_1d<double, 3> normal;
 
        normal[0] = -rVector[1];
        normal[1] = rVector[0];
        normal[2] = 0.0;
 
        return normal;
    }
 
    static inline double AnglePoints(
        const PointType& rPointOrig1,
        const PointType& rPointOrig2,
        const array_1d<double, 3>& rAxis1,
        const array_1d<double, 3>& rAxis2
        )
    {
        array_1d<double, 3> local_edge = rPointOrig2.Coordinates() - rPointOrig1.Coordinates();
        if (norm_2(local_edge) > 0.0)
            local_edge /= norm_2(local_edge);
 
        const double xi  = inner_prod(rAxis1, local_edge);
        const double eta = inner_prod(rAxis2, local_edge);
 
        return (std::atan2(eta, xi));
    }
 
    static inline bool CheckPoints(
        const PointType& rPointOrig,
        const PointType& rPointDest
        )
    {
//         const double tolerance = std::numeric_limits<double>::epsilon(); // NOTE: Giving some problems, too tight
        const double tolerance = 1.0e-15;
        return rPointDest.Distance(rPointOrig) < tolerance;
    }
 
    static inline double FastTriagleCheck2D(
        const PointType& rPointOrig1,
        const PointType& rPointOrig2,
        const PointType& rPointOrig3
        )
    {
        const double x10 = rPointOrig2.X() - rPointOrig1.X();
        const double y10 = rPointOrig2.Y() - rPointOrig1.Y();
 
        const double x20 = rPointOrig3.X() - rPointOrig1.X();
        const double y20 = rPointOrig3.Y() - rPointOrig1.Y();
 
        //Jacobian is calculated:
        //  |dx/dxi  dx/deta|   |x1-x0   x2-x0|
        //J=|               |=  |             |
        //  |dy/dxi  dy/deta|   |y1-y0   y2-y0|
 
        return x10 * y20 - y10 * x20;
    }
 
    template<SizeType TSizeCheck = TNumNodes>
    inline void PushBackPoints(
        VectorPoints& rPointList,
        const array_1d<bool, TSizeCheck>& rAllInside,
        GeometryPointType& rThisGeometry
        )
    {
        for (IndexType i_node = 0; i_node < TSizeCheck; ++i_node) {
            if (rAllInside[i_node]) {
                // We check if the node already exists
                bool add_point = true;
                for (IndexType iter = 0; iter < rPointList.size(); ++iter)
                    if (CheckPoints(rThisGeometry[i_node], rPointList[iter])) add_point = false;
 
                if (add_point)
                    rPointList.push_back(rThisGeometry[i_node]);
            }
        }
    }
 
    template<SizeType TSizeCheck = TNumNodes>
    inline void PushBackPoints(
        VectorPointsBelong& rPointList,
        const array_1d<bool, TSizeCheck>& rAllInside,
        GeometryPointType& rThisGeometry,
        const PointBelongs& rThisBelongs
        )
    {
        for (IndexType i_node = 0; i_node < TSizeCheck; ++i_node) {
            if (rAllInside[i_node]) {
                // We check if the node already exists
                bool add_point = true;
                for (IndexType iter = 0; iter < rPointList.size(); ++iter) {
                    if (CheckPoints(rThisGeometry[i_node], rPointList[iter])) {
                        add_point = false;
                    }
                }
 
                if (add_point) {
                    const IndexType initial_index = (rThisBelongs == PointBelongs::Master) ? TNumNodes : 0;
                    rPointList.push_back(PointBelong<TNumNodes, TNumNodesMaster>(rThisGeometry[i_node].Coordinates(), static_cast<BelongType>(initial_index + i_node)));
                }
            }
        }
    }
 
    template<SizeType TSizeCheck = TNumNodes>
    inline void CheckInside(
        array_1d<bool, TSizeCheck>& rAllInside,
        GeometryPointType& rGeometry1,
        GeometryPointType& rGeometry2,
        const double Tolerance
        )
    {
        for (IndexType i_node = 0; i_node < TSizeCheck; ++i_node) {
            GeometryType::CoordinatesArrayType projected_gp_local;
            rAllInside[i_node] = rGeometry1.IsInside( rGeometry2[i_node].Coordinates( ), projected_gp_local, Tolerance) ;
        }
    }
 
    inline std::vector<IndexType> ComputeAnglesIndexes(
        PointListType& rPointList,
        const array_1d<double, 3>& rNormal
        ) const;
 
    inline void ComputeClippingIntersections(
        PointListType& rPointList,
        const GeometryPointType& rSlaveGeometry,
        const GeometryPointType& rMasterGeometry,
        const PointType& rRefCenter
        );
 
    template <class TGeometryType = GeometryType>
    inline bool TriangleIntersections(
        ConditionArrayListType& rConditionsPointsSlave,
        PointListType& rPointList,
        const TGeometryType& rOriginalSlaveGeometry,
        const GeometryPointType& rSlaveGeometry,
        const GeometryPointType& rMasterGeometry,
        const array_1d<double, 3>& rSlaveTangentXi,
        const array_1d<double, 3>& rSlaveTangentEta,
        const PointType& rRefCenter,
        const bool IsAllInside = false
        );
 
    static inline bool CheckCenterIsInside(
        const array_1d<double, 2>& rAuxiliarCenterLocalCoordinates,
        const SizeType NumNodes = TNumNodes
        );
 
 
 
private:
 
    SizeType mIntegrationOrder;              
    double mDistanceThreshold;               
    SizeType mEchoLevel;                     
    double mZeroToleranceFactor;             
    IntegrationMethod mAuxIntegrationMethod; 
    bool mConsiderDelaunator;                
 
 
 
 
 
 
};  // Class ExactMortarIntegrationUtility
}
#endif /* KRATOS_EXACT_MORTAR_INTEGRATION_UTILITY_H_INCLUDED defined */