triangle_2d_3.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//                   Janosch Stascheit
//                   Felix Nagel
//  contributors:    Hoang Giang Bui
//                   Josep Maria Carbonell
//
 
#if !defined(KRATOS_TRIANGLE_2D_3_H_INCLUDED )
#define  KRATOS_TRIANGLE_2D_3_H_INCLUDED
 
 
// System includes
 
// External includes
 
// Project includes
#include "geometries/line_2d_2.h"
#include "integration/triangle_gauss_legendre_integration_points.h"
#include "integration/triangle_collocation_integration_points.h"
#include "utilities/intersection_utilities.h"
 
namespace Kratos
{
 
 
 
 
 
template<class TPointType> class Triangle2D3
    : public Geometry<TPointType>
{
public:
 
    typedef Geometry<TPointType> BaseType;
 
    typedef Line2D2<TPointType> EdgeType;
 
    KRATOS_CLASS_POINTER_DEFINITION( Triangle2D3 );
 
    typedef GeometryData::IntegrationMethod IntegrationMethod;
 
    typedef typename BaseType::GeometriesArrayType GeometriesArrayType;
 
    typedef TPointType PointType;
 
    typedef typename BaseType::IndexType IndexType;
 
    typedef typename BaseType::SizeType SizeType;
 
    typedef  typename BaseType::PointsArrayType PointsArrayType;
 
    typedef typename BaseType::CoordinatesArrayType CoordinatesArrayType;
 
    typedef typename BaseType::IntegrationPointType IntegrationPointType;
 
    typedef typename BaseType::IntegrationPointsArrayType IntegrationPointsArrayType;
 
    typedef typename BaseType::IntegrationPointsContainerType IntegrationPointsContainerType;
 
    typedef typename BaseType::ShapeFunctionsValuesContainerType ShapeFunctionsValuesContainerType;
 
    typedef typename BaseType::ShapeFunctionsLocalGradientsContainerType ShapeFunctionsLocalGradientsContainerType;
 
    typedef typename BaseType::JacobiansType JacobiansType;
 
    typedef typename BaseType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType;
 
    typedef typename BaseType::ShapeFunctionsSecondDerivativesType
    ShapeFunctionsSecondDerivativesType;
 
    typedef typename BaseType::ShapeFunctionsThirdDerivativesType
    ShapeFunctionsThirdDerivativesType;
 
    typedef typename BaseType::NormalType NormalType;
 
 
//     Triangle2D3( const PointType& FirstPoint,
//                  const PointType& SecondPoint,
//                  const PointType& ThirdPoint )
//         : BaseType( PointsArrayType(), &msGeometryData )
//     {
//         this->Points().push_back( typename PointType::Pointer( new PointType( FirstPoint ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( SecondPoint ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( ThirdPoint ) ) );
//     }
 
    Triangle2D3( typename PointType::Pointer pFirstPoint,
                 typename PointType::Pointer pSecondPoint,
                 typename PointType::Pointer pThirdPoint )
        : BaseType( PointsArrayType(), &msGeometryData )
    {
        this->Points().push_back( pFirstPoint );
        this->Points().push_back( pSecondPoint );
        this->Points().push_back( pThirdPoint );
    }
 
    explicit Triangle2D3( const PointsArrayType& ThisPoints )
        : BaseType( ThisPoints, &msGeometryData )
    {
        if ( this->PointsNumber() != 3 )
            KRATOS_ERROR << "Invalid points number. Expected 3, given " << this->PointsNumber() << std::endl;
    }
 
    explicit Triangle2D3(
        const IndexType GeometryId,
        const PointsArrayType& rThisPoints
        ) : BaseType(GeometryId, rThisPoints, &msGeometryData)
    {
        KRATOS_ERROR_IF( this->PointsNumber() != 3 ) << "Invalid points number. Expected 3, given " << this->PointsNumber() << std::endl;
    }
 
    explicit Triangle2D3(
        const std::string& rGeometryName,
        const PointsArrayType& rThisPoints
        ) : BaseType(rGeometryName, rThisPoints, &msGeometryData)
    {
        KRATOS_ERROR_IF(this->PointsNumber() != 3) << "Invalid points number. Expected 3, given " << this->PointsNumber() << std::endl;
    }
 
    Triangle2D3( Triangle2D3 const& rOther )
        : BaseType( rOther )
    {
    }
 
    template<class TOtherPointType> explicit Triangle2D3( Triangle2D3<TOtherPointType> const& rOther )
        : BaseType( rOther )
    {
    }
 
    ~Triangle2D3() override {}
 
    GeometryData::KratosGeometryFamily GetGeometryFamily() const override
    {
        return GeometryData::KratosGeometryFamily::Kratos_Triangle;
    }
 
    GeometryData::KratosGeometryType GetGeometryType() const override
    {
        return GeometryData::KratosGeometryType::Kratos_Triangle2D3;
    }
 
 
    Triangle2D3& operator=( const Triangle2D3& rOther )
    {
        BaseType::operator=( rOther );
        return *this;
    }
 
    template<class TOtherPointType>
    Triangle2D3& operator=( Triangle2D3<TOtherPointType> const & rOther )
    {
        BaseType::operator=( rOther );
        return *this;
    }
 
 
    typename BaseType::Pointer Create(
        PointsArrayType const& rThisPoints
        ) const override
    {
        return typename BaseType::Pointer( new Triangle2D3( rThisPoints ) );
    }
 
    typename BaseType::Pointer Create(
        const IndexType NewGeometryId,
        PointsArrayType const& rThisPoints
        ) const override
    {
        return typename BaseType::Pointer( new Triangle2D3( NewGeometryId, rThisPoints ) );
    }
 
    typename BaseType::Pointer Create(
        const BaseType& rGeometry
        ) const override
    {
        auto p_geometry = typename BaseType::Pointer( new Triangle2D3( rGeometry.Points() ) );
        p_geometry->SetData(rGeometry.GetData());
        return p_geometry;
    }
 
    typename BaseType::Pointer Create(
        const IndexType NewGeometryId,
        const BaseType& rGeometry
        ) const override
    {
        auto p_geometry = typename BaseType::Pointer( new Triangle2D3( NewGeometryId, rGeometry.Points() ) );
        p_geometry->SetData(rGeometry.GetData());
        return p_geometry;
    }
 
    Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
    {
        rResult = ZeroMatrix( 3, 2 );
        rResult( 0, 0 ) =  0.0;
        rResult( 0, 1 ) =  0.0;
        rResult( 1, 0 ) =  1.0;
        rResult( 1, 1 ) =  0.0;
        rResult( 2, 0 ) =  0.0;
        rResult( 2, 1 ) =  1.0;
        return rResult;
    }
 
    Vector& LumpingFactors(
        Vector& rResult,
        const typename BaseType::LumpingMethods LumpingMethod = BaseType::LumpingMethods::ROW_SUM
        )  const override
    {
        if(rResult.size() != 3)
            rResult.resize( 3, false );
        std::fill( rResult.begin(), rResult.end(), 1.00 / 3.00 );
        return rResult;
    }
 
 
    double Length() const override
    {
        //return sqrt(fabs( DeterminantOfJacobian(PointType()))*0.5);
        double length = 0.000;
        length = 1.1283791670955 * sqrt( fabs( Area() ) );  // (4xA/Pi)^(1/2)
        return length;
    }
 
 
    double Area() const override {
        const PointType& p0 = this->operator [](0);
        const PointType& p1 = this->operator [](1);
        const PointType& p2 = this->operator [](2);
 
        double x10 = p1.X() - p0.X();
        double y10 = p1.Y() - p0.Y();
 
        double x20 = p2.X() - p0.X();
        double y20 = p2.Y() - p0.Y();
 
        double detJ = x10 * y20-y10 * x20;
        return 0.5*detJ;
    }
 
    // TODO: Code activated in June 2023
    // /**
    //  * @brief This method calculates and returns the volume of this geometry.
    //  * @return Error, the volume of a 2D geometry is not defined
    //  * @see Length()
    //  * @see Area()
    //  * @see Volume()
    //  */
    // double Volume() const override
    // {
    //     KRATOS_ERROR << "Triangle2D3:: Method not well defined. Replace with DomainSize() instead" << std::endl;
    //     return 0.0;
    // }
 
    bool HasIntersection( const BaseType& rThisGeometry ) const override
    {
        if (rThisGeometry.LocalSpaceDimension() < this->LocalSpaceDimension()) {
            return IntersectionUtilities::TriangleLineIntersection2D(
                *this, rThisGeometry[0], rThisGeometry[1]);
        }  // Both geometries are 2D
        const BaseType& geom_1 = *this;
        const BaseType& geom_2 = rThisGeometry;
        return  NoDivTriTriIsect(geom_1[0], geom_1[1], geom_1[2], geom_2[0], geom_2[1], geom_2[2]);
    }
 
    bool HasIntersection( const Point& rLowPoint, const Point& rHighPoint ) const override
    {
        Point box_center;
        Point box_half_size;
 
        box_center[0] = 0.50 * (rLowPoint[0] + rHighPoint[0]);
        box_center[1] = 0.50 * (rLowPoint[1] + rHighPoint[1]);
        box_center[2] = 0.00;
 
        box_half_size[0] = 0.50 * std::abs(rHighPoint[0] - rLowPoint[0]);
        box_half_size[1] = 0.50 * std::abs(rHighPoint[1] - rLowPoint[1]);
        box_half_size[2] = 0.00;
 
        return TriBoxOverlap(box_center, box_half_size);
    }
 
    double DomainSize() const override {
      return this->Area();
    }
 
    virtual double Semiperimeter() const {
      return CalculateSemiperimeter(
        MathUtils<double>::Norm3(this->GetPoint(0)-this->GetPoint(1)),
        MathUtils<double>::Norm3(this->GetPoint(1)-this->GetPoint(2)),
        MathUtils<double>::Norm3(this->GetPoint(2)-this->GetPoint(0))
      );
    }
 
    double MinEdgeLength() const override {
      auto a = this->GetPoint(0) - this->GetPoint(1);
      auto b = this->GetPoint(1) - this->GetPoint(2);
      auto c = this->GetPoint(2) - this->GetPoint(0);
 
      auto sa = (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]);
      auto sb = (b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]);
      auto sc = (c[0]*c[0])+(c[1]*c[1])+(c[2]*c[2]);
 
      return CalculateMinEdgeLength(sa, sb, sc);
    }
 
    double MaxEdgeLength() const override {
      auto a = this->GetPoint(0) - this->GetPoint(1);
      auto b = this->GetPoint(1) - this->GetPoint(2);
      auto c = this->GetPoint(2) - this->GetPoint(0);
 
      auto sa = (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]);
      auto sb = (b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]);
      auto sc = (c[0]*c[0])+(c[1]*c[1])+(c[2]*c[2]);
 
      return CalculateMaxEdgeLength(sa, sb, sc);
    }
 
    double AverageEdgeLength() const override {
      return CalculateAvgEdgeLength(
        MathUtils<double>::Norm3(this->GetPoint(0)-this->GetPoint(1)),
        MathUtils<double>::Norm3(this->GetPoint(1)-this->GetPoint(2)),
        MathUtils<double>::Norm3(this->GetPoint(2)-this->GetPoint(0))
      );
    }
 
    double Circumradius() const override {
      return CalculateCircumradius(
        MathUtils<double>::Norm3(this->GetPoint(0)-this->GetPoint(1)),
        MathUtils<double>::Norm3(this->GetPoint(1)-this->GetPoint(2)),
        MathUtils<double>::Norm3(this->GetPoint(2)-this->GetPoint(0))
      );
    }
 
    double Inradius() const override {
      return CalculateInradius(
        MathUtils<double>::Norm3(this->GetPoint(0)-this->GetPoint(1)),
        MathUtils<double>::Norm3(this->GetPoint(1)-this->GetPoint(2)),
        MathUtils<double>::Norm3(this->GetPoint(2)-this->GetPoint(0))
      );
    }
 
 
    double InradiusToCircumradiusQuality() const override {
      constexpr double normFactor = 1.0;
 
      double a = MathUtils<double>::Norm3(this->GetPoint(0)-this->GetPoint(1));
      double b = MathUtils<double>::Norm3(this->GetPoint(1)-this->GetPoint(2));
      double c = MathUtils<double>::Norm3(this->GetPoint(2)-this->GetPoint(0));
 
      return normFactor * CalculateInradius(a,b,c) / CalculateCircumradius(a,b,c);
    };
 
    double InradiusToLongestEdgeQuality() const override {
      constexpr double normFactor = 1.0; // TODO: This normalization coeficient is not correct.
 
      auto a = this->GetPoint(0) - this->GetPoint(1);
      auto b = this->GetPoint(1) - this->GetPoint(2);
      auto c = this->GetPoint(2) - this->GetPoint(0);
 
      double sa = (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]);
      double sb = (b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]);
      double sc = (c[0]*c[0])+(c[1]*c[1])+(c[2]*c[2]);
 
      return normFactor * CalculateInradius(std::sqrt(sa),std::sqrt(sb),std::sqrt(sc)) / CalculateMaxEdgeLength(sa,sb,sc);
    }
 
    double AreaToEdgeLengthRatio() const override {
      constexpr double normFactor = 1.0;
 
      auto a = this->GetPoint(0) - this->GetPoint(1);
      auto b = this->GetPoint(1) - this->GetPoint(2);
      auto c = this->GetPoint(2) - this->GetPoint(0);
 
      double sa = (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]);
      double sb = (b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]);
      double sc = (c[0]*c[0])+(c[1]*c[1])+(c[2]*c[2]);
 
      return normFactor * Area() / (sa+sb+sc);
    }
 
    double ShortestAltitudeToEdgeLengthRatio() const override {
      constexpr double normFactor = 1.0;
 
      auto a = this->GetPoint(0) - this->GetPoint(1);
      auto b = this->GetPoint(1) - this->GetPoint(2);
      auto c = this->GetPoint(2) - this->GetPoint(0);
 
      double sa = (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]);
      double sb = (b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]);
      double sc = (c[0]*c[0])+(c[1]*c[1])+(c[2]*c[2]);
 
      // Shortest altitude is the one intersecting the largest base.
      double base = CalculateMaxEdgeLength(sa,sb,sc);
 
      return normFactor * (Area() * 2 / base ) / std::sqrt(sa+sb+sc);
    }
 
    array_1d<double, 3> Normal(const CoordinatesArrayType& rPointLocalCoordinates) const override
    {
        const array_1d<double, 3> tangent_xi  = this->GetPoint(1) - this->GetPoint(0);
        const array_1d<double, 3> tangent_eta = this->GetPoint(2) - this->GetPoint(0);
 
        array_1d<double, 3> normal;
        MathUtils<double>::CrossProduct(normal, tangent_xi, tangent_eta);
 
        return 0.5 * normal;
    }
 
    bool IsInside(
        const CoordinatesArrayType& rPoint,
        CoordinatesArrayType& rResult,
        const double Tolerance = std::numeric_limits<double>::epsilon()
        ) const override
    {
        this->PointLocalCoordinates( rResult, rPoint );
 
        if ( (rResult[0] >= (0.0-Tolerance)) && (rResult[0] <= (1.0+Tolerance)) )
        {
            if ( (rResult[1] >= (0.0-Tolerance)) && (rResult[1] <= (1.0+Tolerance)) )
            {
                if ( (rResult[0] + rResult[1]) <= (1.0+Tolerance) )
                {
                    return true;
                }
            }
        }
 
        return false;
    }
 
    CoordinatesArrayType& PointLocalCoordinates(
        CoordinatesArrayType& rResult,
        const CoordinatesArrayType& rPoint
        ) const override {
 
        rResult = ZeroVector(3);
 
        const TPointType& point_0 = this->GetPoint(0);
 
        // Compute the Jacobian matrix and its determinant
        BoundedMatrix<double, 2, 2> J;
        J(0,0) = this->GetPoint(1).X() - this->GetPoint(0).X();
        J(0,1) = this->GetPoint(2).X() - this->GetPoint(0).X();
        J(1,0) = this->GetPoint(1).Y() - this->GetPoint(0).Y();
        J(1,1) = this->GetPoint(2).Y() - this->GetPoint(0).Y();
        const double det_J = J(0,0)*J(1,1) - J(0,1)*J(1,0);
 
        // Compute eta and xi
        const double eta = (J(1,0)*(point_0.X() - rPoint(0)) +
                            J(0,0)*(rPoint(1) - point_0.Y())) / det_J;
        const double xi  = (J(1,1)*(rPoint(0) - point_0.X()) +
                            J(0,1)*(point_0.Y() - rPoint(1))) / det_J;
 
        rResult(0) = xi;
        rResult(1) = eta;
        rResult(2) = 0.0;
 
        return rResult;
    }
 
 
    SizeType EdgesNumber() const override
    {
        return 3;
    }
 
    GeometriesArrayType GenerateEdges() const override
    {
        GeometriesArrayType edges = GeometriesArrayType();
        edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 1 ), this->pGetPoint( 2 ) ) );
        edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 2 ), this->pGetPoint( 0 ) ) );
        edges.push_back( Kratos::make_shared<EdgeType>( this->pGetPoint( 0 ), this->pGetPoint( 1 ) ) );
        return edges;
    }
 
    SizeType FacesNumber() const override
    {
        return 3;
    }
 
    //Connectivities of faces required
    void NumberNodesInFaces (DenseVector<unsigned int>& NumberNodesInFaces) const override
    {
        if(NumberNodesInFaces.size() != 3 )
            NumberNodesInFaces.resize(3,false);
        // Linear Triangles have elements of 2 nodes as faces
        NumberNodesInFaces[0]=2;
        NumberNodesInFaces[1]=2;
        NumberNodesInFaces[2]=2;
 
    }
 
    void NodesInFaces (DenseMatrix<unsigned int>& NodesInFaces) const override
    {
        // faces in columns
        if(NodesInFaces.size1() != 3 || NodesInFaces.size2() != 3)
            NodesInFaces.resize(3,3,false);
 
        //face 1
        NodesInFaces(0,0)=0;//contrary node to the face
        NodesInFaces(1,0)=1;
        NodesInFaces(2,0)=2;
        //face 2
        NodesInFaces(0,1)=1;//contrary node to the face
        NodesInFaces(1,1)=2;
        NodesInFaces(2,1)=0;
        //face 3
        NodesInFaces(0,2)=2;//contrary node to the face
        NodesInFaces(1,2)=0;
        NodesInFaces(2,2)=1;
 
    }
 
 
 
 
    double ShapeFunctionValue( IndexType ShapeFunctionIndex,
                                       const CoordinatesArrayType& rPoint ) const override
    {
        switch ( ShapeFunctionIndex )
        {
        case 0:
            return( 1.0 -rPoint[0] - rPoint[1] );
        case 1:
            return( rPoint[0] );
        case 2:
            return( rPoint[1] );
        default:
            KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
        }
 
        return 0;
    }
 
    Vector& ShapeFunctionsValues (Vector &rResult, const CoordinatesArrayType& rCoordinates) const override
    {
      if(rResult.size() != 3) rResult.resize(3,false);
      rResult[0] =  1.0 -rCoordinates[0] - rCoordinates[1];
      rResult[1] =  rCoordinates[0] ;
      rResult[2] =  rCoordinates[1] ;
 
        return rResult;
    }
 
 
    void ShapeFunctionsIntegrationPointsGradients(
        ShapeFunctionsGradientsType& rResult,
        IntegrationMethod ThisMethod ) const override
    {
        const unsigned int integration_points_number =
            msGeometryData.IntegrationPointsNumber( ThisMethod );
 
        BoundedMatrix<double,3,2> DN_DX;
        double x10 = this->Points()[1].X() - this->Points()[0].X();
        double y10 = this->Points()[1].Y() - this->Points()[0].Y();
 
        double x20 = this->Points()[2].X() - this->Points()[0].X();
        double y20 = this->Points()[2].Y() - this->Points()[0].Y();
 
        //Jacobian is calculated:
        //  |dx/dxi  dx/deta|   |x1-x0   x2-x0|
        //J=|               |=  |             |
        //  |dy/dxi  dy/deta|   |y1-y0   y2-y0|
 
 
        double detJ = x10 * y20-y10 * x20;
 
        DN_DX(0,0) = -y20 + y10;
        DN_DX(0,1) =  x20 - x10;
        DN_DX(1,0) =  y20;
        DN_DX(1,1) = -x20;
        DN_DX(2,0) = -y10;
        DN_DX(2,1) =  x10;
 
        DN_DX /= detJ;
 
        //workaround by riccardo
        if(rResult.size() != integration_points_number)
        {
            rResult.resize(integration_points_number,false);
        }
        for(unsigned int i=0; i<integration_points_number; i++)
            rResult[i] = DN_DX;
    }
 
    void ShapeFunctionsIntegrationPointsGradients(
        ShapeFunctionsGradientsType& rResult,
        Vector& determinants_of_jacobian,
        IntegrationMethod ThisMethod ) const override
    {
        const unsigned int integration_points_number =
            msGeometryData.IntegrationPointsNumber( ThisMethod );
 
        BoundedMatrix<double,3,2> DN_DX;
        double x10 = this->Points()[1].X() - this->Points()[0].X();
        double y10 = this->Points()[1].Y() - this->Points()[0].Y();
 
        double x20 = this->Points()[2].X() - this->Points()[0].X();
        double y20 = this->Points()[2].Y() - this->Points()[0].Y();
 
        //Jacobian is calculated:
        //  |dx/dxi  dx/deta|   |x1-x0   x2-x0|
        //J=|               |=  |             |
        //  |dy/dxi  dy/deta|   |y1-y0   y2-y0|
 
 
        double detJ = x10 * y20-y10 * x20;
 
        DN_DX(0,0) = -y20 + y10;
        DN_DX(0,1) =  x20 - x10;
        DN_DX(1,0) =  y20;
        DN_DX(1,1) = -x20;
        DN_DX(2,0) = -y10;
        DN_DX(2,1) =  x10;
 
        DN_DX /= detJ;
 
        //workaround by riccardo
        if(rResult.size() != integration_points_number)
        {
            rResult.resize(integration_points_number,false);
        }
        for(unsigned int i=0; i<integration_points_number; i++)
            rResult[i] = DN_DX;
 
        if(determinants_of_jacobian.size() != integration_points_number )
            determinants_of_jacobian.resize(integration_points_number,false);
 
        for(unsigned int i=0; i<integration_points_number; i++)
            determinants_of_jacobian[i] = detJ;
    }
 
    Vector& DeterminantOfJacobian( Vector& rResult,
                                           IntegrationMethod ThisMethod ) const override
    {
        const unsigned int integration_points_number = msGeometryData.IntegrationPointsNumber( ThisMethod );
        if(rResult.size() != integration_points_number)
        {
            rResult.resize(integration_points_number,false);
        }
 
        const double detJ = 2.0*(this->Area());
 
        for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
        {
            rResult[pnt] = detJ;
        }
        return rResult;
    }
 
    double DeterminantOfJacobian( IndexType IntegrationPointIndex,
                                          IntegrationMethod ThisMethod ) const override
    {
        return 2.0*(this->Area());
    }
 
    double DeterminantOfJacobian( const CoordinatesArrayType& rPoint ) const override
    {
        return 2.0*(this->Area());
    }
 
 
 
    std::string Info() const override
    {
        return "2 dimensional triangle with three nodes in 2D space";
    }
 
    void PrintInfo( std::ostream& rOStream ) const override
    {
        rOStream << "2 dimensional triangle with three nodes in 2D space";
    }
 
    void PrintData( std::ostream& rOStream ) const override
    {
        // Base Geometry class PrintData call
        BaseType::PrintData( rOStream );
        std::cout << std::endl;
 
        // If the geometry has valid points, calculate and output its data
        if (this->AllPointsAreValid()) {
            Matrix jacobian;
            this->Jacobian( jacobian, PointType() );
            rOStream << "    Jacobian in the origin\t : " << jacobian;
        }
    }
 
    virtual ShapeFunctionsGradientsType ShapeFunctionsLocalGradients(
        IntegrationMethod ThisMethod )
    {
        ShapeFunctionsGradientsType localGradients
            = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
        const int integration_points_number
            = msGeometryData.IntegrationPointsNumber( ThisMethod );
        ShapeFunctionsGradientsType Result( integration_points_number );
 
        for ( int pnt = 0; pnt < integration_points_number; pnt++ )
        {
            Result[pnt] = localGradients[pnt];
        }
 
        return Result;
    }
 
    virtual ShapeFunctionsGradientsType ShapeFunctionsLocalGradients()
    {
        IntegrationMethod ThisMethod = msGeometryData.DefaultIntegrationMethod();
        ShapeFunctionsGradientsType localGradients
            = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
        const int integration_points_number
            = msGeometryData.IntegrationPointsNumber( ThisMethod );
        ShapeFunctionsGradientsType Result( integration_points_number );
 
        for ( int pnt = 0; pnt < integration_points_number; pnt++ )
        {
            Result[pnt] = localGradients[pnt];
        }
 
        return Result;
    }
 
    Matrix& ShapeFunctionsLocalGradients( Matrix& rResult,
            const CoordinatesArrayType& rPoint ) const override
    {
        rResult = ZeroMatrix( 3, 2 );
        rResult( 0, 0 ) = -1.0;
        rResult( 0, 1 ) = -1.0;
        rResult( 1, 0 ) =  1.0;
        rResult( 1, 1 ) =  0.0;
        rResult( 2, 0 ) =  0.0;
        rResult( 2, 1 ) =  1.0;
        return rResult;
    }
 
 
 
    virtual Matrix& ShapeFunctionsGradients( Matrix& rResult, PointType& rPoint )
    {
        rResult = ZeroMatrix( 3, 2 );
        rResult( 0, 0 ) = -1.0;
        rResult( 0, 1 ) = -1.0;
        rResult( 1, 0 ) =  1.0;
        rResult( 1, 1 ) =  0.0;
        rResult( 2, 0 ) =  0.0;
        rResult( 2, 1 ) =  1.0;
        return rResult;
    }
 
    ShapeFunctionsSecondDerivativesType& ShapeFunctionsSecondDerivatives( ShapeFunctionsSecondDerivativesType& rResult, const CoordinatesArrayType& rPoint ) const override
    {
        if ( rResult.size() != this->PointsNumber() )
        {
            // KLUDGE: While there is a bug in
            // ublas vector resize, I have to put this beside resizing!!
            ShapeFunctionsGradientsType temp( this->PointsNumber() );
            rResult.swap( temp );
        }
 
        if(rResult[0].size1() != 2 || rResult[0].size2() != 2 )
            rResult[0].resize( 2, 2,false );
 
        if(rResult[1].size1() != 2 || rResult[1].size2() != 2 )
            rResult[1].resize( 2, 2,false );
 
        if(rResult[2].size1() != 2 || rResult[2].size2() != 2 )
            rResult[2].resize( 2, 2,false );
 
        rResult[0]( 0, 0 ) = 0.0;
        rResult[0]( 0, 1 ) = 0.0;
        rResult[0]( 1, 0 ) = 0.0;
        rResult[0]( 1, 1 ) = 0.0;
        rResult[1]( 0, 0 ) = 0.0;
        rResult[1]( 0, 1 ) = 0.0;
        rResult[1]( 1, 0 ) = 0.0;
        rResult[1]( 1, 1 ) = 0.0;
        rResult[2]( 0, 0 ) = 0.0;
        rResult[2]( 0, 1 ) = 0.0;
        rResult[2]( 1, 0 ) = 0.0;
        rResult[2]( 1, 1 ) = 0.0;
        return rResult;
    }
 
    ShapeFunctionsThirdDerivativesType& ShapeFunctionsThirdDerivatives( ShapeFunctionsThirdDerivativesType& rResult, const CoordinatesArrayType& rPoint ) const override
    {
        if ( rResult.size() != this->PointsNumber() )
        {
            // KLUDGE: While there is a bug in
            // ublas vector resize, I have to put this beside resizing!!
//                 ShapeFunctionsGradientsType
            ShapeFunctionsThirdDerivativesType temp( this->PointsNumber() );
            rResult.swap( temp );
        }
 
        for ( IndexType i = 0; i < rResult.size(); i++ )
        {
            DenseVector<Matrix> temp( this->PointsNumber() );
            rResult[i].swap( temp );
        }
 
        rResult[0][0].resize( 2, 2,false );
 
        rResult[0][1].resize( 2, 2,false );
        rResult[1][0].resize( 2, 2,false );
        rResult[1][1].resize( 2, 2,false );
        rResult[2][0].resize( 2, 2,false );
        rResult[2][1].resize( 2, 2,false );
 
        for ( int i = 0; i < 3; i++ )
        {
 
            rResult[i][0]( 0, 0 ) = 0.0;
            rResult[i][0]( 0, 1 ) = 0.0;
            rResult[i][0]( 1, 0 ) = 0.0;
            rResult[i][0]( 1, 1 ) = 0.0;
            rResult[i][1]( 0, 0 ) = 0.0;
            rResult[i][1]( 0, 1 ) = 0.0;
            rResult[i][1]( 1, 0 ) = 0.0;
            rResult[i][1]( 1, 1 ) = 0.0;
        }
 
        return rResult;
    }
 
 
 
protected:
private:
 
    static const GeometryData msGeometryData;
 
    static const GeometryDimension msGeometryDimension;
 
 
    friend class Serializer;
 
    void save( Serializer& rSerializer ) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS( rSerializer, BaseType );
    }
 
    void load( Serializer& rSerializer ) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, BaseType );
    }
 
    Triangle2D3():BaseType( PointsArrayType(), &msGeometryData ) {}
 
 
 
 
 
 
    static Matrix CalculateShapeFunctionsIntegrationPointsValues(
        typename BaseType::IntegrationMethod ThisMethod )
    {
        IntegrationPointsContainerType all_integration_points =
            AllIntegrationPoints();
        IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
        //number of integration points
        const int integration_points_number = integration_points.size();
        //number of nodes in current geometry
        const int points_number = 3;
        //setting up return matrix
        Matrix shape_function_values( integration_points_number, points_number );
        //loop over all integration points
 
        for ( int pnt = 0; pnt < integration_points_number; pnt++ )
        {
            shape_function_values( pnt, 0 ) = 1.0
                                              - integration_points[pnt].X()
                                              - integration_points[pnt].Y();
            shape_function_values( pnt, 1 ) = integration_points[pnt].X();
            shape_function_values( pnt, 2 ) = integration_points[pnt].Y();
        }
 
        return shape_function_values;
    }
 
    static ShapeFunctionsGradientsType
    CalculateShapeFunctionsIntegrationPointsLocalGradients(
        typename BaseType::IntegrationMethod ThisMethod )
    {
        IntegrationPointsContainerType all_integration_points =
            AllIntegrationPoints();
        IntegrationPointsArrayType integration_points = all_integration_points[static_cast<int>(ThisMethod)];
        //number of integration points
        const int integration_points_number = integration_points.size();
        ShapeFunctionsGradientsType d_shape_f_values( integration_points_number );
        //initialising container
        //std::fill(d_shape_f_values.begin(), d_shape_f_values.end(), Matrix(4,2));
        //loop over all integration points
 
        for ( int pnt = 0; pnt < integration_points_number; pnt++ )
        {
            Matrix result( 3, 2 );
            result( 0, 0 ) = -1.0;
            result( 0, 1 ) = -1.0;
            result( 1, 0 ) =  1.0;
            result( 1, 1 ) =  0.0;
            result( 2, 0 ) =  0.0;
            result( 2, 1 ) =  1.0;
            d_shape_f_values[pnt] = result;
        }
 
        return d_shape_f_values;
    }
 
    static const IntegrationPointsContainerType AllIntegrationPoints()
    {
        IntegrationPointsContainerType integration_points =
        {
            {
                Quadrature<TriangleGaussLegendreIntegrationPoints1, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleGaussLegendreIntegrationPoints2, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleGaussLegendreIntegrationPoints3, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleGaussLegendreIntegrationPoints4, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleGaussLegendreIntegrationPoints5, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleCollocationIntegrationPoints1, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleCollocationIntegrationPoints2, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleCollocationIntegrationPoints3, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleCollocationIntegrationPoints4, 2, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<TriangleCollocationIntegrationPoints5, 2, IntegrationPoint<3> >::GenerateIntegrationPoints()
            }
        };
        return integration_points;
    }
 
    static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
    {
        ShapeFunctionsValuesContainerType shape_functions_values =
        {
            {
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_1 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_2 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_3 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_4 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_5 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_1 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_2 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_3 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_4 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_5 ),
            }
        };
        return shape_functions_values;
    }
 
    static const ShapeFunctionsLocalGradientsContainerType
    AllShapeFunctionsLocalGradients()
    {
        ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
        {
            {
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_3 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_4 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_5 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_1 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_2 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_3 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_4 ),
                Triangle2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_EXTENDED_GAUSS_5 ),
            }
        };
        return shape_functions_local_gradients;
    }
 
 
//*************************************************************************************
//*************************************************************************************
 
    /* Triangle/triangle intersection test routine,
    * by Tomas Moller, 1997.
    * See article "A Fast Triangle-Triangle Intersection Test",
    * Journal of Graphics Tools, 2(2), 1997
    *
    * Updated June 1999: removed the divisions -- a little faster now!
    * Updated October 1999: added {} to CROSS and SUB macros
    *
    * int NoDivTriTriIsect(float V0[3],float V1[3],float V2[3],
    *                      float U0[3],float U1[3],float U2[3])
    *
    * parameters: vertices of triangle 1: V0,V1,V2
    *             vertices of triangle 2: U0,U1,U2
    * result    : returns 1 if the triangles intersect, otherwise 0
    *
    */
 
    bool NoDivTriTriIsect( const Point& V0,
                           const Point& V1,
                           const Point& V2,
                           const Point& U0,
                           const Point& U1,
                           const Point& U2) const
    {
        short index;
        double d1,d2;
        double du0,du1,du2,dv0,dv1,dv2;
        double du0du1,du0du2,dv0dv1,dv0dv2;
        double vp0,vp1,vp2;
        double up0,up1,up2;
        double bb,cc,max;
        array_1d<double,2> isect1, isect2;
        array_1d<double,3> D;
        array_1d<double,3> E1,E2;
        array_1d<double,3> N1,N2;
 
 
        const double epsilon =  1E-6;
 
// compute plane equation of triangle(V0,V1,V2) //
        noalias(E1) = V1-V0;
        noalias(E2) = V2-V0;
        MathUtils<double>::CrossProduct(N1, E1, E2);
        d1=-inner_prod(N1,V0);
// plane equation 1: N1.X+d1=0 //
 
// put U0,U1,U2 into plane equation 1 to compute signed distances to the plane//
        du0=inner_prod(N1,U0)+d1;
        du1=inner_prod(N1,U1)+d1;
        du2=inner_prod(N1,U2)+d1;
 
// coplanarity robustness check //
        if(fabs(du0)<epsilon) du0=0.0;
        if(fabs(du1)<epsilon) du1=0.0;
        if(fabs(du2)<epsilon) du2=0.0;
 
        du0du1=du0*du1;
        du0du2=du0*du2;
 
        if(du0du1>0.00 && du0du2>0.00)// same sign on all of them + not equal 0 ? //
            return false;                   // no intersection occurs //
 
// compute plane of triangle (U0,U1,U2) //
        noalias(E1) = U1 - U0;
        noalias(E2) = U2 - U0;
        MathUtils<double>::CrossProduct(N2, E1, E2);
        d2=-inner_prod(N2,U0);
// plane equation 2: N2.X+d2=0 //
 
// put V0,V1,V2 into plane equation 2 //
        dv0=inner_prod(N2,V0)+d2;
        dv1=inner_prod(N2,V1)+d2;
        dv2=inner_prod(N2,V2)+d2;
 
        if(fabs(dv0)<epsilon) dv0=0.0;
        if(fabs(dv1)<epsilon) dv1=0.0;
        if(fabs(dv2)<epsilon) dv2=0.0;
 
        dv0dv1=dv0*dv1;
        dv0dv2=dv0*dv2;
 
        if(dv0dv1>0.00 && dv0dv2>0.00)// same sign on all of them + not equal 0 ? //
            return false;                   // no intersection occurs //
 
// compute direction of intersection line //
        MathUtils<double>::CrossProduct(D, N1, N2);
 
// compute and index to the largest component of D //
        max=(double)fabs(D[0]);
        index=0;
        bb=(double)fabs(D[1]);
        cc=(double)fabs(D[2]);
        if(bb>max)
        {
            max=bb,index=1;
        }
        if(cc>max)
        {
            max=cc,index=2;
        }
 
// this is the simplified projection onto L//
        vp0=V0[index];
        vp1=V1[index];
        vp2=V2[index];
 
        up0=U0[index];
        up1=U1[index];
        up2=U2[index];
 
 
// compute interval for triangle 1 //
        double a,b,c,x0,x1;
        if( New_Compute_Intervals(vp0,vp1,vp2,dv0, dv1,dv2,dv0dv1,dv0dv2, a,b,c,x0,x1)==true )
        {
            return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2);
        }
 
// compute interval for triangle 2 //
        double d,e,f,y0,y1;
        if( New_Compute_Intervals(up0, up1, up2, du0, du1, du2, du0du1, du0du2, d, e, f, y0, y1)==true )
        {
            return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2);
        }
 
 
        double xx,yy,xxyy,tmp;
        xx=x0*x1;
        yy=y0*y1;
        xxyy=xx*yy;
 
        tmp=a*xxyy;
        isect1[0]=tmp+b*x1*yy;
        isect1[1]=tmp+c*x0*yy;
 
        tmp=d*xxyy;
        isect2[0]=tmp+e*xx*y1;
        isect2[1]=tmp+f*xx*y0;
 
 
        Sort(isect1[0],isect1[1]);
        Sort(isect2[0],isect2[1]);
 
        if(isect1[1]<isect2[0] || isect2[1]<isect1[0]) return false;
        return true;
    }
 
 
//*************************************************************************************
//*************************************************************************************
 
 
// sort so that a<=b //
    void Sort(double& a, double& b) const
    {
        if(a>b)
        {
            double c;
            c=a;
            a=b;
            b=c;
        }
    }
 
//*************************************************************************************
//*************************************************************************************
 
    bool New_Compute_Intervals( double& VV0,
                                double& VV1,
                                double& VV2,
                                double& D0,
                                double& D1,
                                double& D2,
                                double& D0D1,
                                double& D0D2,
                                double& A,
                                double& B,
                                double& C,
                                double& X0,
                                double& X1
                              ) const
    {
        if(D0D1>0.00)
        {
            // here we know that D0D2<=0.0 //
            // that is D0, D1 are on the same side, D2 on the other or on the plane //
            A=VV2;
            B=(VV0-VV2)*D2;
            C=(VV1-VV2)*D2;
            X0=D2-D0;
            X1=D2-D1;
        }
        else if(D0D2>0.00)
        {
            // here we know that d0d1<=0.0 //
            A=VV1;
            B=(VV0-VV1)*D1;
            C=(VV2-VV1)*D1;
            X0=D1-D0;
            X1=D1-D2;
        }
        else if(D1*D2>0.00 || D0!=0.00)
        {
            // here we know that d0d1<=0.0 or that D0!=0.0 //
            A=VV0;
            B=(VV1-VV0)*D0;
            C=(VV2-VV0)*D0;
            X0=D0-D1;
            X1=D0-D2;
        }
        else if(D1!=0.00)
        {
            A=VV1;
            B=(VV0-VV1)*D1;
            C=(VV2-VV1)*D1;
            X0=D1-D0;
            X1=D1-D2;
        }
        else if(D2!=0.00)
        {
            A=VV2;
            B=(VV0-VV2)*D2;
            C=(VV1-VV2)*D2;
            X0=D2-D0;
            X1=D2-D1;
        }
        else
        {
            return true;
        }
 
        return false;
 
    }
 
//*************************************************************************************
//*************************************************************************************
 
    bool coplanar_tri_tri( const array_1d<double, 3>& N,
                           const Point& V0,
                           const Point& V1,
                           const Point& V2,
                           const Point& U0,
                           const Point& U1,
                           const Point& U2) const
    {
        array_1d<double, 3 > A;
        short i0,i1;
 
        // first project onto an axis-aligned plane, that maximizes the area //
        // of the triangles, compute indices: i0,i1. //
        A[0]=fabs(N[0]);
        A[1]=fabs(N[1]);
        A[2]=fabs(N[2]);
        if(A[0]>A[1])
        {
            if(A[0]>A[2])
            {
                i0=1;      // A[0] is greatest //
                i1=2;
            }
            else
            {
                i0=0;      // A[2] is greatest //
                i1=1;
            }
        }
        else   // A[0]<=A[1] //
        {
            if(A[2]>A[1])
            {
                i0=0;      // A[2] is greatest //
                i1=1;
            }
            else
            {
                i0=0;      // A[1] is greatest //
                i1=2;
            }
        }
 
        // test all edges of triangle 1 against the edges of triangle 2 //
        //std::cout<< "Proof One " << std::endl;
        if ( Edge_Against_Tri_Edges(i0, i1, V0,V1,U0,U1,U2)==true) return true;
 
        //std::cout<< "Proof Two " << std::endl;
        if ( Edge_Against_Tri_Edges(i0, i1, V1,V2,U0,U1,U2)==true) return true;
 
        //std::cout<< "Proof Three " << std::endl;
        if ( Edge_Against_Tri_Edges(i0, i1, V2,V0,U0,U1,U2)==true) return true;
 
        // finally, test if tri1 is totally contained in tri2 or vice versa //
        if (Point_In_Tri( i0, i1, V0, U0, U1, U2)==true) return true;
        if (Point_In_Tri( i0, i1, U0, V0, V1, V2)==true) return true;
 
        return false;
    }
 
//*************************************************************************************
//*************************************************************************************
 
    bool Edge_Against_Tri_Edges(const short& i0,
                                const short& i1,
                                const Point& V0,
                                const Point& V1,
                                const Point&U0,
                                const Point&U1,
                                const Point&U2) const
    {
 
        double Ax,Ay,Bx,By,Cx,Cy,e,d,f;
        Ax=V1[i0]-V0[i0];
        Ay=V1[i1]-V0[i1];
        // test edge U0,U1 against V0,V1 //
 
        //std::cout<< "Proof One B " << std::endl;
        if(Edge_Edge_Test(Ax, Ay, Bx, By, Cx, Cy, e, d, f, i0, i1, V0, U0, U1)==true) return true;
        // test edge U1,U2 against V0,V1 //
        //std::cout<< "Proof Two B " << std::endl;
        if(Edge_Edge_Test(Ax, Ay, Bx, By, Cx, Cy, e, d, f, i0, i1, V0, U1, U2)==true) return true;
        // test edge U2,U1 against V0,V1 //
        if(Edge_Edge_Test(Ax, Ay, Bx, By, Cx, Cy, e, d, f, i0, i1, V0, U2, U0)==true) return true;
 
        return false;
    }
 
 
//*************************************************************************************
//*************************************************************************************
 
//   this edge to edge test is based on Franlin Antonio's gem:
//   "Faster Line Segment Intersection", in Graphics Gems III,
//   pp. 199-202
    bool Edge_Edge_Test(const double& Ax,
                        const double& Ay,
                        double& Bx,
                        double& By,
                        double& Cx,
                        double& Cy,
                        double& e,
                        double& d,
                        double& f,
                        const short& i0,
                        const short& i1,
                        const Point&V0,
                        const Point&U0,
                        const Point&U1) const
    {
        Bx=U0[i0]-U1[i0];
        By=U0[i1]-U1[i1];
        Cx=V0[i0]-U0[i0];
        Cy=V0[i1]-U0[i1];
        f=Ay*Bx-Ax*By;
        d=By*Cx-Bx*Cy;
 
        if(std::fabs(f)<1E-10) f = 0.00;
        if(std::fabs(d)<1E-10) d = 0.00;
 
 
        if((f>0.00 && d>=0.00 && d<=f) || (f<0.00 && d<=0.00 && d>=f))
        {
            e=Ax*Cy-Ay*Cx;
            //std::cout<< "e =  "<< e << std::endl;
            if(f>0.00)
            {
                if(e>=0.00 && e<=f) return true;
            }
            else
            {
                if(e<=0.00 && e>=f) return true;
            }
        }
        return false;
    }
 
//*************************************************************************************
//*************************************************************************************
 
 
    bool Point_In_Tri(const short& i0,
                      const short& i1,
                      const Point& V0,
                      const Point& U0,
                      const Point& U1,
                      const Point& U2) const
    {
        double a,b,c,d0,d1,d2;
        // is T1 completly inside T2? //
        // check if V0 is inside tri(U0,U1,U2) //
        a=U1[i1]-U0[i1];
        b=-(U1[i0]-U0[i0]);
        c=-a*U0[i0]-b*U0[i1];
        d0=a*V0[i0]+b*V0[i1]+c;
 
        a=U2[i1]-U1[i1];
        b=-(U2[i0]-U1[i0]);
        c=-a*U1[i0]-b*U1[i1];
        d1=a*V0[i0]+b*V0[i1]+c;
 
        a=U0[i1]-U2[i1];
        b=-(U0[i0]-U2[i0]);
        c=-a*U2[i0]-b*U2[i1];
        d2=a*V0[i0]+b*V0[i1]+c;
        if(d0*d1>0.0)
        {
            if(d0*d2>0.0) return true;
        }
 
        return false;
    }
 
 
    inline bool TriBoxOverlap(Point& rBoxCenter, Point& rBoxHalfSize) const
    {
        double abs_ex, abs_ey;
        array_1d<double,3 > vert0, vert1, vert2;
        array_1d<double,3 > edge0, edge1, edge2;
        std::pair<double, double> min_max;
 
        // move everything so that the boxcenter is in (0,0,0)
        noalias(vert0) = this->GetPoint(0) - rBoxCenter;
        noalias(vert1) = this->GetPoint(1) - rBoxCenter;
        noalias(vert2) = this->GetPoint(2) - rBoxCenter;
 
        // compute triangle edges
        noalias(edge0) = vert1 - vert0;
        noalias(edge1) = vert2 - vert1;
        noalias(edge2) = vert0 - vert2;
 
        // Bullet 3:
        //  test the 9 tests first (this was faster)
        //  We are only interested on z-axis test, which defines the {x,y} plane
        //  Note that projections onto crossproduct tri-{x,y} are always 0:
        //    that means there is no separating axis on X,Y-axis tests
        abs_ex = std::abs(edge0[0]);
        abs_ey = std::abs(edge0[1]);
        if (AxisTestZ(edge0[0],edge0[1],abs_ex,abs_ey,vert0,vert2,rBoxHalfSize)) return false;
 
        abs_ex = std::abs(edge1[0]);
        abs_ey = std::abs(edge1[1]);
        if (AxisTestZ(edge1[0],edge1[1],abs_ex,abs_ey,vert1,vert0,rBoxHalfSize)) return false;
 
        abs_ex = std::abs(edge2[0]);
        abs_ey = std::abs(edge2[1]);
        if (AxisTestZ(edge2[0],edge2[1],abs_ex,abs_ey,vert2,vert1,rBoxHalfSize)) return false;
 
        // Bullet 1:
        //  first test overlap in the {x,y,(z)}-directions
        //  find min, max of the triangle each direction, and test for overlap in
        //  that direction -- this is equivalent to testing a minimal AABB around
        //  the triangle against the AABB
 
        // test in X-direction
        min_max = std::minmax({vert0[0],vert1[0],vert2[0]});
        if(min_max.first>rBoxHalfSize[0] || min_max.second<-rBoxHalfSize[0]) return false;
 
        // test in Y-direction
        min_max = std::minmax({vert0[1],vert1[1],vert2[1]});
        if(min_max.first>rBoxHalfSize[1] || min_max.second<-rBoxHalfSize[1]) return false;
 
        // test in Z-direction
        // we do not consider rBoxHalfSize[2] since we are working in 2D
 
        // Bullet 2:
        //  test if the box intersects the plane of the triangle
        //  compute plane equation of triangle: normal*x+d=0
        //  Triangle and box are on the {x-y} plane, so this test won't return false
 
        return true;   // box and triangle overlaps
    }
 
    bool AxisTestZ(double& rEdgeX, double& rEdgeY,
                   double& rAbsEdgeX, double& rAbsEdgeY,
                   array_1d<double,3>& rVertA,
                   array_1d<double,3>& rVertC,
                   Point& rBoxHalfSize) const
    {
        double proj_a, proj_c, rad;
        proj_a = rEdgeX*rVertA[1] - rEdgeY*rVertA[0];
        proj_c = rEdgeX*rVertC[1] - rEdgeY*rVertC[0];
        std::pair<double, double> min_max = std::minmax(proj_a, proj_c);
 
        rad = rAbsEdgeY*rBoxHalfSize[0] + rAbsEdgeX*rBoxHalfSize[1];
 
        if(min_max.first>rad || min_max.second<-rad) return true;
        else return false;
    }
 
    inline double CalculateSemiperimeter(const double a, const double b, const double c) const {
      return 0.5 * (a+b+c);
    }
 
    inline double CalculateMinEdgeLength(const double sa, const double sb, const double sc) const {
      return std::sqrt(std::min({sa, sb, sc}));
    }
 
    inline double CalculateMaxEdgeLength(const double sa, const double sb, const double sc) const {
      return std::sqrt(std::max({sa, sb, sc}));
    }
 
    inline double CalculateAvgEdgeLength(const double a, const double b, const double c) const {
      return 1.0 / 3.0 * (a+b+c);
    }
 
    inline double CalculateInradius(const double a, const double b, const double c) const {
      return 0.5 * std::sqrt((b+c-a) * (c+a-b) * (a+b-c) / (a+b+c));
    }
 
    inline double CalculateCircumradius(const double a, const double b, const double c) const {
      return (a*b*c) / std::sqrt((a+b+c) * (b+c-a) * (c+a-b) * (a+b-c));
    }
 
 
 
 
 
 
    template<class TOtherPointType> friend class Triangle2D3;
 
 
 
}; // Class Geometry
 
 
 
 
template<class TPointType> inline std::istream& operator >> (
    std::istream& rIStream,
    Triangle2D3<TPointType>& rThis );
template<class TPointType> inline std::ostream& operator << (
    std::ostream& rOStream,
    const Triangle2D3<TPointType>& rThis )
{
    rThis.PrintInfo( rOStream );
    rOStream << std::endl;
    rThis.PrintData( rOStream );
    return rOStream;
}
 
 
template<class TPointType> const
GeometryData Triangle2D3<TPointType>::msGeometryData(
    &msGeometryDimension,
    GeometryData::IntegrationMethod::GI_GAUSS_1,
    Triangle2D3<TPointType>::AllIntegrationPoints(),
    Triangle2D3<TPointType>::AllShapeFunctionsValues(),
    AllShapeFunctionsLocalGradients()
);
 
template<class TPointType> const
GeometryDimension Triangle2D3<TPointType>::msGeometryDimension(2, 2);
 
}// namespace Kratos.
 
#endif // KRATOS_TRIANGLE_2D_3_H_INCLUDED  defined