prism_interface_3d_6.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Ignasi de Pouplana
//
 
#if !defined(KRATOS_PRISM_INTERFACE_3D_6_H_INCLUDED )
#define  KRATOS_PRISM_INTERFACE_3D_6_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "geometries/triangle_3d_3.h"
#include "geometries/quadrilateral_3d_4.h"
#include "integration/prism_gauss_lobatto_integration_points.h"
 
namespace Kratos
{
 
 
 
 
 
template<class TPointType> class PrismInterface3D6
    : public Geometry<TPointType>
{
public:
 
    typedef Geometry<TPointType> BaseType;
 
    typedef Line3D2<TPointType> EdgeType;
    typedef Triangle3D3<TPointType> FaceType1;
    typedef Quadrilateral3D4<TPointType> FaceType2;
 
    KRATOS_CLASS_POINTER_DEFINITION( PrismInterface3D6 );
 
    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;
 
 
//     PrismInterface3D6( const PointType& Point1, const PointType& Point2,
//               const PointType& Point3, const PointType& Point4,
//               const PointType& Point5, const PointType& Point6 )
//         : BaseType( PointsArrayType(), &msGeometryData )
//     {
//         this->Points().reserve( 6 );
//         this->Points().push_back( typename PointType::Pointer( new PointType( Point1 ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( Point2 ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( Point3 ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( Point4 ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( Point5 ) ) );
//         this->Points().push_back( typename PointType::Pointer( new PointType( Point6 ) ) );
//     }
 
    PrismInterface3D6( typename PointType::Pointer pPoint1,
              typename PointType::Pointer pPoint2,
              typename PointType::Pointer pPoint3,
              typename PointType::Pointer pPoint4,
              typename PointType::Pointer pPoint5,
              typename PointType::Pointer pPoint6 )
        : BaseType( PointsArrayType(), &msGeometryData )
    {
        this->Points().reserve( 6 );
        this->Points().push_back( pPoint1 );
        this->Points().push_back( pPoint2 );
        this->Points().push_back( pPoint3 );
        this->Points().push_back( pPoint4 );
        this->Points().push_back( pPoint5 );
        this->Points().push_back( pPoint6 );
    }
 
    PrismInterface3D6( const PointsArrayType& ThisPoints )
        : BaseType( ThisPoints, &msGeometryData )
    {
        if ( this->PointsNumber() != 6 )
            KRATOS_ERROR << "Invalid points number. Expected 6, given " << this->PointsNumber() << std::endl;
    }
 
    PrismInterface3D6( PrismInterface3D6 const& rOther )
        : BaseType( rOther )
    {
    }
 
    template<class TOtherPointType> PrismInterface3D6( PrismInterface3D6<TOtherPointType> const& rOther )
        : BaseType( rOther )
    {
    }
 
    virtual ~PrismInterface3D6() {}
 
    GeometryData::KratosGeometryFamily GetGeometryFamily() const override
    {
        return GeometryData::KratosGeometryFamily::Kratos_Prism;
    }
 
    GeometryData::KratosGeometryType GetGeometryType() const override
    {
        return GeometryData::KratosGeometryType::Kratos_Prism3D6;
    }
 
 
    PrismInterface3D6& operator=( const PrismInterface3D6& rOther )
    {
        BaseType::operator=( rOther );
        return *this;
    }
 
    template<class TOtherPointType>
    PrismInterface3D6& operator=( PrismInterface3D6<TOtherPointType> const & rOther )
    {
        BaseType::operator=( rOther );
        return *this;
    }
 
 
    typename BaseType::Pointer Create( PointsArrayType const& ThisPoints ) const override
    {
        return typename BaseType::Pointer( new PrismInterface3D6( ThisPoints ) );
    }
 
    Matrix& PointsLocalCoordinates( Matrix& rResult ) const override
    {
        if ( rResult.size1() != 6 || rResult.size2() != 3 )
            rResult.resize( 6, 3, false );
 
        rResult( 0, 0 ) = 0.0;
        rResult( 0, 1 ) = 0.0;
        rResult( 0, 2 ) = 0.0;
 
        rResult( 1, 0 ) = 1.0;
        rResult( 1, 1 ) = 0.0;
        rResult( 1, 2 ) = 0.0;
 
        rResult( 2, 0 ) = 0.0;
        rResult( 2, 1 ) = 1.0;
        rResult( 2, 2 ) = 0.0;
 
        rResult( 3, 0 ) = 0.0;
        rResult( 3, 1 ) = 0.0;
        rResult( 3, 2 ) = 1.0;
 
        rResult( 4, 0 ) = 1.0;
        rResult( 4, 1 ) = 0.0;
        rResult( 4, 2 ) = 1.0;
 
        rResult( 5, 0 ) = 0.0;
        rResult( 5, 1 ) = 1.0;
        rResult( 5, 2 ) = 1.0;
 
        return rResult;
    }
 
 
    double Length() const override
    {
        return std::sqrt(2.0*Area());
    }
 
    double Area() const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        //heron formula
        Vector side_a = ( p0 - p1 );
        double a = MathUtils<double>::Norm3( side_a );
        Vector side_b = ( p1 - p2 );
        double b = MathUtils<double>::Norm3( side_b );
        Vector side_c = ( p2 - p0 );
        double c = MathUtils<double>::Norm3( side_c );
        double s = ( a + b + c ) / 2;
        return( sqrt( s*( s - a )*( s - b )*( s - c ) ) );
    }
 
    double Volume() const override
    {
        return Area();
    }
 
    double DomainSize() const override
    {
        return Area();
    }
 
    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) )
            if ( rResult[1] >= (0.0 - Tolerance) )
                if ( rResult[2] >= (0.0 - Tolerance) )
                    if( rResult[2] <= (1.0 + Tolerance) )
                        if ( (rResult[0] + rResult[1]) <= (1.0 + Tolerance) )
                            return true;
 
        return false;
    }
 
    CoordinatesArrayType& PointLocalCoordinates(
        CoordinatesArrayType& rResult,
        const CoordinatesArrayType& rPoint
        ) const override
    {
        BoundedMatrix<double,3,3> X;
        BoundedMatrix<double,3,2> DN;
        for(unsigned int i=0; i<3;i++) {
            X(0,i ) = 0.5*(this->GetPoint( i ).X() + this->GetPoint( i+3 ).X());
            X(1,i ) = 0.5*(this->GetPoint( i ).Y() + this->GetPoint( i+3 ).Y());
            X(2,i ) = 0.5*(this->GetPoint( i ).Z() + this->GetPoint( i+3 ).Z());
        }
 
        double tol = 1.0e-8;
        int maxiter = 1000;
 
        Matrix J = ZeroMatrix( 2, 2 );
        Matrix invJ = ZeroMatrix( 2, 2 );
 
        //starting with xi = 0
        rResult = ZeroVector( 3 );
        Vector DeltaXi = ZeroVector( 2 );
        array_1d<double,3> CurrentGlobalCoords;
 
 
        //Newton iteration:
        for ( int k = 0; k < maxiter; k++ )
        {
            noalias(CurrentGlobalCoords) = ZeroVector( 3 );
            this->GlobalCoordinates( CurrentGlobalCoords, rResult );
 
            noalias( CurrentGlobalCoords ) = rPoint - CurrentGlobalCoords;
 
 
            //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);
            Vector res = prod(trans(DN),CurrentGlobalCoords);
 
            //deteminant of Jacobian
            const double det_j = J( 0, 0 ) * J( 1, 1 ) - J( 0, 1 ) * J( 1, 0 );
 
            //filling matrix
            invJ( 0, 0 ) = ( J( 1, 1 ) ) / ( det_j );
            invJ( 1, 0 ) = -( J( 1, 0 ) ) / ( det_j );
            invJ( 0, 1 ) = -( J( 0, 1 ) ) / ( det_j );
            invJ( 1, 1 ) = ( J( 0, 0 ) ) / ( det_j );
 
 
            DeltaXi( 0 ) = invJ( 0, 0 ) * res[0] + invJ( 0, 1 ) * res[1];
            DeltaXi( 1 ) = invJ( 1, 0 ) * res[0] + invJ( 1, 1 ) * res[1];
 
            rResult[0] += DeltaXi[0];
            rResult[1] += DeltaXi[1];
            rResult[2] = 0.0;
 
            if ( k>0 && norm_2( DeltaXi ) > 30 )
            {
                KRATOS_ERROR << "Computation of local coordinates failed at iteration " << k<< std::endl;
            }
 
            if ( norm_2( DeltaXi ) < tol )
            {
                break;
            }
        }
 
        return( rResult );
    }
 
    Matrix& Jacobian( Matrix& rResult,
                      IndexType IntegrationPointIndex,
                      IntegrationMethod ThisMethod ) const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        rResult.resize( 3, 2, false );
        rResult( 0, 0 ) = -( p0[0] ) + ( p1[0] ); //on the Gauss points (J is constant at each element)
        rResult( 1, 0 ) = -( p0[1] ) + ( p1[1] );
        rResult( 2, 0 ) = -( p0[2] ) + ( p1[2] );
        rResult( 0, 1 ) = -( p0[0] ) + ( p2[0] );
        rResult( 1, 1 ) = -( p0[1] ) + ( p2[1] );
        rResult( 2, 1 ) = -( p0[2] ) + ( p2[2] );
        return rResult;
    }
 
    Matrix& Jacobian( Matrix& rResult,
                      IndexType IntegrationPointIndex,
                      IntegrationMethod ThisMethod,
                      const Matrix& DeltaPosition ) const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        Matrix deltaPosMid(3, 3);
        for(unsigned int i = 0; i < 3; i++)
            for(unsigned int j = 0; j < 3; j++)
                deltaPosMid(i, j) = 0.5*( DeltaPosition(i, j) + DeltaPosition(i+3, j) );
 
        if(rResult.size1() != 3 || rResult.size2() != 2)
            rResult.resize(3, 2, false);
 
        rResult( 0, 0 ) = -( p0[0] - deltaPosMid(0,0) ) + ( p1[0] - deltaPosMid(1,0) ); //on the Gauss points (J is constant at each element)
        rResult( 1, 0 ) = -( p0[1] - deltaPosMid(0,1) ) + ( p1[1] - deltaPosMid(1,1) );
        rResult( 2, 0 ) = -( p0[2] - deltaPosMid(0,2) ) + ( p1[2] - deltaPosMid(1,2) );
        rResult( 0, 1 ) = -( p0[0] - deltaPosMid(0,0) ) + ( p2[0] - deltaPosMid(2,0) );
        rResult( 1, 1 ) = -( p0[1] - deltaPosMid(0,1) ) + ( p2[1] - deltaPosMid(2,1) );
        rResult( 2, 1 ) = -( p0[2] - deltaPosMid(0,2) ) + ( p2[2] - deltaPosMid(2,2) );
 
        return rResult;
    }
 
    Matrix& Jacobian( Matrix& rResult, const CoordinatesArrayType& rPoint ) const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        rResult.resize( 3, 2 , false);
        rResult( 0, 0 ) = -( p0[0] ) + ( p1[0] ); //on the Gauss points (J is constant at each element)
        rResult( 1, 0 ) = -( p0[1] ) + ( p1[1] );
        rResult( 2, 0 ) = -( p0[2] ) + ( p1[2] );
        rResult( 0, 1 ) = -( p0[0] ) + ( p2[0] );
        rResult( 1, 1 ) = -( p0[1] ) + ( p2[1] );
        rResult( 2, 1 ) = -( p0[2] ) + ( p2[2] );
        return rResult;
    }
 
    Vector& DeterminantOfJacobian( Vector& rResult,
                                   IntegrationMethod ThisMethod ) const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        array_1d<double, 3> Tangent1( p1 - p0 );
        array_1d<double, 3> Tangent2( p2 - p0 );
 
        array_1d<double, 3> Normal;
        MathUtils<double>::CrossProduct( Normal, Tangent1, Tangent2);
 
        double detJ = norm_2(Normal);
 
        //for all integration points
        if ( rResult.size() != this->IntegrationPointsNumber( ThisMethod ) )
            rResult.resize( this->IntegrationPointsNumber( ThisMethod ), false );
 
        for ( unsigned int pnt = 0; pnt < this->IntegrationPointsNumber( ThisMethod ); pnt++ )
            rResult[pnt] = detJ;
 
        return rResult;
    }
 
    double DeterminantOfJacobian( IndexType IntegrationPointIndex,
                                  IntegrationMethod ThisMethod ) const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        array_1d<double, 3> Tangent1( p1 - p0 );
        array_1d<double, 3> Tangent2( p2 - p0 );
 
        array_1d<double, 3> Normal;
        MathUtils<double>::CrossProduct( Normal, Tangent1, Tangent2);
 
        return norm_2(Normal);
    }
 
    double DeterminantOfJacobian( const CoordinatesArrayType& rPoint ) const override
    {
        array_1d<double, 3> p0 = 0.5 * (BaseType::GetPoint( 0 ) + BaseType::GetPoint( 3 ));
        array_1d<double, 3> p1 = 0.5 * (BaseType::GetPoint( 1 ) + BaseType::GetPoint( 4 ));
        array_1d<double, 3> p2 = 0.5 * (BaseType::GetPoint( 2 ) + BaseType::GetPoint( 5 ));
 
        array_1d<double, 3> Tangent1( p1 - p0 );
        array_1d<double, 3> Tangent2( p2 - p0 );
 
        array_1d<double, 3> Normal;
        MathUtils<double>::CrossProduct( Normal, Tangent1, Tangent2);
 
        return norm_2(Normal);
    }
 
    JacobiansType& InverseOfJacobian( JacobiansType& rResult,
                                      IntegrationMethod ThisMethod ) const override
    {
        KRATOS_ERROR << "Jacobian is not square" << std::endl;;
        return rResult;
    }
 
    Matrix& InverseOfJacobian( Matrix& rResult,
                               IndexType IntegrationPointIndex,
                               IntegrationMethod ThisMethod ) const override
    {
        KRATOS_ERROR << "Jacobian is not square" << std::endl;
        return rResult;
    }
 
    Matrix& InverseOfJacobian( Matrix& rResult,
                               const CoordinatesArrayType& rPoint ) const override
    {
        KRATOS_ERROR << "Jacobian is not square" << std::endl;
        return rResult;
    }
 
 
    SizeType EdgesNumber() const override
    {
        return 9;
    }
 
    GeometriesArrayType GenerateEdges() const override
    {
        GeometriesArrayType edges = GeometriesArrayType();
        typedef typename Geometry<TPointType>::Pointer EdgePointerType;
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 0 ),
                                              this->pGetPoint( 1 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 1 ),
                                              this->pGetPoint( 2 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 2 ),
                                              this->pGetPoint( 0 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 3 ),
                                              this->pGetPoint( 4 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 4 ),
                                              this->pGetPoint( 5 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 5 ),
                                              this->pGetPoint( 3 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 0 ),
                                              this->pGetPoint( 3 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 1 ),
                                              this->pGetPoint( 4 ) ) ) );
        edges.push_back( EdgePointerType( new EdgeType(
                                              this->pGetPoint( 2 ),
                                              this->pGetPoint( 5 ) ) ) );
        return edges;
    }
 
 
    SizeType FacesNumber() const override
    {
        return 5;
    }
 
    GeometriesArrayType GenerateFaces() const override
    {
        GeometriesArrayType faces = GeometriesArrayType();
        typedef typename Geometry<TPointType>::Pointer FacePointerType;
        faces.push_back( FacePointerType( new FaceType1(
                                              this->pGetPoint( 0 ),
                                              this->pGetPoint( 2 ),
                                              this->pGetPoint( 1 ) ) ) );
        faces.push_back( FacePointerType( new FaceType1(
                                              this->pGetPoint( 3 ),
                                              this->pGetPoint( 4 ),
                                              this->pGetPoint( 5 ) ) ) );
        faces.push_back( FacePointerType( new FaceType2(
                                              this->pGetPoint( 1 ),
                                              this->pGetPoint( 2 ),
                                              this->pGetPoint( 5 ),
                                              this->pGetPoint( 4 ) ) ) );
        faces.push_back( FacePointerType( new FaceType2(
                                              this->pGetPoint( 0 ),
                                              this->pGetPoint( 3 ),
                                              this->pGetPoint( 5 ),
                                              this->pGetPoint( 2 ) ) ) );
        faces.push_back( FacePointerType( new FaceType2(
                                              this->pGetPoint( 0 ),
                                              this->pGetPoint( 1 ),
                                              this->pGetPoint( 4 ),
                                              this->pGetPoint( 3 ) ) ) );
        return faces;
    }
 
 
    double ShapeFunctionValue( IndexType ShapeFunctionIndex,
                               const CoordinatesArrayType& rPoint ) const override
    {
        switch ( ShapeFunctionIndex )
        {
        case 0:
            return( 1.0 -( rPoint[0] + rPoint[1] + rPoint[2] - ( rPoint[0]*rPoint[2] ) - ( rPoint[1]*rPoint[2] ) ) );
        case 1:
            return( rPoint[0] - ( rPoint[0]*rPoint[2] ) );
        case 2:
            return( rPoint[1] - ( rPoint[1]*rPoint[2] ) );
        case 3:
            return( rPoint[2] - ( rPoint[0]*rPoint[2] ) - ( rPoint[1]*rPoint[2] ) );
        case 4:
            return( rPoint[0]*rPoint[2] );
        case 5:
            return( rPoint[1]*rPoint[2] );
        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() != 6) rResult.resize(6,false);
      rResult[0] =  1.0 -( rCoordinates[0] + rCoordinates[1] + rCoordinates[2] - ( rCoordinates[0]*rCoordinates[2] ) - ( rCoordinates[1]*rCoordinates[2] ) );
      rResult[1] =  rCoordinates[0] - ( rCoordinates[0]*rCoordinates[2] );
      rResult[2] =  rCoordinates[1] - ( rCoordinates[1]*rCoordinates[2] );
      rResult[3] =  rCoordinates[2] - ( rCoordinates[0]*rCoordinates[2] ) - ( rCoordinates[1]*rCoordinates[2] );
      rResult[4] =  rCoordinates[0]*rCoordinates[2];
      rResult[5] =  rCoordinates[1]*rCoordinates[2];
 
        return rResult;
    }
 
    void ShapeFunctionsIntegrationPointsGradients(
        ShapeFunctionsGradientsType& rResult,
        IntegrationMethod ThisMethod ) const override
    {
        const unsigned int integration_points_number =
            msGeometryData.IntegrationPointsNumber( ThisMethod );
 
        if ( integration_points_number == 0 )
            KRATOS_ERROR << "This integration method is not supported" << *this << std::endl;
 
        //workaround by riccardo
        if ( rResult.size() != integration_points_number )
        {
            // KLUDGE: While there is a bug in ublas
            // vector resize, I have to put this beside resizing!!
            ShapeFunctionsGradientsType temp( integration_points_number );
            rResult.swap( temp );
        }
 
        //calculating the local gradients
        ShapeFunctionsGradientsType locG =
            CalculateShapeFunctionsIntegrationPointsLocalGradients( ThisMethod );
 
        //getting the inverse jacobian matrices
        JacobiansType temp( integration_points_number );
 
        JacobiansType invJ = InverseOfJacobian( temp, ThisMethod );
 
        //loop over all integration points
        for ( unsigned int pnt = 0; pnt < integration_points_number; pnt++ )
        {
            rResult[pnt].resize( 6, 3, false );
 
            for ( int i = 0; i < 6; i++ )
            {
                for ( int j = 0; j < 3; j++ )
                {
                    rResult[pnt]( i, j ) =
                        ( locG[pnt]( i, 0 ) * invJ[pnt]( j, 0 ) )
                        + ( locG[pnt]( i, 1 ) * invJ[pnt]( j, 1 ) )
                        + ( locG[pnt]( i, 2 ) * invJ[pnt]( j, 2 ) );
                }
            }
        }//end of loop over integration points
    }
 
 
    std::string Info() const override
    {
        return "3 dimensional interface Prism with six nodes in 3D space";
    }
 
    void PrintInfo( std::ostream& rOStream ) const override
    {
        rOStream << "3 dimensional interface Prism with six nodes in 3D 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;
        }
    }
 
    Matrix& ShapeFunctionsLocalGradients( Matrix& rResult,
            const CoordinatesArrayType& rPoint ) const override
    {
        rResult.resize( 6, 3, false );
        noalias( rResult ) = ZeroMatrix( 6, 3 );
 
        rResult( 0, 0 ) = -1.0 + rPoint[2];
        rResult( 0, 1 ) = -1.0 + rPoint[2];
        rResult( 0, 2 ) =  -1.0 + rPoint[0] + rPoint[1];
 
        rResult( 1, 0 ) =  1.0 - rPoint[2];
        rResult( 1, 1 ) =  0.0;
        rResult( 1, 2 ) =  -rPoint[0];
 
        rResult( 2, 0 ) =  0.0;
        rResult( 2, 1 ) =  1.0 - rPoint[2];
        rResult( 2, 2 ) =  -rPoint[1];
 
        rResult( 3, 0 ) = -rPoint[2];
        rResult( 3, 1 ) = -rPoint[2];
        rResult( 3, 2 ) =  1.0 - rPoint[0] - rPoint[1];
 
        rResult( 4, 0 ) =  rPoint[2];
        rResult( 4, 1 ) =  0.0;
        rResult( 4, 2 ) =  rPoint[0];
 
        rResult( 5, 0 ) =  0.0;
        rResult( 5, 1 ) =  rPoint[2];
        rResult( 5, 2 ) =  rPoint[1];
 
        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 );
        }
 
        //TODO: this is not correct for a prism
 
        rResult[0].resize( 2, 2 , false);
        rResult[1].resize( 2, 2 , false);
        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 );
        }
 
        //TODO: this is not correct for a prism
 
        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 );
    }
 
    PrismInterface3D6(): 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 = 6;
        //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()
                                                - integration_points[pnt].Z()
                                                + ( integration_points[pnt].X() * integration_points[pnt].Z() )
                                                + ( integration_points[pnt].Y() * integration_points[pnt].Z() ) );
            shape_function_values( pnt, 1 ) = integration_points[pnt].X()
                                              - ( integration_points[pnt].X() * integration_points[pnt].Z() );
            shape_function_values( pnt, 2 ) = integration_points[pnt].Y()
                                              - ( integration_points[pnt].Y() * integration_points[pnt].Z() );
            shape_function_values( pnt, 3 ) = integration_points[pnt].Z()
                                              - ( integration_points[pnt].X() * integration_points[pnt].Z() )
                                              - ( integration_points[pnt].Y() * integration_points[pnt].Z() );
            shape_function_values( pnt, 4 ) = ( integration_points[pnt].X() * integration_points[pnt].Z() );
            shape_function_values( pnt, 5 ) = ( integration_points[pnt].Y() * integration_points[pnt].Z() );
        }
 
        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 = ZeroMatrix( 6, 3 );
            result( 0, 0 ) = -1.0 + integration_points[pnt].Z();
            result( 0, 1 ) = -1.0 + integration_points[pnt].Z();
            result( 0, 2 ) = -1.0 + integration_points[pnt].X() + integration_points[pnt].Y();
            result( 1, 0 ) =  1.0 - integration_points[pnt].Z();
            result( 1, 1 ) =  0.0;
            result( 1, 2 ) =  -integration_points[pnt].X();
            result( 2, 0 ) =  0.0;
            result( 2, 1 ) =  1.0 - integration_points[pnt].Z();
            result( 2, 2 ) =  -integration_points[pnt].Y();
            result( 3, 0 ) =  -integration_points[pnt].Z();
            result( 3, 1 ) =  -integration_points[pnt].Z();
            result( 3, 2 ) =  1.0 - integration_points[pnt].X() - integration_points[pnt].Y();
            result( 4, 0 ) =  integration_points[pnt].Z();
            result( 4, 1 ) =  0.0;
            result( 4, 2 ) =  integration_points[pnt].X();
            result( 5, 0 ) =  0.0;
            result( 5, 1 ) =  integration_points[pnt].Z();
            result( 5, 2 ) =  integration_points[pnt].Y();
            d_shape_f_values[pnt] = result;
        }
 
        return d_shape_f_values;
    }
 
    static const IntegrationPointsContainerType AllIntegrationPoints()
    {
        IntegrationPointsContainerType integration_points =
        {
            {
                Quadrature<PrismGaussLobattoIntegrationPoints1,
                3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                Quadrature<PrismGaussLobattoIntegrationPoints2,
                3, IntegrationPoint<3> >::GenerateIntegrationPoints(),
                IntegrationPointsArrayType(),
                IntegrationPointsArrayType()
            }
        };
 
        return integration_points;
    }
 
    static const ShapeFunctionsValuesContainerType AllShapeFunctionsValues()
    {
        ShapeFunctionsValuesContainerType shape_functions_values =
        {
            {
                PrismInterface3D6<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_1 ),
                PrismInterface3D6<TPointType>::CalculateShapeFunctionsIntegrationPointsValues(
                    GeometryData::IntegrationMethod::GI_GAUSS_2 ),
                Matrix(),
                Matrix()
            }
        };
        return shape_functions_values;
    }
 
    static const ShapeFunctionsLocalGradientsContainerType
    AllShapeFunctionsLocalGradients()
    {
        ShapeFunctionsLocalGradientsContainerType shape_functions_local_gradients =
        {
            {
                PrismInterface3D6<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_1 ),
                PrismInterface3D6<TPointType>::CalculateShapeFunctionsIntegrationPointsLocalGradients( GeometryData::IntegrationMethod::GI_GAUSS_2 ),
                ShapeFunctionsGradientsType(),
                ShapeFunctionsGradientsType(),
            }
        };
        return shape_functions_local_gradients;
    }
 
 
 
 
 
 
    template<class TOtherPointType> friend class PrismInterface3D6;
 
 
 
 
}; // Class Geometry
 
 
 
 
template<class TPointType> inline std::istream& operator >> (
    std::istream& rIStream,
    PrismInterface3D6<TPointType>& rThis );
template<class TPointType> inline std::ostream& operator << (
    std::ostream& rOStream,
    const PrismInterface3D6<TPointType>& rThis )
{
    rThis.PrintInfo( rOStream );
    rOStream << std::endl;
    rThis.PrintData( rOStream );
    return rOStream;
}
 
 
template<class TPointType> const
GeometryData PrismInterface3D6<TPointType>::msGeometryData(
    &msGeometryDimension,
    GeometryData::IntegrationMethod::GI_GAUSS_2,
    PrismInterface3D6<TPointType>::AllIntegrationPoints(),
    PrismInterface3D6<TPointType>::AllShapeFunctionsValues(),
    AllShapeFunctionsLocalGradients()
);
 
template<class TPointType> const
GeometryDimension PrismInterface3D6<TPointType>::msGeometryDimension(3, 3);
 
}// namespace Kratos.
 
#endif // KRATOS_PRISM_INTERFACE_3D_6_H_INCLUDED  defined