trigen_pfem_refine_vms.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
 
#if !defined(KRATOS_TRIGEN_PFEM_MODELER_VMS_H_INCLUDED )
#define  KRATOS_TRIGEN_PFEM_MODELER_VMS_H_INCLUDED
 
// System includes
 
// External includes
#include "triangle.h"
 
// Project includes
#include "includes/define.h"
#include "includes/model_part.h"
#include "includes/global_pointer_variables.h"
#include "utilities/timer.h"
#include "geometries/triangle_2d_3.h"
#include "meshing_application_variables.h"
#include "processes/entity_erase_process.h"
#include "spatial_containers/spatial_containers.h"
#include "trigen_pfem_refine.h"
 
namespace Kratos
{
extern "C" {
    void triangulate(char *, struct triangulateio *, struct triangulateio *,struct triangulateio *);
    //void trifree();
}
 
 
 
 
 
 
 
 
class TriGenPFEMModelerVMS : public TriGenPFEMModeler
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(TriGenPFEMModelerVMS);
    typedef Node PointType;
    typedef Node::Pointer PointPointerType;
    typedef std::vector<PointType::Pointer>           PointVector;
    typedef PointVector::iterator PointIterator;
    typedef std::vector<double>               DistanceVector;
    typedef std::vector<double>::iterator     DistanceIterator;
    typedef Bucket<3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType;
    typedef Tree< KDTreePartition<BucketType> > KdtreeType; //Kdtree;
 
 
    TriGenPFEMModelerVMS() :
        mJ(ZeroMatrix(2,2)), //local jacobian
        mJinv(ZeroMatrix(2,2)), //inverse jacobian
        mC(ZeroVector(2)), //dimension = number of nodes
        mRhs(ZeroVector(2)) {}
    //mpNodeEraseProcess(NULL){} //dimension = number of nodes
 
    ~TriGenPFEMModelerVMS() override = default;
 
 
 
 
 
 
    //*******************************************************************************************
    //*******************************************************************************************
    void ReGenerateMesh(
        ModelPart& ThisModelPart ,
        Element const& rReferenceElement,
        Condition const& rReferenceBoundaryCondition,
        EntitiesEraseProcess<Node>& node_erase, bool rem_nodes = true, bool add_nodes=true,
        double my_alpha = 1.4, double h_factor=0.5)
    {
 
        KRATOS_TRY
        if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FREE_SURFACE)==false )
            KRATOS_THROW_ERROR(std::logic_error,"Add  ----IS_FREE_SURFACE---- variable!!!!!! ERROR","");
        if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_STRUCTURE)==false )
            KRATOS_THROW_ERROR(std::logic_error,"Add  ----IS_STRUCTURE---- variable!!!!!! ERROR","");
        if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_BOUNDARY)==false )
            KRATOS_THROW_ERROR(std::logic_error,"Add  ----IS_BOUNDARY---- variable!!!!!! ERROR","");
        if (ThisModelPart.NodesBegin()->SolutionStepsDataHas(IS_FLUID)==false )
            KRATOS_THROW_ERROR(std::logic_error,"Add  ----IS_FLUID---- variable!!!!!! ERROR","");
 
        KRATOS_WATCH("Trigen PFEM Refining Mesher")
        const auto inital_time = std::chrono::steady_clock::now();
 
 
//clearing elements
 
        ThisModelPart.Elements().clear();
        ThisModelPart.Conditions().clear();
 
 
 
        int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
 
        // bucket types
        //typedef Bucket<3, PointType, ModelPart::NodesContainerType, PointPointerType, PointIterator, DistanceIterator > BucketType;
        //typedef Bins< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > StaticBins;
        // bucket types
 
        //typedef Tree< StaticBins > Bin;                //Binstree;
        unsigned int bucket_size = 20;
 
        //performing the interpolation - all of the nodes in this list will be preserved
        unsigned int max_results = 100;
        //PointerVector<PointType> res(max_results);
        //NodeIterator res(max_results);
        PointVector res(max_results);
        DistanceVector res_distances(max_results);
        Node work_point(0,0.0,0.0,0.0);
        //if the remove_node switch is activated, we check if the nodes got too close
 
        if (rem_nodes==true)
        {
            PointVector list_of_nodes;
            list_of_nodes.reserve(ThisModelPart.Nodes().size());
            for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
            {
                (list_of_nodes).push_back(*(i_node.base()));
            }
 
            KdtreeType nodes_tree1(list_of_nodes.begin(),list_of_nodes.end(), bucket_size);
 
            RemoveCloseNodes(ThisModelPart, nodes_tree1, node_erase, h_factor);
        }
 
        //                                                  //
        //      Now we shall pass the Alpha Shape for the second time, having the "bad nodes" removed   //
        //creating the containers for the input and output
        struct triangulateio in_mid;
        struct triangulateio out_mid;
        struct triangulateio vorout_mid;
 
        initialize_triangulateio(in_mid);
        initialize_triangulateio(out_mid);
        initialize_triangulateio(vorout_mid);
 
        //assigning the size of the input containers
 
        in_mid.numberofpoints = ThisModelPart.Nodes().size();
        in_mid.pointlist = (REAL *) malloc(in_mid.numberofpoints * 2 * sizeof(REAL));
 
        //reorder the id's and give the coordinates to the mesher
        ModelPart::NodesContainerType::iterator nodes_begin = ThisModelPart.NodesBegin();
        for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
        {
            int base = i*2;
            //int base = ((nodes_begin + i)->Id()   -  1 ) * 2;
 
            //from now on it is consecutive
            (nodes_begin + i)->SetId(i+1);
//              (nodes_begin + i)->Id() = i+1;
 
            in_mid.pointlist[base] = (nodes_begin + i)->X();
            in_mid.pointlist[base+1] = (nodes_begin + i)->Y();
 
            auto& node_dofs = (nodes_begin + i)->GetDofs();
            for(auto iii = node_dofs.begin();    iii != node_dofs.end(); iii++)
            {
                (**iii).SetEquationId(i+1);
            }
        }
        //in_mid.numberoftriangles = ThisModelPart.Elements().size();
        //in_mid.trianglelist = (int*) malloc(in_mid.numberoftriangles * 3 * sizeof(int));
 
        // "P" suppresses the output .poly file. Saves disk space, but you
        //lose the ability to maintain constraining segments  on later refinements of the mesh.
        // "B" Suppresses boundary markers in the output .node, .poly, and .edge output files
        // "n" outputs a list of triangles neighboring each triangle.
        // "e" outputs edge list (i.e. all the "connectivities")
        char options1[] = "Pne";
        triangulate(options1, &in_mid, &out_mid, &vorout_mid);
        //print out the mesh generation time
        std::cout<<"mesh generation time = " << Timer::ElapsedSeconds(inital_time) << std::endl;
        //number of newly generated triangles
        unsigned int el_number=out_mid.numberoftriangles;
 
        //PASSING THE CREATED ELEMENTS TO THE ALPHA-SHAPE
        std::vector<int> preserved_list1(el_number);
        preserved_list1.resize(el_number);
 
        int number_of_preserved_elems= PassAlphaShape(ThisModelPart, preserved_list1, el_number, my_alpha, out_mid);
 
        //freeing memory
 
        clean_triangulateio(in_mid);
        clean_triangulateio(vorout_mid);
        KRATOS_WATCH("ln367");
        //NOW WE SHALL PERFORM ADAPTIVE REMESHING, i.e. insert and remove nodes based upon mesh quality
        // and prescribed sizes
        struct triangulateio in2;
        struct triangulateio out2;
        struct triangulateio vorout2;
 
        initialize_triangulateio(in2);
        initialize_triangulateio(out2);
        initialize_triangulateio(vorout2);
 
//          in2.firstnumber = 1;
        in2.numberofpoints = ThisModelPart.Nodes().size();
        in2.pointlist = (REAL *) malloc(in2.numberofpoints * 2 * sizeof(REAL));
 
        //writing the input point list
        for(unsigned int i = 0; i<ThisModelPart.Nodes().size(); i++)
        {
            int base = i*2;
            in2.pointlist[base] = (nodes_begin + i)->X();
            in2.pointlist[base+1] = (nodes_begin + i)->Y();
        }
        in2.numberoftriangles=number_of_preserved_elems;
 
        in2.trianglelist = (int *) malloc(in2.numberoftriangles * 3 * sizeof(int));
        in2.trianglearealist = (REAL *) malloc(in2.numberoftriangles * sizeof(REAL));
//          in2.trianglelist = new int[in2.numberoftriangles * 3];
//          in2.trianglearealist = new double[in2.numberoftriangles];
 
        KRATOS_WATCH(el_number);
        int counter = 0;
        //here I will assign a huge number of NODAL_H to the free surface nodes, so that there no nodes will be added
        for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
        {
            if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
            {
 
                double& val=i_node->FastGetSolutionStepValue(NODAL_H);
                val*=2.0;
                //i_node->FastGetSolutionStepValue(NODAL_H,1)=val;
                //KRATOS_WATCH("AAAAAAAAAAAAAAAAAAAAAAAA!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
            }
        }
 
 
        for (unsigned int el=0; el<el_number; el++)
        {
            if( static_cast<bool>(preserved_list1[el]) == true )
            {
                //saving the list of ONLY preserved triangles, the ones that passed alpha-shape check
                int new_base = counter*3;
                int old_base = el*3;
                //copying in case it is preserved
                in2.trianglelist[new_base] = out_mid.trianglelist[old_base];
                in2.trianglelist[new_base+1] = out_mid.trianglelist[old_base+1];
                in2.trianglelist[new_base+2] = out_mid.trianglelist[old_base+2];
 
                //calculate the prescribed h
                double prescribed_h = (nodes_begin + out_mid.trianglelist[old_base]-1)->FastGetSolutionStepValue(NODAL_H);
                prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+1]-1)->FastGetSolutionStepValue(NODAL_H);
                prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+2]-1)->FastGetSolutionStepValue(NODAL_H);
                prescribed_h *= 0.3333;
                //if h is the height of a equilateral triangle, the area is 1/2*h*h
                in2.trianglearealist[counter] = 0.5*(1.5*prescribed_h*1.5*prescribed_h);
                counter += 1;
            }
 
        }
        //now I set back the nodal_h
        for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; i_node++)
        {
            if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
            {
                double& nodal_h=i_node->FastGetSolutionStepValue(NODAL_H);
                nodal_h/=2.0;
            }
        }
 
 
        clean_triangulateio(out_mid);
        KRATOS_WATCH("ln420");
        //here we generate a new mesh adding/removing nodes, by initializing "q"-quality mesh and "a"-area constraint switches
        //
        // MOST IMPORTANT IS "r" switch, that refines previously generated mesh!!!!!!!!!!(that is the one given inside in2)
        //char mesh_regen_opts[] = "YYJaqrn";
        //char mesh_regen_opts[] = "YJq1.4arn";
        if (add_nodes==true)
        {
            char mesh_regen_opts[] = "YJq1.4arn";
            triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
            KRATOS_WATCH("Adaptive remeshing executed")
        }
        else
        {
            char mesh_regen_opts[] = "YJrn";
            triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
            KRATOS_WATCH("Non-Adaptive remeshing executed")
        }
 
        //and now we shall find out where the new nodes belong to
        //defintions for spatial search
        typedef Node PointType;
        typedef Node::Pointer PointPointerType;
        typedef std::vector<PointType::Pointer>           PointVector;
        //typedef std::vector<PointType::Pointer>::iterator PointIterator;
        typedef PointVector::iterator PointIterator;
        typedef std::vector<double>               DistanceVector;
        typedef std::vector<double>::iterator     DistanceIterator;
        KRATOS_WATCH("ln449");
 
        typedef Bucket<3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType;
 
        typedef Tree< KDTreePartition<BucketType> > kd_tree; //Kdtree;
 
        //int step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
 
        //creating an auxiliary list for the new nodes
        PointVector list_of_new_nodes;
 
        //node to get the DOFs from
        Node::DofsContainerType& reference_dofs = (ThisModelPart.NodesBegin())->GetDofs();
 
        double z = 0.0;
        int n_points_before_refinement = in2.numberofpoints;
        //if points were added, we add them as nodes to the ModelPart
        if (out2.numberofpoints > n_points_before_refinement )
        {
            for(int i = n_points_before_refinement; i<out2.numberofpoints; i++)
            {
                int id=i+1;
                int base = i*2;
                double& x= out2.pointlist[base];
                double& y= out2.pointlist[base+1];
 
                Node::Pointer pnode = ThisModelPart.CreateNewNode(id,x,y,z);
 
                pnode->SetBufferSize(ThisModelPart.NodesBegin()->GetBufferSize() );
 
                list_of_new_nodes.push_back( pnode );
 
                //std::cout << "new node id = " << pnode->Id() << std::endl;
                //generating the dofs
                for(Node::DofsContainerType::iterator iii = reference_dofs.begin();    iii != reference_dofs.end(); iii++)
                {
                    Node::DofType& rDof = **iii;
                    Node::DofType::Pointer p_new_dof = pnode->pAddDof( rDof );
 
                    (p_new_dof)->FreeDof();
//                                                (p_new_dof)->EquationId() = -1;
 
                }
 
            }
        }
 
        std::cout << "During refinement we added " << out2.numberofpoints-n_points_before_refinement<< " nodes " <<std::endl;
        //unsigned int bucket_size = 20;
        //performing the interpolation - all of the nodes in this list will be preserved
        //unsigned int max_results = 50;
        //PointVector results(max_results);
        //DistanceVector results_distances(max_results);
        array_1d<double,3> N;
 
        array_1d<double,3> x1,x2,x3,xc;
 
        //int number_of_preserved_elems=0;
 
        int point_base;
        //WHAT ARE THOSE????
//          Node work_point(0,0.0,0.0,0.0);
        unsigned int MaximumNumberOfResults = list_of_new_nodes.size();
        PointVector Results(MaximumNumberOfResults);
        DistanceVector ResultsDistances(MaximumNumberOfResults);
 
        step_data_size = ThisModelPart.GetNodalSolutionStepDataSize();
 
        if(out2.numberofpoints-n_points_before_refinement > 0) //if we added points
        {
 
            kd_tree  nodes_tree2(list_of_new_nodes.begin(),list_of_new_nodes.end(),bucket_size);
            nodes_begin = ThisModelPart.NodesBegin();
 
            for(int el = 0; el< in2.numberoftriangles; el++)
            {
                int base = el * 3;
                //coordinates
                point_base = (in2.trianglelist[base] - 1)*2;
                x1[0] = in2.pointlist[point_base];
                x1[1] = in2.pointlist[point_base+1];
 
                point_base = (in2.trianglelist[base+1] - 1)*2;
                x2[0] = in2.pointlist[point_base];
                x2[1] = in2.pointlist[point_base+1];
 
                point_base = (in2.trianglelist[base+2] - 1)*2;
                x3[0] = in2.pointlist[point_base];
                x3[1] = in2.pointlist[point_base+1];
 
 
                //find the center and "radius" of the element
                double xc,  yc, radius;
                CalculateCenterAndSearchRadius( x1[0], x1[1],
                                                x2[0], x2[1],
                                                x3[0], x3[1],
                                                xc,yc,radius);
 
                //find all of the new nodes within the radius
                int number_of_points_in_radius;
                work_point.X() = xc;
                work_point.Y() = yc;
                work_point.Z() = 0.0;
 
                number_of_points_in_radius = nodes_tree2.SearchInRadius(work_point, radius*1.01, Results.begin(),
                                             ResultsDistances.begin(),  MaximumNumberOfResults);
 
                Triangle2D3<Node > geom(
                    *( (nodes_begin +  in2.trianglelist[base]-1).base()     ),
                    *( (nodes_begin +  in2.trianglelist[base+1]-1).base()   ),
                    *( (nodes_begin +  in2.trianglelist[base+2]-1).base()   )
                );
 
                //check if inside and eventually interpolate
                for( auto it_found = Results.begin(); it_found != Results.begin() + number_of_points_in_radius; it_found++)
                {
                    bool is_inside = false;
                    is_inside = CalculatePosition(x1[0], x1[1],
                                                  x2[0], x2[1],
                                                  x3[0], x3[1],
                                                  (*it_found)->X(),(*it_found)->Y(),N);
 
 
                    if(is_inside == true)
                    {
                        Interpolate(  geom,  N, step_data_size, *(it_found ) );
 
                    }
                }
            }
        }
 
        ThisModelPart.Elements().clear();
 
        //set the coordinates to the original value
        for(auto & list_of_new_node : list_of_new_nodes)
        {
            const array_1d<double,3>& disp = list_of_new_node->FastGetSolutionStepValue(DISPLACEMENT);
            list_of_new_node->X0() = list_of_new_node->X() - disp[0];
            list_of_new_node->Y0() = list_of_new_node->Y() - disp[1];
            list_of_new_node->Z0() = 0.0;
        }
        //cleaning unnecessary data
        //in2.deinitialize();
        //in2.initialize();
        //free( in2.pointlist );
 
        //add preserved elements to the kratos
        Properties::Pointer properties = ThisModelPart.GetMesh().pGetProperties(1);
        nodes_begin = ThisModelPart.NodesBegin();
        (ThisModelPart.Elements()).reserve(out2.numberoftriangles);
 
        for(int iii = 0; iii< out2.numberoftriangles; iii++)
        {
            int id = iii + 1;
            int base = iii * 3;
            Triangle2D3<Node > geom(
                *( (nodes_begin +  out2.trianglelist[base]-1).base()    ),
                *( (nodes_begin +  out2.trianglelist[base+1]-1).base()  ),
                *( (nodes_begin +  out2.trianglelist[base+2]-1).base()  )
            );
 
 
#ifdef _DEBUG
            ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
            if( *(ModelNodes).find( out2.trianglelist[base]).base() == *(ThisModelPart.Nodes().end()).base() )
                KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
            if( *(ModelNodes).find( out2.trianglelist[base+1]).base() == *(ThisModelPart.Nodes().end()).base() )
                KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
            if( *(ModelNodes).find( out2.trianglelist[base+2]).base() == *(ThisModelPart.Nodes().end()).base() )
                KRATOS_THROW_ERROR(std::logic_error,"trying to use an inexisting node","");
 
#endif
 
            Element::Pointer p_element = rReferenceElement.Create(id, geom, properties);
            (ThisModelPart.Elements()).push_back(p_element);
 
        }
        ThisModelPart.Elements().Sort();
 
        //filling the neighbour list
        ModelPart::ElementsContainerType::const_iterator el_begin = ThisModelPart.ElementsBegin();
        for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.ElementsBegin();
                iii != ThisModelPart.ElementsEnd(); iii++)
        {
            //Geometry< Node >& geom = iii->GetGeometry();
            int base = ( iii->Id() - 1 )*3;
 
            (iii->GetValue(NEIGHBOUR_ELEMENTS)).resize(3);
            GlobalPointersVector< Element >& neighb = iii->GetValue(NEIGHBOUR_ELEMENTS);
            for(int i = 0; i<3; i++)
            {
                int index = out2.neighborlist[base+i];
                if(index > 0)
                    neighb(i) = GlobalPointer<Element>(&*(el_begin + index-1));
                else
                    neighb(i) = Element::WeakPointer();
            }
        }
        //identifying boundary, creating skin
        IdentifyBoundary(ThisModelPart, rReferenceBoundaryCondition, properties, out2);
 
        KRATOS_WATCH("ln749");
 
        clean_triangulateio(in2);
        KRATOS_WATCH("ln752");
        clean_triangulateio(out2);
        KRATOS_WATCH("ln754");
        clean_triangulateio(vorout2);
        KRATOS_WATCH("Finished remeshing with Trigen_PFEM_Refine")
 
        KRATOS_CATCH("")
    }
 
 
 
 
 
 
 
    std::string Info() const override
    {
        return "";
    }
 
    void PrintInfo(std::ostream& rOStream) const override {}
 
    void PrintData(std::ostream& rOStream) const override {}
 
 
 
 
 
protected:
 
 
    //KdtreeType mKdtree;
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
    boost::numeric::ublas::bounded_matrix<double,2,2> mJ; //local jacobian
    boost::numeric::ublas::bounded_matrix<double,2,2> mJinv; //inverse jacobian
    array_1d<double,2> mC; //center pos
    array_1d<double,2> mRhs; //center pos
    //EntitiesEraseProcess<Node>* mpNodeEraseProcess;
 
 
    void RemoveCloseNodes(ModelPart& ThisModelPart, KdtreeType& nodes_tree1, EntitiesEraseProcess<Node>& node_erase, double& h_factor)
    {
        //unsigned int bucket_size = 20;
 
        //performing the interpolation - all of the nodes in this list will be preserved
        unsigned int max_results = 100;
        //PointerVector<PointType> res(max_results);
        //NodeIterator res(max_results);
        PointVector res(max_results);
        DistanceVector res_distances(max_results);
        Node work_point(0,0.0,0.0,0.0);
 
        unsigned int n_points_in_radius;
        //radius means the distance, closer than which no node shall be allowd. if closer -> mark for erasing
        double radius;
 
        for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in != ThisModelPart.NodesEnd(); in++)
        {
            radius=h_factor*in->FastGetSolutionStepValue(NODAL_H);
 
            work_point[0]=in->X();
            work_point[1]=in->Y();
            work_point[2]=in->Z();
 
            n_points_in_radius = nodes_tree1.SearchInRadius(work_point, radius, res.begin(),res_distances.begin(), max_results);
            if (n_points_in_radius>1)
            {
                if (in->FastGetSolutionStepValue(IS_BOUNDARY)==0.0 && in->FastGetSolutionStepValue(IS_STRUCTURE)==0.0)
                {
                    //look if we are already erasing any of the other nodes
                    double erased_nodes = 0;
                    for(auto i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
                        erased_nodes += in->Is(TO_ERASE);
 
                    if( erased_nodes < 1) //we cancel the node if no other nodes are being erased
                        in->Set(TO_ERASE,true);
 
                }
                else if ( (in)->FastGetSolutionStepValue(IS_STRUCTURE)!=1.0) //boundary nodes will be removed if they get REALLY close to another boundary node (0.2 * h_factor)
                {
                    //here we loop over the neighbouring nodes and if there are nodes
                    //with IS_BOUNDARY=1 which are closer than 0.2*nodal_h from our we remove the node we are considering
                    unsigned int k = 0;
                    unsigned int counter = 0;
                    for(auto i=res.begin(); i!=res.begin() + n_points_in_radius ; i++)
                    {
                        if ( (*i)->FastGetSolutionStepValue(IS_BOUNDARY,1)==1.0 && res_distances[k] < 0.2*radius && res_distances[k] > 0.0 )
                        {
                            //                                      KRATOS_WATCH( res_distances[k] );
                            counter += 1;
                        }
                        k++;
                    }
                    if(counter > 0)
                        in->Set(TO_ERASE,true);
                }
            }
 
        }
 
 
        node_erase.Execute();
 
        KRATOS_WATCH("Number of nodes after erasing")
        KRATOS_WATCH(ThisModelPart.Nodes().size())
    }
 
 
    int PassAlphaShape(ModelPart& ThisModelPart, std::vector<int>& preserved_list1, unsigned int & el_number, double& my_alpha, struct triangulateio& out_mid)
    {
        //NOTE THAT preserved_list1 will be overwritten, only the elements that passed alpha-shaoe check will enter it
 
        //prepairing for alpha shape passing : creating necessary arrays
        //list of preserved elements is created: at max el_number can be preserved (all elements)
 
 
        array_1d<double,3> x1,x2,x3,xc;
 
        //int number_of_preserved_elems=0;
        int number_of_preserved_elems=0;
        int point_base;
        //loop for passing alpha shape
        for(unsigned int el = 0; el< el_number; el++)
        {
            int base = el * 3;
 
            //coordinates
            point_base = (out_mid.trianglelist[base] - 1)*2;
            x1[0] = out_mid.pointlist[point_base];
            x1[1] = out_mid.pointlist[point_base+1];
 
            point_base = (out_mid.trianglelist[base+1] - 1)*2;
            x2[0] = out_mid.pointlist[point_base];
            x2[1] = out_mid.pointlist[point_base+1];
 
            point_base = (out_mid.trianglelist[base+2] - 1)*2;
            x3[0] = out_mid.pointlist[point_base];
            x3[1] = out_mid.pointlist[point_base+1];
 
            //here we shall temporarily save the elements and pass them afterwards to the alpha shape
            Geometry<Node > temp;
 
            temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base]).base()  ) );
            temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+1]).base()    ) );
            temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+2]).base()    ) );
 
            int number_of_structure_nodes = int( temp[0].FastGetSolutionStepValue(IS_STRUCTURE) );
            number_of_structure_nodes += int( temp[1].FastGetSolutionStepValue(IS_STRUCTURE) );
            number_of_structure_nodes += int( temp[2].FastGetSolutionStepValue(IS_STRUCTURE) );
 
            //check the number of nodes of boundary
            int nfs = int( temp[0].FastGetSolutionStepValue(IS_FREE_SURFACE) );
            nfs += int( temp[1].FastGetSolutionStepValue(IS_FREE_SURFACE) );
            nfs += int( temp[2].FastGetSolutionStepValue(IS_FREE_SURFACE) );
 
            //check the number of nodes of boundary
            int nfluid = int( temp[0].FastGetSolutionStepValue(IS_FLUID) );
            nfluid += int( temp[1].FastGetSolutionStepValue(IS_FLUID) );
            nfluid += int( temp[2].FastGetSolutionStepValue(IS_FLUID) );
 
            //check the number of nodes of boundary
            // int nboundary = int( temp[0].FastGetSolutionStepValue(IS_BOUNDARY) );
            // nboundary += int( temp[1].FastGetSolutionStepValue(IS_BOUNDARY) );
            // nboundary += int( temp[2].FastGetSolutionStepValue(IS_BOUNDARY) );
            //first check that we are working with fluid elements, otherwise throw an error
            //if (nfluid!=3)
            //  KRATOS_THROW_ERROR(std::logic_error,"THATS NOT FLUID or NOT TRIANGLE!!!!!! ERROR","");
            //otherwise perform alpha shape check
 
 
            if(number_of_structure_nodes!=3) //if it is = 3 it is a completely fixed element -> do not add it
            {
                preserved_list1[el] = false;
                if (nfs != 0 || nfluid != 3 )  //in this case it is close to the surface so i should use alpha shape
                {
 
                    if( AlphaShape(my_alpha,temp) && number_of_structure_nodes!=3) //if alpha shape says to preserve
                    {
//                          if(nboundary==3 && number_of_structure_nodes > 1 && nfs > 0) //if it is = 3 pressure problems -> do not add it
//                          {
//                                preserved_list1[el] = false;
//                          }
//                          else
//                          {
                        preserved_list1[el] = true;
                        number_of_preserved_elems += 1;
//                          }
                    }
                }
                else //internal triangle --- should be ALWAYS preserved
                {
                    double bigger_alpha = my_alpha*10.0;
                    if( AlphaShape(bigger_alpha,temp) && number_of_structure_nodes!=3)
                    {
//                          if(nboundary==3 && number_of_structure_nodes > 1 && nfs > 0) //if it is = 3 pressure problems -> do not add it
//                          {
//                                preserved_list1[el] = false;
//                          }
//                          else
//                          {
                        preserved_list1[el] = true;
                        number_of_preserved_elems += 1;
//                          }
                    }
                }
//qua inizia la parte cambiata
 
                if ( preserved_list1[el] == false && number_of_structure_nodes == 1 )
                {
                    int noFSEdge = 0;
                    int Vertexiter[4] = {0,1,2,0};
                    for (int iVertex = 0 ; iVertex < 3 ; iVertex++ )
                    {
                        if (    ( temp[Vertexiter[iVertex]].FastGetSolutionStepValue(IS_FREE_SURFACE) && temp[Vertexiter[iVertex + 1]].FastGetSolutionStepValue(IS_FREE_SURFACE) )
                             || ( temp[Vertexiter[iVertex]].FastGetSolutionStepValue(IS_FREE_SURFACE) && temp[Vertexiter[iVertex + 1]].FastGetSolutionStepValue(IS_STRUCTURE) )
                             || ( temp[Vertexiter[iVertex]].FastGetSolutionStepValue(IS_STRUCTURE) && temp[Vertexiter[iVertex + 1]].FastGetSolutionStepValue(IS_FREE_SURFACE) )
                           )
                        {
                            noFSEdge++;
                        }
                    }
                    if (noFSEdge == 3)
                    {
                        double bigger_alpha = my_alpha*1.23;
                        if( AlphaShape(bigger_alpha,temp) && number_of_structure_nodes!=3)
                        {
    //                          if(nboundary==3 && number_of_structure_nodes > 1 && nfs > 0) //if it is = 3 pressure problems -> do not add it
    //                          {
    //                                preserved_list1[el] = false;
    //                          }
    //                          else
    //                          {
                            preserved_list1[el] = true;
                            number_of_preserved_elems += 1;
    //                          }
                        }
                    }
             }
 
// fine della parte cambiata
 
              }
            else
                preserved_list1[el] = false;
 
        }
        return number_of_preserved_elems;
    }
 
 
 
 
 
    void IdentifyBoundary(ModelPart& ThisModelPart, Condition const& rReferenceBoundaryCondition,  Properties::Pointer& properties, struct triangulateio& out2)
    {
 
        //reset the boundary flag
        for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in!=ThisModelPart.NodesEnd(); in++)
        {
            in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
        }
        //filling the elemental neighbours list (from now on the elements list can not change)
        ModelPart::ElementsContainerType::iterator elements_end = ThisModelPart.Elements().end();
 
        ThisModelPart.Elements().Unique();
 
        //now the boundary faces
        for(ModelPart::ElementsContainerType::iterator iii = ThisModelPart.ElementsBegin(); iii != ThisModelPart.ElementsEnd(); iii++)
        {
            int base = ( iii->Id() - 1 )*3;
 
            ModelPart::ElementsContainerType::iterator el_neighb;
            /*each face is opposite to the corresponding node number so
             0 ----- 1 2
             1 ----- 2 0
             2 ----- 0 1
            */
 
            //
            //********************************************************************
            //first face
            el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base] );
            if( el_neighb == elements_end )
            {
                //std::cout << "node0" << std::endl;
                //if no neighnour is present => the face is free surface
                iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
                iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
 
                //Generate condition
                Condition::NodesArrayType temp1;
                temp1.reserve(2);
                temp1.push_back(iii->GetGeometry()(1));
                temp1.push_back(iii->GetGeometry()(2));
 
                Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
                int id = (iii->Id()-1)*3;
                Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
                ThisModelPart.Conditions().push_back(p_cond);
 
            }
 
            //********************************************************************
            //second face
            el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+1] );
            //if( el != ThisModelPart.Elements().end() )
            if( el_neighb == elements_end )
            {
                //if no neighnour is present => the face is free surface
                iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
                iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
 
                //Generate condition
                Condition::NodesArrayType temp1;
                temp1.reserve(2);
                temp1.push_back(iii->GetGeometry()(2));
                temp1.push_back(iii->GetGeometry()(0));
 
                Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
                int id = (iii->Id()-1)*3+1;
                Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
                ThisModelPart.Conditions().push_back(p_cond);
 
 
            }
 
            //********************************************************************
            //third face
            el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+2] );
            if( el_neighb == elements_end )
            {
                //if no neighnour is present => the face is free surface
                iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
                iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
 
//                  Generate condition
                Condition::NodesArrayType temp1;
                temp1.reserve(2);
                temp1.push_back(iii->GetGeometry()(0));
                temp1.push_back(iii->GetGeometry()(1));
 
                Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(new Geometry< Node >(temp1) );
                int id = (iii->Id()-1)*3+2;
                Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(id, temp1, properties);
                ThisModelPart.Conditions().push_back(p_cond);
 
            }
 
 
        }
    }
 
 
    //returns false if it should be removed
    bool AlphaShape(double alpha_param, Geometry<Node >& pgeom)
    //bool AlphaShape(double alpha_param, Triangle2D<Node >& pgeom)
    {
        KRATOS_TRY
 
 
        double x0 = pgeom[0].X();
        double x1 = pgeom[1].X();
        double x2 = pgeom[2].X();
 
        double y0 = pgeom[0].Y();
        double y1 = pgeom[1].Y();
        double y2 = pgeom[2].Y();
 
        mJ(0,0)=2.0*(x1-x0);
        mJ(0,1)=2.0*(y1-y0);
        mJ(1,0)=2.0*(x2-x0);
        mJ(1,1)=2.0*(y2-y0);
 
 
        double detJ = mJ(0,0)*mJ(1,1)-mJ(0,1)*mJ(1,0);
 
        mJinv(0,0) =  mJ(1,1);
        mJinv(0,1) = -mJ(0,1);
        mJinv(1,0) = -mJ(1,0);
        mJinv(1,1) =  mJ(0,0);
 
        bounded_matrix<double,2,2> check;
 
//calculate average h
        double h;
        h =  pgeom[0].FastGetSolutionStepValue(NODAL_H);
        h += pgeom[1].FastGetSolutionStepValue(NODAL_H);
        h += pgeom[2].FastGetSolutionStepValue(NODAL_H);
        h *= 0.333333333;
 
 
        if(detJ < 5e-3*h*h)
        {
            //std::cout << "detJ = " << detJ << std::endl;
            pgeom[0].GetSolutionStepValue(IS_BOUNDARY) = 1;
            pgeom[1].GetSolutionStepValue(IS_BOUNDARY) = 1;
            pgeom[2].GetSolutionStepValue(IS_BOUNDARY) = 1;
            return false;
        }
 
        else
        {
 
            double x0_2 = x0*x0 + y0*y0;
            double x1_2 = x1*x1 + y1*y1;
            double x2_2 = x2*x2 + y2*y2;
 
            //finalizing the calculation of the inverted matrix
            //std::cout<<"MATR INV"<<MatrixInverse(mJ)<<std::endl;
            mJinv /= detJ;
            //calculating the RHS
            mRhs[0] = (x1_2 - x0_2);
            mRhs[1] = (x2_2 - x0_2);
 
            //calculate position of the center
            noalias(mC) = prod(mJinv,mRhs);
 
            double radius = sqrt(pow(mC[0]-x0,2)+pow(mC[1]-y0,2));
 
 
            if (radius < h*alpha_param)
            {
                return true;
            }
            else
            {
                return false;
            }
        }
 
 
        KRATOS_CATCH("")
    }
    //AUXILLIARY FCTNS
    inline void CalculateCenterAndSearchRadius(const double x0, const double y0,
            const double x1, const double y1,
            const double x2, const double y2,
            double& xc, double& yc, double& R
                                              )
    {
        xc = 0.3333333333333333333*(x0+x1+x2);
        yc = 0.3333333333333333333*(y0+y1+y2);
 
        double R1 = (xc-x0)*(xc-x0) + (yc-y0)*(yc-y0);
        double R2 = (xc-x1)*(xc-x1) + (yc-y1)*(yc-y1);
        double R3 = (xc-x2)*(xc-x2) + (yc-y2)*(yc-y2);
 
        R = R1;
        if(R2 > R) R = R2;
        if(R3 > R) R = R3;
 
        R = sqrt(R);
    }
 
    inline double CalculateVol( const double x0, const double y0,
                                const double x1, const double y1,
                                const double x2, const double y2
                              )
    {
        return 0.5*( (x1-x0)*(y2-y0)- (y1-y0)*(x2-x0) );
    }
 
    inline bool CalculatePosition(  const double x0, const double y0,
                                    const double x1, const double y1,
                                    const double x2, const double y2,
                                    const double xc, const double yc,
                                    array_1d<double,3>& N
                                 )
    {
        double area = CalculateVol(x0,y0,x1,y1,x2,y2);
 
        if(area < 0.000000000001)
        {
            KRATOS_THROW_ERROR(std::logic_error,"element with zero area found","");
        }
 
 
 
        N[0] = CalculateVol(x1,y1,x2,y2,xc,yc)  / area;
        N[1] = CalculateVol(x2,y2,x0,y0,xc,yc)  / area;
        N[2] = CalculateVol(x0,y0,x1,y1,xc,yc)  / area;
 
        /*            N[0] = CalculateVol(x0,y0,x1,y1,xc,yc) * inv_area;
                    N[1] = CalculateVol(x1,y1,x2,y2,xc,yc) * inv_area;
                    N[2] = CalculateVol(x2,y2,x0,y0,xc,yc) * inv_area;*/
        double tol = 1e-4;
        double upper_limit = 1.0+tol;
        double lower_limit = -tol;
 
        if(N[0] >= lower_limit && N[1] >= lower_limit && N[2] >= lower_limit && N[0] <= upper_limit && N[1] <= upper_limit && N[2] <= upper_limit) //if the xc yc is inside the triangle
            //if(N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0) //if the xc yc is inside the triangle return true
            return true;
 
        return false;
    }
 
    void Interpolate( Triangle2D3<Node >& geom, const array_1d<double,3>& N,
                      unsigned int step_data_size,
                      Node::Pointer pnode)
    {
        unsigned int buffer_size = pnode->GetBufferSize();
        //KRATOS_WATCH("Buffer size")
        //KRATOS_WATCH(buffer_size)
 
        for(unsigned int step = 0; step<buffer_size; step++)
        {
 
            //getting the data of the solution step
            double* step_data = (pnode)->SolutionStepData().Data(step);
 
 
            double* node0_data = geom[0].SolutionStepData().Data(step);
            double* node1_data = geom[1].SolutionStepData().Data(step);
            double* node2_data = geom[2].SolutionStepData().Data(step);
 
            //copying this data in the position of the vector we are interested in
            for(unsigned int j= 0; j<step_data_size; j++)
            {
 
                step_data[j] = N[0]*node0_data[j] + N[1]*node1_data[j] + N[2]*node2_data[j];
 
 
            }
        }
        if (N[0]==0.0 && N[1]==0.0 && N[2]==0.0)
            KRATOS_THROW_ERROR(std::logic_error,"SOMETHING's wrong with the added nodes!!!!!! ERROR","");
 
        //if ( pnode->FastGetSolutionStepValue(BULK_MODULUS)==0.0)
        //      KRATOS_THROW_ERROR(std::logic_error,"SOMETHING's wrong with the added nodes!!!!!! ERROR","");
 
        //now we assure that the flag variables are set coorect!! since we add nodes inside of the fluid volume only
        //we manually reset the IS_BOUNDARY, IS_FLUID, IS_STRUCTURE, IS_FREE_SURFACE values in a right way
        //not to have values, like 0.33 0.66 resulting if we would have been interpolating them in the same way
        //as the normal variables, like Velocity etc
 
 
        pnode->FastGetSolutionStepValue(IS_BOUNDARY)=0.0;
        pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
        pnode->FastGetSolutionStepValue(IS_INTERFACE)=0.0;
        pnode->Set(TO_ERASE,false);
        pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=0.0;
        pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
    }
 
    void initialize_triangulateio( triangulateio& tr )
    {
        tr.pointlist                  = (REAL*) nullptr;
        tr.pointattributelist         = (REAL*) nullptr;
        tr.pointmarkerlist            = (int*) nullptr;
        tr.numberofpoints             = 0;
        tr.numberofpointattributes    = 0;
        tr.trianglelist               = (int*) nullptr;
        tr.triangleattributelist      = (REAL*) nullptr;
        tr.trianglearealist           = (REAL*) nullptr;
        tr.neighborlist               = (int*) nullptr;
        tr.numberoftriangles          = 0;
        tr.numberofcorners            = 3;
        tr.numberoftriangleattributes = 0;
        tr.segmentlist                = (int*) nullptr;
        tr.segmentmarkerlist          = (int*) nullptr;
        tr.numberofsegments           = 0;
        tr.holelist                   = (REAL*) nullptr;
        tr.numberofholes              = 0;
        tr.regionlist                 = (REAL*) nullptr;
        tr.numberofregions            = 0;
        tr.edgelist                   = (int*) nullptr;
        tr.edgemarkerlist             = (int*) nullptr;
        tr.normlist                   = (REAL*) nullptr;
        tr.numberofedges              = 0;
    };
 
    void clean_triangulateio( triangulateio& tr )
    {
        if(tr.pointlist != nullptr) free(tr.pointlist );
        if(tr.pointattributelist != nullptr) free(tr.pointattributelist );
        if(tr.pointmarkerlist != nullptr) free(tr.pointmarkerlist   );
        if(tr.trianglelist != nullptr) free(tr.trianglelist  );
        if(tr.triangleattributelist != nullptr) free(tr.triangleattributelist );
        if(tr.trianglearealist != nullptr) free(tr.trianglearealist );
        if(tr.neighborlist != nullptr) free(tr.neighborlist   );
        if(tr.segmentlist != nullptr) free(tr.segmentlist    );
        if(tr.segmentmarkerlist != nullptr) free(tr.segmentmarkerlist   );
        if(tr.holelist != nullptr) free(tr.holelist      );
        if(tr.regionlist != nullptr) free(tr.regionlist  );
        if(tr.edgelist != nullptr) free(tr.edgelist   );
        if(tr.edgemarkerlist != nullptr) free(tr.edgemarkerlist   );
        if(tr.normlist != nullptr) free(tr.normlist  );
    };
 
 
 
 
 
 
 
    TriGenPFEMModelerVMS & operator=(TriGenPFEMModelerVMS const& rOther);
 
 
 
}; // Class TriGenPFEMModelerVMS
 
 
 
 
 
 
inline std::istream& operator >> (std::istream& rIStream,
                                  TriGenPFEMModelerVMS& rThis);
 
inline std::ostream& operator << (std::ostream& rOStream,
                                  const TriGenPFEMModelerVMS& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_TRIGEN_PFEM_MODELER_VMS_H_INCLUDED  defined