line_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
//
 
#pragma once
 
// System includes
 
// External includes
 
// Project includes
#include "geometries/geometry.h"
#include "integration/line_gauss_legendre_integration_points.h"
#include "utilities/integration_utilities.h"
 
namespace Kratos
{
 
 
 
 
 
 
template<class TPointType>
class Line2D3
  : public Geometry<TPointType>
{
public:
 
    typedef Geometry<TPointType> BaseType;
 
    KRATOS_CLASS_POINTER_DEFINITION( Line2D3 );
 
    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::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::NormalType NormalType;
 
    typedef typename BaseType::CoordinatesArrayType CoordinatesArrayType;
 
 
    Line2D3( const PointType& FirstPoint,
             const PointType& SecondPoint,
             const PointType& ThirdPoint )
        : BaseType( PointsArrayType(), &msGeometryData )
    {
        BaseType::Points().push_back( typename PointType::Pointer( new PointType( FirstPoint ) ) );
        BaseType::Points().push_back( typename PointType::Pointer( new PointType( SecondPoint ) ) );
        BaseType::Points().push_back( typename PointType::Pointer( new PointType( ThirdPoint ) ) );
    }
 
    Line2D3( typename PointType::Pointer pFirstPoint,
             typename PointType::Pointer pSecondPoint,
             typename PointType::Pointer pThirdPoint )
        : BaseType( PointsArrayType(), &msGeometryData )
    {
        BaseType::Points().push_back( pFirstPoint );
        BaseType::Points().push_back( pSecondPoint );
        BaseType::Points().push_back( pThirdPoint );
    }
 
    explicit Line2D3( const PointsArrayType& rThisPoints )
        : BaseType( rThisPoints, &msGeometryData )
    {
        KRATOS_ERROR_IF( BaseType::PointsNumber() != 3 ) << "Invalid points number. Expected 3, given " << BaseType::PointsNumber() << std::endl;
    }
 
    explicit Line2D3(
        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 Line2D3(
        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;
    }
 
    Line2D3( Line2D3 const& rOther )
        : BaseType( rOther )
    {
    }
 
    template<class TOtherPointType> explicit Line2D3( Line2D3<TOtherPointType> const& rOther )
        : BaseType( rOther )
    {
    }
 
    ~Line2D3() override {}
 
    GeometryData::KratosGeometryFamily GetGeometryFamily() const override
    {
        return GeometryData::KratosGeometryFamily::Kratos_Linear;
    }
 
    GeometryData::KratosGeometryType GetGeometryType() const override
    {
        return GeometryData::KratosGeometryType::Kratos_Line2D3;
    }
 
 
    Line2D3& operator=( const Line2D3& rOther )
    {
        BaseType::operator=( rOther );
        return *this;
    }
 
    template<class TOtherPointType>
    Line2D3& operator=( Line2D3<TOtherPointType> const & rOther )
    {
        BaseType::operator=( rOther );
        return *this;
    }
 
 
    typename BaseType::Pointer Create(
        const IndexType NewGeometryId,
        PointsArrayType const& rThisPoints
        ) const override
    {
        return typename BaseType::Pointer( new Line2D3( NewGeometryId, rThisPoints ) );
    }
 
    typename BaseType::Pointer Create(
        const IndexType NewGeometryId,
        const BaseType& rGeometry
        ) const override
    {
        auto p_geometry = typename BaseType::Pointer( new Line2D3( NewGeometryId, rGeometry.Points() ) );
        p_geometry->SetData(rGeometry.GetData());
        return p_geometry;
    }
 
    Vector& LumpingFactors(
        Vector& rResult,
        const typename BaseType::LumpingMethods LumpingMethod = BaseType::LumpingMethods::ROW_SUM
        )  const override
    {
    if(rResult.size() != 3)
           rResult.resize( 3, false );
        rResult[0] = 0.25;
        rResult[2] = 0.5;
        rResult[1] = 0.25;
        return rResult;
    }
 
 
    double Length() const override
    {
        const IntegrationMethod integration_method = IntegrationUtilities::GetIntegrationMethodForExactMassMatrixEvaluation(*this);
        return IntegrationUtilities::ComputeDomainSize(*this, integration_method);
    }
 
    double Area() const override
    {
        return Length();
    }
 
 
    double DomainSize() const override
    {
        return Length();
    }
 
    bool IsInside(
        const CoordinatesArrayType& rPoint,
        CoordinatesArrayType& rResult,
        const double Tolerance = std::numeric_limits<double>::epsilon()
        ) const override
    {
        PointLocalCoordinates( rResult, rPoint );
 
        if ( std::abs( rResult[0] ) <= (1.0 + Tolerance) ) {
            return true;
        }
 
        return false;
    }
 
    CoordinatesArrayType& PointLocalCoordinates(
        CoordinatesArrayType& rResult,
        const CoordinatesArrayType& rPoint
        ) const override
    {
        BoundedMatrix<double,3,3> X;
        BoundedMatrix<double,3,1> DN;
        for(IndexType i=0; i<this->size(); ++i) {
            const auto& r_node = this->GetPoint(i);
            X(0, i) = r_node.X();
            X(1, i) = r_node.Y();
            X(2, i) = r_node.Z();
        }
 
        static constexpr double MaxNormPointLocalCoordinates = 300.0;
        static constexpr std::size_t MaxIteratioNumberPointLocalCoordinates = 500;
        static constexpr double MaxTolerancePointLocalCoordinates = 1.0e-8;
 
        Matrix J = ZeroMatrix( 1, 1 );
        Matrix invJ = ZeroMatrix( 1, 1 );
 
        // Starting with xi = 0
        if (rResult.size() != 3)
            rResult.resize(3, false);
        rResult = ZeroVector( 3 );
        double delta_xi = 0.0;
        const array_1d<double, 3> zero_array = ZeroVector(3);
        array_1d<double, 3> current_global_coords;
 
        //Newton iteration:
        for ( IndexType k = 0; k < MaxIteratioNumberPointLocalCoordinates; k++ ) {
            noalias(current_global_coords) = zero_array;
            this->GlobalCoordinates( current_global_coords, rResult );
 
            noalias( current_global_coords ) = rPoint - current_global_coords;
 
            // Derivatives of shape functions
            Matrix shape_functions_gradients;
            shape_functions_gradients = ShapeFunctionsLocalGradients(shape_functions_gradients, rResult );
            noalias(DN) = prod(X, shape_functions_gradients);
 
            noalias(J) = prod(trans(DN), DN);
            const array_1d<double, 1> res = prod(trans(DN), current_global_coords);
 
            // The inverted jacobian matrix
            invJ(0, 0) = 1.0/J( 0, 0 );
 
            delta_xi = invJ(0, 0) * res[0];
 
            rResult[0] += delta_xi;
 
            if ( delta_xi > MaxNormPointLocalCoordinates ) {
                KRATOS_WARNING_IF("Line2D3", k > 0) << "detJ =\t" << J( 0, 0 ) << " DeltaX =\t" << delta_xi << " stopping calculation. Iteration:\t" << k << std::endl;
                break;
            }
 
            if ( delta_xi < MaxTolerancePointLocalCoordinates )
                break;
        }
 
        return rResult;
    }
 
 
    JacobiansType& Jacobian( JacobiansType& rResult, IntegrationMethod ThisMethod ) const override
    {
        // Getting derivatives of shape functions
        const ShapeFunctionsGradientsType shape_functions_gradients = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
 
        const std::size_t number_of_integration_points = this->IntegrationPointsNumber( ThisMethod );
        if (rResult.size() !=  number_of_integration_points) {
            JacobiansType temp( number_of_integration_points );
            rResult.swap( temp );
        }
 
        // Loop over all integration points
        for (std::size_t pnt = 0; pnt < number_of_integration_points; ++pnt) {
            // Initializing jacobian matrix
            noalias(rResult[pnt]) = ZeroMatrix( 2, 1 );
 
            // Loop over all nodes
            for (std::size_t i = 0; i < this->PointsNumber(); ++i) {
                const auto& r_node = this->GetPoint(i);
                rResult[pnt](0, 0) += r_node.X() * shape_functions_gradients[pnt](i, 0);
                rResult[pnt](1, 0) += r_node.Y() * shape_functions_gradients[pnt](i, 0);
            }
        } // End of loop over all integration points
 
        return rResult;
    }
 
    JacobiansType& Jacobian( JacobiansType& rResult, IntegrationMethod ThisMethod, Matrix& rDeltaPosition ) const override
    {
        // Getting derivatives of shape functions
        ShapeFunctionsGradientsType shape_functions_gradients =  CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
        const std::size_t number_of_integration_points = this->IntegrationPointsNumber( ThisMethod );
 
        // Getting values of shape functions
        Matrix shape_functions_values = CalculateShapeFunctionsIntegrationPointsValues( ThisMethod );
 
        if ( rResult.size() != number_of_integration_points ) {
            JacobiansType temp( number_of_integration_points );
            rResult.swap( temp );
        }
 
        // Loop over all integration points
        for (std::size_t pnt = 0; pnt < number_of_integration_points; ++pnt ) {
            // Initializing jacobian matrix
            noalias(rResult[pnt]) = ZeroMatrix( 2, 1 );
 
            // Loop over all nodes
            for (std::size_t i = 0; i < this->PointsNumber(); ++i ) {
                const auto& r_node = this->GetPoint(i);
                rResult[pnt](0, 0) += (r_node.X() - rDeltaPosition(i,0)) * shape_functions_gradients[pnt](i, 0);
                rResult[pnt](1, 0) += (r_node.Y() - rDeltaPosition(i,1)) * shape_functions_gradients[pnt](i, 0);
            }
        }// End of loop over all integration points
 
        return rResult;
    }
 
    Matrix& Jacobian( Matrix& rResult, IndexType IntegrationPointIndex, IntegrationMethod ThisMethod ) const override
    {
        // Setting up size of jacobian matrix
        rResult.resize(2, 1, false);
        noalias(rResult) = ZeroMatrix( 2, 1 );
 
        // Derivatives of shape functions
        ShapeFunctionsGradientsType shape_functions_gradients = CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
        Matrix shape_function_gradient_in_integration_point = shape_functions_gradients( IntegrationPointIndex );
 
        // Values of shape functions in integration points
        DenseVector<double> ShapeFunctionsValuesInIntegrationPoint = ZeroVector( 3 );
        ShapeFunctionsValuesInIntegrationPoint = row( CalculateShapeFunctionsIntegrationPointsValues( ThisMethod ), IntegrationPointIndex );
 
        // Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
        // Loop over all nodes
        for (std::size_t i = 0; i < this->PointsNumber(); ++i ) {
            const auto& r_node = this->GetPoint(i);
            rResult(0, 0) += r_node.X() * shape_function_gradient_in_integration_point(i, 0);
            rResult(1, 0) += r_node.Y() * shape_function_gradient_in_integration_point(i, 0);
        }
 
        return rResult;
    }
 
    Matrix& Jacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
    {
        // Setting up size of jacobian matrix
        rResult.resize( 2, 1, false );
        noalias(rResult) = ZeroMatrix( 2, 1 );
 
        // Derivatives of shape functions
        Matrix shape_functions_gradients;
        shape_functions_gradients = ShapeFunctionsLocalGradients( shape_functions_gradients, rPoint );
 
        // Elements of jacobian matrix (e.g. J(1,1) = dX1/dXi1)
        // Loop over all nodes
        for (std::size_t i = 0; i < this->PointsNumber(); ++i) {
            const auto& r_node = this->GetPoint(i);
            rResult(0, 0) += r_node.X() * shape_functions_gradients(i, 0);
            rResult(1, 0) += r_node.Y() * shape_functions_gradients(i, 0);
        }
 
        return rResult;
    }
 
    Vector& DeterminantOfJacobian( Vector& rResult, IntegrationMethod ThisMethod ) const override
    {
        const std::size_t number_of_integration_points = this->IntegrationPointsNumber( ThisMethod );
        if( rResult.size() != number_of_integration_points)
            rResult.resize( number_of_integration_points, false );
 
        Matrix J(2, 1);
        for (std::size_t pnt = 0; pnt < number_of_integration_points; ++pnt) {
            this->Jacobian( J, pnt, ThisMethod);
            rResult[pnt] = std::sqrt(std::pow(J(0,0), 2) + std::pow(J(1,0), 2));
        }
        return rResult;
    }
 
    double DeterminantOfJacobian( IndexType IntegrationPointIndex, IntegrationMethod ThisMethod ) const override
    {
        Matrix J(2, 1);
        this->Jacobian( J, IntegrationPointIndex, ThisMethod);
        return std::sqrt(std::pow(J(0,0), 2) + std::pow(J(1,0), 2));
    }
 
    double DeterminantOfJacobian( const CoordinatesArrayType& rPoint ) const override
    {
        Matrix J(2, 1);
        this->Jacobian( J, rPoint);
        return std::sqrt(std::pow(J(0,0), 2) + std::pow(J(1,0), 2));
    }
 
 
    SizeType EdgesNumber() const override
    {
        return 2;
    }
 
 
    SizeType FacesNumber() const override
    {
      return 0;
    }
 
 
     Vector& ShapeFunctionsValues(
        Vector& rResult,
        const CoordinatesArrayType& rCoordinates
        ) const override
    {
        if(rResult.size() != 3) {
            rResult.resize(3, false);
        }
 
        rResult[0] = 0.5 * (rCoordinates[0] - 1.0) * rCoordinates[0];
        rResult[1] = 0.5 * (rCoordinates[0] + 1.0) * rCoordinates[0];
        rResult[2] = 1.0 - rCoordinates[0] * rCoordinates[0];
 
        return rResult;
    }
 
    double ShapeFunctionValue(
        IndexType ShapeFunctionIndex,
        const CoordinatesArrayType& rPoint
        ) const override
    {
        switch ( ShapeFunctionIndex )
        {
        case 0:
            return( 0.5*( rPoint[0] - 1.0 )*rPoint[0] );
        case 1:
            return( 0.5*( rPoint[0] + 1.0 )*rPoint[0] );
        case 2:
            return( 1.0 -rPoint[0]*rPoint[0] );
        default:
            KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
        }
 
        return 0;
    }
 
 
    std::string Info() const override
    {
        return "1 dimensional line with 3 nodes in 2D space";
    }
 
    void PrintInfo( std::ostream& rOStream ) const override
    {
        rOStream << "1 dimensional line with 3 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\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
    {
        // Setting up result matrix
        if(rResult.size1() != 3 || rResult.size2() != 1) {
            rResult.resize( 3, 1, false );
        }
 
        noalias( rResult ) = ZeroMatrix( 3, 1 );
        rResult( 0, 0 ) =  rPoint[0] - 0.5;
        rResult( 1, 0 ) =  rPoint[0] + 0.5;
        rResult( 2, 0 ) = -rPoint[0] * 2.0;
        return( rResult );
    }
 
    Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
    {
        if(rResult.size1() != 3 || rResult.size2() != 1)
        {
            rResult.resize( 3, 1, false );
        }
        noalias( rResult ) = ZeroMatrix( 3, 1 );
        rResult( 0, 0 ) = -1.0;
        rResult( 1, 0 ) =  1.0;
        rResult( 2, 0 ) =  0.0;
        return rResult;
    }
 
    virtual Matrix& ShapeFunctionsGradients( Matrix& rResult, CoordinatesArrayType& rPoint )
    {
        if(rResult.size1() != 3 || rResult.size2() != 1)
        {
            rResult.resize( 3, 1, false );
        }
        noalias( rResult ) = ZeroMatrix( 3, 1 );
 
        rResult( 0, 0 ) =  rPoint[0] - 0.5;
        rResult( 1, 0 ) =  rPoint[0] + 0.5;
        rResult( 2, 0 ) = -rPoint[0] * 2.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 );
    }
 
    Line2D3(): BaseType( PointsArrayType(), &msGeometryData ) {}
 
 
 
 
    static Matrix CalculateShapeFunctionsIntegrationPointsValues( typename BaseType::IntegrationMethod ThisMethod )
    {
        const IntegrationPointsContainerType& all_integration_points = AllIntegrationPoints();
        const IntegrationPointsArrayType& IntegrationPoints = all_integration_points[static_cast<int>(ThisMethod)];
        int integration_points_number = IntegrationPoints.size();
        Matrix N( integration_points_number, 3 );
 
        for ( int it_gp = 0; it_gp < integration_points_number; it_gp++ )
        {
            double e = IntegrationPoints[it_gp].X();
            N( it_gp, 0 ) = 0.5 * ( e - 1 ) * e;
            N( it_gp, 2 ) = 1.0 - e * e;
            N( it_gp, 1 ) = 0.5 * ( 1 + e ) * e;
        }
 
        return N;
    }
 
    static ShapeFunctionsGradientsType CalculateShapeFunctionsIntegrationPointsLocalGradients(
        typename BaseType::IntegrationMethod ThisMethod )
    {
        const IntegrationPointsContainerType& all_integration_points = AllIntegrationPoints();
        const IntegrationPointsArrayType& IntegrationPoints = all_integration_points[static_cast<int>(ThisMethod)];
        ShapeFunctionsGradientsType DN_De( IntegrationPoints.size() );
        std::fill( DN_De.begin(), DN_De.end(), Matrix( 3, 1 ) );
 
        for ( unsigned int it_gp = 0; it_gp < IntegrationPoints.size(); it_gp++ )
        {
            double e = IntegrationPoints[it_gp].X();
            DN_De[it_gp]( 0, 0 ) = e - 0.5;
            DN_De[it_gp]( 2, 0 ) = -2.0 * e;
            DN_De[it_gp]( 1, 0 ) = e + 0.5;
        }
 
        return DN_De;
    }
 
    static const IntegrationPointsContainerType AllIntegrationPoints()
    {
        IntegrationPointsContainerType integration_points = {{
                Quadrature<LineGaussLegendreIntegrationPoints1, 1, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<LineGaussLegendreIntegrationPoints2, 1, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<LineGaussLegendreIntegrationPoints3, 1, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<LineGaussLegendreIntegrationPoints4, 1, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<LineGaussLegendreIntegrationPoints5, 1, IntegrationPoint<3> >::GenerateIntegrationPoints()
            }
        };
        return integration_points;
    }
 
    static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
    {
        ShapeFunctionsValuesContainerType shape_functions_values = {{
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_3 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_4 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsValues( GeometryData::IntegrationMethod::GI_GAUSS_5 )
            }
        };
        return shape_functions_values;
    }
 
    static const ShapeFunctionsLocalGradientsContainerType AllShapeFunctionsLocalGradients()
    {
        ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients = {{
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_3 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_4 ),
                Line2D3<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_5 )
            }
        };
        return shape_functions_local_gradients;
    }
 
 
 
 
 
 
    template<class TOtherPointType> friend class Line2D3;
 
 
 
 
 
}; // Class Geometry
 
 
 
 
 
 
template<class TPointType>
inline std::istream& operator >> ( std::istream& rIStream,
                                   Line2D3<TPointType>& rThis );
 
template<class TPointType>
inline std::ostream& operator << ( std::ostream& rOStream,
                                   const Line2D3<TPointType>& rThis )
{
    rThis.PrintInfo( rOStream );
    rOStream << std::endl;
    rThis.PrintData( rOStream );
 
    return rOStream;
}
 
 
template<class TPointType>
const GeometryData Line2D3<TPointType>::msGeometryData(
        &msGeometryDimension,
        GeometryData::IntegrationMethod::GI_GAUSS_2,
        Line2D3<TPointType>::AllIntegrationPoints(),
        Line2D3<TPointType>::AllShapeFunctionsValues(),
        AllShapeFunctionsLocalGradients() );
 
template<class TPointType>
const GeometryDimension Line2D3<TPointType>::msGeometryDimension(2, 1);
 
}  // namespace Kratos.