line_2d_5.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Mohamed Nabi
//
//
//  Contributors:
//
//
 
#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 Line2D5 : public Geometry<TPointType>
    {
    public:
 
        typedef Geometry<TPointType> BaseType;
 
        KRATOS_CLASS_POINTER_DEFINITION(Line2D5);
 
        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;
 
 
        Line2D5(const PointType& Point01, const PointType& Point02, const PointType& Point03,
                const PointType& Point04, const PointType& Point05)
                : BaseType(PointsArrayType(), &msGeometryData)
        {
            BaseType::Points().push_back(typename PointType::Pointer(new PointType(Point01)));
            BaseType::Points().push_back(typename PointType::Pointer(new PointType(Point02)));
            BaseType::Points().push_back(typename PointType::Pointer(new PointType(Point03)));
            BaseType::Points().push_back(typename PointType::Pointer(new PointType(Point04)));
            BaseType::Points().push_back(typename PointType::Pointer(new PointType(Point05)));
        }
 
        Line2D5(typename PointType::Pointer pPoint01, typename PointType::Pointer pPoint02,
                typename PointType::Pointer pPoint03, typename PointType::Pointer pPoint04,
                typename PointType::Pointer pPoint05) : BaseType(PointsArrayType(), &msGeometryData)
        {
            BaseType::Points().push_back(pPoint01);
            BaseType::Points().push_back(pPoint02);
            BaseType::Points().push_back(pPoint03);
            BaseType::Points().push_back(pPoint04);
            BaseType::Points().push_back(pPoint05);
        }
 
        explicit Line2D5(const PointsArrayType& rThisPoints) : BaseType(rThisPoints, &msGeometryData)
        {
            KRATOS_ERROR_IF(BaseType::PointsNumber() != 5) << "Invalid points number. Expected 5, given "
                << BaseType::PointsNumber() << std::endl;
        }
 
        explicit Line2D5(const IndexType GeometryId, const PointsArrayType& rThisPoints)
            : BaseType(GeometryId, rThisPoints, &msGeometryData)
        {
            KRATOS_ERROR_IF(this->PointsNumber() != 5) << "Invalid points number. Expected 5, given "
                << this->PointsNumber() << std::endl;
        }
 
        explicit Line2D5(const std::string& rGeometryName, const PointsArrayType& rThisPoints)
            : BaseType(rGeometryName, rThisPoints, &msGeometryData)
        {
            KRATOS_ERROR_IF(this->PointsNumber() != 5) << "Invalid points number. Expected 5, given "
                << this->PointsNumber() << std::endl;
        }
 
        Line2D5(Line2D5 const& rOther) : BaseType(rOther)
        {
        }
 
        template<class TOtherPointType> explicit Line2D5(Line2D5<TOtherPointType> const& rOther)
            : BaseType(rOther)
        {
        }
 
        ~Line2D5() override {}
 
        GeometryData::KratosGeometryFamily GetGeometryFamily() const override
        {
            return GeometryData::KratosGeometryFamily::Kratos_Linear;
        }
 
        GeometryData::KratosGeometryType GetGeometryType() const override
        {
            return GeometryData::KratosGeometryType::Kratos_Line2D5;
        }
 
 
        Line2D5& operator=(const Line2D5& rOther)
        {
            BaseType::operator=(rOther);
            return *this;
        }
 
        template<class TOtherPointType>
        Line2D5& operator=(Line2D5<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 Line2D5(NewGeometryId, rThisPoints));
        }
 
        typename BaseType::Pointer Create(const IndexType NewGeometryId, const BaseType& rGeometry) const override
        {
            auto p_geometry = typename BaseType::Pointer(new Line2D5(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() != 5) rResult.resize(5, false);
            rResult[0] =  7.0 / 90.0;
            rResult[1] =  7.0 / 90.0;
            rResult[2] = 32.0 / 90.0;
            rResult[3] = 12.0 / 90.0;
            rResult[4] = 32.0 / 90.0;
            return rResult;
        }
 
 
        double Length() const override
        {
            Vector temp;
            const IntegrationMethod integration_method = IntegrationUtilities::GetIntegrationMethodForExactMassMatrixEvaluation(*this);
            this->DeterminantOfJacobian(temp, integration_method);
            const IntegrationPointsArrayType& r_integration_points = this->IntegrationPoints(integration_method);
            double length = 0.0;
 
            for (std::size_t i = 0; i < r_integration_points.size(); ++i) {
                length += temp[i] * r_integration_points[i].Weight();
            }
 
            return length;
        }
 
        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, 5> 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);
            noalias(rResult) = ZeroVector(3);
            double delta_xi = 0.0;
            const array_1d<double, 3> zero_array = ZeroVector(3);
            array_1d<double, 3> current_global_coords;
            array_1d<double, 1> res;
 
            //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);
                noalias(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("Line2D5", 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 EdgesNumber();
        }
 
 
        Vector& ShapeFunctionsValues(Vector& rResult, const CoordinatesArrayType& rCoordinates) const override
        {
            if (rResult.size() != 5) rResult.resize(5, false);
            //
            const double xi = rCoordinates[0];
            const double fx1 = xi - 1.0;
            const double fx2 = xi + 1.0;
            const double fx3 = fx1 * fx2;
            const double gx1 = 2.0 * xi - 1.0;
            const double gx2 = 2.0 * xi + 1.0;
            const double gx3 = gx1 * gx2;
            const double cof1 = 1.0 / 6.0;
            const double cof2 = 8.0 * cof1;
            //
            rResult[0] = xi * fx1 * gx3 * cof1;
            rResult[1] = xi * fx2 * gx3 * cof1;
            rResult[2] = -xi * fx3 * gx1 * cof2;
            rResult[3] = fx3 * gx3;
            rResult[4] = -xi * fx3 * gx2 * cof2;
            //
            return rResult;
        }
 
        double ShapeFunctionValue(const IndexType ShapeFunctionIndex, const CoordinatesArrayType& rPoint) const override
        {
            const double xi = rPoint[0];
            const double fx1 = xi - 1.0;
            const double fx2 = xi + 1.0;
            const double fx3 = fx1 * fx2;
            const double gx1 = 2.0 * xi - 1.0;
            const double gx2 = 2.0 * xi + 1.0;
            const double gx3 = gx1 * gx2;
            const double cof1 = 1.0 / 6.0;
            const double cof2 = 8.0 * cof1;
            double shape = 0.0;
            //
            switch (ShapeFunctionIndex)
            {
            case 0:
                shape = xi * fx1 * gx3 * cof1;
                 break;
            case 1:
                shape = xi * fx2 * gx3 * cof1;
                break;
            case 2:
                shape = -xi * fx3 * gx1 * cof2;
                break;
            case 3:
                shape = fx3 * gx3;
                break;
            case 4:
                shape = -xi * fx3 * gx2 * cof2;
                break;
            default:
                KRATOS_ERROR << "Wrong index of shape function!" << *this << std::endl;
                break;
            }
            //
            return shape;
        }
 
 
        void ShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType& rResult,
            IntegrationMethod ThisMethod) const override
        {
            KRATOS_ERROR << "Jacobian is not square" << std::endl;
        }
 
        void ShapeFunctionsIntegrationPointsGradients(ShapeFunctionsGradientsType& rResult,
            Vector& rDeterminantsOfJacobian, IntegrationMethod ThisMethod) const override
        {
            KRATOS_ERROR << "Jacobian is not square" << std::endl;
        }
 
 
        std::string Info() const override
        {
            return "1 dimensional line with 5 nodes in 2D space";
        }
 
        void PrintInfo(std::ostream& rOStream) const override
        {
            rOStream << "1 dimensional line with 5 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() != 5 || rResult.size2() != 1) rResult.resize(5, 1, false);
            noalias(rResult) = ZeroMatrix(5, 1);
            //
            const double xi = rPoint[0];
            const double fx1 = xi - 1.0;
            const double fx2 = xi + 1.0;
            const double fx3 = fx1 * fx2;
            const double gx1 = 2.0 * xi - 1.0;
            const double gx2 = 2.0 * xi + 1.0;
            const double gx3 = gx1 * gx2;
            const double cof1 = 1.0 / 6.0;
            const double cof2 = 8.0 * cof1;
            //
            rResult(0, 0) = gx3 * gx1 * cof1 + xi * xi * fx1 * cof2;
            rResult(1, 0) = gx3 * gx2 * cof1 + xi * xi * fx2 * cof2;
            rResult(2, 0) = -(gx1 * (gx3 - fx3 - 1.0) + 2.0 * xi * fx3) * cof2;
            rResult(3, 0) = 2.0 * xi * (gx3 + 4.0 * fx3);
            rResult(4, 0) = -(gx2 * (gx3 - fx3 - 1.0) + 2.0 * xi * fx3) * cof2;
            return rResult;
        }
 
        Matrix& PointsLocalCoordinates(Matrix& rResult) const override
        {
            if (rResult.size1() != 5 || rResult.size2() != 1) rResult.resize(5, 1, false);
            noalias(rResult) = ZeroMatrix(5, 1);
            //
            rResult(0, 0) = -1.0;
            rResult(1, 0) = 1.0;
            rResult(2, 0) = -0.5;
            rResult(3, 0) = 0.0;
            rResult(4, 0) = 0.5;
            //
            return rResult;
        }
 
        virtual Matrix& ShapeFunctionsGradients(Matrix& rResult, const CoordinatesArrayType& rPoint)
        {
            if (rResult.size1() != 5 || rResult.size2() != 1) rResult.resize(5, 1, false);
            noalias(rResult) = ZeroMatrix(5, 1);
            //
            const double xi = rPoint[0];
            const double fx1 = xi - 1.0;
            const double fx2 = xi + 1.0;
            const double fx3 = fx1 * fx2;
            const double gx1 = 2.0 * xi - 1.0;
            const double gx2 = 2.0 * xi + 1.0;
            const double gx3 = gx1 * gx2;
            const double cof1 = 1.0 / 6.0;
            const double cof2 = 8.0 * cof1;
            //
            rResult(0, 0) = gx3 * gx1 * cof1 + xi * xi * fx1 * cof2;
            rResult(1, 0) = gx3 * gx2 * cof1 + xi * xi * fx2 * cof2;
            rResult(2, 0) = -(gx1 * (gx3 - fx3 - 1.0) + 2.0 * xi * fx3) * cof2;
            rResult(3, 0) = 2.0 * xi * (gx3 + 4.0 * fx3);
            rResult(4, 0) = -(gx2 * (gx3 - fx3 - 1.0) + 2.0 * xi * fx3) * cof2;
            //
            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);
        }
 
        Line2D5() : 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)];
            const int integration_points_number = IntegrationPoints.size();
            Matrix shape_function_values(integration_points_number, 5);
            //
            const double cof1 = 1.0 / 6.0;
            const double cof2 = 8.0 * cof1;
            //
            for (int pnt = 0; pnt < integration_points_number; ++pnt)
            {
                const double xi = IntegrationPoints[pnt].X();
                const double fx1 = xi - 1.0;
                const double fx2 = xi + 1.0;
                const double fx3 = fx1 * fx2;
                const double gx1 = 2.0 * xi - 1.0;
                const double gx2 = 2.0 * xi + 1.0;
                const double gx3 = gx1 * gx2;
                //
                shape_function_values(pnt, 0) = xi * fx1 * gx3 * cof1;
                shape_function_values(pnt, 1) = xi * fx2 * gx3 * cof1;
                shape_function_values(pnt, 2) = -xi * fx3 * gx1 * cof2;
                shape_function_values(pnt, 3) = fx3 * gx3;
                shape_function_values(pnt, 4) = -xi * fx3 * gx2 * cof2;
            }
            return shape_function_values;
        }
 
        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(5, 1));
            //
            const double cof1 = 1.0 / 6.0;
            const double cof2 = 8.0 * cof1;
            //
            for (unsigned int pnt = 0; pnt < IntegrationPoints.size(); ++pnt)
            {
                const double xi = IntegrationPoints[pnt].X();
                const double fx1 = xi - 1.0;
                const double fx2 = xi + 1.0;
                const double fx3 = fx1 * fx2;
                const double gx1 = 2.0 * xi - 1.0;
                const double gx2 = 2.0 * xi + 1.0;
                const double gx3 = gx1 * gx2;
                //
                DN_De[pnt](0, 0) = gx3 * gx1 * cof1 + xi * xi * fx1 * cof2;
                DN_De[pnt](1, 0) = gx3 * gx2 * cof1 + xi * xi * fx2 * cof2;
                DN_De[pnt](2, 0) = -(gx1 * (gx3 - fx3 - 1.0) + 2.0 * xi * fx3) * cof2;
                DN_De[pnt](3, 0) = 2.0 * xi * (gx3 + 4.0 * fx3);
                DN_De[pnt](4, 0) = -(gx2 * (gx3 - fx3 - 1.0) + 2.0 * xi * fx3) * cof2;
            }
            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 = {{
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryData::IntegrationMethod::GI_GAUSS_1),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryData::IntegrationMethod::GI_GAUSS_2),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryData::IntegrationMethod::GI_GAUSS_3),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryData::IntegrationMethod::GI_GAUSS_4),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(GeometryData::IntegrationMethod::GI_GAUSS_5)
                }
            };
            return shape_functions_values;
        }
 
        static const ShapeFunctionsLocalGradientsContainerType AllShapeFunctionsLocalGradients()
        {
            ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients = {{
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryData::IntegrationMethod::GI_GAUSS_1),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryData::IntegrationMethod::GI_GAUSS_2),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryData::IntegrationMethod::GI_GAUSS_3),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryData::IntegrationMethod::GI_GAUSS_4),
                    Line2D5<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients(GeometryData::IntegrationMethod::GI_GAUSS_5)
                }
            };
            return shape_functions_local_gradients;
        }
 
 
 
 
        template<class TOtherPointType> friend class Line2D5;
 
 
 
    }; // Class Geometry
 
 
 
 
    template<class TPointType>
    inline std::istream& operator >> (std::istream& rIStream, Line2D5<TPointType>& rThis);
 
    template<class TPointType>
    inline std::ostream& operator << (std::ostream& rOStream, const Line2D5<TPointType>& rThis)
    {
        rThis.PrintInfo(rOStream);
        rOStream << std::endl;
        rThis.PrintData(rOStream);
        return rOStream;
    }
 
 
    template<class TPointType>
    const GeometryData Line2D5<TPointType>::msGeometryData(
        &msGeometryDimension,
        GeometryData::IntegrationMethod::GI_GAUSS_4,
        Line2D5<TPointType>::AllIntegrationPoints(),
        Line2D5<TPointType>::AllShapeFunctionsValues(),
        AllShapeFunctionsLocalGradients());
 
    template<class TPointType>
    const GeometryDimension Line2D5<TPointType>::msGeometryDimension(2, 1);
 
}  // namespace Kratos.