calculate_signed_distance_to_3d_skin_process.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Daniel Baumgaertner
//                   Johannes Wolf
//
 
 
#if !defined(KRATOS_CALCULATE_DISTANCE_PROCESS_H_INCLUDED )
#define  KRATOS_CALCULATE_DISTANCE_PROCESS_H_INCLUDED
 
 
 
// System includes
#include <string>
#include <iostream> 
#include <ctime>
 
// External includes 
 
 
// Project includes
#include "includes/define.h"
#include "processes/process.h"
#include "includes/model_part.h"
#include "includes/deprecated_variables.h"
 
#include "spatial_containers/octree_binary.h"
#include "utilities/spatial_containers_configure.h"
#include "utilities/timer.h"
#include "utilities/math_utils.h"
#include "utilities/geometry_utilities.h"
#include "geometries/triangle_3d_3.h"
#include "geometries/quadrilateral_3d_4.h"
#include "utilities/body_normal_calculation_utils.h"
#include "includes/kratos_flags.h"
#include "utilities/binbased_fast_point_locator.h"
#include "utilities/binbased_nodes_in_element_locator.h"
#include "processes/calculate_distance_to_skin_process.h"
 
#ifdef _OPENMP
#include "omp.h"
#endif
 
using namespace boost::numeric::ublas;
 
 
namespace Kratos
{
 
class DistanceSpatialContainersConfigure
{
public:
    class CellNodeData
    {
        double mDistance;
        double mCoordinates[3];
        std::size_t mId;
    public:
        double& Distance(){return mDistance;}
        double& X() {return mCoordinates[0];}
        double& Y() {return mCoordinates[1];}
        double& Z() {return mCoordinates[2];}
        double& operator[](int i) {return mCoordinates[i];}
        std::size_t& Id(){return mId;}
    };
 
 
 
 
 
    enum { Dimension = 3,
           DIMENSION = 3,
           MAX_LEVEL = 12,
           MIN_LEVEL = 2    // this cannot be less than 2!!!
         };
 
    typedef Point                                           PointType;  
    typedef std::vector<double>::iterator                   DistanceIteratorType;
    typedef ModelPart::ElementsContainerType::ContainerType ContainerType;
    typedef ContainerType::value_type                       PointerType;
    typedef ContainerType::iterator                         IteratorType;
    typedef ModelPart::ElementsContainerType::ContainerType ResultContainerType;
    typedef ResultContainerType::value_type                 ResultPointerType;
    typedef ResultContainerType::iterator                   ResultIteratorType;
 
    typedef Element::Pointer                                        pointer_type;
    typedef CellNodeData                cell_node_data_type;
    typedef std::vector<CellNodeData*> data_type;
 
    typedef std::vector<PointerType>::iterator             PointerTypeIterator;
 
 
 
 
    KRATOS_CLASS_POINTER_DEFINITION(DistanceSpatialContainersConfigure);
 
 
    DistanceSpatialContainersConfigure() {}
 
    virtual ~DistanceSpatialContainersConfigure() {}
 
 
 
 
 
    static data_type* AllocateData() {
        return new data_type(27, (CellNodeData*)NULL);
    }
 
    static void CopyData(data_type* source, data_type* destination) {
        *destination = *source;
    }
 
    static void DeleteData(data_type* data) {
        delete data;
    }
 
    static inline void CalculateBoundingBox(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint)
    {
        rHighPoint = rObject->GetGeometry().GetPoint(0);
        rLowPoint  = rObject->GetGeometry().GetPoint(0);
        
        for (unsigned int point = 0; point<rObject->GetGeometry().PointsNumber(); point++)
        {
            for(std::size_t i = 0; i<3; i++)
            {
                rLowPoint[i]  =  (rLowPoint[i]  >  rObject->GetGeometry().GetPoint(point)[i] ) ?  rObject->GetGeometry().GetPoint(point)[i] : rLowPoint[i];
                rHighPoint[i] =  (rHighPoint[i] <  rObject->GetGeometry().GetPoint(point)[i] ) ?  rObject->GetGeometry().GetPoint(point)[i] : rHighPoint[i];
            }
        }
    }
 
    static inline void GetBoundingBox(const PointerType rObject, double* rLowPoint, double* rHighPoint)
    {
 
        for(std::size_t i = 0; i<3; i++)
        {
            rLowPoint[i]  =  rObject->GetGeometry().GetPoint(0)[i];
            rHighPoint[i] =  rObject->GetGeometry().GetPoint(0)[i];
        }
 
        for (unsigned int point = 0; point<rObject->GetGeometry().PointsNumber(); point++)
        {
            for(std::size_t i = 0; i<3; i++)
            {
                rLowPoint[i]  =  (rLowPoint[i]  >  rObject->GetGeometry().GetPoint(point)[i] ) ?  rObject->GetGeometry().GetPoint(point)[i] : rLowPoint[i];
                rHighPoint[i] =  (rHighPoint[i] <  rObject->GetGeometry().GetPoint(point)[i] ) ?  rObject->GetGeometry().GetPoint(point)[i] : rHighPoint[i];
            }
        }
    }
 
    static inline bool Intersection(const PointerType& rObj_1, const PointerType& rObj_2)
    {
        Element::GeometryType& geom_1 = rObj_1->GetGeometry();
        Element::GeometryType& geom_2 = rObj_2->GetGeometry();
        return  geom_1.HasIntersection(geom_2);
 
    }
 
 
    static inline bool  IntersectionBox(const PointerType& rObject,  const PointType& rLowPoint, const PointType& rHighPoint)
    {
        return rObject->GetGeometry().HasIntersection(rLowPoint, rHighPoint);
    }
 
 
    static  inline bool  IsIntersected(const Element::Pointer rObject, double Tolerance, const double* rLowPoint, const double* rHighPoint)
    {
        Point low_point(rLowPoint[0] - Tolerance, rLowPoint[1] - Tolerance, rLowPoint[2] - Tolerance);
        Point high_point(rHighPoint[0] + Tolerance, rHighPoint[1] + Tolerance, rHighPoint[2] + Tolerance);
 
        KRATOS_THROW_ERROR(std::logic_error, "Not Implemented method", "")
                //return HasIntersection(rObject->GetGeometry(), low_point, high_point);
    }
 
 
 
 
 
 
 
    virtual std::string Info() const
    {
        return " Spatial Containers Configure";
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const {}
 
    virtual void PrintData(std::ostream& rOStream) const {}
 
 
 
protected:
 
private:
 
    DistanceSpatialContainersConfigure& operator=(DistanceSpatialContainersConfigure const& rOther);
 
    DistanceSpatialContainersConfigure(DistanceSpatialContainersConfigure const& rOther);
 
 
}; // Class DistanceSpatialContainersConfigure
 
 
 
 
 
 
 
 
class CalculateSignedDistanceTo3DSkinProcess
        : public Process
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(CalculateSignedDistanceTo3DSkinProcess);
 
    typedef DistanceSpatialContainersConfigure ConfigurationType;
    typedef OctreeBinaryCell<ConfigurationType> CellType;
    typedef OctreeBinary<CellType> OctreeType;
    typedef ConfigurationType::cell_node_data_type CellNodeDataType;
    typedef Point PointType;  
    typedef OctreeType::cell_type::object_container_type object_container_type;
    typedef struct{
        array_1d<double,3>  Coordinates;
        array_1d<double,3>  StructElemNormal;
        unsigned int EdgeNode1;
        unsigned int EdgeNode2;
    }IntersectionNodeStruct;
    typedef struct{
        std::vector<IntersectionNodeStruct> IntNodes;
    }TetEdgeStruct;
 
 
 
    CalculateSignedDistanceTo3DSkinProcess(ModelPart& rThisModelPartStruc, ModelPart& rThisModelPartFluid)
        : mrSkinModelPart(rThisModelPartStruc), mrBodyModelPart(rThisModelPartStruc), mrFluidModelPart(rThisModelPartFluid)
    {
    }
 
    ~CalculateSignedDistanceTo3DSkinProcess() override
    {
    }
 
 
 
    void operator()()
    {
        Execute();
    }
 
 
 
 
    void Execute() override
    {
        KRATOS_TRY;
 
        GenerateOctree();
 
        //DistanceFluidStructure();
 
        CalculateDistanceToSkinProcess<3> distance_process(mrFluidModelPart, mrBodyModelPart);
        distance_process.Execute();
 
        //          ------------------------------------------------------------------
        //          GenerateNodes();
        CalculateDistance2(); // I have to change this. Pooyan.
        //mrSkinModelPart.GetCommunicator().AssembleCurrentData(DISTANCE);
        //          std::ofstream mesh_file1("octree1.post.msh");
        //          std::ofstream res_file("octree1.post.res");
        //          Timer::Start("Writing Gid conform Mesh");
        //          PrintGiDMesh(mesh_file1);
        //          PrintGiDResults(res_file);
        //          octree.PrintGiDMeshNew(mesh_file2);
        //          Timer::Stop("Writing Gid conform Mesh");
        //          delete octree. TODO: Carlos
        //          ------------------------------------------------------------------
 
        KRATOS_CATCH("");
    }
 
 
    void MappingPressureToStructure(BinBasedFastPointLocator<3>& node_locator)
    {
        //loop over nodes and find the tetra in which it falls, than do interpolation
        Vector N;
        const int max_results = 10000;
        BinBasedFastPointLocator<3>::ResultContainerType results(max_results);
        const int n_structure_nodes = mrSkinModelPart.Nodes().size();
 
#pragma omp parallel for firstprivate(results,N)
    //MY NEW LOOP: reset the viisted flaf
    for (int i = 0; i < n_structure_nodes; i++)
        {
            ModelPart::NodesContainerType::iterator iparticle = mrSkinModelPart.NodesBegin() + i;
            Node ::Pointer p_structure_node = *(iparticle.base());
        p_structure_node->Set(VISITED, false);
    }
        for (int i = 0; i < n_structure_nodes; i++)
        {
            ModelPart::NodesContainerType::iterator iparticle = mrSkinModelPart.NodesBegin() + i;
            Node ::Pointer p_structure_node = *(iparticle.base());
            BinBasedFastPointLocator<3>::ResultIteratorType result_begin = results.begin();
            Element::Pointer pElement;
 
            bool is_found = node_locator.FindPointOnMesh(p_structure_node->Coordinates(), N, pElement, result_begin, max_results);
 
            if (is_found == true)
            {
                array_1d<double,4> nodalPressures;
                const Vector& ElementalDistances = pElement->GetValue(ELEMENTAL_DISTANCES);
 
                Geometry<Node >& geom = pElement->GetGeometry();
 
                for(unsigned int j=0; j<geom.size(); j++)
                {
                    nodalPressures[j] = geom[j].FastGetSolutionStepValue(PRESSURE);
                }
 
                if(pElement->GetValue(SPLIT_ELEMENT)==true)
                {
                    array_1d<double,4> Npos,Nneg;
 
                    // Do mapping
                    ComputeDiscontinuousInterpolation((*p_structure_node),pElement->GetGeometry(),ElementalDistances,Npos,Nneg);
 
                    // Compute face pressure
                    double p_positive_structure = inner_prod(nodalPressures,Npos);
                    double p_negative_structure = inner_prod(nodalPressures,Nneg);
 
                    // Assign ModelPart::ElementIteratorface pressure to structure node
                    p_structure_node->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE) = p_positive_structure;
                    p_structure_node->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE) = p_negative_structure;
            p_structure_node->Set(VISITED);
                }
                else
                {
                    double p = inner_prod(nodalPressures,N);
                    p_structure_node->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE) = p;
                    p_structure_node->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE) = p;
                p_structure_node->Set(VISITED);
                }
            }
        }
    //AND NOW WE "TREAT" the bad nodes, the ones that belong to the structural faces that by some chance did not cross the fluid elements
    //to such nodes we simply extrapolate the pressure from the neighbors
    int n_bad_nodes=0;
    for (int i = 0; i < n_structure_nodes; i++)
        {
            ModelPart::NodesContainerType::iterator iparticle = mrSkinModelPart.NodesBegin() + i;
            Node ::Pointer p_structure_node = *(iparticle.base());
        if (p_structure_node->IsNot(VISITED))
        n_bad_nodes++;
    }
    //KRATOS_WATCH("THERE WERE THIS MANY BAD NODES ORIGINALLY")
    //KRATOS_WATCH(n_bad_nodes)
    while (n_bad_nodes >= 1.0) {
                int n_bad_nodes_backup = n_bad_nodes;
 
                for (int i = 0; i < n_structure_nodes; i++) {
                    ModelPart::NodesContainerType::iterator iparticle = mrSkinModelPart.NodesBegin() + i;
                    Node ::Pointer p_structure_node = *(iparticle.base());
 
                    //here we store the number of neigbor nodes that were given the pressure in the previous loop (i.e. were found)
                    if (p_structure_node->IsNot(VISITED)) {
                        int n_good_neighbors = 0;
                        double pos_pres = 0.0;
                        double neg_pres = 0.0;
                        GlobalPointersVector< Node >& neighours = p_structure_node->GetValue(NEIGHBOUR_NODES);
                        
                        for (GlobalPointersVector< Node >::iterator j = neighours.begin(); j != neighours.end(); j++) {
                            if (j->Is(VISITED)) {
                                n_good_neighbors++;
                                pos_pres += j->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE);
                                neg_pres += j->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE);
                                //KRATOS_WATCH("Good neighbor found")
                            }
                        }
                        if (n_good_neighbors != 0) {
                            pos_pres /= n_good_neighbors;
                            neg_pres /= n_good_neighbors;
                            p_structure_node->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE) = pos_pres;
                            p_structure_node->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE) = neg_pres;
                            p_structure_node->Set(VISITED);
                            n_bad_nodes--;
                        }
                        //KRATOS_WATCH(pos_pres)
                        //KRATOS_WATCH(neg_pres)          
                    }
                }
                
                if(n_bad_nodes == n_bad_nodes_backup) break; //WE BREAK THE WHILE HERE, OTHERWISE THE CODE HANGS (it was not able to remove any other node)
 
                /*int n_bad_nodes=0;
                for (int i = 0; i < n_structure_nodes; i++)
                {
                   ModelPart::NodesContainerType::iterator iparticle = mrSkinModelPart.NodesBegin() + i;
                   Node ::Pointer p_structure_node = *(iparticle.base());
                   if (p_structure_node->IsNot(VISITED))
                       n_bad_nodes++;         
                }
                 */
                //KRATOS_WATCH(n_bad_nodes)
 
            }
        //THE BELOW ONE IS A "CHEAT".. THERE IS A PROBLEM OF INCORRECT PROJECTION BETWEEN THE MESHES AT SOME POINTS
        //FOR NODES WITH PRESSURE VERY DIFFERENT FROM THAT OF THE NEIGHBORS, I JUST TAKE THE NEIGHBOR PRESSURE AVERAGED
        for (int i = 0; i < n_structure_nodes; i++)
        {
            ModelPart::NodesContainerType::iterator iparticle = mrSkinModelPart.NodesBegin() + i;
            Node ::Pointer p_structure_node = *(iparticle.base());
 
            double pos_pressure=p_structure_node->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE);
            double neg_pressure=p_structure_node->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE);
                    
                    GlobalPointersVector< Node >& neighours = p_structure_node->GetValue(NEIGHBOUR_NODES);
            
            if (neighours.size()>=1.0)
            {               
                double av_pos_pres=0.0;
                double av_neg_pres=0.0;
                for( GlobalPointersVector< Node >::iterator j = neighours.begin();
                        j != neighours.end(); j++)
                {
                
                    av_pos_pres+=j->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE);
                    av_neg_pres+=j->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE);
                                        
                }
                av_pos_pres/=neighours.size();
                av_neg_pres/=neighours.size();
                
                //IF the average pressure of the neighbors is 10 times lower than of the given node, something is bad and we reset its value
                if (fabs(pos_pressure)>3.0*fabs(av_pos_pres))
                    {
                    p_structure_node->FastGetSolutionStepValue(POSITIVE_FACE_PRESSURE) = av_pos_pres;
                    //KRATOS_WATCH("BAD NODE")
                    }
                if (fabs(neg_pressure)>3.0*fabs(av_neg_pres))
                    {
                    p_structure_node->FastGetSolutionStepValue(NEGATIVE_FACE_PRESSURE) = av_neg_pres;
                    //KRATOS_WATCH("BAD NODE")
                    }                       
             
            }
         }
        
 
 
    }
 
 
    void ComputeDiscontinuousInterpolation( const Node& pNode,
                                            Geometry< Node >& geom,
                                            const array_1d<double,4>& distances,
                                            array_1d<double,4>& Npos,
                                            array_1d<double,4>& Nneg)
    {
        //count positives
        int n_positives = 0;
        for(unsigned int i=0; i<distances.size(); i++)
            if(distances[i]>0) n_positives++;
 
        //generate the points on the edges at the zero of the distance function
        //generate "father nodes", defined as the end nodes of the edge on which the local point is located
        std::vector< Point > edge_points;
        edge_points.reserve(4);
        array_1d<unsigned int, 4> positive_fathers,  negative_fathers;  //there are at most 4 cut edges
        unsigned int k=0;
        unsigned int l=0;
 
        for(unsigned int i=0; i<3; i++)
        {
            for(unsigned int j=i+1; j<4; j++) // go through the edges 01, 02, 03, 12, 13, 23
            {
                double di = distances[i];
                double dj = distances[j];
 
                if(di*dj < 0) //edge is cut
                {
                    //generate point on edge by linear interpolation
                    double Ni = fabs(dj) / ( fabs(di) + fabs(dj) );
                    double Nj = 1.0 - Ni;
                    Point edge_point(Ni * geom[i] + Nj * geom[j]);
                    edge_points.push_back(edge_point);
 
                    //store the id of the positive and negative fathers
                    if(di > 0.0)
                    {
                        positive_fathers[k++] = i;
                        negative_fathers[l++] = j;
                    }
                    else
                    {
                        positive_fathers[k++] = j;
                        negative_fathers[l++] = i;
                    }
                }
            }
        }
 
        if(edge_points.size() == 3)
        {
            //compute local shape functions (tell how to interpolate from the edge nodes)
            Vector Nlocal(3);
 
            //form a triangle with the edge nodes
            Triangle3D3< Point > triangle(Point::Pointer(new Point(edge_points[0])), 
                         Point::Pointer(new Point(edge_points[1])), 
                         Point::Pointer(new Point(edge_points[2]))
                         );
 
            array_1d<double,3> local_coords;
            local_coords = triangle.PointLocalCoordinates(local_coords, pNode);
 
            for(unsigned int i=0; i<3;i++)
                Nlocal[i] = triangle.ShapeFunctionValue(i, local_coords );
 
            noalias(Npos) = ZeroVector(4);
            noalias(Nneg) = ZeroVector(4);
            for(unsigned int i=0; i<3; i++)
            {
                Npos[ positive_fathers[i] ] += Nlocal[i];
                Nneg[ negative_fathers[i] ] += Nlocal[i];
            }
        }
 
        if(edge_points.size() == 4)
        {
            //compute local shape functions (tell how to interpolate from the edge nodes)
            Vector Nlocal(4);
 
            //form a quadrilatera with the 4 cut nodes
            array_1d<double,3> x21 = edge_points[1] - edge_points[0];
            array_1d<double,3> x31 = edge_points[2] - edge_points[0];
            array_1d<double,3> x41 = edge_points[3] - edge_points[0];
 
            //define a vector oriented as x21
            array_1d<double,3> v1 = x21 / norm_2(x21);
 
            BoundedMatrix<double,4,3> DN_DX;
            array_1d<double,4> msN;
            double Area;
            GeometryUtils::CalculateGeometryData( geom, DN_DX, msN, Area );
 
            array_1d<double,3> n = prod(trans(DN_DX),distances);
            n /= norm_2(n);
 
            array_1d<double,3> v2;
            MathUtils<double>::CrossProduct(v2,v1,n); // v2 = v1 x n
 
            array_1d<double,3> angles;
            angles[0] = 0.0; //angle between x21 and v1
            angles[1] = atan2( inner_prod(x31,v2), inner_prod(x31,v1) ); //angle between x31 and v1
            angles[2] = atan2( inner_prod(x41,v2), inner_prod(x41,v1) ); //angle between x31 and v1
 
            double max_angle = 0.0;
            double min_angle = 0.0;
            unsigned int min_pos = 1;
            unsigned int max_pos = 1;
            for(unsigned int i=1; i<3; i++)
            {
                if(angles[i] < min_angle)
                {
                    min_pos = i+1; //this is the local index of the edge point which forms the minimal angle
                    min_angle = angles[i];
                }
                else if(angles[i] > max_angle)
                {
                    max_pos = i+1; //this is the local index of the edge point which forms the maximal angle
                    max_angle = angles[i];
                }
            }
 
            //find the pos of the center node
            unsigned int center_pos = 0;
            for(unsigned int i=1; i<4; i++)
            {
                if((i!= min_pos) && (i!=max_pos))
                { center_pos = i; }
            }
 
            //form a quadrilateral with the edge nodes
            Quadrilateral3D4< Point > quad = Quadrilateral3D4< Point >(
            Point::Pointer(new Point(edge_points[0])),
            Point::Pointer(new Point(edge_points[min_pos])),
            Point::Pointer(new Point(edge_points[center_pos])), 
            Point::Pointer(new Point(edge_points[max_pos]))
            );
 
            array_1d<double,3> local_coords;
            local_coords = quad.PointLocalCoordinates(local_coords, pNode);
 
            array_1d<unsigned int, 4> indices;
            indices[0] = 0;
            indices[1] = min_pos;
            indices[2] = center_pos;
            indices[3] = max_pos;
 
            for(unsigned int i=0; i<4;i++)
                Nlocal[ i ]  = quad.ShapeFunctionValue(i, local_coords );
 
            noalias(Npos) = ZeroVector(4);
            noalias(Nneg) = ZeroVector(4);
            for(unsigned int i=0; i<4; i++)
            {
                Npos[ positive_fathers[i] ] += Nlocal[indices[i]];
                Nneg[ negative_fathers[i] ] += Nlocal[indices[i]];
            }
        }
    }
 
 
    void AveragePressureToNode(BinBasedFastPointLocator<3>& node_locator,
                               Node& node)
    {
        //loop over nodes and find the tetra in which it falls, than do interpolation
        Vector N;
        const int max_results = 10000;
        BinBasedFastPointLocator<3>::ResultContainerType results(max_results);
        BinBasedFastPointLocator<3>::ResultIteratorType result_begin = results.begin();
        Element::Pointer pElement;
 
        bool is_found = node_locator.FindPointOnMesh(node.Coordinates(), N, pElement, result_begin, max_results);
 
        if (is_found == true)
        {
            array_1d<double,4> nodalPressures;
            const Vector& ElementalDistances = pElement->GetValue(ELEMENTAL_DISTANCES);
            Geometry<Node >& geom = pElement->GetGeometry();
 
            for(unsigned int i=0; i<4; i++)
                nodalPressures[i] = geom[i].GetSolutionStepValue(PRESSURE);
 
            if(pElement->GetValue(SPLIT_ELEMENT)==true)
            {
                // Compute average of all positive and all negative values
                double positiveAverage = 0;
                double negativeAverage = 0;
                unsigned int nPos = 0;
                unsigned int nNeg = 0;
 
                for(unsigned int i=0 ; i<4 ; i++)
                {
                    if(ElementalDistances[i]>=0)
                    {
                        positiveAverage += nodalPressures[i];
                        nPos++;
                    }
                    else
                    {
                        negativeAverage += nodalPressures[i];
                        nNeg++;
                    }
                }
 
                positiveAverage /= nPos;
                negativeAverage /= nNeg;
 
                // Assign Pressures
                node.GetSolutionStepValue(POSITIVE_FACE_PRESSURE,0) = positiveAverage;
                node.GetSolutionStepValue(NEGATIVE_FACE_PRESSURE,0) = negativeAverage;
            }
            else
            {
                // Compute average of all positive and all negative values
                double Average = 0;
 
                // for output of
                for(unsigned int i = 0 ; i<4 ; i++)
                    Average += nodalPressures[i];
 
                Average /= 4;
 
                // Assign Pressures
                node.GetSolutionStepValue(POSITIVE_FACE_PRESSURE,0) = Average;
                node.GetSolutionStepValue(NEGATIVE_FACE_PRESSURE,0) = Average;
            }
        }
    }
 
 
    void DistanceFluidStructure()
    {
      //std::cout << "Start calculating Elemental distances..." << std::endl;
 
        // Initialize Elemental distances in the domain
        Initialize();
 
        // Initialize index table that defines line Edges of fluid Element
        BoundedMatrix<unsigned int,6,2> TetEdgeIndexTable;
        SetIndexTable(TetEdgeIndexTable);
 
        // loop over all fluid Elements
        // this loop is parallelized using openmp
#ifdef _OPENMP
        int number_of_threads = omp_get_max_threads();
#else
        int number_of_threads = 1;
#endif
 
        ModelPart::ElementsContainerType& pElements = mrFluidModelPart.Elements();
 
        DenseVector<unsigned int> Element_partition;
        CreatePartition(number_of_threads, pElements.size(), Element_partition);
 
#pragma omp parallel for
        for (int k = 0; k < number_of_threads; k++)
        {
            ModelPart::ElementsContainerType::iterator it_begin = pElements.ptr_begin() + Element_partition[k];
            ModelPart::ElementsContainerType::iterator it_end = pElements.ptr_begin() + Element_partition[k+1];
 
            // assemble all Elements
            for (ModelPart::ElementIterator it = it_begin; it != it_end; ++it)
            {
                CalcElementDistances( it , TetEdgeIndexTable );
            }
        }        
 
        // Finally, each tetrahedral Element has 4 distance values. But each node belongs to
        // several Elements, such that it is assigned several distance values
        // --> now synchronize these values by finding the minimal distance and assign to each node a minimal nodal distance
        AssignMinimalNodalDistance();
 
        //std::cout << "Finished calculating Elemental distances..." << std::endl;
    }
 
 
    void Initialize()
    {
        const double initial_distance = 1.0;
 
        ModelPart::NodesContainerType::ContainerType& nodes = mrFluidModelPart.NodesArray();
 
        // reset the node distance to 1.0 which is the maximum distance in our normalized space.
        int nodesSize = nodes.size();
 
#pragma omp parallel for firstprivate(nodesSize)
        for(int i = 0 ; i < nodesSize ; i++)
            nodes[i]->GetSolutionStepValue(DISTANCE) = initial_distance;
 
        ModelPart::ElementsContainerType::ContainerType& fluid_Elements = mrFluidModelPart.ElementsArray();
 
        array_1d<double,4> ElementalDistances;
        ElementalDistances[0] = initial_distance;
        ElementalDistances[1] = initial_distance;
        ElementalDistances[2] = initial_distance;
        ElementalDistances[3] = initial_distance;
 
        // reset the Elemental distance to 1.0 which is the maximum distance in our normalized space.
        // also initialize the embedded velocity of the fluid Element
        int ElementsSize = fluid_Elements.size();
 
#pragma omp parallel for firstprivate(ElementsSize)
        for(int i = 0 ; i < ElementsSize ; i++)
        {
            fluid_Elements[i]->GetValue(ELEMENTAL_DISTANCES) = ElementalDistances;
            fluid_Elements[i]->GetValue(SPLIT_ELEMENT) = false;
            fluid_Elements[i]->GetValue(EMBEDDED_VELOCITY)=ZeroVector(3);
        }
 
    }
 
 
    void SetIndexTable( BoundedMatrix<unsigned int,6,2>& TetEdgeIndexTable )
    {
        // Initialize index table to define line Edges of fluid Element
        TetEdgeIndexTable(0,0) = 0;
        TetEdgeIndexTable(0,1) = 1;
        TetEdgeIndexTable(1,0) = 0;
        TetEdgeIndexTable(1,1) = 2;
        TetEdgeIndexTable(2,0) = 0;
        TetEdgeIndexTable(2,1) = 3;
        TetEdgeIndexTable(3,0) = 1;
        TetEdgeIndexTable(3,1) = 2;
        TetEdgeIndexTable(4,0) = 1;
        TetEdgeIndexTable(4,1) = 3;
        TetEdgeIndexTable(5,0) = 2;
        TetEdgeIndexTable(5,1) = 3;
    }
 
 
    void CalcElementDistances( ModelPart::ElementsContainerType::iterator& i_fluidElement,
                               BoundedMatrix<unsigned int,6,2>            TetEdgeIndexTable )
    {
        std::vector<OctreeType::cell_type*> leaves;
        std::vector<TetEdgeStruct>          IntersectedTetEdges;
        unsigned int NumberIntersectionsOnTetCorner = 0;
 
        // Get leaves of octree intersecting with fluid Element
        mpOctree->GetIntersectedLeaves(*(i_fluidElement).base(),leaves);
 
        int intersection_counter = 0;
 
        // Loop over all 6 line Edges of the tetrahedra
        for(unsigned int i_tetEdge = 0;
            i_tetEdge < 6;
            i_tetEdge++)
        {
            IdentifyIntersectionNodes( i_fluidElement, i_tetEdge, leaves, IntersectedTetEdges, NumberIntersectionsOnTetCorner, TetEdgeIndexTable, intersection_counter );
        }
 
        if (intersection_counter!=0)
        {
            i_fluidElement->GetValue(EMBEDDED_VELOCITY) /= intersection_counter;
        }
 
        if(IntersectedTetEdges.size() > 0)
            CalcDistanceTo3DSkin( IntersectedTetEdges , i_fluidElement , NumberIntersectionsOnTetCorner );
    }
 
 
    void IdentifyIntersectionNodes( ModelPart::ElementsContainerType::iterator&   i_fluidElement,
                                    unsigned int                                  i_tetEdge,
                                    std::vector<OctreeType::cell_type*>&          leaves,
                                    std::vector<TetEdgeStruct>&                   IntersectedTetEdges,
                                    unsigned int&                                 NumberIntersectionsOnTetCorner,
                                    BoundedMatrix<unsigned int,6,2>              TetEdgeIndexTable,
                                    int&                                          intersection_counter )
    {
        std::vector<unsigned int> IntersectingStructElemID;
        TetEdgeStruct             NewTetEdge;
        unsigned int              NumberIntersectionsOnTetCornerCurrentEdge = 0;
 
        // Get nodes of line Edge
        unsigned int EdgeStartIndex = TetEdgeIndexTable(i_tetEdge,0);
        unsigned int EdgeEndIndex   = TetEdgeIndexTable(i_tetEdge,1);
 
        PointType& P1 = i_fluidElement->GetGeometry()[EdgeStartIndex];
        PointType& P2 = i_fluidElement->GetGeometry()[EdgeEndIndex];
 
        double EdgeNode1[3] = {P1.X() , P1.Y() , P1.Z()};
        double EdgeNode2[3] = {P2.X() , P2.Y() , P2.Z()};
    
        // loop over all octree cells which are intersected by the fluid Element
        for(unsigned int i_cell = 0 ; i_cell < leaves.size() ; i_cell++)
        {
            // Structural Element contained in one cell of the octree
            object_container_type* struct_elem = (leaves[i_cell]->pGetObjects());
 
            // loop over all structural Elements within each octree cell
            for(object_container_type::iterator i_StructElement = struct_elem->begin(); i_StructElement != struct_elem->end(); i_StructElement++)
            {
 
                if( StructuralElementNotYetConsidered( (*i_StructElement)->Id() , IntersectingStructElemID ) )
                {
 
                    // Calculate and associate intersection point to the current fluid Element
                    double IntersectionPoint[3] = {0.0 , 0.0 , 0.0};
                    int TetEdgeHasIntersections = IntersectionTriangleSegment( (*i_StructElement)->GetGeometry() , EdgeNode1 , EdgeNode2 , IntersectionPoint );
 
                    if( TetEdgeHasIntersections == 1 )
                    {
                        IntersectionNodeStruct NewIntersectionNode;
 
                        // Assign information to the intersection node
                        NewIntersectionNode.Coordinates[0] = IntersectionPoint[0];
                        NewIntersectionNode.Coordinates[1] = IntersectionPoint[1];
                        NewIntersectionNode.Coordinates[2] = IntersectionPoint[2];
             
                        if( IsIntersectionNodeOnTetEdge( IntersectionPoint , EdgeNode1 , EdgeNode2 ) )
                        {
                            if ( IsNewIntersectionNode( NewIntersectionNode , IntersectedTetEdges ) )
                            {
                                // Calculate normal of the structural Element at the position of the intersection point
                                CalculateNormal3D((*i_StructElement)->GetGeometry()[0],
                                                  (*i_StructElement)->GetGeometry()[1],
                                                  (*i_StructElement)->GetGeometry()[2],
                                                  NewIntersectionNode.StructElemNormal);
 
                                // check, how many intersection nodes are located on corner points of the tetrahedra
                                if ( IsIntersectionOnCorner( NewIntersectionNode , EdgeNode1 , EdgeNode2) )
                                {
                                    NumberIntersectionsOnTetCornerCurrentEdge++;
 
                                    // only allow one intersection node on a tet edge
                                    if(NumberIntersectionsOnTetCornerCurrentEdge < 2)
                                    {
                                        // add the new intersection point to the list of intersection points of the fluid Element
                                        NewIntersectionNode.EdgeNode1 = EdgeStartIndex;
                                        NewIntersectionNode.EdgeNode2 = EdgeEndIndex;
                                        NewTetEdge.IntNodes.push_back(NewIntersectionNode);
 
                                        // if tet edge belonging to this intersection point is not already marked as "IntersectedTetEdge" --> put it into the respective container
                                        // when a second intersection node is found, then it is not necessary to push_back again
                                        if( NewTetEdge.IntNodes.size() == 1 )
                                            IntersectedTetEdges.push_back(NewTetEdge);
                                    }
 
                                    // this corner intersection node is only considered once for each tet edge
                                    if(NumberIntersectionsOnTetCornerCurrentEdge==1)
                                    {
                                        NumberIntersectionsOnTetCorner++;
                                    }
                                }
                                else
                                {
                                    // add the new intersection point to the list of intersection points of the fluid Element
                                    NewIntersectionNode.EdgeNode1 = EdgeStartIndex;
                                    NewIntersectionNode.EdgeNode2 = EdgeEndIndex;
                                    NewTetEdge.IntNodes.push_back(NewIntersectionNode);
 
                                    // velocity mapping structure --> fluid
                                    array_1d<double,3> emb_vel = (*i_StructElement)->GetGeometry()[0].GetSolutionStepValue(VELOCITY);
                                    emb_vel += (*i_StructElement)->GetGeometry()[1].GetSolutionStepValue(VELOCITY);
                                    emb_vel += (*i_StructElement)->GetGeometry()[2].GetSolutionStepValue(VELOCITY);
 
                                    i_fluidElement->GetValue(EMBEDDED_VELOCITY) += emb_vel/3;
                                    intersection_counter++;
                                }
                            }
                        }
                    }
                }
            }
        }
 
        // Finally put the found intersection nodes into the container
        if( NewTetEdge.IntNodes.size() > 0 )
    {
      if(NumberIntersectionsOnTetCornerCurrentEdge == 0)
            IntersectedTetEdges.push_back(NewTetEdge);
    }
    }
 
 
    bool StructuralElementNotYetConsidered( unsigned int                IDCurrentStructElem,
                                            std::vector<unsigned int>&  IntersectingStructElemID )
    {
        // check if the structural Element was already considered as intersecting Element
        for(unsigned int k = 0 ; k < IntersectingStructElemID.size() ; k++)
        {
            if( IDCurrentStructElem == IntersectingStructElemID[k] )
                return false;
        }
 
        // if structural Element has not been considered in another octree, which also intersects the fluid Element
        // add the new object ID to the vector
        IntersectingStructElemID.push_back( IDCurrentStructElem );
        return true;
    }
 
 
    bool IsIntersectionNodeOnTetEdge( double* IntersectionPoint,
                                      double* EdgeNode1,
                                      double* EdgeNode2 )
    {
        // check, if intersection point is located on any edge of the fluid Element
        array_1d<double,3> ConnectVectTetNodeIntNode1;
        array_1d<double,3> ConnectVectTetNodeIntNode2;
        array_1d<double,3> EdgeVector;
 
        ConnectVectTetNodeIntNode1[0] = IntersectionPoint[0] - EdgeNode1[0];
        ConnectVectTetNodeIntNode1[1] = IntersectionPoint[1] - EdgeNode1[1];
        ConnectVectTetNodeIntNode1[2] = IntersectionPoint[2] - EdgeNode1[2];
 
        ConnectVectTetNodeIntNode2[0] = IntersectionPoint[0] - EdgeNode2[0];
        ConnectVectTetNodeIntNode2[1] = IntersectionPoint[1] - EdgeNode2[1];
        ConnectVectTetNodeIntNode2[2] = IntersectionPoint[2] - EdgeNode2[2];
 
        double LengthConnectVect1 = norm_2( ConnectVectTetNodeIntNode1 );
        double LengthConnectVect2 = norm_2( ConnectVectTetNodeIntNode2 );
 
        EdgeVector[0] = EdgeNode2[0] - EdgeNode1[0];
        EdgeVector[1] = EdgeNode2[1] - EdgeNode1[1];
        EdgeVector[2] = EdgeNode2[2] - EdgeNode1[2];
 
        double MaxEdgeLength = norm_2( EdgeVector );
 
        // if both connection vectors (corner point --> intersection point)
        // are smaller or equal to the edge length of tetrahedra,
        // then intersection point is located on the edge
        if( (LengthConnectVect1 <= (MaxEdgeLength)) && (LengthConnectVect2 <= (MaxEdgeLength)) )
            return true;
        else
            return false;
    }
 
 
    bool IsNewIntersectionNode( IntersectionNodeStruct&     NewIntersectionNode,
                                std::vector<TetEdgeStruct>& IntersectedTetEdges )
    {
        array_1d<double,3> DiffVector;
        double NormDiffVector = 0;
        unsigned int NumberIntNodes = 0;
 
        for( unsigned int i_TetEdge = 0 ; i_TetEdge < IntersectedTetEdges.size() ; i_TetEdge++ )
        {
            NumberIntNodes = IntersectedTetEdges[i_TetEdge].IntNodes.size();
            for( unsigned int i_IntNode = 0 ; i_IntNode < NumberIntNodes ; i_IntNode++ )
            {
                DiffVector[0] = NewIntersectionNode.Coordinates[0] - IntersectedTetEdges[i_TetEdge].IntNodes[i_IntNode].Coordinates[0];
                DiffVector[1] = NewIntersectionNode.Coordinates[1] - IntersectedTetEdges[i_TetEdge].IntNodes[i_IntNode].Coordinates[1];
                DiffVector[2] = NewIntersectionNode.Coordinates[2] - IntersectedTetEdges[i_TetEdge].IntNodes[i_IntNode].Coordinates[2];
 
                NormDiffVector = norm_2(DiffVector);
 
                if( NormDiffVector < epsilon )
                    return false;
            }
        }
 
        // if the new intersection node is not existing (as intersection with a corner point), then return false
        return true;
    }
 
 
    bool IsIntersectionOnCorner( IntersectionNodeStruct& NewIntersectionNode,
                                 double*                 EdgeNode1,
                                 double*                 EdgeNode2 )
    {
        array_1d<double,3> DiffVector;
        double NormDiffVector;
 
        DiffVector[0] = EdgeNode1[0] - NewIntersectionNode.Coordinates[0];
        DiffVector[1] = EdgeNode1[1] - NewIntersectionNode.Coordinates[1];
        DiffVector[2] = EdgeNode1[2] - NewIntersectionNode.Coordinates[2];
        NormDiffVector = norm_2(DiffVector);
 
        if( NormDiffVector < epsilon )
            return true;
 
        DiffVector[0] = EdgeNode2[0] - NewIntersectionNode.Coordinates[0];
        DiffVector[1] = EdgeNode2[1] - NewIntersectionNode.Coordinates[1];
        DiffVector[2] = EdgeNode2[2] - NewIntersectionNode.Coordinates[2];
        NormDiffVector = norm_2(DiffVector);
 
        if( NormDiffVector < epsilon )
            return true;
        else
            return false;
    }
 
 
    void CalculateNormal3D( Point&       Point1,
                            Point&       Point2,
                            Point&       Point3,
                            array_1d<double,3>&   rResultNormal )
    {
        array_1d<double,3> v1 = Point2 - Point1;
        array_1d<double,3> v2 = Point3 - Point1;
 
        MathUtils<double>::CrossProduct(rResultNormal,v1,v2);
        rResultNormal *= 0.5;
    }
 
 
    void CalcDistanceTo3DSkin( std::vector<TetEdgeStruct>&                 IntersectedTetEdges,
                               ModelPart::ElementsContainerType::iterator& i_fluid_Element,
                               unsigned int                                NumberIntersectionsOnTetCorner )
    {
        std::vector<IntersectionNodeStruct> NodesOfApproximatedStructure;
        array_1d<double,4> ElementalDistances;
 
        FillIntNodesContainer(IntersectedTetEdges,NodesOfApproximatedStructure);
    
        // Intersection with one corner point
        if( NodesOfApproximatedStructure.size() == 1 && NumberIntersectionsOnTetCorner == 1 )
        {
            CalcSignedDistancesToOneIntNode(i_fluid_Element,NodesOfApproximatedStructure,ElementalDistances);
            i_fluid_Element->GetValue(SPLIT_ELEMENT) = true;
        }
 
        // Intersection with two corner points / one tetrahedra edge
        if( NodesOfApproximatedStructure.size() == 2 && NumberIntersectionsOnTetCorner == 2 )
        {
            CalcSignedDistancesToTwoIntNodes(i_fluid_Element,NodesOfApproximatedStructure,ElementalDistances);
            i_fluid_Element->GetValue(SPLIT_ELEMENT) = true;
        }
 
        // Intersection with three tetrahedra edges
        if( NodesOfApproximatedStructure.size() == 3 )
        {
            CalcSignedDistancesToThreeIntNodes(i_fluid_Element,NodesOfApproximatedStructure,ElementalDistances);
            i_fluid_Element->GetValue(SPLIT_ELEMENT) = true;
        }
 
        // Intersection with more than three tetrahedra edges
        if( NodesOfApproximatedStructure.size() > 3 )
        {
            CalcSignedDistancesToMoreThanThreeIntNodes(i_fluid_Element,NodesOfApproximatedStructure,ElementalDistances,IntersectedTetEdges);
            i_fluid_Element->GetValue(SPLIT_ELEMENT) = true;
        }
 
        // Postprocessing treatment of Elemental distances
        if( i_fluid_Element->GetValue(SPLIT_ELEMENT) == true )
            AvoidZeroDistances(i_fluid_Element, ElementalDistances);
 
        // In case there is intersection with fluid Element: assign distances to the Element
        if( i_fluid_Element->GetValue(SPLIT_ELEMENT) == true )
            i_fluid_Element->GetValue(ELEMENTAL_DISTANCES) = ElementalDistances;
    }
 
 
    void FillIntNodesContainer( std::vector<TetEdgeStruct>&          IntersectedTetEdges,
                                std::vector<IntersectionNodeStruct>& NodesOfApproximatedStructure )
    {
        const unsigned int NumberCutEdges = IntersectedTetEdges.size();
 
        for(unsigned int i_TetEdge = 0 ; i_TetEdge < NumberCutEdges ; i_TetEdge++)
        {
            unsigned int NumberIntNodes = IntersectedTetEdges[i_TetEdge].IntNodes.size();
 
            for( unsigned int i_IntNode = 0 ; i_IntNode < NumberIntNodes ; i_IntNode++ )
            {
                NodesOfApproximatedStructure.push_back(IntersectedTetEdges[i_TetEdge].IntNodes[i_IntNode]);
            }
        }
    }
 
 
    void CalcSignedDistancesToOneIntNode( ModelPart::ElementsContainerType::iterator& i_fluid_Element,
                                          std::vector<IntersectionNodeStruct>         NodesOfApproximatedStructure,
                                          array_1d<double,4>&                         ElementalDistances )
    {
        Geometry< Node >& rFluidGeom = i_fluid_Element->GetGeometry();
 
        Point  P1;
        P1.Coordinates() = NodesOfApproximatedStructure[0].Coordinates;
 
        array_1d<double,3>&  Normal = NodesOfApproximatedStructure[0].StructElemNormal;
 
        // Compute distance values for all tet-nodes
        for(unsigned int i_TetNode = 0 ; i_TetNode < 4 ; i_TetNode++)
        {
            ElementalDistances[i_TetNode] = PointDistanceToPlane(P1, Normal, rFluidGeom[i_TetNode]);
        }
    }
 
 
    void CalcSignedDistancesToTwoIntNodes( ModelPart::ElementsContainerType::iterator& i_fluid_Element,
                                           std::vector<IntersectionNodeStruct>         NodesOfApproximatedStructure,
                                           array_1d<double,4>&                         ElementalDistances )
    {
        Geometry< Node >& rFluidGeom = i_fluid_Element->GetGeometry();
 
        Point  P1;
        P1.Coordinates() = NodesOfApproximatedStructure[0].Coordinates;
 
        // Get normal at intersections, average them and check direction of distances
        array_1d<double,3> NormalAtIntersectionNode1 = NodesOfApproximatedStructure[0].StructElemNormal;
        array_1d<double,3> NormalAtIntersectionNode2 = NodesOfApproximatedStructure[1].StructElemNormal;
 
        // Compute normal of surface plane
        array_1d<double,3> Normal;
        Normal[0] = 0.5*(NormalAtIntersectionNode1[0] + NormalAtIntersectionNode2[0]);
        Normal[1] = 0.5*(NormalAtIntersectionNode1[1] + NormalAtIntersectionNode2[1]);
        Normal[2] = 0.5*(NormalAtIntersectionNode1[2] + NormalAtIntersectionNode2[2]);
 
        // Check whether orientation of normal is in direction of the normal of the intersecting structure
        // Note: The normal of the approx. surface can be max. 90deg to every surrounding normal of the structure at the intersection nodes
        const array_1d<double,3> NormalAtOneIntersectionNode = NodesOfApproximatedStructure[0].StructElemNormal;
 
        bool NormalWrongOriented = false;
        if(inner_prod(NormalAtOneIntersectionNode,Normal)<0)
            NormalWrongOriented = true;
 
        // switch direction of normal
        if(NormalWrongOriented)
            Normal *=-1;    
    
        // Compute distance values for all tet-nodes
        for(unsigned int i_TetNode = 0 ; i_TetNode < 4 ; i_TetNode++)
        {
            ElementalDistances[i_TetNode] = PointDistanceToPlane(P1, Normal, rFluidGeom[i_TetNode]);
        }
    }
 
 
    void CalcSignedDistancesToThreeIntNodes( ModelPart::ElementsContainerType::iterator& i_fluid_Element,
                                             std::vector<IntersectionNodeStruct>&        NodesOfApproximatedStructure,
                                             array_1d<double,4>&                         ElementalDistances )
    {
        Geometry< Node >& rFluidGeom = i_fluid_Element->GetGeometry();
 
        Point P1;
        Point P2;
        Point P3;
 
        P1.Coordinates() = NodesOfApproximatedStructure[0].Coordinates;
        P2.Coordinates() = NodesOfApproximatedStructure[1].Coordinates;
        P3.Coordinates() = NodesOfApproximatedStructure[2].Coordinates;
 
        array_1d<double,3> Normal;
        CalculateNormal3D(P1,P2,P3,Normal);
 
        // Check whether orientation of normal is in direction of the normal of the intersecting structure
        // Note: The normal of the approx. surface can be max. 90deg to every surrounding normal of the structure at the intersection nodes
        const array_1d<double,3> NormalAtOneIntersectionNode = NodesOfApproximatedStructure[0].StructElemNormal;
 
        bool NormalWrongOriented = false;
        if(inner_prod(NormalAtOneIntersectionNode,Normal)<0)
            NormalWrongOriented = true;
 
        // switch direction of normal
        if(NormalWrongOriented)
            Normal *=-1;
 
        // Compute distance values for all tet-nodes
        for(unsigned int i_TetNode = 0 ; i_TetNode < 4 ; i_TetNode++)
        {
            ElementalDistances[i_TetNode] = PointDistanceToPlane(P1, Normal, rFluidGeom[i_TetNode] );
        }
    }
 
 
    void CalcSignedDistancesToMoreThanThreeIntNodes(  ModelPart::ElementsContainerType::iterator& i_fluid_Element,
                                                      std::vector<IntersectionNodeStruct>         NodesOfApproximatedStructure,
                                                      array_1d<double,4>&                         ElementalDistances,
                                                      std::vector<TetEdgeStruct>&                 IntersectedTetEdges )
    {
        unsigned int numberCutEdges = NodesOfApproximatedStructure.size();
 
        // Compute average of the intersection nodes which is a node on the plane we look for
        Point P_mean;
        for(unsigned int k=0; k<numberCutEdges; k++)
            for(unsigned int i=0; i<3; i++)
                P_mean.Coordinates()[i] += NodesOfApproximatedStructure[k].Coordinates[i];
 
        for(unsigned int i=0; i<3; i++)
            P_mean.Coordinates()[i] /= numberCutEdges;
 
        // Compute normal for the best-fitted plane
        array_1d<double,3> N_mean;
 
        Matrix coordinates(numberCutEdges,3);
        for(unsigned int i=0; i<numberCutEdges; i++)
            for(unsigned int j=0; j<3; j++)
                coordinates(i,j) = NodesOfApproximatedStructure[i].Coordinates[j] - P_mean[j];
 
        Matrix A = prod(trans(coordinates),coordinates);
        Matrix V(3,3);
        Vector lambda(3);
 
        // Calculate the eigenvectors V and the corresponding eigenvalues lambda
        EigenVectors(A, V, lambda);
 
        // Look for the minimal eigenvalue all lambdas
        unsigned int min_pos = 0;
        double min_lambda = lambda[min_pos];
        for(unsigned int i=1;i<3; i++)
            if(min_lambda > lambda[i])
            {
                min_lambda = lambda[i];
                min_pos = i;
            }
 
        // the normal equals to the eigenvector which corresponds to the minimal eigenvalue
        for(unsigned int i=0;i<3; i++) N_mean[i] = V(min_pos,i);
        N_mean /= norm_2(N_mean);
 
        // Check whether orientation of normal is in direction of the normal of the intersecting structure
        // Note: The normal of the approx. surface can be max. 90deg to every surrounding normal of the structure at the intersection nodes
        array_1d<double,3> NormalAtOneIntersectionNode;
        NormalAtOneIntersectionNode = NodesOfApproximatedStructure[0].StructElemNormal;
 
        bool NormalWrongOriented = false;
        if(inner_prod(NormalAtOneIntersectionNode,N_mean)<0)
            NormalWrongOriented = true;
 
        // switch direction of normal
        if(NormalWrongOriented)
            N_mean *=-1;
 
        // Determine about the minimal distance by considering the distances to both triangles
        for(unsigned int i_TetNode = 0 ; i_TetNode < 4 ; i_TetNode++)
        {
            ElementalDistances[i_TetNode] = PointDistanceToPlane(P_mean, N_mean, i_fluid_Element->GetGeometry()[i_TetNode] );
        }
 
        // #################################################
        unsigned int numberDoubleCutEdges = 0;
        unsigned int indexDoubleCutEdge = 0;
 
        // figure out the edges which are cut more than once
        for(unsigned int i_TetEdge = 0 ; i_TetEdge < IntersectedTetEdges.size() ; i_TetEdge++)
        {
            unsigned int NumberIntNodes = IntersectedTetEdges[i_TetEdge].IntNodes.size();
            if(NumberIntNodes == 2)
            {
                numberDoubleCutEdges++;
                indexDoubleCutEdge = i_TetEdge;
            }
        }
 
        if((numberDoubleCutEdges >= 1))
        {
            array_1d<double,3> normal_1 = IntersectedTetEdges[indexDoubleCutEdge].IntNodes[0].StructElemNormal;
            array_1d<double,3> normal_2 = IntersectedTetEdges[indexDoubleCutEdge].IntNodes[1].StructElemNormal;
 
            // normalize normals
            normal_1 /= norm_2(normal_1);
            normal_2 /= norm_2(normal_2);
 
            const double pi = 3.1415926;
 
            // compute angle between normals
            double angle_n1n2 = acos( inner_prod(normal_1,normal_2) );
            // rad --> degree
            angle_n1n2 *= 180 / pi;
 
            // if angle between -60º and 120º, take the mean
            if( (angle_n1n2 > -60) && (angle_n1n2 < 120) )
            {
                // take the mean of the normals
                N_mean = 0.5 * (normal_1 + normal_2);
            }
            else
            {
                N_mean = 0.5 * (normal_1 - normal_2);
            }
 
            // Based on N_mean and P_mean compute the distances to that plane
            for(unsigned int i_TetNode = 0 ; i_TetNode < 4 ; i_TetNode++)
            {
                ElementalDistances[i_TetNode] = PointDistanceToPlane(P_mean, N_mean, i_fluid_Element->GetGeometry()[i_TetNode] );
            }
        }
    }
 
 
    double PointDistanceToPlane( Point&            planeBasePoint,
                                 array_1d<double, 3>& planeNormal,
                                 Point&            ToPoint)
    {
        // calculate vector pointing from a node in the plane (e.g. triangle point 1) to the considered node ToPoint
        array_1d<double,3> planeToPointVec = ToPoint - planeBasePoint;
 
        // projection of node on the plane
        const double sn = inner_prod(planeToPointVec,planeNormal);
        const double sd = inner_prod(planeNormal,planeNormal);
        double DistanceToPlane = sn / sqrt(sd);
 
        if( fabs(DistanceToPlane) < epsilon )
            DistanceToPlane = 0;
 
        return DistanceToPlane;
    }
 
 
    void AssignMinimalNodalDistance()
    {
        // loop over all fluid Elements
        for( ModelPart::ElementIterator i_fluid_Element = mrFluidModelPart.ElementsBegin();
             i_fluid_Element != mrFluidModelPart.ElementsEnd();
             i_fluid_Element++)
        {
            Geometry< Node >& geom = i_fluid_Element->GetGeometry();
            const Vector& ElementalDistances = i_fluid_Element->GetValue(ELEMENTAL_DISTANCES);
 
            // Assign distances to the single nodes, if a smaller value is found
            for( unsigned int i_TetNode = 0; i_TetNode < 4; i_TetNode++ )
            {
                double currentNodeDist = geom[i_TetNode].GetSolutionStepValue(DISTANCE);
                double nodeInElemDist  = ElementalDistances[i_TetNode];
 
                if( fabs( nodeInElemDist ) < fabs( currentNodeDist ) )
                    geom[i_TetNode].GetSolutionStepValue(DISTANCE) = nodeInElemDist; // overwrite nodal distance (which is global)
            } // loop i_TetNode
        } // loop i_fluidElement
    }
 
 
    void AvoidZeroDistances( ModelPart::ElementsContainerType::iterator& Element,
                             array_1d<double,4>&                         ElementalDistances)
    {
        // Assign a distance limit
        double dist_limit = 1e-5;
//         bool distChangedToLimit = false; //variable to indicate that a distance value < tolerance is set to a limit distance = tolerance
// 
//         for(unsigned int i_node = 0; i_node < 4; i_node++)
//         {
//             if(fabs(ElementalDistances[i_node]) < dist_limit)
//             {
//                 ElementalDistances[i_node] = dist_limit;
//                 distChangedToLimit = true;
//             }
//         }
// 
//         // Check, if this approach changes the split-flag (might be, that Element is not cut anymore if node with zero distance gets a positive limit distance value
//         unsigned int numberNodesPositiveDistance = 0;
//         for(unsigned int i_node = 0; i_node < 4; i_node++)
//         {
//             if((ElementalDistances[i_node]) > 0)
//                 numberNodesPositiveDistance++;
//         }
        
        for(unsigned int i_node = 0; i_node < 4; i_node++)
        {
            double & di = ElementalDistances[i_node];
            if(fabs(di) < dist_limit)
            {
                if(di >= 0) di = dist_limit;
                else di = -dist_limit;
            }
        }
 
        // Element is not set
//         if(numberNodesPositiveDistance == 4 && distChangedToLimit == true)
//             Element->GetValue(SPLIT_ELEMENT) = false;
    }
 
 
    void GenerateSkinModelPart( ModelPart& mrNewSkinModelPart )
    {
        unsigned int id_node = mrFluidModelPart.NumberOfNodes() + 1;
        unsigned int id_condition = mrFluidModelPart.NumberOfConditions() + 1;
 
        mrNewSkinModelPart.Nodes().reserve(mrFluidModelPart.Nodes().size());
        mrNewSkinModelPart.Conditions().reserve(mrFluidModelPart.Elements().size());
 
        for(ModelPart::ElementIterator i_fluid_element = mrFluidModelPart.ElementsBegin();
            i_fluid_element != mrFluidModelPart.ElementsEnd();
            i_fluid_element++)
        {
            bool is_split = i_fluid_element->Is(TO_SPLIT);
            if(is_split == true)
            {
                const Vector& distances = i_fluid_element->GetValue(ELEMENTAL_DISTANCES);
                Geometry< Node >& geom = i_fluid_element->GetGeometry();
 
                // generate the points on the edges at the zero of the distance function
                std::vector< Point > edge_points;
                edge_points.reserve(4);
 
                // loop over all 6 edges of the tetrahedra
                for(unsigned int i=0; i<3; i++)
                {
                    for(unsigned int j=i+1; j<4; j++) // go through the edges 01, 02, 03, 12, 13, 23
                    {
                        double di = distances[i];
                        double dj = distances[j];
 
                        if(di*dj < 0) //edge is cut
                        {
                            // generate point on edge by linear interpolation
                            double Ni = fabs(dj) / ( fabs(di) + fabs(dj) );
                            double Nj = 1.0 - Ni;
                            Point edge_point(Ni * geom[i] + Nj * geom[j]);
                            edge_points.push_back(edge_point);
                        }
                    }
                }
 
                // three intersection nodes
                if(edge_points.size() == 3)
                {
                    // ######## ADDING NEW NODE #########
                    Node::Pointer pnode1 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[0].X(),edge_points[0].Y(),edge_points[0].Z());
                    Node::Pointer pnode2 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[1].X(),edge_points[1].Y(),edge_points[1].Z());
                    Node::Pointer pnode3 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[2].X(),edge_points[2].Y(),edge_points[2].Z());
 
                    // ######## ADDING NEW CONDITION #########
                    //form a triangle
                    Triangle3D3< Node > triangle(pnode1, pnode2, pnode3);
 
                    Condition const& rReferenceCondition = KratosComponents<Condition>::Get("SurfaceCondition3D");
                    Properties::Pointer properties = mrNewSkinModelPart.rProperties()(0);
                    Condition::Pointer p_condition = rReferenceCondition.Create(id_condition++, triangle, properties);
 
                    mrNewSkinModelPart.Conditions().push_back(p_condition);
                }
 
                // four intersection nodes
                if(edge_points.size() == 4)
                {
                    //form a quadrilatera with the 4 cut nodes
                    array_1d<double,3> x21 = edge_points[1] - edge_points[0];
                    array_1d<double,3> x31 = edge_points[2] - edge_points[0];
                    array_1d<double,3> x41 = edge_points[3] - edge_points[0];
 
                    //define a vector oriented as x21
                    array_1d<double,3> v1 = x21 / norm_2(x21);
 
                    BoundedMatrix<double,4,3> DN_DX;
                    array_1d<double,4> msN;
                    double Area;
                    GeometryUtils::CalculateGeometryData( geom, DN_DX, msN, Area );
 
                    array_1d<double,3> n = prod(trans(DN_DX),distances);
                    n /= norm_2(n);
 
                    array_1d<double,3> v2;
                    MathUtils<double>::CrossProduct(v2,v1,n); // v2 = v1 x n
 
                    array_1d<double,3> angles;
                    angles[0] = 0.0; //angle between x21 and v1
                    angles[1] = atan2( inner_prod(x31,v2), inner_prod(x31,v1) ); //angle between x31 and v1
                    angles[2] = atan2( inner_prod(x41,v2), inner_prod(x41,v1) ); //angle between x31 and v1
 
                    double max_angle = 0.0;
                    double min_angle = 0.0;
                    unsigned int min_pos = 1;
                    unsigned int max_pos = 1;
                    for(unsigned int i=1; i<3; i++)
                    {
                        if(angles[i] < min_angle)
                        {
                            min_pos = i+1; //this is the local index of the edge point which forms the minimal angle
                            min_angle = angles[i];
                        }
                        else if(angles[i] > max_angle)
                        {
                            max_pos = i+1; //this is the local index of the edge point which forms the maximal angle
                            max_angle = angles[i];
                        }
                    }
 
                    //find the pos of the center node
                    unsigned int center_pos = 0;
                    for(unsigned int i=1; i<4; i++)
                    {
                        if((i!= min_pos) && (i!=max_pos))
                        { center_pos = i; }
                    }
 
                    // ######## ADDING NEW NODE #########
                    Node::Pointer pnode1 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[0].X(),edge_points[0].Y(),edge_points[0].Z());
                    Node::Pointer pnode2 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[min_pos].X(),edge_points[min_pos].Y(),edge_points[min_pos].Z());
                    Node::Pointer pnode3 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[center_pos].X(),edge_points[center_pos].Y(),edge_points[center_pos].Z());
                    Node::Pointer pnode4 = mrNewSkinModelPart.CreateNewNode(id_node++,edge_points[max_pos].X(),edge_points[max_pos].Y(),edge_points[max_pos].Z());
 
                    // ######## ADDING NEW CONDITION #########
                    //form two triangles
                    Triangle3D3< Node > triangle1(pnode1, pnode2, pnode3);
                    Triangle3D3< Node > triangle2(pnode1, pnode3, pnode4);
 
                    Condition const& rReferenceCondition = KratosComponents<Condition>::Get("SurfaceCondition3D");
                 
                    Properties::Pointer properties = mrNewSkinModelPart.rProperties()(0);
 
                    Condition::Pointer p_condition1 = rReferenceCondition.Create(id_condition++, triangle1, properties);
                    Condition::Pointer p_condition2 = rReferenceCondition.Create(id_condition++, triangle2, properties);
 
                    mrNewSkinModelPart.Conditions().push_back(p_condition1);
                    mrNewSkinModelPart.Conditions().push_back(p_condition2);
                    
                }
 
            }
        }
 
    }
 
 
    void GenerateOctree()
    {
        Timer::Start("Generating Octree");
        //std::cout << "Generating the Octree..." << std::endl;
        auto temp_octree =  Kratos::make_shared<OctreeType>();
        //OctreeType::Pointer temp_octree = OctreeType::Pointer(new OctreeType() );
        mpOctree.swap(temp_octree);
        
        double low[3];
        double high[3];
        
        for (int i = 0 ; i < 3; i++)
        {
            low[i] = high[i] = mrFluidModelPart.NodesBegin()->Coordinates()[i];
        }
        
        // loop over all nodes in the bounding box
        for(ModelPart::NodeIterator i_node = mrFluidModelPart.NodesBegin();
            i_node != mrFluidModelPart.NodesEnd();
            i_node++)
        {
            const array_1d<double,3>& r_coordinates = i_node->Coordinates();
            for (int i = 0 ; i < 3; i++)
            {
                low[i]  = r_coordinates[i] < low[i]  ? r_coordinates[i] : low[i];
                high[i] = r_coordinates[i] > high[i] ? r_coordinates[i] : high[i];
            }
        }
        
        // loop over all skin nodes
        for(ModelPart::NodeIterator i_node = mrSkinModelPart.NodesBegin();
            i_node != mrSkinModelPart.NodesEnd();
            i_node++)
        {
            const array_1d<double,3>& r_coordinates = i_node->Coordinates();
            for (int i = 0 ; i < 3; i++)
            {
                low[i]  = r_coordinates[i] < low[i]  ? r_coordinates[i] : low[i];
                high[i] = r_coordinates[i] > high[i] ? r_coordinates[i] : high[i];
            }
        }
                
        mpOctree->SetBoundingBox(low,high);
 
        //mpOctree->RefineWithUniformSize(0.0625);
 
        // loop over all structure nodes
        for(ModelPart::NodeIterator i_node = mrSkinModelPart.NodesBegin();
            i_node != mrSkinModelPart.NodesEnd();
            i_node++)
        {
            double temp_point[3];
            temp_point[0] = i_node->X();
            temp_point[1] = i_node->Y();
            temp_point[2] = i_node->Z();
            mpOctree->Insert(temp_point);
        }
 
        //mpOctree->Constrain2To1(); // To be removed. Pooyan.
 
        // loop over all structure elements
        for(ModelPart::ElementIterator i_element = mrSkinModelPart.ElementsBegin();
            i_element != mrSkinModelPart.ElementsEnd();
            i_element++)
        {
            mpOctree->Insert(*(i_element).base());
        }
 
        Timer::Stop("Generating Octree");
 
//        KRATOS_WATCH(mpOctree);
 
//        std::cout << "######## WRITING OCTREE MESH #########" << std::endl;
//        std::ofstream myfile;
//        myfile.open ("octree.post.msh");
//        mpOctree.PrintGiDMesh(myfile);
//        myfile.close();
 
        //std::cout << "Generating the Octree finished" << std::endl;
    }
 
 
    void GenerateNodes()
    {
        Timer::Start("Generating Nodes");
        std::vector<OctreeType::cell_type*> all_leaves;
        mpOctree->GetAllLeavesVector(all_leaves);
 
        int leaves_size = all_leaves.size();
 
#pragma omp parallel for
        for (int i = 0; i < leaves_size; i++)
        {
            *(all_leaves[i]->pGetDataPointer()) = ConfigurationType::AllocateData();
        }
 
 
        std::size_t last_id = mrBodyModelPart.NumberOfNodes() + 1;
 
        for (std::size_t i = 0; i < all_leaves.size(); i++)
        {
                    CellType* cell = all_leaves[i];
            GenerateCellNode(cell, last_id);
        }
 
        Timer::Stop("Generating Nodes");
 
    }
 
 
    void GenerateCellNode(CellType* pCell, std::size_t& LastId)
    {
        for (int i_pos=0; i_pos < 8; i_pos++) // position 8 is for center
        {
            DistanceSpatialContainersConfigure::cell_node_data_type* p_node = (*(pCell->pGetData()))[i_pos];
            if(p_node == 0)
            {
                (*(pCell->pGetData()))[i_pos] = new DistanceSpatialContainersConfigure::cell_node_data_type;
 
                (*(pCell->pGetData()))[i_pos]->Id() = LastId++;
 
                        mOctreeNodes.push_back((*(pCell->pGetData()))[i_pos]);
 
                SetNodeInNeighbours(pCell,i_pos,(*(pCell->pGetData()))[i_pos]);
            }
 
        }
    }
 
 
    void SetNodeInNeighbours(CellType* pCell, int Position, CellNodeDataType* pNode)
    {
        CellType::key_type point_key[3];
        pCell->GetKey(Position, point_key);
 
        for (std::size_t i_direction = 0; i_direction < 8; i_direction++) {
            CellType::key_type neighbour_key[3];
            if (pCell->GetNeighbourKey(Position, i_direction, neighbour_key)) {
                CellType* neighbour_cell = mpOctree->pGetCell(neighbour_key);
                if (!neighbour_cell || (neighbour_cell == pCell))
                    continue;
 
                std::size_t position = neighbour_cell->GetLocalPosition(point_key);
                if((*neighbour_cell->pGetData())[position])
                {
                    std::cout << "ERROR!! Bad Position calculated!!!!!!!!!!! position :" << position << std::endl;
                    continue;
                }
 
                (*neighbour_cell->pGetData())[position] = pNode;
            }
        }
    }
 
 
    void CalculateDistance2()
    {
        Timer::Start("Calculate Distances2");
        ModelPart::NodesContainerType::ContainerType& nodes = mrFluidModelPart.NodesArray();
        int nodes_size = nodes.size();
        //         // first of all we reset the node distance to 1.00 which is the maximum distnace in our normalized space.
        //#pragma omp parallel for firstprivate(nodes_size)
        //         for(int i = 0 ; i < nodes_size ; i++)
        //             nodes[i]->GetSolutionStepValue(DISTANCE) = 1.00;
 
        std::vector<CellType*> leaves;
 
        mpOctree->GetAllLeavesVector(leaves);
        //int leaves_size = leaves.size();
 
        //         for(int i = 0 ; i < leaves_size ; i++)
        //             CalculateNotEmptyLeavesDistance(leaves[i]);
 
#pragma omp parallel for firstprivate(nodes_size)
        for(int i = 0 ; i < nodes_size ; i++)
        {
            CalculateNodeDistance(*(nodes[i]));
        }
        Timer::Stop("Calculate Distances2");
 
    }
 
 
    //     void CalculateDistance3()
    //     {
    //         Timer::Start("Calculate Distances2");
    //         ModelPart::NodesContainerType::ContainerType& nodes = mrFluidModelPart.NodesArray();
    //         int nodes_size = nodes.size();
    //#pragma omp parallel for firstprivate(nodes_size)
    //         for(int i = 0 ; i < nodes_size ; i++)
    //             nodes[i]->GetSolutionStepValue(DISTANCE) = 1.00;
 
    //           std::vector<CellType*> leaves;
 
    //         mpOctree->GetAllLeavesVector(leaves);
    //         int leaves_size = leaves.size();
 
    //         for(int i = 0 ; i < leaves_size ; i++)
    //             CalculateNotEmptyLeavesDistance(leaves[i]);
 
    //#pragma omp parallel for firstprivate(nodes_size)
    //         for(int i = 0 ; i < nodes_size ; i++)
    //         {
    //             CalculateNodeDistance(*(nodes[i]));
    //         }
    //         Timer::Stop("Calculate Distances2");
 
    //     }
    //     void CalculateDistance4()
    //     {
    //         Timer::Start("Calculate Distances3");
    //         ModelPart::NodesContainerType::ContainerType& nodes = mrFluidModelPart.NodesArray();
    //         int nodes_size = nodes.size();
    //           std::vector<CellType*> leaves;
 
    //         mpOctree->GetAllLeavesVector(leaves);
    //         int leaves_size = leaves.size();
 
    //#pragma omp parallel for firstprivate(nodes_size)
    //         for(int i = 0 ; i < nodes_size ; i++)
    //         {
    //             CalculateNodeDistanceFromCell(*(nodes[i]));
    //         }
    //         Timer::Stop("Calculate Distances3");
 
    //     }
 
    void CalculateDistance()
    {
        Timer::Start("Calculate Distances");
        DistanceSpatialContainersConfigure::data_type& nodes = mOctreeNodes;
        int nodes_size = nodes.size();
        // first of all we reste the node distance to 1.00 which is the maximum distnace in our normalized space.
#pragma omp parallel for firstprivate(nodes_size)
        for(int i = 0 ; i < nodes_size ; i++)
            nodes[i]->Distance() = 1.00;
 
 
        std::vector<CellType*> leaves;
 
        mpOctree->GetAllLeavesVector(leaves);
        int leaves_size = leaves.size();
 
        for(int i = 0 ; i < leaves_size ; i++)
            CalculateNotEmptyLeavesDistance(leaves[i]);
 
        for(int i_direction = 0 ; i_direction < 1 ; i_direction++)
        {
 
            //#pragma omp parallel for firstprivate(nodes_size)
            for(int i = 0 ; i < nodes_size ; i++)
            {
                if(nodes[i]->X() < 1.00 && nodes[i]->Y() < 1.00 && nodes[i]->Z() < 1.00)
                    // if((*nodes[i])[i_direction] == 0.00)
                    CalculateDistance(*(nodes[i]), i_direction);
            }
        }
        Timer::Stop("Calculate Distances");
 
    }
 
    void CalculateDistance(CellNodeDataType& rNode, int i_direction)
    {
        double coords[3] = {rNode.X(), rNode.Y(), rNode.Z()};
        // KRATOS_WATCH_3(coords);
 
        //This function must color the positions in space defined by 'coords'.
        //coords is of dimension (3) normalized in (0,1)^3 space
 
        typedef Element::GeometryType triangle_type;
        typedef std::vector<std::pair<double, triangle_type*> > intersections_container_type;
 
        intersections_container_type intersections;
        DistanceSpatialContainersConfigure::data_type nodes_array;
 
 
        const double epsilon = 1e-12;
 
        double distance = 1.0;
 
        // Creating the ray
        double ray[3] = {coords[0], coords[1], coords[2]};
        
        mpOctree->NormalizeCoordinates(ray);
        ray[i_direction] = 0; // starting from the lower extreme
 
        //            KRATOS_WATCH_3(ray)
        GetIntersectionsAndNodes(ray, i_direction, intersections, nodes_array);
        //            KRATOS_WATCH(nodes_array.size())
        for (std::size_t i_node = 0; i_node < nodes_array.size() ; i_node++)
        {
            double coord = (*nodes_array[i_node])[i_direction];
            //             KRATOS_WATCH(intersections.size());
 
            int ray_color= 1;
            std::vector<std::pair<double, Element::GeometryType*> >::iterator i_intersection = intersections.begin();
            while (i_intersection != intersections.end()) {
                double d = coord - i_intersection->first;
                if (d > epsilon) {
 
                    ray_color = -ray_color;
                    distance = d;
                } else if (d > -epsilon) {//interface
                    distance = 0.00;
                    break;
                } else {
                    if(distance > -d)
                        distance = -d;
                    break;
                }
 
                i_intersection++;
            }
 
            distance *= ray_color;
 
            double& node_distance = nodes_array[i_node]->Distance();
            if(fabs(distance) < fabs(node_distance))
                node_distance = distance;
            else if (distance*node_distance < 0.00) // assigning the correct sign
                node_distance = -node_distance;
 
 
        }
    }
 
    void CalculateNotEmptyLeavesDistance(CellType* pCell)
    {
      //typedef Element::GeometryType triangle_type;
        typedef OctreeType::cell_type::object_container_type object_container_type;
 
        object_container_type* objects = (pCell->pGetObjects());
 
        // There are no intersection in empty cells
        if (objects->empty())
            return;
 
 
        for (int i_pos=0; i_pos < 8; i_pos++) // position 8 is for center
        {
            double distance = 1.00; // maximum distance is 1.00
 
            for(object_container_type::iterator i_object = objects->begin(); i_object != objects->end(); i_object++)
            {
                CellType::key_type keys[3];
                pCell->GetKey(i_pos,keys);
 
                double cell_point[3];
                mpOctree->CalculateCoordinates(keys,cell_point);
 
                double d = GeometryUtils::PointDistanceToTriangle3D((*i_object)->GetGeometry()[0], (*i_object)->GetGeometry()[1], (*i_object)->GetGeometry()[2], Point(cell_point[0], cell_point[1], cell_point[2]));
 
                if(d < distance)
                    distance = d;
            }
 
            double& node_distance = (*(pCell->pGetData()))[i_pos]->Distance();
            if(distance < node_distance)
                node_distance = distance;
 
        }
 
    }
 
    void CalculateNodeDistance(Node& rNode)
    {
        double coord[3] = {rNode.X(), rNode.Y(), rNode.Z()};
        double distance = DistancePositionInSpace(coord);
        double& node_distance =  rNode.GetSolutionStepValue(DISTANCE);
 
        //const double epsilon = 1.00e-12;
        //if(fabs(node_distance) > fabs(distance))
        //    node_distance = distance;
        /*else*/ if (distance*node_distance < 0.00) // assigning the correct sign
            node_distance = -node_distance;
    }
 
    //      void CalculateNodeDistanceFromCell(Node& rNode)
    //      {
    //          OctreeType::key_type node_key[3] = {octree->CalcKeyNormalized(rNode.X()), octree->CalcKeyNormalized(rNode.Y()), octree->CalcKeyNormalized(rNode.Z())};
    //          OctreeType::cell_type* pcell = octree->pGetCell(node_key);
 
    //          object_container_type* objects = (pCell->pGetObjects());
 
    //          // We interpolate the cell distances for the node in empty cells
    //          if (objects->empty())
    //          {
 
    //          }
 
    //          double distance = DistancePositionInSpace(coord);
    //          double& node_distance =  rNode.GetSolutionStepValue(DISTANCE);
 
    //          //const double epsilon = 1.00e-12;
    //          if(fabs(node_distance) > fabs(distance))
    //            node_distance = distance;
    //          else if (distance*node_distance < 0.00) // assigning the correct sign
    //              node_distance = -node_distance;
    //      }
 
    double DistancePositionInSpace(double* coords)
    {
        //This function must color the positions in space defined by 'coords'.
        //coords is of dimension (3) normalized in (0,1)^3 space
 
        typedef Element::GeometryType triangle_type;
        typedef std::vector<std::pair<double, triangle_type*> > intersections_container_type;
 
        intersections_container_type intersections;
 
        const int dimension = 3;
        const double epsilon = 1e-12;
 
        double distances[3] = {1.0, 1.0, 1.0};
 
        for (int i_direction = 0; i_direction < dimension; i_direction++)
        {
            // Creating the ray
            double ray[3] = {coords[0], coords[1], coords[2]};
            
            mpOctree->NormalizeCoordinates(ray);
            ray[i_direction] = 0; // starting from the lower extreme
 
            GetIntersections(ray, i_direction, intersections);
 
            //                if(intersections.size() == 1)
            //                    KRATOS_WATCH_3(ray)
 
            //             KRATOS_WATCH(intersections.size());
 
            int ray_color= 1;
            std::vector<std::pair<double, Element::GeometryType*> >::iterator i_intersection = intersections.begin();
            while (i_intersection != intersections.end()) {
                double d = coords[i_direction] - i_intersection->first;
                if (d > epsilon) {
 
                    ray_color = -ray_color;
                    distances[i_direction] = d;
                    //                        if(distances[i_direction] > d) // I think this is redundunt. Pooyan.
                    //                        {
                    //                            if(ray_color > 0.00)
                    //                                    distances[i_direction] = d;
                    //                            else
                    //                                distances[i_direction] = -d;
                    //                        }
                } else if (d > -epsilon) {//interface
                    distances[i_direction] = 0.00;
                    break;
                } else {
                    if(distances[i_direction] > -d)
                        distances[i_direction] = -d;
                    break;
                }
 
                i_intersection++;
            }
 
            distances[i_direction] *= ray_color;
        }
 
        //            if(distances[0]*distances[1] < 0.00 || distances[2]*distances[1] < 0.00)
        //                KRATOS_WATCH_3(distances);
 
//#ifdef _DEBUG
//        std::cout << "colors : " << colors[0] << ", " << colors[1] << ", " << colors[2] << std::endl;
//#endif
        double distance = (fabs(distances[0]) > fabs(distances[1])) ? distances[1] : distances[0];
        distance = (fabs(distance) > fabs(distances[2])) ? distances[2] : distance;
 
        return distance;
 
    }
 
    void GetIntersectionsAndNodes(double* ray, int direction, std::vector<std::pair<double,Element::GeometryType*> >& intersections, DistanceSpatialContainersConfigure::data_type& rNodesArray)
    {
        //This function passes the ray through the model and gives the hit point to all objects in its way
        //ray is of dimension (3) normalized in (0,1)^3 space
        // direction can be 0,1,2 which are x,y and z respectively
 
        const double epsilon = 1.00e-12;
 
        // first clearing the intersections points vector
        intersections.clear();
 
        //OctreeType* octree = &mOctree;
        OctreeType* octree = mpOctree.get();
 
        OctreeType::key_type ray_key[3] = {octree->CalcKeyNormalized(ray[0]), octree->CalcKeyNormalized(ray[1]), octree->CalcKeyNormalized(ray[2])};
        OctreeType::key_type cell_key[3];
 
        // getting the entrance cell from lower extreme
        ray_key[direction] = 0;
        OctreeType::cell_type* cell = octree->pGetCell(ray_key);
 
        while (cell) {
            std::size_t position = cell->GetLocalPosition(ray_key); // Is this the local position!?!?!?!
            OctreeType::key_type node_key[3];
            cell->GetKey(position, node_key);
            if((node_key[0] == ray_key[0]) && (node_key[1] == ray_key[1]) && (node_key[2] == ray_key[2]))
            {
                if(cell->pGetData())
                {
                    if(cell->pGetData()->size() > position)
                    {
                        CellNodeDataType* p_node = (*cell->pGetData())[position];
                        if(p_node)
                        {
                            //KRATOS_WATCH(p_node->Id())
                            rNodesArray.push_back(p_node);
                        }
                    }
                    else
                        KRATOS_WATCH(cell->pGetData()->size())
                }
            }
 
 
            //        std::cout << ".";
            GetCellIntersections(cell, ray, ray_key, direction, intersections);
 
            // Add the cell's middle node if existed
            //      cell->GetKey(8, cell_key); // 8 is the central position
            //      ray_key[direction]=cell_key[direction]; // positioning the ray in the middle of cell in its direction
 
            //      position = cell->GetLocalPosition(ray_key);
            //      if(position < 27) // principal nodes
            //      {
            //          if(cell->pGetData())
            //          {
            //              if(cell->pGetData()->size() > position)
            //              {
            //                  Node* p_node = (*cell->pGetData())[position];
            //                  if(p_node)
            //                  {
            //                      //KRATOS_WATCH(p_node->Id())
            //                      rNodesArray.push_back(p_node);
            //                  }
            //              }
            //              else
            //                  KRATOS_WATCH(cell->pGetData()->size())
            //          }
            //      }
            //      else
            //      {
            //          KRATOS_WATCH(position);
            //          KRATOS_WATCH(*cell);
            //      }
 
 
            // go to the next cell
            if (cell->GetNeighbourKey(1 + direction * 2, cell_key)) {
                ray_key[direction] = cell_key[direction];
                cell = octree->pGetCell(ray_key);
                ray_key[direction] -= 1 ;//the key returned by GetNeighbourKey is inside the cell (minkey +1), to ensure that the corresponding
                //cell get in pGetCell is the right one.
//#ifdef _DEBUG
//                Octree_Pooyan::key_type min_key[3];
//                cell->GetMinKey(min_key[0],min_key[1],min_key[2]);
//                Octree_Pooyan::key_type tmp;
//                tmp= min_key[direction];
//                assert(ray_key[direction]==tmp);
//#endif
            } else
                cell = NULL;
        }
 
 
 
        //   KRATOS_WATCH(rNodesArray.size());
        // now eliminating the repeated objects
        if (!intersections.empty()) {
            //sort
            std::sort(intersections.begin(), intersections.end());
            // unique
            std::vector<std::pair<double, Element::GeometryType*> >::iterator i_begin = intersections.begin();
            std::vector<std::pair<double, Element::GeometryType*> >::iterator i_intersection = intersections.begin();
            while (++i_begin != intersections.end()) {
                // considering the very near points as the same points
                if (fabs(i_begin->first - i_intersection->first) > epsilon) // if the hit points are far enough they are not the same
                    *(++i_intersection) = *i_begin;
            }
            intersections.resize((++i_intersection) - intersections.begin());
 
        }
    }
 
    void GetIntersections(double* ray, int direction, std::vector<std::pair<double,Element::GeometryType*> >& intersections)
    {
        //This function passes the ray through the model and gives the hit point to all objects in its way
        //ray is of dimension (3) normalized in (0,1)^3 space
        // direction can be 0,1,2 which are x,y and z respectively
 
        const double epsilon = 1.00e-12;
 
        // first clearing the intersections points vector
        intersections.clear();
 
        //OctreeType* octree = &mOctree;
        OctreeType* octree = mpOctree.get();
 
        OctreeType::key_type ray_key[3] = {octree->CalcKeyNormalized(ray[0]), octree->CalcKeyNormalized(ray[1]), octree->CalcKeyNormalized(ray[2])};
        OctreeType::key_type cell_key[3];
 
        // getting the entrance cell from lower extreme
        OctreeType::cell_type* cell = octree->pGetCell(ray_key);
 
        while (cell) {
            //        std::cout << ".";
            GetCellIntersections(cell, ray, ray_key, direction, intersections);
            // go to the next cell
            if (cell->GetNeighbourKey(1 + direction * 2, cell_key)) {
                ray_key[direction] = cell_key[direction];
                cell = octree->pGetCell(ray_key);
                ray_key[direction] -= 1 ;//the key returned by GetNeighbourKey is inside the cell (minkey +1), to ensure that the corresponding
                //cell get in pGetCell is the right one.
//#ifdef _DEBUG
//                Octree_Pooyan::key_type min_key[3];
//                cell->GetMinKey(min_key[0],min_key[1],min_key[2]);
//                Octree_Pooyan::key_type tmp;
//                tmp= min_key[direction];
//                assert(ray_key[direction]==tmp);
//#endif
            } else
                cell = NULL;
        }
 
 
        // now eliminating the repeated objects
        if (!intersections.empty()) {
            //sort
            std::sort(intersections.begin(), intersections.end());
            // unique
            std::vector<std::pair<double, Element::GeometryType*> >::iterator i_begin = intersections.begin();
            std::vector<std::pair<double, Element::GeometryType*> >::iterator i_intersection = intersections.begin();
            while (++i_begin != intersections.end()) {
                // considering the very near points as the same points
                if (fabs(i_begin->first - i_intersection->first) > epsilon) // if the hit points are far enough they are not the same
                    *(++i_intersection) = *i_begin;
            }
            intersections.resize((++i_intersection) - intersections.begin());
 
        }
    }
 
    int GetCellIntersections(OctreeType::cell_type* cell, double* ray,
                             OctreeType::key_type* ray_key, int direction,
                             std::vector<std::pair<double, Element::GeometryType*> >& intersections)  {
        //This function passes the ray through the cell and gives the hit point to all objects in its way
        //ray is of dimension (3) normalized in (0,1)^3 space
        // direction can be 0,1,2 which are x,y and z respectively
 
        //typedef Element::GeometryType triangle_type;
        typedef OctreeType::cell_type::object_container_type object_container_type;
 
        object_container_type* objects = (cell->pGetObjects());
 
        // There are no intersection in empty cells
        if (objects->empty())
            return 0;
 
        //      std::cout << "X";
        // calculating the two extreme of the ray segment inside the cell
        double ray_point1[3] = {ray[0], ray[1], ray[2]};
        double ray_point2[3] = {ray[0], ray[1], ray[2]};
        double normalized_coordinate;
        mpOctree->CalculateCoordinateNormalized(ray_key[direction], normalized_coordinate);
        ray_point1[direction] = normalized_coordinate;
        ray_point2[direction] = ray_point1[direction] + mpOctree->CalcSizeNormalized(cell);
 
        mpOctree->ScaleBackToOriginalCoordinate(ray_point1);
        mpOctree->ScaleBackToOriginalCoordinate(ray_point2);
 
        for (object_container_type::iterator i_object = objects->begin(); i_object != objects->end(); i_object++) {
            double intersection[3]={0.00,0.00,0.00};
 
            int is_intersected = IntersectionTriangleSegment((*i_object)->GetGeometry(), ray_point1, ray_point2, intersection); // This intersection has to be optimized for axis aligned rays
 
            if (is_intersected == 1) // There is an intersection but not coplanar
                intersections.push_back(std::pair<double, Element::GeometryType*>(intersection[direction], &((*i_object)->GetGeometry())));
            //else if(is_intersected == 2) // coplanar case
        }
 
        return 0;
    }
 
    int IntersectionTriangleSegment(Element::GeometryType& rGeometry, double* RayPoint1, double* RayPoint2, double* IntersectionPoint)
    {
        // This is the adaption of the implemnetation provided in:
        // http://www.softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm#intersect_RayTriangle()
 
        const double epsilon = 1.00e-12;
 
        array_1d<double,3>    u, v, n;             // triangle vectors
        array_1d<double,3>    dir, w0, w;          // ray vectors
        double     r, a, b;             // params to calc ray-plane intersect
 
 
        // get triangle edge vectors and plane normal
        u = rGeometry[1] - rGeometry[0];
        v = rGeometry[2] - rGeometry[0];
 
        MathUtils<double>::CrossProduct(n, u, v);             // cross product
 
        if (norm_2(n) == 0)            // triangle is degenerate
            return -1;                 // do not deal with this case
 
        double triangle_origin_distance = -inner_prod(n, rGeometry[0]);
        Point ray_point_1, ray_point_2;
        
        for(int i = 0 ; i < 3 ; i++)
        {
            dir[i] = RayPoint2[i] - RayPoint1[i];             // ray direction vector
            w0[i] = RayPoint1[i] - rGeometry[0][i];
            ray_point_1[i] = RayPoint1[i];
            ray_point_2[i] = RayPoint2[i];
        }
 
        double sign_distance_1 = inner_prod(n, ray_point_1) + triangle_origin_distance;
        double sign_distance_2 = inner_prod(n, ray_point_2) + triangle_origin_distance;
 
        if (sign_distance_1*sign_distance_2 > epsilon) // segment line point on the same side of plane
            return 0;
        a = -inner_prod(n,w0);
        b = inner_prod(n,dir);
 
        if (fabs(b) < epsilon) {     // ray is parallel to triangle plane
            if (a == 0)                // ray lies in triangle plane
                return 2;
            else return 0;             // ray disjoint from plane
        }
 
        // get intersect point of ray with triangle plane
        r = a / b;
        if (r < 0.0)                   // ray goes away from triangle
            return 0;                  // => no intersect
        // for a segment, also test if (r > 1.0) => no intersect
 
        for(int i = 0 ; i < 3 ; i++)
            IntersectionPoint[i]  = RayPoint1[i] + r * dir[i];           // intersect point of ray and plane
 
        // is I inside T?
        double    uu, uv, vv, wu, wv, D;
        uu = inner_prod(u,u);
        uv = inner_prod(u,v);
        vv = inner_prod(v,v);
 
 
        for(int i = 0 ; i < 3 ; i++)
            w[i] = IntersectionPoint[i] - rGeometry[0][i];
 
 
        wu = inner_prod(w,u);
        wv = inner_prod(w,v);
        D = uv * uv - uu * vv;
 
        // get and test parametric coords
        double s, t;
        s = (uv * wv - vv * wu) / D;
        if (s < 0.0 - epsilon || s > 1.0 + epsilon)        // I is outside T
            return 0;
        t = (uv * wu - uu * wv) / D;
        if (t < 0.0 - epsilon || (s + t) > 1.0 + epsilon)  // I is outside T
            return 0;
 
        return 1;                      // I is in T
 
    }
 
 
 
 
 
 
    std::string Info() const override
    {
        return "CalculateSignedDistanceTo3DSkinProcess";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "CalculateSignedDistanceTo3DSkinProcess";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
    }
 
    void PrintGiDMesh(std::ostream & rOStream) const {
        std::vector<CellType*> leaves;
 
        mpOctree->GetAllLeavesVector(leaves);
 
        std::cout << "writing " << leaves.size() << " leaves" << std::endl;
        rOStream << "MESH \"leaves\" dimension 3 ElemType Hexahedra Nnode 8" << std::endl;
        rOStream << "# color 96 96 96" << std::endl;
        rOStream << "Coordinates" << std::endl;
        rOStream << "# node number coordinate_x coordinate_y coordinate_z  " << std::endl;
 
        for(DistanceSpatialContainersConfigure::data_type::const_iterator i_node = mOctreeNodes.begin() ; i_node != mOctreeNodes.end() ; i_node++)
        {
            rOStream << (*i_node)->Id() << "  " << (*i_node)->X() << "  " << (*i_node)->Y() << "  " << (*i_node)->Z() << std::endl;
            //mpOctree->Insert(temp_point);
        }
        std::cout << "Nodes written..." << std::endl;
        rOStream << "end coordinates" << std::endl;
        rOStream << "Elements" << std::endl;
        rOStream << "# Element node_1 node_2 node_3 material_number" << std::endl;
 
        for (std::size_t i = 0; i < leaves.size(); i++) {
            if ((leaves[i]->pGetData()))
            {
                DistanceSpatialContainersConfigure::data_type& nodes = (*(leaves[i]->pGetData()));
 
                rOStream << i + 1;
                for(int j = 0 ; j < 8 ; j++)
                    rOStream << "  " << nodes[j]->Id();
                rOStream << std::endl;
            }
        }
        rOStream << "end Elements" << std::endl;
 
    }
 
    void PrintGiDResults(std::ostream & rOStream) const {
        std::vector<CellType*> leaves;
 
        mpOctree->GetAllLeavesVector(leaves);
 
        rOStream << "GiD Post Results File 1.0" << std::endl << std::endl;
 
        rOStream << "Result \"Distance\" \"Kratos\" 1 Scalar OnNodes" << std::endl;
 
        rOStream << "Values" << std::endl;
 
        for(DistanceSpatialContainersConfigure::data_type::const_iterator i_node = mOctreeNodes.begin() ; i_node != mOctreeNodes.end() ; i_node++)
        {
            rOStream << (*i_node)->Id() << "  " << (*i_node)->Distance() << std::endl;
        }
        rOStream << "End Values" << std::endl;
 
    }
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
    ModelPart& mrSkinModelPart;
    ModelPart& mrBodyModelPart;
    ModelPart& mrFluidModelPart;
 
    DistanceSpatialContainersConfigure::data_type mOctreeNodes;
 
    Kratos::shared_ptr<OctreeType> mpOctree;
 
    static const double epsilon;
 
    static inline void EigenVectors(const Matrix& A, Matrix& vectors, Vector& lambda, double zero_tolerance =1e-9, int max_iterations = 10)
    {
        Matrix Help= A;
 
        for(int i=0; i<3; i++)
            for(int j=0; j<3; j++)
                Help(i,j)= Help(i,j);
 
 
        vectors.resize(Help.size1(),Help.size2(),false);
 
        lambda.resize(Help.size1(),false);
 
        Matrix HelpDummy(Help.size1(),Help.size2());
 
        bool is_converged = false;
 
        Matrix unity=ZeroMatrix(Help.size1(),Help.size2());
 
        for(unsigned int i=0; i< Help.size1(); i++)
            unity(i,i)= 1.0;
 
        Matrix V= unity;
 
        Matrix VDummy(Help.size1(),Help.size2());
 
        Matrix Rotation(Help.size1(),Help.size2());
 
 
        for(int iterations=0; iterations<max_iterations; iterations++)
        {
 
            is_converged= true;
 
            double a= 0.0;
 
            unsigned int index1= 0;
 
            unsigned int index2= 1;
 
            for(unsigned int i=0; i< Help.size1(); i++)
            {
                for(unsigned int j=(i+1); j< Help.size2(); j++)
                {
                    if((fabs(Help(i,j)) > a ) && (fabs(Help(i,j)) > zero_tolerance))
                    {
                        a= fabs(Help(i,j));
 
                        index1= i;
                        index2= j;
 
                        is_converged= false;
                    }
                }
            }
 
            //                 KRATOS_WATCH(Help);
 
            if(is_converged)
                break;
 
            //Calculation of Rotationangle
 
            double gamma= (Help(index2,index2)-Help(index1,index1))/(2*Help(index1,index2));
 
            double u=1.0;
 
            if(fabs(gamma) > zero_tolerance && fabs(gamma)< (1/zero_tolerance))
            {
                u= gamma/fabs(gamma)*1.0/(fabs(gamma)+sqrt(1.0+gamma*gamma));
            }
            else
            {
                if  (fabs(gamma)>= (1.0/zero_tolerance))
                    u= 0.5/gamma;
            }
 
            double c= 1.0/(sqrt(1.0+u*u));
 
            double s= c*u;
 
            double teta= s/(1.0+c);
 
            //Ratotion of the Matrix
            HelpDummy= Help;
 
            HelpDummy(index2,index2)= Help(index2,index2)+u*Help(index1,index2);
            HelpDummy(index1,index1)= Help(index1,index1)-u*Help(index1,index2);
            HelpDummy(index1,index2)= 0.0;
            HelpDummy(index2,index1)= 0.0;
 
            for(unsigned int i=0; i<Help.size1(); i++)
            {
                if((i!= index1) && (i!= index2))
                {
                    HelpDummy(index2,i)=Help(index2,i)+s*(Help(index1,i)- teta*Help(index2,i));
                    HelpDummy(i,index2)=Help(index2,i)+s*(Help(index1,i)- teta*Help(index2,i));
 
                    HelpDummy(index1,i)=Help(index1,i)-s*(Help(index2,i)+ teta*Help(index1,i));
                    HelpDummy(i,index1)=Help(index1,i)-s*(Help(index2,i)+ teta*Help(index1,i));
                }
            }
 
 
            Help= HelpDummy;
 
            //Calculation of the eigenvectors V
            Rotation =unity;
            Rotation(index2,index1)=-s;
            Rotation(index1,index2)=s;
            Rotation(index1,index1)=c;
            Rotation(index2,index2)=c;
 
            //                 Help=ZeroMatrix(A.size1(),A.size1());
 
            VDummy = ZeroMatrix(Help.size1(), Help.size2());
 
            for(unsigned int i=0; i< Help.size1(); i++)
            {
                for(unsigned int j=0; j< Help.size1(); j++)
                {
                    for(unsigned int k=0; k< Help.size1(); k++)
                    {
                        VDummy(i,j) += V(i,k)*Rotation(k,j);
                    }
                }
            }
            V= VDummy;
        }
 
        if(!(is_converged))
        {
            std::cout<<"########################################################"<<std::endl;
            std::cout<<"Max_Iterations exceed in Jacobi-Seidel-Iteration (eigenvectors)"<<std::endl;
            std::cout<<"########################################################"<<std::endl;
        }
 
        for(unsigned int i=0; i< Help.size1(); i++)
        {
            for(unsigned int j=0; j< Help.size1(); j++)
            {
                vectors(i,j)= V(j,i);
            }
        }
 
        for(unsigned int i=0; i<Help.size1(); i++)
            lambda(i)= Help(i,i);
 
        return;
    }
 
 
    inline void CreatePartition(unsigned int number_of_threads, const int number_of_rows, DenseVector<unsigned int>& partitions)
    {
        partitions.resize(number_of_threads + 1);
        int partition_size = number_of_rows / number_of_threads;
        partitions[0] = 0;
        partitions[number_of_threads] = number_of_rows;
        for (unsigned int i = 1; i < number_of_threads; i++)
            partitions[i] = partitions[i - 1] + partition_size;
    }
 
 
 
 
 
 
 
 
 
 
    CalculateSignedDistanceTo3DSkinProcess& operator=(CalculateSignedDistanceTo3DSkinProcess const& rOther);
 
    //CalculateSignedDistanceTo3DSkinProcess(CalculateSignedDistanceTo3DSkinProcess const& rOther);
 
 
 
}; // Class CalculateSignedDistanceTo3DSkinProcess
 
 
 
 
 
 
inline std::istream& operator >> (std::istream& rIStream,
                                  CalculateSignedDistanceTo3DSkinProcess& rThis);
 
inline std::ostream& operator << (std::ostream& rOStream,
                                  const CalculateSignedDistanceTo3DSkinProcess& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
const double CalculateSignedDistanceTo3DSkinProcess::epsilon = 1e-18;
 
 
}  // namespace Kratos.
 
#endif // KRATOS_CALCULATE_DISTANCE_PROCESS_H_INCLUDED  defined