CAD4FE_PolyCurve.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       CAD4FE_PolyCurve.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:54 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "gestionversion.h"
#include "mg_geometrie.h"
#include "mg_gestionnaire.h"
#include "mg_arete.h"
#include "lc_point.h"
#include "tpl_fonction.h"
#include "ot_algorithme_geometrique.h"
#include "ot_decalage_parametre.h"
#include "CAD4FE_MG_ARETE_ClosestPointOn.h"
 
 
#include <map>
#include <ostream>
#include <string>
#include <vector>
#include <sstream>
 
 
#pragma hdrstop
 
#include "CAD4FE_PolyCurve.h"
 
 
#pragma package(smart_init)
 
#include <math.h>
 
#ifdef __BORLANDC__
#pragma warn -8012      
#pragma warn -8037
#endif
 
using namespace CAD4FE;
 
PolyCurve::PolyCurve()
{
 
}
 
PolyCurve::PolyCurve(MG_ARETE * __edge)
{
        t_min = (0);
        t_max = (-1);
        InsertCurve( __edge );
}
 
PolyCurve::PolyCurve(MG_SOMMET * __vertex)
{
    t_min = (0);
    t_max = (0);
    lst_vertices.push_back(__vertex);
    lst_length.push_back(0);
}       
 
double
PolyCurve::GetLength(MG_ARETE * __edge)
{
        return __edge->get_longueur(__edge->get_tmin(), __edge->get_tmax());
}
 
double
PolyCurve::GetLength(unsigned __index)
{
        return lst_length[__index];
}
             
double equation_longueur(MG_ARETE & __edge,double t)
{
        PolyCurve::VerifyRefEdgeT(&__edge, t);
        double dxyz[3];
        __edge.deriver(t,dxyz);
        return sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
}
 
 
MG_ARETE * PolyCurve::GetRefEdge(unsigned int __index)
{
    return lst_ref_edges[__index];
}
 
MG_SOMMET * PolyCurve::GetRefVertex(unsigned int __index)
{
    return lst_vertices[__index];
}
 
unsigned PolyCurve::GetRefEdgeCount() 
{
    return lst_ref_edges.size();
}
 
unsigned PolyCurve::GetRefVertexCount() 
{
    return lst_vertices.size();
}
 
 
bool PolyCurve::ContainsRefEdge(MG_ARETE * __refEdge)
{
    return ( std::find(lst_ref_edges.begin(),lst_ref_edges.end(),__refEdge) != lst_ref_edges.end());
}
 
bool PolyCurve::ContainsRefVertex(MG_SOMMET * __v)
{
    return ( std::find(lst_vertices.begin(),lst_vertices.end(),__v) != lst_vertices.end());
}
 
bool PolyCurve::Contains(MG_ELEMENT_TOPOLOGIQUE * __topo)
{
    if (__topo->get_dimension() == 0)
        return ContainsRefVertex((MG_SOMMET *)__topo);
    if (__topo->get_dimension() == 1)
        return ContainsRefEdge((MG_ARETE *)__topo);
    return false;
}       
 
void
PolyCurve::VerifyRefEdgeT(MG_ARETE * __edge, double & __t)
{
    double t1 = __edge->get_tmin();
    double t2 = __edge->get_tmax();
 
    if ( __edge->get_courbe()->est_periodique() )
    {
            if (__t > t2 )
                    __t -= __edge->get_courbe()->get_periode();
            if (__t < t1 )
                    __t += __edge->get_courbe()->get_periode();
    }
    else
    {
        double tStart=(t1<t2)?t1:t2;
        double tEnd=(t1>t2)?t1:t2;
            if (__t > tEnd )
                    __t = tEnd;
            if (__t < tStart )
                    __t = tStart;
    }
}
 
void
PolyCurve::inverser(double & __t, double __point[3], double precision, bool __curvilinearLength)
{
        int j;
        double s;
 
        double edge_length = t_max;
 
        // Test wether the requested point is coincident with a vertex of the polyline
        for (std::vector<MG_SOMMET *>::iterator it_vertex = lst_vertices.begin();
                it_vertex != lst_vertices.end();
                it_vertex++)
        {
                double vertex_point [3];
                MG_SOMMET * vertex = *it_vertex;
                vertex->get_point()->evaluer(vertex_point);
 
                for (j = 0; j < 3; j++)
                        if (fabs(vertex_point[j] - __point[j]) > edge_length*precision)
                                break;
 
                if ( j == 3 )
                {
                        unsigned iVertex = std::distance ( lst_vertices.begin(), it_vertex );
                        s = 0;
                        for (unsigned k = 0; k < iVertex; k++)
                               s += lst_length[k];
 
                        __t = s;
 
                        return;
                }
        }
 
        // If the requested point is not on a vertex
        // then test the reference edges
        unsigned iClosestEdge;
        double tClosestEdge = 0;
        double min_distance = 1E99;
        MG_ARETE * edge;
        for (std::vector <MG_ARETE *>::iterator it_edge = lst_ref_edges.begin();
                it_edge != lst_ref_edges.end();
                it_edge++)
        {
                edge = *it_edge;
                double t_edge;
                //edge->inverser ( t_edge, __point );
                MG_ARETE_ClosestPointOn func(__point, edge);
                t_edge = func.Find();
                double xEdge = (t_edge-edge->get_tmin())/(edge->get_tmax()-edge->get_tmin());
                if (xEdge>1) t_edge = edge->get_tmax();
                if (xEdge<0) t_edge = edge->get_tmin();
                double point[3];
                edge->evaluer (t_edge, point );
                double distance = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(point, __point);
                if (distance < min_distance)
                {
                        min_distance = distance;
                        iClosestEdge = std::distance (lst_ref_edges.begin(), it_edge );
                        tClosestEdge = t_edge;
                }
        }
 
        Parameter_RefEdgeTToS ( tClosestEdge, lst_ref_edges[iClosestEdge], &s, __curvilinearLength);
 
        __t = s;
 
        return;
}
 
bool PolyCurve::est_sur_courbe(double* xyz, double precision)
{
    std::cout <<" *** PolyCurve::est_sur_courbe : FONCTION NON IMPLEMENTE ***" << std::endl;
    return false;
}
 
 
void
PolyCurve::VerifyS(double & __s)
{
                if ( est_periodique() )
                {
                        while (__s > get_tmax() )
                                __s -= get_periode();
                }
                else
                {
                        if (__s >  get_tmax() )
                                __s = get_tmax();
                }
}
 
double 
PolyCurve::RefVertex_GetS(MG_SOMMET * __refVertex)
{
    int N = RefVertex_GetIndex(__refVertex);
    double s;
    s=0;
    for (int i=0;i<N;i++)
        s += GetLength(i);
 
    return s;
}
 
void
PolyCurve::RefEdge_GetT(unsigned __iEdge,
                double __s,
                double *__t, double *__dt, bool __curvilinearLength
                )
{
        double t1,t2;
        MG_SOMMET * vertex = lst_vertices[ __iEdge ];
        MG_ARETE * edge = lst_ref_edges[ __iEdge ];
 
        if (vertex == edge->get_cosommet1()->get_sommet())
        {
            t1=edge->get_tmin();
            t2=edge->get_tmax();
        }
        else
        {
            t1=edge->get_tmax();
            t2=edge->get_tmin();
        }
        
        if ( __curvilinearLength == false )
        {
                *__dt = (t2  - t1) / GetLength( __iEdge );
                *__t = t1  + __s * (*__dt);
        }
        else
        {
                TPL_FONCTION1<double,MG_ARETE,double> longueur_bsp(*edge, equation_longueur);
                double s = __s;
                double increment = s/t_max*(t2-t1)/40;
                if (increment < 1E-4 * (t2-t1))
                {
                    *__dt = (t2  - t1) / GetLength( __iEdge );
                    *__t = t1  + __s * (*__dt);
 
                    return;
                }
                longueur_bsp.integrer_jusqua_gauss_2(t1, increment, s, __t);
                *__dt = (*__t - t1) / s;
        }
}
 
void
PolyCurve::RefEdge_GetS(unsigned __iEdge,
                double __t,
                double *__s, bool __curvilinearLength
                )
{
        double t1,t2;
        MG_SOMMET * vertex = lst_vertices[ __iEdge ];
        MG_ARETE * edge = lst_ref_edges[ __iEdge ];
        double period=edge->get_courbe()->get_periode();
 
        t1=edge->get_tmin();
        t2=edge->get_tmax();
 
        if (edge->get_courbe()->est_periodique())
        {
            if (__t < t1 - period*1E-6)
                    __t += period;
        }
 
        if ( __curvilinearLength == false)
        {
            if (vertex == edge->get_cosommet1()->get_sommet())
                *__s = lst_length[__iEdge] * (__t-t1)/(t2-t1);
            else
                *__s = lst_length[__iEdge] * (__t-t2)/(t1-t2);
        }
        else
        {
            if (vertex == edge->get_cosommet1()->get_sommet())
                *__s = edge->get_longueur(t1,__t);
            else
                *__s = edge->get_longueur(__t,t2);
        }
}
 
void
PolyCurve::InsertCurve(MG_ARETE *__edge)
{                
        if (lst_ref_edges.size() == 0)
        {
                lst_ref_edges.push_back ( __edge );
                lst_length.push_back ( GetLength(__edge) );
                lst_vertices.push_back ( __edge->get_cosommet1()->get_sommet());
                lst_vertices.push_back ( __edge->get_cosommet2()->get_sommet());
                t_min = 0;
                t_max = lst_length[0];
                return;
        }
        
        // Get the two vertices of this edge
        MG_SOMMET * V[2];
        int iV[2], // index of V[i] in lst_vertices[]  ie. V [ i ] == lst_vertices[ iV[ i ] ]
        iiVc=-1, // index < 2 of the common vertex in V[0,1]  ie. if (iVc != -1) V [ iiVc ] == lst_vertices[ iV[ iiVc ] ]
        iiVo=-1; // index < 2 of the other  vertex in V[0,1]  ie. if (iVo != -1) V [ iiVo ] == lst_vertices[ iV[ iiVo ] ]
        V[0] = __edge->get_cosommet1()->get_sommet();
        V[1] = __edge->get_cosommet2()->get_sommet();
 
        // Verify wether the vertices are already referenced in
        // the polycurve
        std::vector<MG_SOMMET*>::iterator itV[2];
        for (int i=0; i<2; i++)
        {
                itV[i] = std::find (lst_vertices.begin(), lst_vertices.end(), V[i]);
                if (itV[i] != lst_vertices.end())
                {  // V[i] is already in a curve of the polycurve
                        iV[i] = std::distance(lst_vertices.begin(), itV[i]);
                        if (iiVc == -1)
                                iiVc = i;
                        else
                        {
                            if (iV[iiVc] == 0)
                            {  // case of a periodic polycurve : two vertices are common, the
                               // common vertex should be the last vertex of this polycurve (iV[iiVc]=N-1),
                               // and the other vertex the first vertex of this polycurve (iV[iiVo]=0)
                                iiVo = iiVc;
                                iiVc = i;
                            }
                            else
                            {
                                iiVo = i;
                            }
                        }
 
                        if ( iV[i] != 0 && iV[i]+1 != (int)lst_vertices.size() )
                        {
                                printf (" Error in PolyCurve::InsertCurve : can't insert a MG_ARETE (extremities = %lu, %lu) not adjacent to the extremities of this PolyCurve (extremities = %lu, %lu)\n", __edge->get_cosommet1()->get_sommet()->get_id(), __edge->get_cosommet2()->get_sommet()->get_id(), get_sommet1()->get_id(), get_sommet2()->get_id());
                                return;
                        }
                }
                else               // vertex V[i] not found
                {
                        iV[i] = -1;
                        iiVo = i;
                }
        }
 
 
        if (iiVc == -1)
        {
                printf (" Error in PolyCurve::InsertCurve : the curve requested to append is not connected to the existing polycurve\n");
                return;
        }
 
        if ( iV [iiVc] == 0 ) // the inserted curve is adjacent to the first vertex of the polycurve
        {                                          
                std::reverse(lst_ref_edges.begin(), lst_ref_edges.end());
                std::reverse(lst_vertices.begin(), lst_vertices.end());
                std::reverse(lst_length.begin(), lst_length.end());
                InsertCurve(__edge);
        }
        else if ( iV [iiVc]+1 == lst_vertices.size() )  // the inserted curve is adjacent to the last vertex of the polycurve
        {
                lst_vertices.push_back ( V [ iiVo ] );
                lst_ref_edges.push_back ( __edge );
                lst_length.push_back ( GetLength(__edge) );
        }
        else
        {
            printf("PolyCurve::InsertCurve Error : the edge %lu is neither adjacent to end vertex nor to start vertex of polycurve %lu\n", __edge->get_id(), get_id());
            printf("Polycurve Start vertex = %lu, End vertex = %lu\n", get_sommet1()->get_id(),get_sommet2()->get_id());
            printf("Edge Start vertex = %lu, End vertex = %lu\n", V[0]->get_id(), V[1]->get_id());
            return;               
        }
 
        // update t_max        
        t_max = 0;
        for (std::vector<double>::const_iterator it_length = lst_length.begin();
                it_length != lst_length.end();
                it_length++)
                t_max += (*it_length);
 
}
 
void
PolyCurve::Parameter_SToRefEdgeT (double __s, unsigned * __iEdge, double *__t, double *__dt, bool __curvilinearLength)
{
        double s;
        unsigned i;
 
        VerifyS(__s);
 
        i=0;
        s=0;
        while (s + GetLength(i) < __s && i+1 < lst_ref_edges.size() )
        {
                s += GetLength(i);
                i++;
        }
 
        RefEdge_GetT(i, __s - s, __t, __dt, __curvilinearLength);
        *__iEdge = i;
}
 
void
PolyCurve::Parameter_SToRefEdgeT (double __s, MG_ARETE ** __edge, double *__t, double *__dt, bool __curvilinearLength)
{
    unsigned iEdge;
    Parameter_SToRefEdgeT(__s, &iEdge, __t, __dt, __curvilinearLength);
    *__edge = GetRefEdge(iEdge);
}
 
unsigned 
PolyCurve::RefEdge_GetIndex(MG_ARETE * __edge)
{
    unsigned i = std::distance ( lst_ref_edges.begin(), std::find ( lst_ref_edges.begin(), lst_ref_edges.end(), __edge) );
    return i;
}
 
unsigned 
PolyCurve::RefVertex_GetIndex(MG_SOMMET * __vertex)
{
    unsigned i = std::distance ( lst_vertices.begin(), std::find ( lst_vertices.begin(), lst_vertices.end(), __vertex) );
    return i;
}
 
void
PolyCurve::Parameter_RefEdgeTToS (double __t, MG_ARETE * __edge, double *__s, bool __curvilinearLength)
{                
        double s, ds;
        unsigned i;
 
        i = RefEdge_GetIndex(__edge);
 
        s=0;
        for ( unsigned k = 0; k < i; k++)
        {
                s += GetLength(k);
        }
        RefEdge_GetS(i, __t, &ds, __curvilinearLength);
        s += ds;
 
        *__s = s;
}
                                         
void
PolyCurve::evaluer (double __s, double __X[3])
{
        bool curvilinearLength = false;
        evaluer(__s, __X, curvilinearLength);
}
 
void
PolyCurve::evaluer (double __s, double __X[3], bool __curvilinearLength) // curvilinear abcissa
{
        if (IsPoint())
            {
                lst_vertices[0]->get_point()->evaluer(__X);
                return;
            }
 
        VerifyS (__s);
        unsigned iEdge;
        MG_ARETE * edge;
        double t, dt;
        Parameter_SToRefEdgeT ( __s, & iEdge , &t, &dt, __curvilinearLength );
        edge = lst_ref_edges [ iEdge ]; 
        PolyCurve::VerifyRefEdgeT(edge, t); 
        edge->evaluer ( t, __X );
}
                                                                                           
void
PolyCurve::deriver (double __s, double __X[3]) // curvilinear abcissa
{
        bool curvilinearLength = false;
        deriver(__s, __X, curvilinearLength);
}
 
void
PolyCurve::deriver (double __s, double __X[3], bool __curvilinearLength) // curvilinear abcissa
{                       
        if (IsPoint())
            {
                for (int i=0;i<3;i++) __X[i] = 1;
                return;
            }
             
        VerifyS (__s);
        unsigned iEdge;
        MG_ARETE * edge;
        double t, dt;
        Parameter_SToRefEdgeT ( __s, & iEdge , &t, &dt, __curvilinearLength );
        edge = lst_ref_edges [ iEdge ];   
        PolyCurve::VerifyRefEdgeT(edge, t);
        edge->deriver ( t, __X );
        
        // dC/ds = dC/dt * dt/ds
        for (int i = 0; i<3; i++)
                __X[i] *= dt;
}
 
                                                                                                                                          
void
PolyCurve::deriver_seconde (double __s, double __ddxyz[3], double *__dxyz, double *__xyz) // curvilinear abcissa
{
        bool curvilinearLength = false;
        deriver_seconde(__s, __ddxyz, __dxyz, __xyz, curvilinearLength);
 
}
void
PolyCurve::deriver_seconde (double __s, double __ddxyz[3], double *__dxyz, double *__xyz, bool __curvilinearLength) // curvilinear abcissa
{                   
        unsigned iEdge;
        MG_ARETE * edge;
        double t, dt;
      
        if (IsPoint())
            {
                if (__ddxyz) for (int i=0;i<3;i++) __ddxyz[i]=1;
                if (__dxyz) for (int i=0;i<3;i++) __dxyz[i]=1;
                evaluer(0,__xyz);
                return;
            }
 
        VerifyS (__s);
 
        Parameter_SToRefEdgeT ( __s, & iEdge , &t, &dt, __curvilinearLength );
        edge = lst_ref_edges [ iEdge ];
 
        for (int i=0; i<3; i++)
                __ddxyz[i] = 0;
 
        if (__dxyz)
                edge->deriver(t, __dxyz);
        // dC/ds = dC/dt * dt/ds
        for (int i = 0; i<3; i++)
                __dxyz[i] *= dt;
 
        if (__xyz)
                edge->evaluer(t, __xyz);
}
 
 
 
void
PolyCurve::Merge( PolyCurve & __polycurve)
{
        if (lst_vertices.size() == 0)
        {
                for (std::vector < MG_ARETE * >::const_iterator it = __polycurve.lst_ref_edges.begin();
                        it != __polycurve.lst_ref_edges.end();
                        it++ )
                                InsertCurve ( *it );
        }
        else if (get_sommet2() == __polycurve.get_sommet2() || get_sommet1() == __polycurve.get_sommet2() )
        {;
                for (std::vector < MG_ARETE * >::reverse_iterator it = __polycurve.lst_ref_edges.rbegin();
                        it != __polycurve.lst_ref_edges.rend();
                        it++)
                                InsertCurve ( *it );
        }
        else if ( get_sommet1() == __polycurve.get_sommet1() || get_sommet2() == __polycurve.get_sommet1() )
        {
                for (std::vector < MG_ARETE * >::const_iterator it = __polycurve.lst_ref_edges.begin();
                        it != __polycurve.lst_ref_edges.end();
                        it++ )
                                InsertCurve ( *it );
        }
        else
        {
            printf("PolyCurve::Merge Error : the polycurve %lu is neither adjacent to end vertex nor to start vertex of polycurve %lu\n", __polycurve.get_id(), get_id());
            printf("Polycurve Start vertex = %lu, End vertex = %lu\n", get_sommet1()->get_id(),get_sommet2()->get_id());
            printf("Other __polycurve Start vertex = %lu, End vertex = %lu\n", __polycurve.get_sommet1()->get_id(), __polycurve.get_sommet2()->get_id());
            return;  
 
        }
}
 
void
PolyCurve::SetPeriodicPoleRefVertex(MG_SOMMET * __v)
{
    if ( est_periodique() == false)
        return;
        
    std::vector<MG_ARETE*> e;
    std::vector<MG_SOMMET*> v;
    std::vector<double> l;
    unsigned k=RefVertex_GetIndex(__v);
    unsigned N=lst_vertices.size();
    for (unsigned i=0; i<N; i++)
    {
        unsigned m = (i+k)%N;  
        unsigned n = (i+k+1)%N;           
        MG_SOMMET * tempv = lst_vertices[ m ];
        MG_SOMMET * lastV = (v.size()) ? v[v.size()-1] : 0;
        if ( lastV && lastV == tempv )
            tempv = lst_vertices[n];
        v.push_back(tempv);
    }
    for (unsigned i=0; i+1<N; i++)
        e.push_back(lst_ref_edges[ (i+k)%(N-1) ]);
    for (unsigned i=0; i+1<N; i++)
        l.push_back(lst_length[ (i+k)%(N-1) ]);
    lst_vertices = v;
    lst_ref_edges = e;
    lst_length = l;
}
 
double
PolyCurve::get_longueur (double __s_min, double __s_max, double precision) 
{
        double result;
        if ( __s_min >= 0 && __s_max >= 0 )
                return fabs (__s_min - __s_max);
        else
        {
                return t_max;
        }
}
 
double PolyCurve::get_tmin()
{
        return t_min;
}
 
double PolyCurve::get_tmax()
{
        return t_max;
}
 
MG_SOMMET *
PolyCurve::get_sommet1() 
{
        return lst_vertices[0];
}
 
 
MG_SOMMET *
PolyCurve::get_sommet2() 
{
        return lst_vertices[lst_vertices.size()-1];
}
 
double
PolyCurve::get_sommet1_s()
{
        return get_tmin();
}
 
 
double
PolyCurve::get_sommet2_s() 
{
        return get_tmax();
}
 
void
PolyCurve::inverser(double & __t, double __point[3], double precision)
{
        bool curvilinearLength = false;
 
        inverser(__t, __point, precision, curvilinearLength);
}
 
int PolyCurve::est_periodique(void)
{
        if ( get_sommet1 () == get_sommet2 () )
                return 1;
        else
                return 0;
}
 
double PolyCurve::get_periode(void)
{
    if ( est_periodique() )
        return get_longueur();
    else
        return 0;
}
 
 
int PolyCurve::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
{
        return 4343; // TYPE = Poly curve
}
 
void PolyCurve::enregistrer(std::ostream& o,double version)
{
    unsigned int i;
 
    if (IsPoint() == false)
        {
    o << "%" << get_id()
      << "=CAD4FE_POLYCURVE(" << GetRefEdgeCount() << ",(";
 
    for (i=0; i < GetRefEdgeCount(); i++)
    {
        if ( i ) o <<",";
      o << "$" << GetRefEdge(i)->get_id();
    }                    
 
    o << "));" << std::endl;
        }
    else // The polycurve is a point !
        {
          o << "%"   << get_id() << "=CAD4FE_POLYCURVE(0,($"<<GetRefVertex(0)->get_id()<<"));"<< std::endl;
        }
}
 
 std::vector<MG_ARETE*> &
PolyCurve::GetRefEdges() 
{
        return lst_ref_edges;
}
 
void
CAD4FE::SplitRefEdge(MG_ARETE * __refEdge, MG_SOMMET * __refVertex1, MG_SOMMET * __refVertex2, double __xyz[3], MG_VOLUME * __refBody, MG_GEOMETRIE * __geom, MG_ARETE * edges[2], MG_SOMMET ** __splitVertex)
{
    int i,j;
    double tSplit; 
    MG_ARETE * splitEdge = __refEdge;
    splitEdge->inverser(tSplit, __xyz);
 
    // Create the split vertex at position xyz
    MG_POINT * point = new LC_POINT(__xyz);
    __geom->ajouter_mg_point(point);
    std::ostringstream idoriginal;
    idoriginal << "Vertex of "<< splitEdge->get_idoriginal() <<" split at t="<<tSplit;
    MG_SOMMET * splitVertex = new MG_SOMMET(idoriginal.str(), point);
    *__splitVertex = splitVertex;
    __geom->ajouter_mg_sommet(splitVertex);
 
    // split the edge in two edges
    MG_COURBE * curves[2];
    MG_SOMMET * vertices[2][2];
    MG_COSOMMET * cov[2][2];
                                      
    vertices[0][0] = splitEdge->get_cosommet1()->get_sommet();
    vertices[0][1] = splitVertex;
    vertices[1][0] = splitVertex;
    vertices[1][1] = splitEdge->get_cosommet2()->get_sommet();
                                                   
    curves[0] = curves[1] = splitEdge->get_courbe();
 
    double verticesT[2][2];
    verticesT[0][0] = splitEdge->get_tmin(); 
    verticesT[1][1] = splitEdge->get_tmax();
    verticesT[0][1] = verticesT[1][0] = tSplit;
 
    //  create two edges
    for (i=0; i<2; i++)
    {
        // orientation
         int orientation = splitEdge->get_orientation();
        // Create an edge using the start vertex and the split vertex
        std::ostringstream idoriginal;
        idoriginal << __refEdge->get_idoriginal();
        idoriginal << " split with t0 " << verticesT[i][0] << " t1 " << verticesT[i][1];
        edges[i] = new MG_ARETE (idoriginal.str(), 0, curves[i], orientation);
        __geom->ajouter_mg_arete (edges[i]);
        // create covertices
        for (j=0;j<2;j++)
        {
            cov[i][j] = __geom->ajouter_mg_cosommet (edges[i], vertices[i][j]);
        }
 
        edges[i]->changer_cosommet1(cov[i][0]);
        edges[i]->changer_cosommet2(cov[i][1]);
    }
 
    MG_BOUCLE * split_edge_loops[10];
    int nb_split_edge_loops = 0;
    MG_COARETE * split_edge_coedges[10], * coedges[20];
    int nb_split_edge_coedges = 0, nb_coedges = 0;
    // create four coedges (in manifold body case)
    // replace each coedges belong to initial edge in both loops of the ref faces
    //
        // for each shell of the ref body
        unsigned nb_shells = __refBody->get_nb_mg_coquille ();
        for (unsigned it_shells = 0; it_shells < nb_shells; it_shells++)
        {
                MG_COQUILLE * shell = __refBody->get_mg_coquille(it_shells);
        // for each coface->face F of the shell
                unsigned nb_coface = shell->get_nb_mg_coface();
                for (unsigned it_coface = 0; it_coface < nb_coface; it_coface++)
                {
                        MG_FACE * face = shell->get_mg_coface(it_coface)->get_face();
 
        // for each loop L of this face
                        unsigned nb_loop = face->get_nb_mg_boucle();
                        for (unsigned it_loop = 0; it_loop < nb_loop; it_loop++)
                        {
                                MG_BOUCLE * loop = face->get_mg_boucle(it_loop);
        // for each coedge coe of L
                                unsigned nb_edge = loop->get_nb_mg_coarete();
                                for (unsigned it_edge = 0; it_edge < nb_edge; it_edge++)
                                {
                                        MG_COARETE * coedge = loop->get_mg_coarete(it_edge);
                                        MG_ARETE * edge = coedge->get_arete();
                                          
        // if coe contains split_edge
                                        if (edge == splitEdge)
                                        {
                                            split_edge_loops[ nb_split_edge_loops++ ] = loop;        
                                            split_edge_coedges[ nb_split_edge_coedges++ ] = coedge;
        // create coe[0] ( edges[0], F, coe->sense )
        // create coe[1] ( edges[1], F, coe->sense )
        // delete coe in L
        // add coe[0] in L
        // add coe[1] in L                                            
                                            for (i=0; i<2; i++)
                                            {
                                                coedges[nb_coedges++] = __geom->ajouter_mg_coarete (edges[i], loop, coedge->get_orientation());
                                                loop->ajouter_mg_coarete(coedges[nb_coedges-1]);
                                            }
 
                                            loop->supprimer_mg_coarete(coedge);
                                            __geom->supprimer_mg_coarete(coedge);
                                        }
        // fi
                                }
                        }
                }
        }                          
 
}
 
void
CAD4FE::SplitPolyCurve(PolyCurve * __polyCurve, double __xyz[3], MG_VOLUME * __refBody, MG_GEOMETRIE * __geom, PolyCurve * __result[2], MG_ARETE **__origRefEdge,  MG_SOMMET ** __splitRefVertex, MG_ARETE * __splitRefEdges[2])
{
    int i,j;
    *__splitRefVertex = NULL;
    __splitRefEdges[0]=__splitRefEdges[1]=NULL;
 
    // get the parameter of __splitVertex
    double sSplitVertex;
    __polyCurve->inverser(sSplitVertex,__xyz,1E-6);
 
    unsigned index; // index of the reference edge split
    double tEdge, dtEdge; // real parameter of the reference edge split
    __polyCurve->Parameter_SToRefEdgeT(sSplitVertex,&index,&tEdge,&dtEdge, false);
    *__origRefEdge = __polyCurve->GetRefEdge(index);
    MG_SOMMET * v[2];
    v[0] = __polyCurve->GetRefVertex(index);
    v[1] = __polyCurve->GetRefVertex(index+1);
 
    // Test wether __xyz is located on one extremity of the ref edge
    for (i=0; i<2; i++)
    {
        double v_xyz[3];
 
        v[i]->get_point()->evaluer(v_xyz);
 
        if (OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__xyz, v_xyz) / __polyCurve->GetLength(index) < 1E-6)
        {
            if (index+i == 0 || index+i+1 == __polyCurve->GetRefVertexCount())
            {
                printf("Can't split a polycurve on vertex ID %lu because its extremities are vertices : %lu and %lu\n", v[i]->get_id(), __polyCurve->GetRefVertex(0)->get_id(), __polyCurve->GetRefVertex(__polyCurve->GetRefVertexCount()-1)->get_id());
                return;
            }
 
            *__origRefEdge = NULL;
            *__splitRefVertex = v[i];
            __splitRefEdges[0] = __polyCurve->GetRefEdge(index-1+i);
            __splitRefEdges[1] = __polyCurve->GetRefEdge(index+i);
                      
            // create the two polyCurves
            for (int k=0; k<2; k++)
            {
                // create polycurve
                __result[k] = new PolyCurve ( __splitRefEdges[k] );
                __geom->ajouter_mg_courbe(__result[k]);
                for (j = (k==0)?index-2+i:index+i+1; j >= 0 && j < __polyCurve->GetRefEdgeCount(); j += -1+2*k)
                    __result[k]->InsertCurve(__polyCurve->GetRefEdge(j));
            }
        
            break;
        }
    }
    if (*__splitRefVertex == NULL)
    {
        CAD4FE::SplitRefEdge(*__origRefEdge, v[0], v[1], __xyz, __refBody, __geom, __splitRefEdges, __splitRefVertex);
 
        if ( __splitRefEdges[0]->get_cosommet1()->get_sommet() != v[0] && __splitRefEdges[0]->get_cosommet2()->get_sommet() != v[0] )
            std::swap(__splitRefEdges[0], __splitRefEdges[1] );  
 
        // create the two polyCurves
        for (i=0; i<2; i++)
        {
            // create polycurve
            __result[i] = new PolyCurve ( __splitRefEdges[i] );
            __geom->ajouter_mg_courbe(__result[i]);
            for (j = index-1+2*i; j >= 0 && j < __polyCurve->GetRefEdgeCount(); j += -1+2*i)
                __result[i]->InsertCurve(__polyCurve->GetRefEdge(j));
        }
    }
 
    if (__result[0]->get_sommet1() !=__result[1]->get_sommet1() &&__result[0]->get_sommet1() !=__result[1]->get_sommet2() &&
    __result[0]->get_sommet2() !=__result[1]->get_sommet1() &&__result[0]->get_sommet2() !=__result[1]->get_sommet2())
        printf("Error");
 
}
 
bool PolyCurve::IsPoint() const
{
    return lst_vertices.size() == 1 && lst_ref_edges.size() == 0;
}
 
void
PolyCurve::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
{
return;
}