CAD4FE_MCSegment.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_MCSegment.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:56 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
#include <fstream>
 
#pragma hdrstop
 
#include "gestionversion.h"
#include "CAD4FE_MCSegment.h"      
#include "CAD4FE_MCFace.h"         
#include "CAD4FE_MCEdge.h"
 
#include "ot_algorithme_geometrique.h" 
#include "ot_decalage_parametre.h"
#include "CAD4FE_MCNode.h"  
#include "CAD4FE_PolySurface.h"        
#include "CAD4FE_PolyCurve.h"
#include "CAD4FE_MCVertex.h"    
#include "CAD4FE_MCEdge.h"
#include "CAD4FE_Intersection_Plane_PolySurface.h"
#include "CAD4FE_ShortestPathByClosestPointOnEdge.h"     
#include "CAD4FE_MCNodePolyline.h"                     
#include "CAD4FE_Geometric_Tools.h"
#include "CAD4FE_InventorText_MG_MAILLAGE.h"
#include "ot_root_find.h"
#include "CAD4FE_mailleur2d.h"
 
 
#pragma package(smart_init)
 
using namespace CAD4FE;
 
MCSegment::MCSegment (unsigned long num, MG_ELEMENT_TOPOLOGIQUE *topo, MCNode *mgnoeud1, MCNode *mgnoeud2, double longue)
: MG_SEGMENT(num,topo,mgnoeud1,mgnoeud2,longue),_polylineEvaluator(0),_sense(1), _saveFormat(0)
{
    _cstrMCFaceScheme = CstrMCFaceByPlaneIntr;
    ConstructGeometry();
}     
MCSegment::MCSegment (MG_ELEMENT_TOPOLOGIQUE *topo, MCNode *mgnoeud1, MCNode *mgnoeud2, double longue, long __cstrMCFaceScheme)
: MG_SEGMENT(topo,mgnoeud1,mgnoeud2,longue),_polylineEvaluator(0),_sense(1),_cstrMCFaceScheme(__cstrMCFaceScheme), _saveFormat(0)
{
    ConstructGeometry();
}
MCSegment::MCSegment( MCSegment & __s)
: MG_SEGMENT(0,__s.liaison_topologique,__s.get_noeud1(),__s.get_noeud2(),__s.get_longueur()), _saveFormat(0)
{
    unsigned i;
    for (FMapCIterator itF = __s.F.begin(); itF != __s.F.end(); itF++)
    {
        MG_FACE * face = itF->first;
        F.insert(std::make_pair(face,itF->second));
    }
    for (EMapCIterator itE = __s.E.begin(); itE != __s.E.end(); itE++)
    {
        MG_ARETE * edge = itE->first;
        E.insert(std::make_pair(edge,itE->second));
    }
    //V = __s.V;
    _polylineEvaluator=new MCNodePolyline(*__s._polylineEvaluator);
    _sense = __s._sense;
}                     
MCSegment::FMap & MCSegment::GetRefFaceMapping()
{
    return F;
}
MCSegment::EMap & MCSegment::GetRefEdgeMapping()
{
    return E;
}
MCSegment::~MCSegment()
{
    unsigned i,j;       
    if (_polylineEvaluator) delete _polylineEvaluator;
}
std::vector <MCNode*> MCSegment::GetPolylineNodes()
{
    return _polylineEvaluator->GetPolylineNodes();
}
MCNode* MCSegment::GetPolylineNode(unsigned __index)
{
    return _polylineEvaluator->GetPolylineNode(__index);
}
std::vector <MG_ELEMENT_TOPOLOGIQUE*> MCSegment::GetPolylineTopos()
{
    return _polylineEvaluator->GetPolylineTopos();
}
MG_ELEMENT_TOPOLOGIQUE* MCSegment::GetPolylineTopo(unsigned __index)
{
    return _polylineEvaluator->GetPolylineTopo(__index);
}
unsigned MCSegment::GetPolylineNodeCount()
{
    return _polylineEvaluator->GetPolylineNodes().size();
}
double MCSegment::get_longueur_geo()
{
    double length = _polylineEvaluator->GetLength();
    return length;
}                                                                   
void MCSegment::evaluer_geo(double __t, MCNode * __result, double * tangent, double * curvature)
{
    double t,L;
    if (_sense == -1)
        t = (1-__t);
    else
        t = __t;
 
    L=get_longueur_geo();
    * __result = _polylineEvaluator->Evaluate( t * L, tangent, curvature );
    if (tangent)
        for (int i=0; i<3; i++) tangent[i] *= L;
    if (_sense == -1)
    {
        if (tangent)
            for (int i=0; i<3; i++) tangent[i]*=-1;
    }
}
void MCSegment::inverser_geo(double & __t, MCNode * __node)
{
    _polylineEvaluator->Inverse(__t, __node);
    double length = _polylineEvaluator->GetLength();
    __t = __t/length;
 
    if (_sense == -1)
        __t = (1-__t);
}                                      
int MCSegment::get_orientation()
{
    return _sense;
}                                       
void MCSegment::change_orientation(int __sense)
{
    _sense = __sense;
}                            
void MCSegment::UpdateGeometry ()
{              
    calcule_longueur();
    F.clear();
    E.clear();
    //V.clear();
    delete _polylineEvaluator;
    _polylineEvaluator = 0;
    ConstructGeometry ();
}
void MCSegment::ConstructGeometry ()
{
    MCNode * mcNode1 = (MCNode*)get_noeud1();
    MCNode * mcNode2 = (MCNode*)get_noeud2();
    CAD4FE::MCFace * mcFace = NULL;
    CAD4FE::MCEdge * mcEdge = NULL;
    if (get_lien_topologie()->get_dimension() == 2)
    {
        mcFace = (MCFace*)get_lien_topologie();
    }
    if (get_lien_topologie()->get_dimension() == 1)
    {
        mcEdge = ((CAD4FE::MCEdge*)get_lien_topologie());
    }
    if (mcEdge)
    {
        Construct_MCEdge(mcEdge);
        Construct_MergedVertices(mcNode1);
        Construct_MergedVertices(mcNode2);
    }
    else if (mcFace)
    {
        int res;
        if (_cstrMCFaceScheme == CstrMCFaceByPlaneIntr )
        {
            res = Construct_MCFace(mcFace,mcNode1,mcNode2);
            if (res != 1 )
            {
                res = Construct_MCFaceByShortestPath(mcFace,mcNode1,mcNode2);
                if (res == 1) _cstrMCFaceScheme = CstrMCFaceByShortestPath;
            }
        }
        else if (_cstrMCFaceScheme == CstrMCFaceByPlaneIntrNormalByShortestPath)
        {
            res = Construct_MCFace(mcFace,mcNode1,mcNode2);
            if (res != 1)
            {                                
                res = Construct_MCFaceByShortestPath(mcFace,mcNode1,mcNode2);
                if (res == 1)
                    _cstrMCFaceScheme = CstrMCFaceByShortestPath;
            }
        }
        else if (_cstrMCFaceScheme == CstrMCFaceByShortestPath )
        {
            res = Construct_MCFaceByShortestPath(mcFace,mcNode1,mcNode2);
            if (res != 1)
            {  
                res = Construct_MCFace(mcFace,mcNode1,mcNode2);
                if (res == 1)
                    _cstrMCFaceScheme = CstrMCFaceByPlaneIntr;
            }
        }
        else if (_cstrMCFaceScheme == CstrNoRefTopologyMapping)
        {
            res = Construct_NoTopology(mcNode1,mcNode2);
        }
    }
 
    _sense = 1;
    ConstructBoundaryBox();
}
 
void MCSegment::ConstructBoundaryBox()
{
    for (unsigned j=0; j<_polylineEvaluator->GetPolylineNodes().size();j++)
    {
        MCNode * n = _polylineEvaluator->GetPolylineNode(j);
        double * x = n->get_coord();
        if (j == 0)
        {
            for (int i=0;i<3;i++)
                _bbox[i] = _bbox[i+3] = x[i];
        }
        else
        {
            for (int i=0;i<3;i++)
            {
                if (_bbox[i]>x[i])
                    _bbox[i] = x[i];
                else if (_bbox[i+3]<x[i])
                    _bbox[i+3] = x[i];
            }
        }
    }
    for (double t=.25; t<=.75; t+=.25)
    {
        MCNode n = _polylineEvaluator->Evaluate(t);
        double * x = n.get_coord();
 
            for (int i=0;i<3;i++)
            {
                if (_bbox[i]>x[i])
                    _bbox[i] = x[i];
                else if (_bbox[i+3]<x[i])
                    _bbox[i+3] = x[i];
            }
    }
}
 
extern MG_MAILLAGE * debugMesh;
 
int MCSegment::Construct_MCFaceByShortestPath(MCFace * __mcFace, MCNode* __n1, MCNode* __n2)
{
    unsigned i,j;
    double dN1N2=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__n1->get_coord(),__n2->get_coord());
    int shortestPathOptions = 0;
    ShortestPathByClosestPointOnEdge shortestPath(__mcFace,__n1,__n2,dN1N2, shortestPathOptions);
    ShortestPathByClosestPointOnEdge *shortestPathDoubleDist = NULL;
    std::vector<MCNode*> vecPts;
    std::vector<MG_ELEMENT_TOPOLOGIQUE*> vecTopo;
    double distance = shortestPath.Find(&vecPts,&vecTopo);
    
    if (distance > 100) // failure, second try with 2*dist(n1,n2)
    {
        unsigned maxNbTries = 4;
        unsigned maxDistFactor = 1;
        for (i=0; i<maxNbTries; i++)
        {
            maxDistFactor *= 2;
            vecPts.clear();vecTopo.clear();
            shortestPathDoubleDist = new ShortestPathByClosestPointOnEdge (__mcFace,__n1,__n2,maxDistFactor*dN1N2, shortestPathOptions);
            distance = shortestPathDoubleDist->Find(&vecPts,&vecTopo);
            if (distance > 100)
            {
                delete shortestPathDoubleDist;     
                vecPts.clear();vecTopo.clear();
                shortestPathOptions=ShortestPathByClosestPointOnEdge::NoCheckFaceIntr;
                shortestPathDoubleDist = new ShortestPathByClosestPointOnEdge (__mcFace,__n1,__n2,maxDistFactor*dN1N2, shortestPathOptions);
                distance = shortestPathDoubleDist->Find(&vecPts,&vecTopo);
                if (distance > 100)
                {
                    delete shortestPathDoubleDist;
                    shortestPathDoubleDist = NULL;
                }
            }
            else
            {
                break;
            }
        }
    }
 
    if (distance > 100)
    {
        printf("Warning : Construct_MCFaceByShortestPath failed !\n");
        distance = shortestPath.Find(&vecPts,&vecTopo);
        Construct_MCFace(__mcFace,__n1,__n2);
        return 0;
    }
    if (distance < .8*dN1N2)
    {                
        printf("Warning : %f < .8*%f in Construct_MCFaceByShortestPath !\n", distance, dN1N2);
    }
 
    if (_polylineEvaluator == NULL) _polylineEvaluator = new MCNodePolyline(liaison_topologique);
 
    for (i=0;i+1<vecPts.size();i++)
    {
        MCNode * xi=vecPts[i], *xii=vecPts[i+1];
        MG_ELEMENT_TOPOLOGIQUE * topo = vecTopo[i];
 
        if (topo->get_dimension() == 2)
        {
                MG_FACE * face = (MG_FACE *) topo;
                std::vector<MCNode*> Ci;
                OT_VECTEUR_3D normalXi,normalXii,normalAvg,dir,v,xi_xyz(xi->get_coord()),xii_xyz(xii->get_coord());
                xi->NormalMCFace(__mcFace, normalXi);
                xii->NormalMCFace(__mcFace, normalXii);
                normalAvg = .5*normalXi+.5*normalXii;
                dir=xii_xyz-xi_xyz;
                v=normalAvg&dir;
                double L=dir.get_longueur();
                Intersection_Plane_MG_FACE * intr = NULL;
                if (L > 1E-6 && ((shortestPathOptions&ShortestPathByClosestPointOnEdge::NoCheckFaceIntr)==0))
                {
                    intr = new Intersection_Plane_MG_FACE(xi,v,face,__mcFace);
                    intr->MakeSegment(xi,xii,L/50,L/3000,3,dir);
                    Ci = intr->Get()[0];
                    if (Ci.size()<2||Ci[Ci.size()-1]!=xii||Ci[0]!=xi)
                    {
                        delete intr;
                        intr = 0;
                        Ci.clear();
                        Ci.push_back(xi);
                        Ci.push_back(xii);
                    }
                }
                else
                {
                    Ci.push_back(xi);
                    Ci.push_back(xii);
                }
                F.insert( std::make_pair(face, Ci) );
 
                for (j=0; j<Ci.size(); j++)
                {
                    _polylineEvaluator->Add(Ci[j],face);
                }
                if (intr != NULL)
                    delete intr;
        }
        if (topo->get_dimension() == 1)
        {
            MG_ARETE * edge = (MG_ARETE *) topo;
 
            //double dist  = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xi->get_coord(), xii->get_coord());
            //double epsilon = dist*1E-2;
            //int insertedNodeCount = RefineInRefEdge(__mcFace, edge, xi, xii, epsilon, 1E300, refinedNodes, refinedEdges );
 
            VPts Ci;
            Ci.push_back(xi);
            Ci.push_back(xii);
 
            E.insert (std::make_pair(edge, Ci));
            for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
            {
                MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
                F.insert(std::make_pair(face, Ci));
            }
 
            for (j=0; j<Ci.size(); j++)
            {
                _polylineEvaluator->Add(Ci[j],edge);
            }
        }  
        if (topo->get_dimension() == 0)
        {
                
        }
    }
 
    if (shortestPathDoubleDist != NULL)
        delete shortestPathDoubleDist;
 
    return 1;
}
 
int MCSegment::Construct_NoTopology(MCNode * __n1, MCNode * __n2)
{
    _polylineEvaluator = new MCNodePolyline(liaison_topologique);
    _polylineEvaluator->Add(__n1,NULL);
    _polylineEvaluator->Add(__n2,NULL);
    return 1;
}
 
int MCSegment::Construct_MCFace(MCFace * __mcFace, MCNode* __n1, MCNode* __n2)
{
    unsigned i,j;
    PolySurface * polysurface = __mcFace->GetPolySurface();
                                                  
    OT_VECTEUR_3D normal1, normal2, xyz1, xyz2, v1, v2;    
    if ( _cstrMCFaceScheme != CstrMCFaceByPlaneIntrNormalByShortestPath )
    {    
            int nbRefFaces;
            polysurface->calcul_normale_unitaire(__n1->GetRefFaceMapping(),normal1,&nbRefFaces);
            polysurface->calcul_normale_unitaire(__n2->GetRefFaceMapping(),normal2,&nbRefFaces);
            v2 = .5*(normal1+normal2);
            if (v2.get_longueur2() < 1E-100)
            {
                printf("Warning: surface's normal at node 1 is opposite to surface's normal at node 2\n");
                printf("\t trying to evaluate surface's normal with the shortest path algorithm\n"); 
                _cstrMCFaceScheme = CstrMCFaceByPlaneIntrNormalByShortestPath;
            }
    }
    if (_cstrMCFaceScheme==CstrMCFaceByPlaneIntrNormalByShortestPath)
    {
        double maxDistBetween3DSegmentAndShortestPath=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__n1->get_coord(),__n2->get_coord());
        ShortestPathByClosestPointOnEdge * shortestPath = new ShortestPathByClosestPointOnEdge(__mcFace,__n1,__n2,maxDistBetween3DSegmentAndShortestPath);
        std::vector<MCNode*> vecPts;
        std::vector<MG_ELEMENT_TOPOLOGIQUE*> vecTopo;
        double distance=shortestPath->Find(&vecPts, &vecTopo);
        if (distance > 100)
        {
            double *x1=__n1->get_coord(),*x2=__n2->get_coord();
            printf("%f %f %f %f %f %f ",x1[0],x1[1],x1[2],x2[0],x2[1],x2[2]);
            if (distance > 100)
                printf("Error in average normal\n");
 
        }
        if (distance == 0.0)
        {            
            delete shortestPath;
            return 0;
        }
        // estimation of the face normal by the average of the normals
        // in the polyline
        OT_VECTEUR_3D avgNormal(0,0,0);
        MCNodePolyline polyline(liaison_topologique,vecPts,vecTopo);
        unsigned nb_points=10, real_nb_points=0;
        double s_max=polyline.GetLength();
        for (double s=0; real_nb_points < 1 || s < s_max; s+=.9999*(s_max/nb_points))
        {
            real_nb_points++;
            MCNode mcNode=polyline.Evaluate(s);
            OT_VECTEUR_3D tmpNormal;
            int nbRefFaces;
            polysurface->calcul_normale_unitaire(mcNode.GetRefFaceMapping(),tmpNormal,&nbRefFaces);
            avgNormal += tmpNormal;
        }
        avgNormal /= real_nb_points;
        v2=avgNormal;
        delete shortestPath;
    }
 
    xyz1=OT_VECTEUR_3D(__n1->get_coord());
    xyz2=OT_VECTEUR_3D(__n2->get_coord());
    v1=xyz2-xyz1;
 
    // Make a special treatment to manage nodes of merged vertices
    MCNode * polysurfaceNodes[2]={__n1,__n2};
    MCNode * n[2]={__n1,__n2};
    VPts MergedVertexCurveNodes[2];
    std::vector < MG_ELEMENT_TOPOLOGIQUE * > MergedVertexCurveTopo[2];
    for (int i=0;i<2;i++)
    if (n[i]->GetRefVertexMapping().size() > 1) // if node is in a merged vertex
    {
        MCVertex * mcVertex = (MCVertex*) n[i]->get_lien_topologie();
        MG_SOMMET * refVertex = (MG_SOMMET*) n[i]->get_lien_topologie_reference();
        MG_SOMMET * mergedVertex = NULL;
        if (GeometricTools::PolySurface_Contains_RefVertex(polysurface,refVertex) == false)
        {
            for (MCNode::VMapIterator itV = n[i]->GetRefVertexMapping().begin();
                itV != n[i]->GetRefVertexMapping().end();
                itV++)
                if ( GeometricTools::PolySurface_Contains_RefVertex(polysurface,*itV) )
                {
                    mergedVertex = *itV;
                    break;
                }
 
            double xyz[3];mergedVertex->get_point()->evaluer(xyz);
            double d=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xyz,n[i]->get_coord());
            MCNode * mergedVertexNode = new MCNode(mcVertex,mergedVertex,xyz[0],xyz[1],xyz[2]);
            polysurfaceNodes[i] = mergedVertexNode;
 
            MCFace * mcFacePath = NULL;
            std::set<MG_FACE*> mergedVertexMCFaces;
            GeometricTools::MG_SOMMET_GetAdjacent_MG_FACE(mcVertex,mergedVertexMCFaces);
            for (std::set<MG_FACE*>::iterator itF=mergedVertexMCFaces.begin();
                itF != mergedVertexMCFaces.end(); itF++)
            {
                MCFace* adjacentMCFace = (MCFace*)(*itF);
                if (GeometricTools::PolySurface_Contains_RefVertex(adjacentMCFace->GetPolySurface(),mergedVertex)
                && GeometricTools::PolySurface_Contains_RefVertex(adjacentMCFace->GetPolySurface(),refVertex)
                )
                {
                    mcFacePath = adjacentMCFace;
                    break;
                }
            }
            ShortestPathByClosestPointOnEdge shortestPath(mcFacePath,n[i],mergedVertexNode,d);
            double distance = shortestPath.Find(&(MergedVertexCurveNodes[i]),&(MergedVertexCurveTopo[i]));
 
            if (distance < 100)
                for (unsigned j=0; j<MergedVertexCurveNodes[i].size();j++)
                    MergedVertexCurveNodes[i][j]->incrementer();
        }
    }
 
    Intersection_Plane_PolySurface intrPolySurface (polysurfaceNodes[0], v2&v1, __mcFace, polysurfaceNodes[1]);
    MCNodePolyline * intersectionPlaneMCFace = intrPolySurface.MakeSegment(polysurfaceNodes[0], polysurfaceNodes[1], v1);
 
    if (intersectionPlaneMCFace != NULL)
    {
        _polylineEvaluator = intersectionPlaneMCFace;
        _polylineEvaluator->incrementer();
 
        Intersection_Plane_PolySurface::FMapIterator it;
        Intersection_Plane_MG_FACE * intr = NULL;
        for (intr = intrPolySurface.GetFirstCurve(it); intr; intr = intrPolySurface.GetNextCurve(it) )
        {
            MG_FACE * face = intr->GetFace();
            VPts nodes = intr->Get()[0];
            F.insert(std::make_pair(face, nodes));
        }
 
        for ( unsigned i=0;1+i<GetPolylineNodeCount();i++)
        {
            MG_ELEMENT_TOPOLOGIQUE * topo = GetPolylineTopo(i);
            if ( topo != NULL && topo->get_dimension()==1)
            {
                MCNode * xi = GetPolylineNode(i);
                MCNode * xii = GetPolylineNode(i+1);
                MG_ARETE * edge = (MG_ARETE *) topo;
                VPts Ci;
                Ci.push_back(xi);
                Ci.push_back(xii);
 
                E.insert (std::make_pair(edge, Ci));
                for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
                {
                    MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
                    F.insert(std::make_pair(face, Ci));
                }
            }
        }
 
        return 1;
    }    
    else
    {
        return 0;
    }
}
 
int MCSegment::RefineInRefEdge(MCEdge * __mcEdge, MG_ARETE * __refEdge, MCNode * __n1, MCNode * __n2, double __epsilon, double __lastDistanceToEdge, std::map < double , MCNode * > & __mapMCEdgeNodesBySParameter, std::map < double , MG_ARETE * > & __mapMCEdgeRefEdgesBySParameter, int __iRecursiveCalls)
{
    static const int RefineInRefEdge_MAX_RECURSIVE_CALL_COUNT = 100;
    if (__iRecursiveCalls > RefineInRefEdge_MAX_RECURSIVE_CALL_COUNT)
    {
        printf("Warning : RefineInRefEdge : The number of recursive calls has been reached (%d)\n",RefineInRefEdge_MAX_RECURSIVE_CALL_COUNT);
        return 0;
    }
 
    int insertedNodeCount = 0;
    MCNode * n[2] = { __n1, __n2};
    for (int i=0; i<2; i++)
    {
        MG_SOMMET * v = 0;
        if (n[i]->get_lien_topologie_reference()->get_dimension() == 0)
            v = (MG_SOMMET * ) n[i]->get_lien_topologie_reference();
        if (v)
        if (__mcEdge->GetPolyCurve()->Contains( v ) == false )
            return 0;
        if (n[i]->IsInEdge(__refEdge) == false)
            return 0;
        // nodes have merged ref vertices
        if (n[(i+1)%2]->IsInVertex(v))
            return 0;
    }                               
    __n1->incrementer();__n2->incrementer();
    double t1 = __n1->T(__refEdge);
    double t2 = __n2->T(__refEdge);
    if (fabs(t2-t1) < 1E-6*(__refEdge->get_tmax()-__refEdge->get_tmin()))
        return 0;
 
    double period = __refEdge->get_courbe()->get_periode();
    if (period != 0.0 && fabs(t2 - t1) > .5*period)
        if (t2<t1)
            t2 += period;
        else t1 += period;
    double t3 = t1 + (t2-t1)*.5;
    double s3;
    double xyz[3]; __refEdge->evaluer(t3,xyz);
    MCNode edgeNode (__mcEdge, __refEdge, t3, xyz);
    for (int i=0; i<__refEdge->get_nb_mg_coarete(); i++)
    {
        MG_FACE * face = __refEdge->get_mg_coarete(i)->get_boucle()->get_mg_face();
 
        MCNodePolyline faceLine(liaison_topologique);
        faceLine.Add(__n1,face);
        faceLine.Add(__n2,face);
        double L =  faceLine.GetLength();
        MCNode faceNode = faceLine.Evaluate(L*.5);
        double distanceToEdge = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(faceNode.get_coord(),edgeNode.get_coord());
        double ratioConvergence = distanceToEdge / __lastDistanceToEdge;
        double distN1N2 = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__n2->get_coord(),__n1->get_coord());
        bool testConvergence = ( fabs(ratioConvergence - 1) < .001 || distN1N2 < __epsilon );
        if ( (distanceToEdge > __epsilon  )
            && testConvergence == false )
        {
            insertedNodeCount++;
            __mcEdge->GetPolyCurve()->Parameter_RefEdgeTToS ( t3, __refEdge, &s3, false);
            if (__mcEdge->GetPolyCurve()->est_periodique()
            && __mcEdge->get_tmin() > s3+1E-5*(__mcEdge->GetPolyCurve()->get_periode())
               )
                s3  += __mcEdge->GetPolyCurve()->get_periode();
            MCNode * newMCNode = new MCNode ( edgeNode );
            __mapMCEdgeNodesBySParameter.insert(std::make_pair(s3,newMCNode));
            __mapMCEdgeRefEdgesBySParameter.insert(std::make_pair(s3, __refEdge));
            newMCNode->incrementer();
            insertedNodeCount += RefineInRefEdge(__mcEdge, __refEdge, __n1, newMCNode, __epsilon, distanceToEdge, __mapMCEdgeNodesBySParameter,__mapMCEdgeRefEdgesBySParameter, __iRecursiveCalls+1);
            insertedNodeCount += RefineInRefEdge(__mcEdge, __refEdge, newMCNode, __n2, __epsilon, distanceToEdge, __mapMCEdgeNodesBySParameter,__mapMCEdgeRefEdgesBySParameter, __iRecursiveCalls+1);
            newMCNode->decrementer();
            break;
        }
    }
    __n1->decrementer();__n2->decrementer();
    return insertedNodeCount;
}
 
int MCSegment::Construct_MCEdge(MCEdge * __mcEdge)
{
    double s[2];        
    double t,dt;
    int i,j;
    unsigned iEdge;
    MG_ARETE * edge;
    PolyCurve * polycurve = __mcEdge->GetPolyCurve();
    double period = polycurve->get_periode();
    double smin = __mcEdge->get_tmin();
    double smax = __mcEdge->get_tmax();
 
    VPts CurveNodes;
    std::vector < MG_ELEMENT_TOPOLOGIQUE * > CurveTopo;
 
    MCNode * n[2];
    n[0]=(MCNode *)get_noeud1();
    n[1]=(MCNode *)get_noeud2();
 
    std::map < double , MCNode * > nodes;
    std::map < double , MG_ARETE * > edges;
 
    // get the parameter of extremities n1 n2 of the segment
    for (i=0; i<2; i++)
    {
        if (n[i]->GetRefVertexMapping().size() >= 1) // MC node is on a ref vertex
        {
            MG_SOMMET * refVertex = (MG_SOMMET*)n[i]->get_lien_topologie_reference();
            if (polycurve->ContainsRefVertex(refVertex) == false)
            {
                for (std::set<MG_SOMMET*>::iterator itV = n[i]->GetRefVertexMapping().begin();
                    itV !=  n[i]->GetRefVertexMapping().end();
                    itV++)
                {
                    MG_SOMMET * v = *itV;
                    if (polycurve->ContainsRefVertex(v))
                    {
                        s[i] = polycurve->RefVertex_GetS(v);
                        if (n[i]->get_lien_topologie_reference() != v)
                        {
                            double xyzV[3]; v->get_point()->evaluer(xyzV);
                            n[i] = new MCNode (__mcEdge,v,xyzV[0],xyzV[1],xyzV[2]);
                            //n[i]->incrementer();
                        }
                        break;
                    }
                }
            }
            else
            {
                s[i] = polycurve->RefVertex_GetS(refVertex);
            }
        }
        else if (n[i]->GetRefEdgeMapping().size() == 1) // MC node is in a ref edge
        {
            polycurve->inverser(s[i], n[i]->get_coord(), 1E-6);
        }
    }
                        
    if (period != 0.0)
    {                    
        for (i=0;i<2;i++){
        if (s[i] > smax + fabs(smax-smin)*1E-6)
            s[i] -= period;
        if (s[i] < smin - fabs(smax-smin)*1E-6)
            s[i] += period; }
    }
 
    if (period != 0.0 && fabs(s[1] - s[0]) > .5*period)
        if (s[1]<s[0])
            s[1] += period;
        else s[1] += period;
    
    for (i=0; i<2; i++)
    {
        nodes[s[i]] = n[i];
    }
 
    {
        MG_ARETE * edge; double t,dt;
        double si = std::min(s[0],s[1]);
        polycurve->Parameter_SToRefEdgeT ( si+1E-6*fabs(s[1]-s[0]), &edge,&t,&dt,false);
        edges [ si ] = edge;
    }
           
 
    for (int i=0; i<2; i++)
    {
        OT_VECTEUR_3D v, v_node(n[i]->get_coord());
        polycurve->evaluer(s[i], v);
        double dist = (v_node-v).get_longueur();
        if (dist > .001*__mcEdge->GetPolyCurve()->get_longueur())
        {
            printf("Distance between projection and node = %f\n", dist);
            polycurve->inverser(s[i], n[i]->get_coord(), 1E-6);     
            polycurve->evaluer(s[i], v);
        }
    }
 
    // Create nodes in vertices between n1 and n2
    double si=0;
    for (i=0;si<std::max(s[1],s[0]);i++)
    {
        int iEdge = i%polycurve->GetRefEdgeCount();
        int iVertex = (period != 0.0)? i%(polycurve->GetRefVertexCount()-1): i%polycurve->GetRefVertexCount();
        if ((si-s[0])*(s[1]-si) > 0)
        {
            MG_SOMMET * vertex = polycurve->GetRefVertex(iVertex);
            double xyz[3];vertex->get_point()->evaluer(xyz);
            MCNode * ni = new MCNode ( __mcEdge, vertex, xyz[0], xyz[1], xyz[2] );
            nodes [ si ] = ni;
            edges [ si ] = polycurve->GetRefEdge(iEdge);
            //ni->incrementer();
        }
        si += polycurve->GetLength(iEdge);
    }
 
    std::map < double , MCNode * > refinedNodes = nodes;
    std::map < double , MG_ARETE * > refinedEdges = edges;
 
    for (std::map<double, MCNode * >::iterator itN = nodes.begin();
        itN != nodes.end();
        itN++)
    {
        std::map<double, MCNode * >::iterator itNii = itN;
        itNii++;
        if (itNii == nodes.end())
            break;
        MCNode * ni = itN->second;
        MCNode * nii = itNii->second;
        double si = itN->first;
        double sii = itNii->first;
        double smid = .5*(si+sii);
 
        if (ni->get_lien_topologie_reference()->get_dimension() == 0)
        {
            MG_SOMMET * vi = (MG_SOMMET*)ni->get_lien_topologie_reference();
            if (nii->IsInVertex(vi))
                continue;
        }
 
        if (nii->get_lien_topologie_reference()->get_dimension() == 0)
        {
            MG_SOMMET * vii = (MG_SOMMET*)nii->get_lien_topologie_reference();
            if (ni->IsInVertex(vii))
                continue;
        }
 
        MG_ARETE * refEdge=0; double tMid; double dtMid;
        polycurve->Parameter_SToRefEdgeT(smid,&refEdge,&tMid,&dtMid,false);
        // and if the edge is periodic and tii-ti > period/2 then insert a node in the middle
        {
            if (refEdge->get_courbe()->est_periodique() != 0 )
            {
                MG_ARETE * refEdgei=0, * refEdgeii=0;
                double tii, dtii, ti, dti;
                polycurve->Parameter_SToRefEdgeT(si,&refEdgei,&ti,&dti,false);
                polycurve->Parameter_SToRefEdgeT(sii,&refEdgeii,&tii,&dtii,false);
                
                if (fabs(tii-ti) > .5*refEdge->get_courbe()->get_periode()) 
                {
                    double xyz[3];
                    refEdge->evaluer(tMid,xyz);
                    MCNode * newMCNode = new MCNode(__mcEdge, refEdge, tMid, xyz);
                    //newMCNode->incrementer();
                    refinedNodes[smid] = newMCNode;
                    refinedEdges[smid] = refEdge;
                }
            }
        }
    }
    nodes = refinedNodes;
    edges = refinedEdges;
 
    for (std::map<double, MCNode * >::iterator itN = nodes.begin();
        itN != nodes.end();
        itN++)
    {
        std::map<double, MCNode * >::iterator itNii = itN;
        itNii++;
        if (itNii == nodes.end())
            break;
        MCNode * ni = itN->second;
        MCNode * nii = itNii->second;
        double si = itN->first;
 
        MG_ARETE * refEdge=edges[si];
        
        // Refine the mc segment discretization with discretization error in faces < segment length / 1000
        double deltaMax = get_longueur()*1E-2;
        if (deltaMax < 1E-7)
            deltaMax = 1E-7;
        RefineInRefEdge(__mcEdge, refEdge, ni, nii, deltaMax, 1E300, refinedNodes, refinedEdges );
    }
    nodes=refinedNodes;
    edges=refinedEdges;
 
 
    // detect edges which contain nodes, and make segments of this edge on these nodes
    for (std::map<double, MCNode * >::iterator itN = nodes.begin();
        itN != nodes.end();
        itN++)
    {
        std::map<double, MCNode * >::iterator itNii = itN;
        itNii++;
        if (itNii == nodes.end())
            break;
                           
        MCNode * ni = itN->second;
        MCNode * nii = itNii->second;
        double si = itN->first;
        std::map<double, MG_ARETE*>::iterator itRefEdge = edges.find(si);
        if (itRefEdge != edges.end())
            edge = itRefEdge->second;
        else
            edge = 0;
        if (ni->IsInEdge(edge) == false ||
            nii->IsInEdge(edge) == false )
        {
            for (MCNode::EMapIterator itNiRefEdge = ni->GetRefEdgeMapping().begin();
                itNiRefEdge != ni->GetRefEdgeMapping().end();
                itNiRefEdge++)
            {
                edge = itNiRefEdge->first;
                if ( nii->IsInEdge(edge) )
                {
                    break;
                }
            }
        }
 
        double params[2];
        params[0]=itN->first;
        params[1]=itNii->first;
 
        for (int i=0;i<2;i++)
        if ((s[1]-params[i])*(s[0]-params[i])>0)
            printf("pb\n");
 
        std::vector < MCNode * > segEdge;
        segEdge.push_back (ni);
        segEdge.push_back (nii);
 
        E.insert(std::make_pair(edge, segEdge));
        // copy the curve in adjacent faces
        for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
        {
            MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
            F.insert(std::make_pair(face, segEdge));
        }      
 
        if (CurveNodes.size () == 0 || CurveNodes[CurveNodes.size()-1] != ni)
            CurveNodes.push_back(ni);
        CurveNodes.push_back(nii);
        CurveTopo.push_back(edge);
    }
 
    if (_polylineEvaluator == NULL) _polylineEvaluator = new MCNodePolyline(liaison_topologique,CurveNodes,CurveTopo);
    else
            for (unsigned i=0;i<CurveNodes.size();i++)
                _polylineEvaluator->Add(CurveNodes[i],CurveTopo[i]);
                
    return 1;
}
 
MCFace  * MCSegment::Find_MCFace(MCEdge * __mcEdge, MG_SOMMET * __refVertex1, MG_SOMMET * __refVertex2)
{
    for (int itF = 0; itF != __mcEdge->get_nb_mg_coarete(); itF++)
    {
        MCFace * F = (MCFace*)__mcEdge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
        bool v1Found = false, v2Found = false;
        for (int itFLoop = 0; itFLoop != F->get_nb_mg_boucle(); itFLoop ++)
        {
            MG_BOUCLE * loop = F->get_mg_boucle(itFLoop);
            for (int itFLoopEdge = 0; itFLoopEdge != loop->get_nb_mg_coarete(); itFLoopEdge++)
            {
                MCEdge * mcEdge = (MCEdge *)loop->get_mg_coarete( itFLoopEdge )->get_arete();
                if ( mcEdge->GetPolyCurve()->ContainsRefVertex(__refVertex1) )
                    v1Found = true;                                
                if ( mcEdge->GetPolyCurve()->ContainsRefVertex(__refVertex2) )
                    v2Found = true;   
                if (v1Found && v2Found)
                    return F;
            }
        }
    }
    return NULL;
}
 
int MCSegment::Construct_MergedVertices(MCNode * __mcNode)
{
    if (get_lien_topologie()->get_dimension() != 1)
    {
        printf("Warning: MCSegment topology is not a MCEdge ! then it's not possible to construct a mapping on a merged vertex\n");
        return 0;
    }
    MCEdge * mcEdge = (MCEdge*)get_lien_topologie();
    MG_SOMMET * polycurveRefVertex = NULL;
    MG_SOMMET * mcVertexRefVertex = (MG_SOMMET*)__mcNode->get_lien_topologie();
    MCFace * mcFace;
    for (std::set<MG_SOMMET*>::iterator itV = __mcNode->GetRefVertexMapping().begin();
        itV != __mcNode->GetRefVertexMapping().end();
        itV++ )
    {
        MG_SOMMET * v = *itV;
        if (v != mcVertexRefVertex)
        {
            polycurveRefVertex = v;
            mcFace = Find_MCFace(mcEdge, polycurveRefVertex, mcVertexRefVertex);
            if (mcFace)
            {
                MCNode * polycurve_mcNode = NULL;
                double v_xyz[3];
                v->get_point()->evaluer(v_xyz);
 
                polycurve_mcNode = new MCNode (get_lien_topologie(), v, v_xyz[0],v_xyz[1],v_xyz[2]);
 
                polycurve_mcNode->incrementer();
                if (polycurve_mcNode)
                    Construct_MCFaceByShortestPath(mcFace, __mcNode, polycurve_mcNode);
                    
                break;
            }
        }
    }
 
    return 1;
}
 
BOITE_3D MCSegment::get_boite_3D(void)
{
    return BOITE_3D(_bbox[0],_bbox[1],_bbox[2],_bbox[3],_bbox[4],_bbox[5]);
}
 
void MCSegment::SetSaveFormat(char __format)
{
    _saveFormat=__format;
    if (__format == 2)
    {
        MG_ELEMENT_TOPOLOGIQUE * refTopo = 0;
        switch (liaison_topologique->get_dimension())
        {
        case 1:
        {
            refTopo = ((MCEdge*)liaison_topologique)->GetPolyCurve()->GetRefEdge(0);
            break;
        }
        case 2:
        {
            refTopo = ((MCFace*)liaison_topologique)->GetPolySurface()->GetRefFace(0);
            break;
        }
        }
        change_lien_topologie(refTopo);
    }
}
 
void MCSegment::enregistrer(std::ostream& o,double version)
{
    if (_saveFormat == 0) // Save as CAD4FE_MC_SEGMENT in magic file
    {
        if(liaison_topologique->get_dimension()==1) o << "%" << get_id() << "=CAD4FE_MCSEGMENT($"<< get_lien_topologie()->get_id() << ",$" << noeud1->get_id() << ",$" << noeud2->get_id() <<  ");" << std::endl;
    }
    else if (_saveFormat == 1) // Save as MG_SEGMENT in magic file
    {
        MG_SEGMENT::enregistrer(o,version);
    }
    else if (_saveFormat == 2)
    {                      
        MG_SEGMENT::enregistrer(o,version);
    }
}                       
int MCSegment::get_type_entite()
{
    return IDMCSEGMENT;
}