CAD4FE/src/CAD4FE_MCAA.cpp

//---------------------------------------------------------------------------

#include <sstream>
#include <fstream>
#pragma hdrstop

//---------------------------------------------------------------------------
// include MAGIC Headers
#include "gestionversion.h"
#include <mg_geometrie.h>
#include <mg_maillage.h>
#include <mg_maillage_outils.h>
#include <mg_gestionnaire.h>
#include <mailleur0d.h>
#include <mailleur1d.h>
#include <mailleur2d.h>
#include <fct_generateur_3d.h>
#include <fct_generateur_fichier.h>    
#include <fct_generateur_constante.h>
#include <fct_taille.h>
#include <ot_algorithme_geometrique.h>

#include "CAD4FE_MCBody.h"
#include "CAD4FE_Geometric_Tools.h"
#include "CAD4FE_PolyCurve.h"
#include "CAD4FE_PolySurface.h"
#include "CAD4FE_MCFace.h"
#include "CAD4FE_MCEdge.h"
#include "CAD4FE_MCVertex.h"
#include "CAD4FE_GlobalEdgeCriteria.h"
#include "CAD4FE_LoopCriteria.h"
#include "CAD4FE_LocalEdgeCriteria.h"
#include "CAD4FE_EdgeCollapseCriteria.h"
#include "CAD4FE_VertexCriteria.h"
#include "CAD4FE_mg_utils.h"
#include "CAD4FE_InventorText_MCAA.h"
#include "CAD4FE_InventorText_MG_MAILLAGE.h"
#include "CAD4FE_Intersection_Plane_MG_MAILLAGE.h"
#include "CAD4FE_Intersection_Plane_MG_SEGMENT.h"
#include "CAD4FE_ColorMap.h"
#include "CAD4FE_Criteria.h"        
#include "CAD4FE_ShortestPathByClosestPointOnEdge.h"  
#include "CAD4FE_MCNode.h"

#include "CAD4FE_MCAA.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)


using namespace CAD4FE;

MCAA::MCAA(MG_VOLUME * __refBody, MG_MAILLAGE * __tessellation, MG_GESTIONNAIRE * __gest, MG_GEOMETRIE * __geom)
: _refBody(__refBody), _refTess(__tessellation), _refGest(__gest), _geom(__geom), _sizeMap(NULL), _mcTessSegGrid(NULL), _feMesh(NULL)
{

        _mcBody = new MCBody(_geom, __refBody);

        std::ofstream file;
        char filename [5000]="";    
        strcat (filename, getenv("TEMP"));
        strcat (filename, "\\RefTessellation.iv");
        file.open(filename);
        file<<InventorText_MG_MAILLAGE(__tessellation).GetText();
        file<<"}";
        file.close();
        
        _InitializeMCBodyTessellation();

        // Setting default meshing parameters :
        _maxOverdensity = 2;
        _limitAngle = Criteria::PlanarCurve_BetaEpsilon(.22);
        SetMinMCEdgeCount(0);
        SetMinMCVertexCount(1);

        // set the time variable, representing the
        // simplification step number to zero
        time = 0;

        _nbA=0; // edge deletion
        _nbB=0; // vertex deletion
        _nbC=0; // edge collapse
        _nbD=0; // constricted section
        _nbE=0; // loops

        //_mcBody->G2ToNeatoFile(_mcBody->G10(), "c:\\gilles\\g10.neato", "c:\\gilles\\g10.neato");

                  /*
        // Try to do a 2D mesh
        //_DebugTraditionnalMesh();
        if (__sizeMap != 0)
            _sizeMap = __sizeMap;
        else
            UseDefaultMeshSize();
    */
}
//---------------------------------------------------------------------------
void
MCAA::Update()
{
        // Initialize the mesh size
        InitializeFEMesh();

        // Initialize the MC Tessellation Segments Grid
        _InitializeMCTessSegGrid();

        // Initialize Edge local properties
        EC_Init();

        // Initialize Vertex properties
        VC_Init();
}
//---------------------------------------------------------------------------

MCAA::~MCAA()
{
    delete _mcTessSegGrid;
    for (ECMapIterator it_ec = _ecMap.begin();
        it_ec != _ecMap.end();
        it_ec++)
    {
        GlobalEdgeCriteria * gec = it_ec->second;
        delete gec;
    }
    for (VCMapIterator it_vc = _vcMap.begin();
        it_vc != _vcMap.end();
        it_vc++)
    {
        VertexCriteria * vc = it_vc->second;
        delete vc;
    }
    delete _mcBody;
}
                   
//---------------------------------------------------------------------------
// Input Meshing criteria
// Limit deviation angle
double MCAA::GetLimitAngle()
{
    return _limitAngle;
}
//---------------------------------------------------------------------------
void MCAA::SetLimitAngle(double __limitAngle)
{
    _limitAngle = __limitAngle;
}
//---------------------------------------------------------------------------  
int MCAA::GetMinMCEdgeCount()
{
    return _minMCEdgeCount;
}
//---------------------------------------------------------------------------  
void MCAA::SetMinMCEdgeCount(int __minMCEdgeCount)
{
    _minMCEdgeCount =  __minMCEdgeCount;
}
//---------------------------------------------------------------------------  
int MCAA::GetMinMCVertexCount()
{
    return _minMCVertexCount;
}
//---------------------------------------------------------------------------  
void MCAA::SetMinMCVertexCount(int __minMCVertexCount)
{
    _minMCVertexCount=__minMCVertexCount;
}
// Input Meshing criteria
// Relative Sag
void MCAA::SetRelativeSag(double __relativeSag )
{
    _relativeSag = __relativeSag;

    // Browse the edge criterions
    // Set the relative sag on each edge suppression criterion
    ECMapCIterator it_ec;
    for (it_ec = _ecMap.begin();
            it_ec != _ecMap.end() ;
            it_ec++ )
    {
        // update local edge criteria
        // while keeping the normal offset polyline
        it_ec->second->Update(false);
    }

                                       
    // Set the relative sag on each vertex suppression criterion
    VCMapIterator itV;
    for (itV = _vcMap.begin(); itV != _vcMap.end(); itV++)
    {                                   
        // update local edge criteria
        // while keeping the normal offset polyline
        it_ec->second->Update(false);
    }
}
//---------------------------------------------------------------------------
double MCAA::GetRelativeSag()
{
    return _relativeSag;
}   
//---------------------------------------------------------------------------
// Overdensity factor
void MCAA::SetMaxOverdensity(double __maxOverdensity)
{
    _maxOverdensity = __maxOverdensity;    
    ECMapCIterator it_ec;
    for (it_ec = _ecMap.begin();
            it_ec != _ecMap.end() ;
            it_ec++ )
    {                             
        // update local edge criteria
        // while keeping the normal offset polyline
        it_ec->second->Update(false);
    }                          
    
    VCMapIterator itV;
    for (itV = _vcMap.begin(); itV != _vcMap.end(); itV++)
    {                            
        // update local edge criteria
        // while keeping the normal offset polyline
        it_ec->second->Update(false);
    }
}
double MCAA::GetMaxOverdensity()
{
    return _maxOverdensity;
}                        
//---------------------------------------------------------------------------
bool MCAA::CheckIfMeshIsReferencedInTopo(MG_ELEMENT_MAILLAGE * __elem)
{
    if ( __elem->get_dimension() == 1 && __elem->get_lien_topologie()->get_dimension() == 2 )
    {
        MG_FACE * face =  (MG_FACE*)__elem->get_lien_topologie();
        MG_SEGMENT* seg = (MG_SEGMENT*) __elem;
        return CheckIfMeshIsReferencedInTopo(seg, face);
    }

    if ( ! __elem->get_lien_topologie()->get_lien_maillage()->contient(__elem) )
    {
        printf("Mesh element ID \"%ld\" is not referenced in the topology ! ", __elem->get_id());
        return false;
    }
    else
        return true;
}              
//---------------------------------------------------------------------------
bool MCAA::CheckIfMeshIsReferencedInTopo(MG_SEGMENT * __seg, MG_FACE * __face)
{
    int i;
    for (i=0; i<__seg->get_lien_triangle()->get_nb(); i++)
    {
        MG_TRIANGLE * triang = __seg->get_lien_triangle()->get(i);

        if ( ! CheckIfMeshIsReferencedInTopo(triang) )
            return false;

    }                  

    return true;
}  
//---------------------------------------------------------------------------
bool MCAA::CheckIfTopoExists(MG_ELEMENT_TOPOLOGIQUE * topo)
{          
    if (topo)
    {
                switch (topo->get_dimension())
                {
                case 0:
                {
                        MCVertex * mctopo = _mcBody->GetMCVertex(topo->get_id());


                        if (mctopo != topo)
                        {
                                printf("MC Vertex %ld does not exist !\n", topo->get_id());
                                return false;
                        }
                        else
                                return true;


                }
                case 1:
                {
                        MCEdge * mctopo = _mcBody->GetMCEdge(topo->get_id());
                        MCEdge * topo2 = (MCEdge*) _geom->get_mg_areteid(topo->get_id());

                        if (mctopo != topo || topo2 != topo)
                        {
                                printf("MC Edge %ld does not exist !\n", topo->get_id());
                                return false;
                        }
                        else
                                return true;
                }
                case 2:
                {
                        MCFace * mctopo = _mcBody->GetMCFace(topo->get_id());

                        if (mctopo != topo)
                        {
                            printf("MC Face %ld does not exist !\n", topo->get_id());
                            return false;
                        }
                        else
                            return true;
                }
                default:
                printf("The topological link claims a dimension \"%d\" > 2 or < 0 !\n", topo->get_dimension());
                }
        }
        else
        {
                printf("Topo = NULL\n");
                return false;
        }
        return false;
}
//---------------------------------------------------------------------------

void MCAA::CheckMCMesh(MG_MAILLAGE * __mesh, std::set<MG_SEGMENT*> & __setSegmentBadTopology, std::set<MG_NOEUD*> & __setNodeBadTopology )
{
        LISTE_MG_TRIANGLE::iterator it_triang;
        for (MG_TRIANGLE* triang = __mesh->get_premier_triangle(it_triang);
                triang; triang =__mesh->get_suivant_triangle(it_triang))
        {
                MG_ELEMENT_TOPOLOGIQUE * topo = triang->get_lien_topologie();
                CheckIfTopoExists(topo);
                CheckIfMeshIsReferencedInTopo(triang);
        }

        LISTE_MG_SEGMENT::iterator it_seg;
        for (MG_SEGMENT* seg = __mesh->get_premier_segment(it_seg);
                seg; seg =__mesh->get_suivant_segment(it_seg))
        {
                MG_ELEMENT_TOPOLOGIQUE * topo = seg->get_lien_topologie();
                if (!CheckIfTopoExists(topo) || !CheckIfMeshIsReferencedInTopo(seg))
                    __setSegmentBadTopology.insert(seg);
        }

        LISTE_MG_NOEUD::iterator it_no;
        for (MG_NOEUD* no = __mesh->get_premier_noeud(it_no);
                no; no =__mesh->get_suivant_noeud(it_no))
        {
                MG_ELEMENT_TOPOLOGIQUE * topo = no->get_lien_topologie();
                if (!CheckIfTopoExists(topo))
                    __setNodeBadTopology.insert(no);
        }

        {MCBODY_FOR_EACH_MCFACE(_mcBody,face)
        {
                TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage = face->get_lien_maillage();

                TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
                MG_ELEMENT_MAILLAGE* element;
                int nb = lien_maillage->get_nb();
                int i=0;
                for (element = lien_maillage->get_premier(it); i++ < nb && element ; element = lien_maillage->get_suivant(it) )
                {
                    if (face != element->get_lien_topologie())
                        printf("face %ld : element id %ld NOT ok\n", face->get_id(), element->get_id());
                }
        }
        }
        {MCBODY_FOR_EACH_MCEDGE(_mcBody,edge)
        {
                TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage = edge->get_lien_maillage();

                TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
                MG_ELEMENT_MAILLAGE* element;
                int nb = lien_maillage->get_nb();
                int i=0;
                for (element = lien_maillage->get_premier(it); i++ < nb && element ; element = lien_maillage->get_suivant(it) )
                {                        
                    if (edge != element->get_lien_topologie())
                    {
                        printf("edge %ld : element id %ld NOT ok\n", edge->get_id(), element->get_id());
                    }
                }


                int nbOfNodesDim0 = 0;
                i=0;
                for (element = lien_maillage->get_premier(it); i++ < nb && element ; element = lien_maillage->get_suivant(it) )
                {
                if ( ! __mesh->contient(element) ) continue;
                MG_SEGMENT * seg = (MG_SEGMENT *) element;
                MG_NOEUD * V[2];
                V[0] = seg->get_noeud1();
                V[1] = seg->get_noeud2();
                for (int i=0; i<2; i++)
                {
                        if ( V[i]->get_lien_topologie() && V[i]->get_lien_topologie()->get_dimension() == 0)
                                nbOfNodesDim0++;
                }
                }
                if (nbOfNodesDim0 != 2)
                    printf("Edge %ld Problem : there has %d nodes having a vertex attached!\n", edge->get_id(), nbOfNodesDim0);

        }
        }
        {MCBODY_FOR_EACH_MCVERTEX(_mcBody,vertex)
        {
                TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage = vertex->get_lien_maillage();
                                         
                TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
                MG_ELEMENT_MAILLAGE* element;
                int nb = lien_maillage->get_nb();
                int i=0;
                for (element = lien_maillage->get_premier(it); i++ < nb && element ; element = lien_maillage->get_suivant(it) )
                {        
                    if (vertex != element->get_lien_topologie())
                        printf("vertex %ld : element id %ld NOT ok\n", vertex->get_id(), element->get_id());
                }
        }
        } 
}                   
                                                                               
//---------------------------------------------------------------------------
void MCAA::MergeVertices(MCVertex * __deleteVertex, MCVertex * __targetVertex)
{
    MCVertex * V1 =  __deleteVertex;      
    MCVertex * V2 =  __targetVertex;
    std::set<MCEdge*>edges;
    _mcBody->MergeVertices(V1,V2,edges);

    std::set<MCEdge*>::iterator itEdge;
    for (itEdge = edges.begin(); itEdge != edges.end(); itEdge++)
    {
        EC_Delete(*itEdge);
        DeleteFEMesh(*itEdge);
    }
    TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = V1->get_lien_maillage();
    TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
    MG_ELEMENT_MAILLAGE* element2;
    for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
        element2->change_lien_topologie(NULL);

    int N=V1->get_nb_mg_cosommet();
    for (int i=0; i<N; i++)
        _geom->supprimer_mg_cosommet(V1->get_mg_cosommet(0));
    _mcBody->DeleteMCVertex(V1);
    VC_Delete(V1);

    for (itEdge = edges.begin(); itEdge != edges.end(); itEdge++)
        InitializeFEMesh(*itEdge);

    for (itEdge = edges.begin(); itEdge != edges.end(); itEdge++)
        EC_Init(*itEdge);
}
//---------------------------------------------------------------------------
void
MCAA::_InitializeMCBodyTessellation()
{       
        // Duplicate the reference geometry's mesh !
        _mcTess = _refTess->dupliquer ( _refGest );

        long id_ref_to_mc = _mcTess->get_id() - _refTess->get_id();

        // update the link of tessellation entities of the MC Body Tessellation
        {
        MCBODY_FOR_EACH_MCFACE(_mcBody, mcFace)
        {
            const std::set<MG_FACE *> & lst_ref_faces = mcFace->GetPolySurface()->GetRefFaces();
            for (std::set<MG_FACE *>::const_iterator it_ref_face = lst_ref_faces.begin();
                   it_ref_face != lst_ref_faces.end();
                   it_ref_face ++ )
            {
                MG_FACE * refFace = *it_ref_face;

                TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = refFace->get_lien_maillage();
                TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
                MG_ELEMENT_MAILLAGE* element2;
                for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
                {
                    if (_refTess->contient(element2))
                    {
                        MG_TRIANGLE * ref_triang = (MG_TRIANGLE *) element2;
                        MG_TRIANGLE * mc_triang = _mcTess->get_mg_triangleid ( ref_triang->get_id() + id_ref_to_mc );

                        mc_triang->change_lien_topologie ( mcFace );
                            
                        MG_SEGMENT * segs[3];
                        MG_NOEUD * nos[3];
                        MG_TRIANGLE_GET_SEGMENTS(mc_triang,segs);
                        MG_TRIANGLE_GET_NOEUDS(mc_triang,nos);
                        for (int i=0;i<3;i++)
                        {
                            if (!segs[i]->get_lien_topologie() || segs[i]->get_lien_topologie()->get_dimension() == 2)
                                segs[i]->change_lien_topologie2(mcFace);//change_lien_topologie_mct(segs[i],mapRefToMCFace,mapRefToMCEdge,mapRefToMCVertex);
                            if (! nos[i]->get_lien_topologie() || nos[i]->get_lien_topologie()->get_dimension() == 2)
                                nos[i]->change_lien_topologie2(mcFace);//change_lien_topologie_mct(nos[i],mapRefToMCFace,mapRefToMCEdge,mapRefToMCVertex);
                        }
                    }

                }
            }
        }}

        {
        MCBODY_FOR_EACH_MCEDGE(_mcBody, edge)
        {
                const std::vector<MG_ARETE *> & lst_ref_edges = edge->GetPolyCurve()->GetRefEdges();
                for (std::vector<MG_ARETE *>::const_iterator it_ref_edge = lst_ref_edges.begin();
                       it_ref_edge != lst_ref_edges.end();
                       it_ref_edge ++ )
                {
                        MG_ARETE * ref_edge = *it_ref_edge;

                        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = ref_edge->get_lien_maillage();
                        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
                        MG_ELEMENT_MAILLAGE* element2;
                        int counter = 0;
                        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
                        {
                            if (_refTess->contient(element2))
                            {
                                MG_SEGMENT * ref_seg = (MG_SEGMENT *) element2;

                                if ( ! ref_seg ) continue;

                                MG_SEGMENT * mc_seg = _mcTess->get_mg_segmentid ( ref_seg->get_id() + id_ref_to_mc );
                                mc_seg->change_lien_topologie( edge );

                                MG_NOEUD * nos[2];
                                nos[0] = mc_seg->get_noeud1();
                                nos[1] = mc_seg->get_noeud2();
                                for (int i=0;i<2;i++)
                                    if (! nos[i]->get_lien_topologie() || nos[i]->get_lien_topologie()->get_dimension() > 0)
                                        nos[i]->change_lien_topologie2(edge); 

                                counter++;
                            }
                        }
                        if (counter == 0)
                            printf("warning : no segments are connected to ref edge id=(%ld,%s) !\n",edge->get_id(),edge->get_idoriginal().c_str());
                }
        }
        }

        {MCBODY_FOR_EACH_MCVERTEX(_mcBody, mc_vertex)
        {
                MG_SOMMET * ref_vertex = mc_vertex->GetRefVertex();                                

                TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = ref_vertex->get_lien_maillage();
                TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
                MG_ELEMENT_MAILLAGE* element2;
                for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
                {
                    if (_refTess->contient(element2))
                    {
                        MG_NOEUD * ref_node = (MG_NOEUD*) element2;
                        MG_NOEUD * mc_node = _mcTess->get_mg_noeudid ( ref_node->get_id() + id_ref_to_mc );

                        mc_node->change_lien_topologie(mc_vertex);
                        break;
                    }
                }
        }}

        MG_UTILS::MG_MAILLAGE_GetBoundingBox(_mcTess, _bbox);
}

//---------------------------------------------------------------------------
void
MCAA::_InitializeMCTessSegGrid()
{
        if (_mcTessSegGrid)
                delete _mcTessSegGrid;

        _mcTessSegGrid = new TPL_GRILLE < MG_SEGMENT * >;

        int nb_steps[3];
        for (unsigned i = 0; i < 3; i++)
        {
                double cellSize = GetSize(_bbox);

                if ( cellSize != 0 )
                    nb_steps[i] = std::min(floor((_bbox[i+3]-_bbox[i])/cellSize),40.0);
                else
                    nb_steps[i] = 40;
                
                if (nb_steps[i] == 0)
                    nb_steps[i] ++;
        }
        MG_UTILS utils;
        utils.TPL_GRILLE_Construct<MG_SEGMENT*>(_bbox, nb_steps, *_mcTessSegGrid); 

        {
        MCBODY_FOR_EACH_MCEDGE(_mcBody, edge)
        {

                        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = edge->get_lien_maillage();
                        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
                        MG_ELEMENT_MAILLAGE* element2;
                        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
                        {
                                MG_SEGMENT * mc_seg = (MG_SEGMENT *)element2;

                                if (_mcTess->contient(mc_seg) )
                                        _mcTessSegGrid->inserer(mc_seg);
                        }
        }
        }
}
//---------------------------------------------------------------------------
void
MCAA::_UpdateMCTessSegGrid(MCEdge * __mcEdge, MG_SEGMENT * __originalSegment, MG_SEGMENT * __splitSegments[2])
{
    int i;

    if (__originalSegment)
        _mcTessSegGrid->supprimer( __originalSegment );

    for (i=0; i<2; i++)
    if ( __splitSegments[i] )
        _mcTessSegGrid->inserer( __splitSegments[i] );

}
TPL_GRILLE <MG_SEGMENT *> * MCAA::GetMCTessSegGrid()
{
    return _mcTessSegGrid;
}
//---------------------------------------------------------------------------
MG_VOLUME * MCAA::GetRefBody () const
{
        return this->_refBody;
}

//---------------------------------------------------------------------------
MCBody * MCAA::GetMCBody () const
{             
        return this->_mcBody;
}


//---------------------------------------------------------------------------
MG_MAILLAGE * MCAA::GetRefTess () const
{
        return this->_refTess;
}


//---------------------------------------------------------------------------
MG_MAILLAGE * MCAA::GetMCTess () const
{
        return this->_mcTess;
}


//---------------------------------------------------------------------------
MG_GESTIONNAIRE * MCAA::GetGest () const
{
        return this->_refGest;
}

//---------------------------------------------------------------------------
void MCAA::InitializeFEMesh ()
{
        if ( ! _sizeMap )
        {
                printf("InitializeFEMesh() : Can't mesh this MC Body : the size map is NULL \n");
                return;
        }

        if (_feMesh)
            _refGest->supprimer_mg_maillageid(_feMesh->get_id() );
        _feMesh = new MG_MAILLAGE(_geom);
        _refGest->ajouter_mg_maillage(_feMesh);

        printf("MAILLAGE 0D\n");
        {
                MCBODY_FOR_EACH_MCVERTEX(_mcBody, vertex)
                {
                        printf("    MC vertex %ld\n", vertex->get_id() );
                        InitializeFEMesh(vertex);
                }
        }


        printf("MAILLAGE 1D\n");
        {
                MCBODY_FOR_EACH_MCEDGE(_mcBody, edge)
                {
                        printf("    MC edge %ld\n",edge->get_id());
                        InitializeFEMesh(edge);
                }
        }       
}

//---------------------------------------------------------------------------
#undef debug_InitializeFEMesh_Edge
void MCAA::InitializeFEMesh(MCEdge * __mcEdge)
{
        MAILLEUR1D m1d(_feMesh,_geom,_sizeMap,__mcEdge);
        m1d.maille();

#ifdef debug_InitializeFEMesh_Edge
        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = __mcEdge->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
        MG_ELEMENT_MAILLAGE* element2;
        int nb=0;
        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
        {
            if (_feMesh->contient(element2))
                nb++;
        }
        printf("There are %ld FE mesh elements on the MC Edge %ld\n",  nb, __mcEdge->get_id());
#endif
}
//---------------------------------------------------------------------------
#undef debug_DeleteFEMesh_
void MCAA::DeleteFEMesh(MCEdge * __mcEdge)
{
        std::vector <MG_SEGMENT *> vecSegs;
        /*LISTE_MG_SEGMENT::iterator it_seg;
        for (MG_SEGMENT* seg = _feMesh->get_premier_segment(it_seg);
                seg; seg =_feMesh->get_suivant_segment(it_seg))
        {
                MG_ELEMENT_TOPOLOGIQUE * topo = seg->get_lien_topologie();
                if (__mcEdge == topo)
                {
                    vecSegs.push_back(seg);
                }
        } */

        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = __mcEdge->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
        MG_ELEMENT_MAILLAGE* element2;
        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
        {
            if (_feMesh->contient(element2))
                    vecSegs.push_back((MG_SEGMENT*)element2);
        }

        for (std::vector<MG_SEGMENT *>::iterator it = vecSegs.begin();
            it != vecSegs.end();
            it++)
        {
            MG_SEGMENT * seg = *it;

            _feMesh->supprimer_mg_segmentid(seg->get_id());
        }

#ifdef debug_DeleteFEMesh_
        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = __mcEdge->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
        MG_ELEMENT_MAILLAGE* element2;
        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
        {
            if (_feMesh->contient(element2))
                _feMesh->supprimer_mg_segmentid(element2->get_id());
            else if ( ! _mcTess->contient(element2))
                printf("A segment of edge %ld is neither attached to MC Tess nor FE Mesh !\n",element2->get_id(), __mcEdge->get_id());
        }
#endif        
}                                    
//---------------------------------------------------------------------------
int MCAA::GetFEMeshSegmentCount(MCEdge * __mcEdge)
{
    int result = 0;
    
    TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = __mcEdge->get_lien_maillage();
    TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
    MG_ELEMENT_MAILLAGE* element2;
    for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
    {
        if (_feMesh->contient(element2))
            result ++;
    }

    return result;
}
//---------------------------------------------------------------------------
void MCAA::DeleteMesh(MCEdge * __mcEdge)
{           
        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = __mcEdge->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
        MG_ELEMENT_MAILLAGE* element2;
        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
        {
            if (_feMesh->contient(element2))
                _feMesh->supprimer_mg_segmentid(element2->get_id());
            else if (_mcTess->contient(element2))
                _mcTess->supprimer_mg_segmentid(element2->get_id());
            else
            {
                delete element2;
                printf("A segment of edge %ld is neither attached to MC Tess nor FE Mesh !\n",element2->get_id(), __mcEdge->get_id());
            }
        }
}
//---------------------------------------------------------------------------
void MCAA::InitializeFEMesh(MCVertex * __mcVertex, MG_NOEUD ** __n)
{
        MG_NOEUD* n = NULL;

            double coo[3];
            __mcVertex->get_point()->evaluer(coo);
            n = new MG_NOEUD(__mcVertex,coo[0],coo[1],coo[2],MAILLEUR_AUTO);
            _feMesh->ajouter_mg_noeud(n);

        if (__n) *__n = n;
}

//---------------------------------------------------------------------------
void MCAA::DeleteFEMesh(MCVertex * __mcVertex)
{
        TPL_SET<MG_ELEMENT_MAILLAGE*> lien_maillage2 = *__mcVertex->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
        MG_ELEMENT_MAILLAGE* element2;
        for (element2 = lien_maillage2.get_premier(it2); element2; element2 = lien_maillage2.get_suivant(it2) )
        {
            if (_feMesh->contient(element2))
                _feMesh->supprimer_mg_noeudid(element2->get_id());
        }
}
//---------------------------------------------------------------------------
void
MCAA::_DebugTraditionnalMesh()
{

        _feMesh = new MG_MAILLAGE(_geom);
        _refGest->ajouter_mg_maillage(_feMesh);

        printf("MAILLAGE 0D\n");
        {
                for (unsigned it_vertex = 0; it_vertex != _geom->get_nb_mg_sommet(); it_vertex++)
                {
                        MG_SOMMET * vertex = _geom->get_mg_sommet(it_vertex);
                        printf("    MC vertex %ld\n", vertex->get_id() );
                        double coo[3];
                        vertex->get_point()->evaluer(coo);
                        MG_NOEUD* n = new MG_NOEUD(vertex,coo[0],coo[1],coo[2],MAILLEUR_AUTO);
                        _feMesh->ajouter_mg_noeud(n);
                        vertex->get_lien_maillage()->ajouter(n);
                }
        }


        printf("MAILLAGE 1D\n");
        {
                for (unsigned it_edge = 0; it_edge != _geom->get_nb_mg_arete(); it_edge++)
                {
                        MG_ARETE * edge = _geom->get_mg_arete(it_edge);
                        printf("    MC edge %ld\n",edge->get_id());
                        MAILLEUR1D m1d(_feMesh,_geom,_sizeMap,edge);
                        m1d.maille();
                }
        }

                       
        printf("MAILLAGE 2D\n");
        {
                for (unsigned it_face = 0; it_face != _geom->get_nb_mg_face(); it_face++)
                {
                        MG_FACE * face = _geom->get_mg_face(it_face);
                        printf("    MC face %ld\n",face->get_id());
                        MAILLEUR2D m2d(_feMesh,_geom,_sizeMap,face);
                         m2d.change_priorite_metrique(.65);
                        m2d.maille();
                }
        }

}
//---------------------------------------------------------------------------
MG_MAILLAGE *
MCAA::GetFEMesh () const
{
        return _feMesh;
}
//---------------------------------------------------------------------------
void
MCAA::UseDefaultMeshSize()
{
        // Setting the default Size Map
        double s = CalculateDefaultMeshSize();
        printf (" Mesh size set to default value = %f\n", s);
        SetConstantMeshSize ( s );
        Update();
}
//---------------------------------------------------------------------------
void MCAA::ReadSizeMapFile(const char * __filename)
{
    FCT_GENERATEUR_3D<4> * sizeMap = new FCT_GENERATEUR_3D<4>;
    sizeMap->lire((char*)__filename);
    if (_sizeMap)
        delete _sizeMap;
    _sizeMap = sizeMap;

    Update();
}                
//---------------------------------------------------------------------------
void MCAA::SetSizeMap(FCT_TAILLE * __sizeMap)
{
    if (_sizeMap)
        delete _sizeMap;
    _sizeMap = __sizeMap;    
    Update();
}
//---------------------------------------------------------------------------
void MCAA::ReadSizeMapRefinementFile(const char * __filename)
{
    // Setting the default Size Map
    double s = CalculateDefaultMeshSize();
    FCT_GENERATEUR_FICHIER * sizeMap = new FCT_GENERATEUR_FICHIER((char*)__filename, s);
    _sizeMap = sizeMap;
    sizeMap->construit();
    char outfilename[10000]="";
    sprintf(outfilename,"%s.ctt",__filename);
    sizeMap->enregistrer(outfilename);
    sprintf(outfilename,"%s.ctt","c:\\gilles\\mcbody.magic");  
    sizeMap->enregistrer(outfilename);      

    Update();
}                                  
//---------------------------------------------------------------------------
void
MCAA::SetConstantMeshSize(double __meshSize)
{
        if (_sizeMap)
                delete _sizeMap;

    char outputfilename[]="c:\\gilles\\mcbody.magic.ctt";
    _sizeMap = new FCT_GENERATEUR_CONSTANTE ( *_refGest, __meshSize );
    ((FCT_GENERATEUR_CONSTANTE*)_sizeMap)->construit(_bbox[0],_bbox[1],_bbox[2],_bbox[0+3],_bbox[1+3],_bbox[2+3]);
    ((FCT_GENERATEUR_CONSTANTE*)_sizeMap)->enregistrer(outputfilename);
}

//---------------------------------------------------------------------------
double
MCAA::CalculateDefaultMeshSize ()
{
        // Default Size = .3 * Volume / Surface
        double V, S;
        V = MG_UTILS::MG_MAILLAGE_GetVolume(_mcTess);
        S = MG_UTILS::MG_MAILLAGE_GetSurface(_mcTess);
        printf("Volume = %f\n", V );
        printf("Surface = %f\n", S );
        return .3 * V / S;
}

//---------------------------------------------------------------------------
void MCAA::GetBoundingBox(double __bbox[6]) const
{
        for (unsigned i = 0; i<6; i++)
                __bbox[i] = _bbox[i];
}

//---------------------------------------------------------------------------
double
MCAA::GetSize(double __xyz[3])
{
    if (_sizeMap == NULL)
        return 0;

    double density_tensor[9];
    _sizeMap->evaluer(__xyz, density_tensor);
    return pow(density_tensor[0], -0.5);
}
//---------------------------------------------------------------------------
double
MCAA::GetMinSize(double __xyz[3])
{
    return GetSize(__xyz) / GetMaxOverdensity();
}
//---------------------------------------------------------------------------
MG_SEGMENT * MCAA::_FindClosestSeg(TPL_MAP_ENTITE<MG_SEGMENT *> & __lstSegs, MCEdge * __mcEdge, double __P[3], double __N[3])
{
        // minimum distance found between the tess edge and P
        double min_distanceFrom_P = 9E99;
        // vertices of the edge which are closest to edgePnt
        MG_SEGMENT * closest_seg = NULL;
        
 TPL_MAP_ENTITE<MG_SEGMENT *>::ITERATEUR itSeg; 
 for (MG_SEGMENT * seg = __lstSegs.get_premier(itSeg);
      seg;
      seg = __lstSegs.get_suivant(itSeg))
        {
                Intersection_Plane_MG_SEGMENT intersection ( __P, __N, seg, 1E-12);

                if (intersection.Test () )
                {
                        intersection.Find();

                        if (__mcEdge != seg->get_lien_topologie() || intersection.GetSegmentT() < 0 || intersection.GetSegmentT() > 1)
                                continue;        

                        double intersectionPt [3];
                        intersection.GetPoint( intersectionPt );
                        double distance = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3 ( intersectionPt, __P );

                        if (distance < min_distanceFrom_P)
                        {
                                min_distanceFrom_P = distance;
                                closest_seg = seg;
                        }
                }
        }

        return  closest_seg;
}
                                                                             
//---------------------------------------------------------------------------
MG_SEGMENT * MCAA::FindClosestSeg(MCEdge * __mcEdge, double __P[3], double __N[3])
{
    TPL_MAP_ENTITE<MG_SEGMENT *>  lstSegs;
    _mcTessSegGrid->rechercher( __P[0], __P[1], __P[2], GetSize(__P), lstSegs);
    return _FindClosestSeg(lstSegs, __mcEdge, __P, __N);
}
//---------------------------------------------------------------------------
void MCAA::EC_Init(MCEdge * __mcEdge)
{
        // arguments : _mcBody, __mcEdge, _feMesh, __mcTess, _mcTessSegGrid,
        GlobalEdgeCriteria * gec = new GlobalEdgeCriteria(__mcEdge, this);
        gec->Init();
        _ecMap.insert( std::make_pair(__mcEdge,gec) );
}                                                           
//---------------------------------------------------------------------------
void MCAA::EC_ChangeTopo(MCEdge * __old, MCEdge * __new)
{
    GlobalEdgeCriteria * old_gec, * new_gec;
    ECMapIterator it_ec = _ecMap.find(__old);

    if ( it_ec != _ecMap.end() )
    {
        old_gec = it_ec->second;
    }
    else
        return;
                              
    it_ec = _ecMap.find(__new);
    if ( it_ec != _ecMap.end() )
    {
        new_gec = it_ec->second;
    }
    else
    {
        new_gec = new GlobalEdgeCriteria (__new, this);
        _ecMap.insert( std::make_pair(__new,new_gec) );
    }

    new_gec->ChangeTopo(old_gec);
}
//---------------------------------------------------------------------------
void MCAA::EC_Delete(MCEdge * __mcEdge)
{
    ECMapIterator it_ec = _ecMap.find(__mcEdge);

    if ( it_ec != _ecMap.end() )
    {
        GlobalEdgeCriteria * gec = it_ec->second;
        delete gec;
        _ecMap.erase (it_ec);
    }
    else
    {       
        for (it_ec = _ecMap.begin();
            it_ec != _ecMap.end();
            it_ec ++)
        {
                GlobalEdgeCriteria * gec = it_ec->second;
                if (gec->GetEdge() == __mcEdge)
                {
                    delete gec;
                    _ecMap.erase (it_ec);
                    it_ec = _ecMap.begin();
                }
        }

    }
}

//---------------------------------------------------------------------------
void MCAA::EC_Update(MCEdge * __mcEdge)
{
    ECMapIterator it_ec = _ecMap.find(__mcEdge);
    if (it_ec != _ecMap.end())
    {
        GlobalEdgeCriteria * ec = it_ec->second;
        // update local edge criteria
        // while keeping the normal offset polyline
        ec->Update(true);
        _ecMap.erase(it_ec);
        _ecMap.insert( std::make_pair(ec->GetEdge(),ec) );
    }
    else
        printf("_ecMap.end() !\n");
}                 
//---------------------------------------------------------------------------
void MCAA::EC_Init()
{                 
        printf("Initializing edge local properties \n");
        {
                MCBODY_FOR_EACH_MCEDGE(_mcBody, edge)
                {
                        printf("    Edge local properties of MC edge %ld\n",edge->get_id());
                        EC_Delete(edge);
                        EC_Init(edge);
                }
        }
}

//---------------------------------------------------------------------------
std::vector <LocalEdgeCriteria *>
MCAA::EC_Get(MCEdge * __mcEdge)
{
    std::vector <LocalEdgeCriteria *> res;
    ECMapIterator it_ec = _ecMap.find(__mcEdge);
    if (it_ec != _ecMap.end())
    {
        GlobalEdgeCriteria * ec = it_ec->second;
        res = ec->GetLocalEdgeCriteria();
    }
    else
        printf("_ecMap.end() !\n");

    return res;
}                                
//---------------------------------------------------------------------------
GlobalEdgeCriteria * MCAA::GetGlobalEdgeCriteria(MCEdge * __mcEdge)
{                
    ECMapIterator it_ec = _ecMap.find(__mcEdge);
    if (it_ec != _ecMap.end())
    {
        return it_ec->second;
    }
    else
    {
        return NULL;
    }
}
//---------------------------------------------------------------------------
GlobalEdgeCriteria * MCAA::GetHighestEdgeDeletionScore()
{
    GlobalEdgeCriteria * BestDeletionEC = NULL;
    double bestscore = 0;

        ECMapIterator it_ec;
        for (it_ec = _ecMap.begin();
            it_ec != _ecMap.end();
            it_ec++)
        {
            MCEdge * mcEdge = it_ec->first;
            GlobalEdgeCriteria * ec = it_ec->second;
            double score = ec->DeletionScore();

            if (mcEdge != ec->GetEdge())
                continue;

            if ( !BestDeletionEC  ||
                score>bestscore )
            {
                BestDeletionEC = ec;
                bestscore = score;
            }
        }

        if ( bestscore < 1E-3 )
        {
            bestscore = 0;
            BestDeletionEC = NULL;
        }

        return BestDeletionEC;
}
        
//---------------------------------------------------------------------------
VertexCriteria * MCAA::GetHighestVCScore()
{
        VertexCriteria * vc = NULL;
        double highestScore = -1E308;

        for (VCMapCIterator itVC = _vcMap.begin();
                itVC != _vcMap.end();
                itVC++)
        {
                if (vc == NULL)
                {
                        vc = itVC->second;
                }
                else if (itVC->second->GetScore() > highestScore)
                {
                    MCVertex * mcVertex = itVC->second->GetVertex();
                    
                    if (_mapSuppressRefVertexCount.find(mcVertex->GetRefVertex()) ==  _mapSuppressRefVertexCount.end() ||
                        _mapSuppressRefVertexCount[mcVertex->GetRefVertex()] < 4 )
                    {
                        vc = itVC->second;   
                        highestScore = vc->GetScore();
                    }
                }
        }
        
        return vc;
}
       
//---------------------------------------------------------------------------
void
MCAA::VC_Init()
{    
        {
                MCBODY_FOR_EACH_MCVERTEX(_mcBody, vertex)
                {
                        printf("    Vertex local properties of MC vertex %ld\n",vertex->get_id());
                        VC_Delete(vertex);
                        VC_Init(vertex);
                }
        }
}
                                  
//---------------------------------------------------------------------------
void MCAA::VC_Init(MCVertex * __mcVertex)
{
        double xyz[3];
        __mcVertex->get_point()->evaluer(xyz);
        VertexCriteria * vc = new VertexCriteria(__mcVertex, this);
        _vcMap [ __mcVertex ] = vc;
}
                                              
//---------------------------------------------------------------------------
void MCAA::VC_Delete(MCVertex * __mcVertex)
{
        VCMapIterator itV = _vcMap.find(__mcVertex);
        if ( itV != _vcMap.end())
        {
            delete itV->second;
            _vcMap.erase(itV);
        }                        
}                                             
//---------------------------------------------------------------------------

VertexCriteria * MCAA::VC_Get(MCVertex * __mcVertex)
{
        VCMapCIterator it = _vcMap.find(__mcVertex);
        if (it != _vcMap.end() )
        {
                return it->second;
        }
        else
        {
                printf("Warning : did not find the vertex properties of vertex %p id = %ld\n", __mcVertex, __mcVertex->get_id());
                return NULL;
        }
}

//---------------------------------------------------------------------------
bool MCAA::SuppressNextVertex(MCTChanges * __mctChanges)
{
    VertexCriteria * vc = GetHighestVCScore();

    if (vc == NULL)
        return false;
    
    bool bSimplifDone = (vc->GetScore() > 1E-3);
    if (bSimplifDone)
    {
        SuppressMCVertex(vc->GetVertex(), __mctChanges);
        vc = GetHighestVCScore();
    }

    return bSimplifDone;
} 
//---------------------------------------------------------------------------
bool MCAA::Node_Segment_Same_Face(MG_NOEUD * __n, MG_SEGMENT * __seg, std::set <MCFace *> & __commonFaces)
{
    Graph::Arc * segEdge_G21Arc = _mcBody->G21()->GetArc(__seg->get_lien_topologie()->get_id());
    
    for (int i=0; i<__n->get_lien_segment()->get_nb(); i++)
    {
        MG_SEGMENT * seg = __n->get_lien_segment()->get(i);
        MG_ARETE * edge = (MG_ARETE *) seg->get_lien_topologie();

            Graph::Arc *G21Arc = _mcBody->G21()->GetArc(edge->get_id());
        for (Graph::Arc::MultimapNodesById::const_iterator itNode = G21Arc->Nodes().begin(); itNode != G21Arc->Nodes().end(); itNode++)
            if ( segEdge_G21Arc->Nodes().find(itNode->first) != segEdge_G21Arc->Nodes().end() )
                __commonFaces.insert( _mcBody->GetMCFace(itNode->second) );
    }

    return (__commonFaces.size() != 0);
}

Graph::Graph * MCAA::FaceBoundaryMesh(MG_MAILLAGE * __mesh, MCFace * __face, std::set <MG_NOEUD * > & nodes, std::set <MG_SEGMENT *> & segs)
{
    Graph::Graph * faceBoundaryMeshGraph = new Graph::Graph();
    nodes.clear();
    segs.clear();

    Graph::Node * F_G21Node = _mcBody->G21()->GetNode(__face->get_id());
    Graph::Node::MultimapArcsById & F_G21NodeArcs = F_G21Node->IncidentArcs();

    LISTE_MG_SEGMENT::iterator it_seg;
    for (MG_SEGMENT* seg = __mesh->get_premier_segment(it_seg);
            seg; seg =__mesh->get_suivant_segment(it_seg))
    {
        int seg_Topo_Id = seg->get_lien_topologie()->get_id();
        if (F_G21NodeArcs.find(seg_Topo_Id) != F_G21NodeArcs.end())
        {
            nodes.insert(seg->get_noeud1());
            nodes.insert(seg->get_noeud2());
            segs.insert(seg);
        }
    }

    for ( std::set<MG_SEGMENT*>::iterator itSeg = segs.begin();
        itSeg != segs.end();
        itSeg ++)
    {
        MG_SEGMENT * seg = *itSeg;
        Graph::Node * graphNode = faceBoundaryMeshGraph->AddNode(seg->get_id());
        graphNode->SetUserData(seg);
        graphNode->SetUserData(2,__mesh);
    }

    for (std::set<MG_NOEUD*>::iterator itNode = nodes.begin();
        itNode != nodes.end();
        itNode ++)
    {
        MG_NOEUD * n = *itNode;
        std::multimap< int, Graph::Node * > arcNodes;

        TPL_LISTE_ENTITE< class MG_SEGMENT * > * adjSegList = n->get_lien_segment();
        int nbAdjSeg = adjSegList->get_nb();
        for (int i=0; i<nbAdjSeg; i++)
        {
            MG_SEGMENT * seg = adjSegList->get(i);
            Graph::Node * graphNode = faceBoundaryMeshGraph->GetNode(seg->get_id());
            if (graphNode)
                arcNodes.insert(std::make_pair(graphNode->Id(), graphNode));
        }
        Graph::Arc * graphArc = faceBoundaryMeshGraph->AddArc(n->get_id(), arcNodes);
        graphArc->SetUserData(n);
        graphArc->SetUserData(2,__mesh);
    }

    return faceBoundaryMeshGraph;
}

//---------------------------------------------------------------------------
int MCAA::SimplifyConstrictedSections()
{
    int result = 0;
    {MCBODY_FOR_EACH_MCFACE(_mcBody, face)
    {
        bool stop = false;
        while (stop == false)
        {
            bool done = ProcessNextConstrictedSectionInFace(face);
            if (done)
            {
                result++;
            }
            else
            {
                stop = true;
            }
        }
    }
    }
    return result;
}

bool MCAA::ProcessNextConstrictedSectionInFace(MCFace * __mcFace)
{
    std::set <MG_NOEUD * > nodes;
    std::set <MG_SEGMENT *> segs;

        MCFace * face = __mcFace;
        Graph::Graph * G = FaceBoundaryMesh(_feMesh, face, nodes, segs);

        ConstrictedSection cs,cs2,cs3;
        MG_NOEUD *no2,*no3;
        for (std::set <MG_NOEUD * >::iterator itNo = nodes.begin();
                itNo != nodes.end();
                itNo++)
        {
            MG_NOEUD * no = *itNo;

            if (NodeConstrictedSection(face, no, cs, G))
            {
                MCVertex *V1, *V2;

                {
                    MCEdge * edge = 0;
                    if (no->get_lien_topologie()->get_dimension() > 0)
                    {
                        edge = (MCEdge*)no->get_lien_topologie();
                        if (((MCVertex*)edge->get_cosommet1()->get_sommet())->GetMergedRefVertices().size() > 0)
                           continue;
                        if (((MCVertex*)edge->get_cosommet2()->get_sommet())->GetMergedRefVertices().size() > 0)
                           continue;
                    }
                    if (cs.Seg1->get_noeud1()->get_lien_topologie()->get_dimension() > 0 && cs.Seg1->get_noeud2()->get_lien_topologie()->get_dimension() > 0)
                    {                           
                        edge = (MCEdge*)cs.Seg1->get_lien_topologie();
                        if (((MCVertex*)edge->get_cosommet1()->get_sommet())->GetMergedRefVertices().size() > 0)
                            continue;

                        if (((MCVertex*)edge->get_cosommet2()->get_sommet())->GetMergedRefVertices().size() > 0)
                            continue;
                    }
                }

                // if node is not on a vertex, then insert a vertex on node
                if (no->get_lien_topologie()->get_dimension() > 0)
                {
                    MCEdge * edge = (MCEdge*)no->get_lien_topologie();

                    MCTChanges mctChanges;
                    RTChanges rtChanges;
                    SplitEdge(edge, no->get_coord(), &mctChanges, &rtChanges);
                    V1 = (mctChanges.V.begin())->second;
                }
                else // V1 = this vertex
                {
                    V1 = (MCVertex *)no->get_lien_topologie();
                }
                // if the extremities of the segment are not on a vertex, then insert a vertex
                // on the closest point
                if (cs.Seg1->get_noeud1()->get_lien_topologie()->get_dimension() > 0 && cs.Seg1->get_noeud2()->get_lien_topologie()->get_dimension() > 0)
                {                          
                    MCEdge * edge = (MCEdge*)cs.Seg1->get_lien_topologie();
                    MCTChanges mctChanges;
                    RTChanges rtChanges;
                    SplitEdge(edge, cs.Seg1P, &mctChanges, &rtChanges);
                    V2 = (mctChanges.V.begin())->second;
                }
                else // V2 = this vertex
                {
                    if (cs.Seg1->get_noeud1()->get_lien_topologie()->get_dimension() == 0 && cs.Seg1->get_noeud2()->get_lien_topologie()->get_dimension() == 0 )
                    {
                        double dist1 = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(cs.Seg1->get_noeud1()->get_coord(), cs.Seg1P);
                        double dist2 = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(cs.Seg1->get_noeud2()->get_coord(), cs.Seg1P);
                        if(dist1<dist2)
                            V2 = (MCVertex *) cs.Seg1->get_noeud1()->get_lien_topologie();
                        else
                            V2 = (MCVertex *) cs.Seg1->get_noeud2()->get_lien_topologie();
                    }
                    else if (cs.Seg1->get_noeud1()->get_lien_topologie()->get_dimension() == 0)
                        V2 = (MCVertex *) cs.Seg1->get_noeud1()->get_lien_topologie();
                    else if (cs.Seg1->get_noeud2()->get_lien_topologie()->get_dimension() == 0)     
                        V2 = (MCVertex *) cs.Seg1->get_noeud2()->get_lien_topologie();
                }
                // Merge vertices V1 with V2
                std::set<MCEdge*> edges;
                std::set<MCEdge*>::iterator itEdge;
                MergeVertices(V1,V2);

                delete G;    
                return true;
            }
        }
        delete G;
        return false;
}

//---------------------------------------------------------------------------
bool MCAA::NodeConstrictedSection(MCFace * __mcFace, MG_NOEUD * __n, ConstrictedSection & __cs, Graph::Graph * __faceBoundaryMeshGraph)
{
    int i,j,k;
    double LMin;

    if (__n->get_lien_topologie()->get_nb_ccf() != 0)
        return false;

    double meshSize = GetSize(__n->get_coord());
    LMin = .5*meshSize/GetMaxOverdensity();

    std::map<double,ConstrictedSection> mapCSByDistance;

    const Graph::Graph::MapNodesById & graphNodes = __faceBoundaryMeshGraph->GetNodes();
    for (Graph::Graph::MapNodesById::const_iterator itGN = graphNodes.begin();
        itGN != graphNodes.end();
        itGN++)
    {
        Graph::Node * graphNode = itGN->second;
        MG_SEGMENT * seg = (MG_SEGMENT *)graphNode->GetUserData();


        std::set<MCFace*>commonF;
        double dist;
        double *SP0 = seg->get_noeud1()->get_coord();
        double *SP1 = seg->get_noeud2()->get_coord();
        double *P = __n->get_coord();
        double P2[3];
        OT_ALGORITHME_GEOMETRIQUE::Closest_Point_to_Segment_3d(SP0,SP1,P,P2);
        dist=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(P,P2);
        if (dist < LMin)
            {
                if (__n==seg->get_noeud1() || __n==seg->get_noeud2())
                    continue;
                std::set <Graph::Node*> adjacentGNSegs, adjacentGNSegs2;
                graphNode->AdjacentNodes (adjacentGNSegs);
                for (std::set <Graph::Node*>::iterator itN1 = adjacentGNSegs.begin(); itN1 != adjacentGNSegs.end(); itN1++)
                {
                    Graph::Node * gn1 = *itN1;
                    MG_SEGMENT * seg1 = (MG_SEGMENT *)gn1->GetUserData();
                    if (__n==seg1->get_noeud1() || __n==seg1->get_noeud2())
                        continue;
                }

                ConstrictedSection cs;
                std::set<MCFace*>::iterator itF=commonF.begin();
                cs.F = *itF;
                cs.Node = __n;
                cs.Seg1 = seg;
                OT_ALGORITHME_GEOMETRIQUE::VEC3_ASSIGN_VEC3(P2,cs.Seg1P);
                mapCSByDistance[dist] = cs;
            }
    }

    for (std::map<double, ConstrictedSection>::iterator itCS=mapCSByDistance.begin();
        itCS != mapCSByDistance.end();
        itCS++)
    {
        double distance = itCS->first;
        MG_SEGMENT *sa = itCS->second.Seg1;
        OT_VECTEUR_3D point_a = itCS->second.Seg1P;

        for (int i=0; i<__n->get_lien_segment()->get_nb(); i++)
        {
            MG_SEGMENT * sb = __n->get_lien_segment()->get(i);
            OT_VECTEUR_3D point_b(__n->get_coord());

            std::vector <Graph::Node *> pathNodes;
            std::vector <Graph::Arc *> pathArcs;
            Graph::Node * nStart = __faceBoundaryMeshGraph->GetNode(sa->get_id());
            Graph::Node * nEnd = __faceBoundaryMeshGraph->GetNode(sb->get_id());
            if (nEnd == 0 || nStart == 0) // segment adjacent to node is not in the face
                continue;
            
            double distance = Graph::Dijkstra(__faceBoundaryMeshGraph,nStart,nEnd,GraphDistanceBetweenSegs,pathNodes,pathArcs);

            bool testShortestPath = (distance > 3*meshSize);

            if (testShortestPath == false)
                break;

            bool testInteriorDir;
            {
                    OT_VECTEUR_3D direction_b_a;
                    for (int i=0; i<3; i++)
                        direction_b_a[i] = point_a[i]-point_b[i];

                    if ( direction_b_a.get_longueur2() > 1E-32 )
                    {
                        double * xyz_topo[2]={point_b,point_a};
                        MCNode * mcNode[2]={0,0};
                        MG_ELEMENT_TOPOLOGIQUE * mcTopo[2], * refTopo[2];
                        mcTopo[0] = __n->get_lien_topologie();
                        mcTopo[1] = sa->get_lien_topologie();

                        for (int k=0;k<2;k++)
                        {
                            switch(mcTopo[k]->get_dimension())
                            {
                            case 0:
                            {
                                MCVertex * mcVertex = (MCVertex *)mcTopo[k];
                                mcNode[k] = new MCNode(mcVertex,mcVertex->GetRefVertex(),xyz_topo[k][0],xyz_topo[k][1],xyz_topo[k][2]);
                                break;
                            }         
                            case 1:
                            {          
                                MCEdge * mcEdge = (MCEdge*) mcTopo[k];
                                PolyCurve * polyCurve = mcEdge->GetPolyCurve();
                                MG_ARETE * refEdge = 0;
                                double s_polycurve,t_refEdge,dt_refEdge;
                                polyCurve->inverser(s_polycurve,xyz_topo[k],1E-6);
                                polyCurve->Parameter_SToRefEdgeT (s_polycurve, &refEdge, &t_refEdge, &dt_refEdge,false);
                                mcNode[k] = new MCNode(mcEdge, refEdge, t_refEdge, xyz_topo[k]);
                                break;
                            }
                            default:
                            {
                                refTopo[k] = 0;
                            }
                            }
                        }

                        ShortestPathByClosestPointOnEdge operateur_distGeo(__mcFace,mcNode[0],mcNode[1],.25*meshSize);
                        double distanceGeo = operateur_distGeo.Find();
                        if (distanceGeo > .5*meshSize)
                            testInteriorDir = false;
                        else
                            testInteriorDir = true;

                    }
            }

            if (testInteriorDir == false)
                break;

            if (testInteriorDir && testShortestPath)
            {
                __cs = itCS->second;
                return true;
            }
       }
    }

    return false;
}
double MCAA::GraphDistanceBetweenSegs(Graph::Node * __a, Graph::Node * __b, Graph::Arc * __arc )
{
    MG_SEGMENT * seg1 = (MG_SEGMENT *)__a->GetUserData();
    MG_SEGMENT * seg2 = (MG_SEGMENT *)__b->GetUserData();
    return .5*(seg1->get_longueur()+seg2->get_longueur());
}
//---------------------------------------------------------------------------
struct DJSeg {
    int nbVisit;
    double d;
    MG_SEGMENT * seg;
    MG_SEGMENT * previous;
};
//---------------------------------------------------------------------------
// Calcule le chemin le plus court entre les milieux des segments sa et sb
double
MCAA::ShortestPath(MG_MAILLAGE * __mesh, MG_SEGMENT * sa, MG_SEGMENT * sb)
{
    int i,j,k;

    if ( ! __mesh->contient(sa) || ! __mesh->contient(sb))
    {
            printf("ShortestPath : error : the mesh does not contain the querried segment!\n");
            return -1;
    }

     /*fonction Dijkstra (noeuds, fils, distance, debut, fin)*/
    double distance;
    std::map<MG_SEGMENT*, struct DJSeg> nodes;
    std::multimap<double, struct DJSeg*> nodesByD;

    /*
  Pour n parcourant noeuds
    n.parcouru = infini   // Peut �tre impl�ment� avec -1
    n.precedent = 0
  Fin pour*/
        MG_SEGMENT *seg;
        LISTE_MG_SEGMENT::iterator it_seg;
        for (MG_SEGMENT* seg = __mesh->get_premier_segment(it_seg);
                seg; seg =__mesh->get_suivant_segment(it_seg))
        {

            struct DJSeg djseg;
            
            djseg.previous = NULL;
            djseg.seg = seg;
            djseg.nbVisit = 0;
            djseg.d = 1E300;

            nodes[seg] = djseg;
        }

  if (nodes.find(sa) == nodes.end() || nodes.find(sb) == nodes.end())
    return 1E300; 

  /*  
  debut.parcouru = 0
  DejaVu = liste vide
  PasEncoreVu = noeuds
  */
  nodes[sa].d = 0;    

            
  std::set <MG_SEGMENT*> visited, unvisited;
  
            nodesByD.insert( std::make_pair( 0, &(nodes[sa])) );
  
        for (MG_SEGMENT* seg = __mesh->get_premier_segment(it_seg);
                seg; seg =__mesh->get_suivant_segment(it_seg))
        {
            unvisited.insert(seg);

        }
                      //

  // Tant que PasEncoreVu != liste vide
  while (unvisited.size() != 0)
  {
    // n1 = minimum(PasEncoreVu)   // Le noeud dans PasEncoreVu avec parcouru le plus petit
    struct DJSeg closestSeg;

    std::multimap<double, struct DJSeg *>::iterator itDJSegByDist;
    for (itDJSegByDist = nodesByD.begin();
        itDJSegByDist != nodesByD.end(); itDJSegByDist++)
    {
        if ( unvisited.find(itDJSegByDist->second->seg) != unvisited.end() )
        {
            closestSeg = *(itDJSegByDist->second);
            break;
        }
    }
    if (itDJSegByDist == nodesByD.end())
        closestSeg = nodes[*(unvisited.begin())];

    MG_SEGMENT * seg = closestSeg.seg;

    // DejaVu.ajouter(n1)
    //
    visited.insert(seg);
    nodes[seg].nbVisit ++;


    //      PasEncoreVu.enlever(n1)
    unvisited.erase(seg);


  /*
    Pour n2 parcourant fils(n1)   // Les noeuds reli�s � n1 par un arc
      Si n2.parcouru > n1.parcouru + distance(n1, n2)   // distance correspond au poids de l'arc reliant n1 et n2
        n2.parcouru = n1.parcouru + distance(n1, n2)
        n2.precedent = n1   // Dis que pour aller � n2, il faut passer par n1
      Fin si
    Fin pour
  Fin tant que
  */
    std::set <MG_SEGMENT*>children;
    TPL_LISTE_ENTITE<MG_SEGMENT*> *lien_segment1 = seg->get_noeud1()->get_lien_segment();
    for (i=0; i<lien_segment1->get_nb(); i++)
    {
        MG_SEGMENT * child = lien_segment1->get(i);
        if (child != seg )
            children.insert(child);
    }
    TPL_LISTE_ENTITE<MG_SEGMENT*> *lien_segment2 = seg->get_noeud2()->get_lien_segment();
    for (i=0; i<lien_segment2->get_nb(); i++)
    {
        MG_SEGMENT * child = lien_segment2->get(i);
        if (child != seg )
            children.insert(child);
    }

    for ( std::set <MG_SEGMENT*>::iterator itChild = children.begin();
        itChild != children.end();
        itChild++)
    {
        MG_SEGMENT * child = *itChild;
        struct DJSeg * sDJSegChild = &(nodes[child]);
        double dist_child = .5*seg->get_longueur() + .5*child->get_longueur();

        if ( sDJSegChild->d > closestSeg.d + dist_child)
        {
            sDJSegChild->d = closestSeg.d + dist_child;
            sDJSegChild->previous = seg;

            if ( unvisited.find(sDJSegChild->seg) != unvisited.end() )
                nodesByD.insert(std::make_pair(sDJSegChild->d, sDJSegChild));
        }
    }

  }

  /*   
  chemin = liste vide
  n = fin
  Tant que n != debut
    chemin.ajouterAvant(n)
    n = n.precedent
  Fin tant que
Retourner chemin
Fin fonction Dijkstra.
  */

  std::vector<MG_SEGMENT*> path;
  struct DJSeg n = nodes[sb];
  distance = n.d;

  int nb_seg = __mesh->get_nb_mg_segment();

  for (i=0; i<nb_seg && n.seg != sa; i++)
  {
    path.insert(path.begin(), n.seg);
    n = nodes[n.previous];
  }
  if (i == nb_seg)
    distance = 1E300;
  
  return distance;
}
//---------------------------------------------------------------------------
/*bool MCAA::TestEdgeSupprBoundary(MCVertex * __mcVertex)
{
    int i;

    if (_mcBody->G10()->GetArc(__mcVertex->get_id())->Rank() != 2)
        return false;

    // Get the 2 edges adjacent to this vertex
    std::vector < MCEdge * > edges;
    _mcBody->Vertex_GetAdjacentEdges(__mcVertex, edges);

    // Get node of this vertex in the FE mesh
    TPL_SET<MG_ELEMENT_MAILLAGE *> * lien_maillage = __mcVertex->get_lien_maillage();
    TPL_SET<MG_ELEMENT_MAILLAGE *>::ITERATEUR it;
    MG_ELEMENT_MAILLAGE * elem;
    MG_NOEUD * node = NULL;
    for (elem = lien_maillage->get_premier(it); elem ; elem = lien_maillage->get_suivant(it))
    {
        if (_feMesh->contient(elem))
        {
            node = (MG_NOEUD*)elem;
            break;
        }
    }          
    if (!node)
    {
        printf("Warning : no node was found for MCVertex %ld in the FE mesh !\n", __mcVertex->get_id());
        return false;
    }

    // Get the two segments adjacent to this node  
    MG_SEGMENT * segs[2];
    segs[0] = node->get_lien_segment()->get(0);
    segs[1] = node->get_lien_segment()->get(1);

    // Get the two local edge criteria of these segments
    LocalEdgeCriteria * edgeCrit[2];
    for (i=0;i<2;i++)
    {
        MCEdge * edge = (MCEdge*)segs[i]->get_lien_topologie();

        ECMapIterator it_ec = _ecMap.find(edge);

        if (it_ec != _ecMap.end())
        {
            GlobalEdgeCriteria * gec = it_ec->second;
            edgeCrit[i] = gec->GetLocalEdgeCriteria(segs[i]);
            if (!edgeCrit[i])
            {
                printf("Warning : no criteria were found for segment %ld of edge %ld\n", segs[i]->get_id(), edge->get_id());
                return false;
            }
        }
        else
        {
            printf("TestEdgeSupprBoundary : Error : vertex %ld is not found\n");
            return false; 
        }
    }

    return ( edgeCrit[0]->DeletionScore() * edgeCrit[1]->DeletionScore() < 0);
}  */

void V(MCAA * mcaa)
{
        MCBody * _mcBody = mcaa->GetMCBody();
    {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge22){
               if (_mcBody->G21()->GetArcs().size() <= 1)
                   break;
               if (mcaa->CheckIfTopoExists(edge22) == false)
                   printf("MC Edge ID=%lu not OK\n", edge22->get_id());
           }
    }

}
//---------------------------------------------------------------------------
void
MCAA::SuppressMCEdge(MCEdge * __mcEdge, MCTChanges * __mctChanges)
{
        MCFace * F[3];

        // update the simplification time
        time++;

        // Get the list of vertex that are adjacent to the
        // mc edge
        std::vector<MCVertex*> edge_vertices;
        _mcBody->Edge_GetVertices(__mcEdge, edge_vertices);

        // remove the FE mesh of this edge
        DeleteFEMesh(__mcEdge);
        // Delete the criteria of the MC Edge
        EC_Delete(__mcEdge);

        // Suppress the mc edge of the mc body
        _mcBody->SuppressMCEdge(__mcEdge, F);
        //        V(this);

        // the current MCFace containing to the deleted edge
        MCFace * newFace = F[2];

        // move the link of the edge to the current MCFace containing to the deleted edge
        MG_MAILLAGE_OUTILS::change_lien_maillage( __mcEdge, newFace );
        for (int i=0; i<2; i++)
        {
            if (F[i] != NULL)
                {
                    // a face that was deleted from MCT
                    MCFace * oldFace = F[i];
                    // update mesh link of old face to new face
                    MG_MAILLAGE_OUTILS::change_lien_maillage( oldFace, newFace );
                    // delete old face
                    _mcBody->DeleteMCFace(oldFace);
                    //V(this);
                }
        }

        std::set<MCEdge*> adjedge = _mcBody->Face_GetAdjacentEdges(newFace);
        for ( std::set<MCEdge*>::iterator itAdjEdge = adjedge.begin();
               itAdjEdge != adjedge.end(); itAdjEdge++)
        {
                MCEdge * edge = *itAdjEdge;
                GlobalEdgeCriteria * gec = GetGlobalEdgeCriteria(edge);
                gec->Update(__mcEdge, false);
                //V(this);
        }

        // Dereference the MCEdge from the geometry manager, and free the memory
        _mcBody->DeleteMCEdge(__mcEdge);
        //V(this);

        // update the vertex properties on the adjacent vertices
        for (std::vector<MCVertex*>::iterator itVertex = edge_vertices.begin();
                itVertex != edge_vertices.end();
                itVertex++ )
        {
                MCVertex * v = *itVertex;
                VC_Get(v)->Update();
                //V(this);
        }              

        // update the modification time of the MC Edge
        newFace->time = time;
        __mctChanges->E.insert(std::make_pair(__mcEdge, (MCEdge*)NULL));
        if (F[0]) __mctChanges->F.insert(std::make_pair(F[0], F[2]));
        if (F[1]) __mctChanges->F.insert(std::make_pair(F[1], F[2]));

        //V(this);
}
//---------------------------------------------------------------------------
void
MCAA::SuppressMCVertex(MCVertex * __mcVertex, MCTChanges * __mctChanges)
{                                          
        std::multimap<MCEdge*,MCEdge*>::iterator itE;
        
        // update the simplification time
        time++;

        if (_mapSuppressRefVertexCount.find(__mcVertex->GetRefVertex()) ==  _mapSuppressRefVertexCount.end())
            _mapSuppressRefVertexCount[ __mcVertex->GetRefVertex() ] = 1;
        else
        {
            _mapSuppressRefVertexCount[ __mcVertex->GetRefVertex() ]++;
        }

        // Get the edges adjacent to this vertex
        std::set< MCEdge * > vertex_edges;
        _mcBody->Vertex_GetAdjacentEdges(__mcVertex, vertex_edges);

        if ( vertex_edges.size() != 0 && vertex_edges.size() != 2)
        {
          printf("P B vertex %ld has %lx edges\n", __mcVertex->get_id(), (long unsigned int)vertex_edges.size());
        }

        std::set< MCFace * > vertex_faces;
        _mcBody->Vertex_GetAdjacentFaces(__mcVertex, vertex_faces);
        if (vertex_edges.size() == 0) //  case of a vertex isolated in the face
        {
            if (vertex_faces.size() == 1)
            {
                std::set< MCFace * >::iterator vertex_faces_it = vertex_faces.begin();
                // delete the node of the FE mesh
                DeleteFEMesh(__mcVertex);
                // update the link of the deleted vertex in the MC Tessellation
                MG_MAILLAGE_OUTILS::change_lien_maillage(__mcVertex, *vertex_faces_it);
            }
            else
            {
                printf("Problem : vertex has 0 edge and 0 face !!!\n");
            }
        } 

        // Suppress this vertex in the MCBody
        MCEdge * newEdge = NULL;

        _mcBody->SuppressMCVertex(__mcVertex, __mctChanges);

        if ( __mctChanges->E.size() != 0)
        {
            itE = __mctChanges->E.begin();
            newEdge = itE->second;

            for (itE = __mctChanges->E.begin(); itE != __mctChanges->E.end(); itE++)
            {
                MCEdge * edge0 = itE->first;
                MCEdge * edge1 = itE->second;

                printf("Created edge %ld containing %ld %ld", newEdge->get_id(), edge0->get_id(), edge1->get_id());

                // If edge0 was deleted
                if ( edge0 != 0 )
                {
                    // move the link of the egde to the current MCFace containing to the deleted edge
                    MG_MAILLAGE_OUTILS::change_lien_maillage( edge0, edge1, _mcTess);
                    MG_MAILLAGE_OUTILS::change_lien_maillage(__mcVertex, edge1, _mcTess);

                    // update the modification time of the MC Edge
                    newEdge->time = time;
                }
            }
            printf("\n");

        std::set <MCEdge*> faces_edges;
        for (std::set<MCFace*>::iterator itFace = vertex_faces.begin();
            itFace != vertex_faces.end();
            itFace++)
        {
            MCFace * mcFace = *itFace;

            // Get the G21 arcs to G21 node
            // they compose the MC edges adjacent to the MC face    
            Graph::Node * G21Node = _mcBody->G21()->GetNode(mcFace->get_id());

            Graph::Node::MultimapArcsById & arcsOfG21Node = G21Node->IncidentArcs();

            for (Graph::Node::MultimapArcsById::const_iterator it = arcsOfG21Node.begin();
                it != arcsOfG21Node.end();
                    it++)
                faces_edges.insert( _mcBody->GetMCEdge(it->second) );
        }

        for (std::set<MCEdge*>::iterator itFaceEdge = faces_edges.begin();
            itFaceEdge != faces_edges.end();
            itFaceEdge++)
        {
            MCEdge * faceEdge = *itFaceEdge;
            if ( vertex_edges.find(faceEdge) != vertex_edges.end() )
                continue;
            GlobalEdgeCriteria * gec = GetGlobalEdgeCriteria(faceEdge);
            if (gec == NULL ) continue;
            
            for (std::set<MCEdge*>::iterator itDeletedVertexEdge = vertex_edges.begin();
                itDeletedVertexEdge != vertex_edges.end();
                itDeletedVertexEdge++)
            {
                MCEdge * deletedEdge = *itDeletedVertexEdge;
                gec->Update(deletedEdge, false);
            }
        }

            // delete edge criteria of each original mc Edge
            // initialize the edge criteria on the new merged MCEdge
            bool bRemeshMergedEdges = true;
            if ( bRemeshMergedEdges )
            {
                DeleteFEMesh(newEdge); 
                DeleteFEMesh(__mcVertex);
                InitializeFEMesh(newEdge);
                EC_Init(newEdge);

                for (itE = __mctChanges->E.begin(); itE != __mctChanges->E.end(); itE++)
                {
                    MCEdge * oldEdge2 =  itE->first;

                    // delete the edge criteria of the original mc edge
                    EC_Delete(oldEdge2);
                }
            }
            else
            {
                for (itE = __mctChanges->E.begin(); itE != __mctChanges->E.end(); itE++)
                {
                    MCEdge * oldEdge2 =  itE->first;
                    MCEdge * newEdge2 =  itE->second;
                    if (oldEdge2)
                    {


                        // move local edge criterie to new merged mc edge
                        EC_ChangeTopo(oldEdge2, newEdge2);
                        GetGlobalEdgeCriteria(newEdge2)->Update(false);
                
                        // delete the edge criteria of the original mc edge
                        EC_Delete(oldEdge2);
                    }
                }
            }
        }

        // the vertex properties of the deleted mc vertex are deleted
        VC_Delete(__mcVertex);

        // Delete Old vertex
        _mcBody->DeleteMCVertex(__mcVertex);

        // Delete old MC Edges
        for (std::set<MCEdge*>::iterator itE = vertex_edges.begin();
                itE != vertex_edges.end();
                itE++ )
        {
            MCEdge * mcEdge =  *itE;
            _mcBody->DeleteMCEdge(mcEdge);
        }
}

//---------------------------------------------------------------------------
bool MCAA::SplitEdgesAtIncompatibleLocalCriteria(MCTChanges * __mctChanges)
{
    int i;
    double bestSplitScore = 0, bestSplitPoint[3];
    MCEdge * bestSplitEdge = NULL;
    RTChanges rtChanges;

    for (ECMapIterator it_ec = _ecMap.begin();
        it_ec != _ecMap.end();
        it_ec++)
    {
        MCEdge * mcEdge = it_ec->first;
        GlobalEdgeCriteria * gec = it_ec->second;
        double splitScore = 0, splitPoint[3];
        splitScore = gec->SplitScore(splitPoint);
        if (splitScore > bestSplitScore)
        {
            bestSplitScore = splitScore;
            bestSplitEdge = mcEdge;
            for (i=0;i<3;i++)
                bestSplitPoint[i] = splitPoint[i];

            break;
        }
    }

    if (bestSplitEdge)
    {
        SplitEdge(bestSplitEdge, bestSplitPoint, __mctChanges, &rtChanges);
        return true;
    }
    else
        return false;
}
//---------------------------------------------------------------------------
void MCAA::SplitEdge(MCEdge * __mcEdge, double xyz[3], MCTChanges * __mctChanges, RTChanges * __rtChanges, bool __splitTessellation)
{
    int i;

    // Get the list of faces that are adjacent to the
    // mc edge
    std::set < MCFace * > edge_faces = _mcBody->Edge_GetAdjacentFaces(__mcEdge);


    // Split the MC Edge
    MG_ARETE *__origRefEdge, *__splitRefEdges[2];
    MG_SOMMET *__splitRefVertex;
    MCVertex * __splitVertex;
    MCEdge * __splitEdges[2];

    __origRefEdge = __mcEdge->GetPolyCurve()->GetRefEdge(0);

    _mcBody->SplitEdge(__mcEdge, xyz, __mctChanges, __rtChanges);

    if (__mctChanges->E.size() == 0)
        return;

    std::multimap<MCEdge*, MCEdge*>::iterator itMCE;
    std::multimap<MCVertex*, MCVertex*>::iterator itMCV;

    itMCE = __mctChanges->E.begin();
    for (i=0; i<2; i++, itMCE++)
        __splitEdges[i] = itMCE->second;

    itMCV = __mctChanges->V.begin();
    __splitVertex = itMCV->second;

    std::multimap<MG_ARETE*, MG_ARETE*>::iterator itRE;
    std::multimap<MG_SOMMET*, MG_SOMMET*>::iterator itRV;
    
    itRE = __rtChanges->E.begin();
    for (i=0; i<2; i++, itRE++)
    {
        __origRefEdge = itRE->first;
        __splitRefEdges[i] = itRE->second;
    }

    itRV = __rtChanges->V.begin();
    __splitRefVertex = itRV->second;

    printf("Edge %ld has been substituted by edges %ld and %ld\n", __mcEdge->get_id(), __splitEdges[0]->get_id(), __splitEdges[1]->get_id());

    MG_SEGMENT * mcTessOriginalSegment, *refTessOriginalSegment;
    // split the tessellation
    if (__splitTessellation)
    {
        MG_NOEUD * mcTessSplitNode;
        MG_SEGMENT * mcTessSplitSegments[2];
        MG_MAILLAGE_OUTILS::SplitEdgeInMesh(_mcTess, __mcEdge, __splitVertex, __splitEdges[0], __splitEdges[1], &mcTessOriginalSegment, &mcTessSplitNode, mcTessSplitSegments);

        MG_NOEUD * refTessSplitNode;
        MG_SEGMENT * refTessSplitSegments[2];
        if (__origRefEdge) MG_MAILLAGE_OUTILS::SplitEdgeInMesh(_refTess, __origRefEdge, __splitRefVertex, __splitRefEdges[0], __splitRefEdges[1], &refTessOriginalSegment, &refTessSplitNode, refTessSplitSegments);

        // update MC Tessellation Grid
        _UpdateMCTessSegGrid(__mcEdge, mcTessOriginalSegment,mcTessSplitSegments);
    }
    if (__origRefEdge) _geom->supprimer_mg_areteid(__origRefEdge->get_id());
    
    // remesh the split edges
    DeleteFEMesh(__mcEdge);
    DeleteFEMesh(__splitVertex);
    MG_NOEUD * feNode;
    InitializeFEMesh(__splitVertex, &feNode);
    for (i=0; i<2; i++)
    {
        InitializeFEMesh(__splitEdges[i]);
    }

    // update vertex criteria
    VC_Init(__splitVertex);

    // update edge criteria
    EC_Delete(__mcEdge);       
    for (i=0; i<2; i++)   
        EC_Init(__splitEdges[i]);
    _mcBody->DeleteMCEdge(__mcEdge);

    if (mcTessOriginalSegment)
    { MCBODY_FOR_EACH_MCEDGE ( _mcBody, mcEdge )
    {
        GlobalEdgeCriteria * gec = GetGlobalEdgeCriteria(mcEdge);
        gec->ReInit(mcTessOriginalSegment);
        gec->Update(__mcEdge, false);
    }
    }

    
    // delete the original segments in MC and Ref tessellations
    if (mcTessOriginalSegment)
        _mcTess->supprimer_mg_segmentid(mcTessOriginalSegment->get_id());
    if (refTessOriginalSegment)
        _refTess->supprimer_mg_segmentid(refTessOriginalSegment->get_id());

    // update the simplification time
    time++;
    __splitVertex->time = time;  
    for (i=0; i<2; i++)
        __splitEdges[i]->time = time;
        
    Graph::Arc * G21Arc = _mcBody->G21()->GetArc( __splitEdges[0]->get_id() );
    for (Graph::Arc::MultimapNodesById::const_iterator itN = G21Arc->Nodes().begin();
        itN != G21Arc->Nodes().end();
        itN++)
        {
            MCFace * adjacentFace = _mcBody->GetMCFace(itN->second);
            adjacentFace->time = time;
        }
}

//---------------------------------------------------------------------------
void MCAA::Simplify()
{
    MCTChanges mctChanges;
    Simplify2(&mctChanges);
}
//---------------------------------------------------------------------------
void
MCAA::Simplify2( MCTChanges * __mctChanges )
{
    int nbFace0, nbFace1;
    int nbEdge0, nbEdge1;
    int nbVertex0, nbVertex1;

    nbFace0 = _mcBody->G21()->GetNodes().size();
    nbEdge0 = _mcBody->G21()->GetArcs().size();
    nbVertex0=_mcBody->G20()->GetArcs().size();

    std::set<MG_SEGMENT*> setSegmentBadTopology;
    std::set<MG_NOEUD*> setNodeBadTopology;

    char animFilename[] = "";//"c:/gilles/anim.iv";
    std::ofstream animFile;
    double animTime = 15; double animTimeStep = 15.0f/nbEdge0;
    InventorText_MCAA * ivText = NULL;
    if ( strlen(animFilename) )
    {
        animFile.open(animFilename);
        ivText = new InventorText_MCAA(this);
        ivText->ShowFEMesh = false;
        ivText->ShowRefModel = false;
        ivText->ShowEdgeIds = false;
        ivText->ShowCriteria = false;
        ivText->ColorMCEdge = InventorText_MCAA::GlobalSuppressionCriterion;
        animFile << ivText->GetText_BeginAnimation(0,animTime);
    }

    bool result= true;
    while ( result )
    {
        MCTChanges mctChanges;

        result = Next2(&mctChanges);
        __mctChanges->Merge(mctChanges);
        if (ivText)
        {
            animFile << ivText->GetText_StepAnimation(animTime,animTime+animTimeStep);
            animTime+=animTimeStep;
        }
        // Stop simplification if there remains only 1 MC Edge in the MCBody !
        if (_mcBody->G21()->GetArcs().size() == _minMCEdgeCount)
            break;


    }                        
    _nbE += SimplifyFaceLoop();  
    if (ivText)
    {
        animFile << ivText->GetText_StepAnimation(animTime,animTime+animTimeStep);
        animTime+=animTimeStep;
    }
    {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge22){
        if (_mcBody->G21()->GetArcs().size() <= 1)
                break;
        if (CheckIfTopoExists(edge22) == false)
                printf("MC Edge ID=%lu not OK\n", edge22->get_id());
    }
    }
    result= true;
    while ( result )
    {
        MCTChanges mctChanges;

        result = Next2(&mctChanges);    
        if (ivText)
        {
            animFile << ivText->GetText_StepAnimation(animTime,animTime+animTimeStep);
            animTime+=animTimeStep;
        }
        __mctChanges->Merge(mctChanges);
    }
    _nbC += SimplifyEdgeCollapse();

    if (ivText)
    {
        animFile << ivText->GetText_StepAnimation(animTime,animTime+animTimeStep);
        animTime+=animTimeStep;
    }
    _nbD += SimplifyConstrictedSections();


    if (ivText)
    {
        animTimeStep = 1E9;
        animFile << ivText->GetText_StepAnimation(animTime,animTime+animTimeStep);
        animTime+=animTimeStep;
        animFile << ivText->GetText_EndAnimation();
    }


    nbFace1 = _mcBody->G21()->GetNodes().size();
    nbEdge1 = _mcBody->G21()->GetArcs().size();
    nbVertex1=_mcBody->G20()->GetArcs().size();

    printf("The initial CAD model contained :\n%d faces\n%d edges\n%d vertices\n", nbFace0,nbEdge0,nbVertex0);
    printf("%d edge deleted\n", _nbA);
    printf("%d vertex deleted\n", _nbB);
    printf("%d edge collapsed\n", _nbC);
    printf("%d constricted sections collapsed\n", _nbD);
    printf("%d loops deleted\n", _nbE);
    printf("The final CAD model contains :\n%d faces\n%d edges\n%d vertices\n", nbFace1,nbEdge1,nbVertex1);
}

int MCAA::SimplifyFaceLoop()
{
    int nbLoopDeleted = 0;
    int i;
    do
    {
        i = NextSimplifyFaceLoop();
        nbLoopDeleted += i;
    }
    while (i);
    return nbLoopDeleted;
}

int MCAA::NextSimplifyFaceLoop()
{
    int loopCount = 0;
    std::multimap < double , LoopCriteria * > mapLoopCriteria;
    {MCBODY_FOR_EACH_MCFACE(_mcBody,face)
    {
        std::vector < std::vector < MCEdge * > > loops = _mcBody->Face_GetCycles(face);
        std::vector < LoopCriteria * > vecLoopCriteria;
        for (unsigned itLoop = 0; itLoop < loops.size(); itLoop++)
        {
            LoopCriteria * crit = new LoopCriteria(this, face, loops[itLoop], 5);
            double score = crit->DeletionScore();
            mapLoopCriteria.insert ( std::make_pair( score, crit ) );
        }
    }}


    std::multimap < double , LoopCriteria * >::reverse_iterator itHighestScore = mapLoopCriteria.rbegin();
    if (itHighestScore != mapLoopCriteria.rend() && itHighestScore->first > 0)
    {
        LoopCriteria * highestLoopCriteria = itHighestScore->second;

            if (_mcBody->Contains(highestLoopCriteria->GetFace()) &&
                highestLoopCriteria->DeletionScore() > 0)
            {
                std::set < MCVertex * > loopVertices;
                std::vector < MCEdge* > loop = highestLoopCriteria->GetLoop();   
                loopCount++;

#define SUPPRESS_ONE_EDGE_OF_SMALL_LOOP
#ifdef SUPPRESS_ONE_EDGE_OF_SMALL_LOOP
                MCEdge * mcEdgeMax; double maxL = -1;   
                MCEdge * mcEdgeMin; double minL = 1E300;
                for (unsigned itEdge = 0; itEdge < loop.size(); itEdge ++ )
                {                                   
                    std::vector<MCVertex *> extremities;
                    MCEdge * mcEdge = loop[itEdge];
                    _mcBody->Edge_GetVertices(mcEdge, extremities);
                    for (unsigned i=0; i<extremities.size(); i++)
                        loopVertices.insert ( extremities [i] );
                        
                    double L = _ecMap[mcEdge]->GetLength();
                    if (L > maxL)
                    {
                        maxL = L;
                        mcEdgeMax = loop[itEdge];
                    }   
                    if (L < minL)
                    {
                        minL = L;
                        mcEdgeMin = loop[itEdge];
                    }
                }

                MCTChanges mctChanges3,mctChanges1;
                SuppressMCEdge(mcEdgeMin, &mctChanges3);
                mctChanges1.Merge(mctChanges3);
#endif

#ifdef SUPPRESS_ALL_EDGES_OF_SMALL_LOOP
                for (unsigned itEdge = 0; itEdge < loop.size(); itEdge ++ )
                {
                    MCEdge * mcEdge = loop[itEdge];
                    std::vector<MCVertex *> extremities;
                    _mcBody->Edge_GetVertices(mcEdge, extremities);
                    for (unsigned i=0; i<extremities.size(); i++)
                        loopVertices.insert ( extremities [i] );

                    double splitScore, splitPoint[3];
                    splitScore = highestLoopCriteria->EdgeSplitScore(mcEdge, splitPoint);
                                          
                    MCTChanges mctChanges1;

                    MCTChanges mctChanges3;
                    SuppressMCEdge(mcEdge, &mctChanges3);
                    mctChanges1.Merge(mctChanges3);
                }

#endif //#ifdef SUPPRESS_ALL_EDGES_OF_SMALL_LOOP
                for (std::set<MCVertex*>::iterator itVertex = loopVertices.begin();
                    itVertex != loopVertices.end();
                    itVertex++)
                {
                    MCVertex * mcVertex = *itVertex;
                    VertexCriteria * vCrit = _vcMap[mcVertex];
                    if ( vCrit->GetScore() > 0)
                    {
                        MCTChanges mctChanges1;
                        SuppressMCVertex(mcVertex, &mctChanges1);
                        _nbB++;
                    }
                }
            }
    }

    for ( std::multimap < double , LoopCriteria * >::iterator itLoopCriteria = mapLoopCriteria.begin();
        itLoopCriteria != mapLoopCriteria.end();
        itLoopCriteria++)
    {
        LoopCriteria * crit = itLoopCriteria->second;
        delete crit;
    }

    return loopCount;
}

int MCAA::SimplifyEdgeCollapse()
{                 
    bool result;
    int counter;

    {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge22){
               if (_mcBody->G21()->GetArcs().size() <= 1)
                   break;
               if (CheckIfTopoExists(edge22) == false)
                   printf("MC Edge ID=%lu not OK\n", edge22->get_id());
           }
    }

    counter=0;
    do
    {
        counter++;
        MCEdge * collapseEdge = NULL;
        MCVertex * deleteVertex = NULL;
        double bestScore = 0;    

        if (_mcBody->G21()->GetArcs().size() <= 1)
            break;

        {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge)
        {
            if (_mcBody->G21()->GetArcs().size() <= 1)
                break;
                
            MCVertex * vertex[2];
            vertex[0] = (MCVertex * ) edge->get_cosommet1()->get_sommet();
            vertex[1] = (MCVertex * ) edge->get_cosommet2()->get_sommet();

            for (int i=0; i<2; i++)
            {
                double score = EdgeCollapseCriteria(edge, vertex[i], this).GetScore();
                if (score > bestScore)
                {                   
                    bestScore = score;
                    collapseEdge = edge;
                    deleteVertex = vertex[i];
                }
            }


        }
        }
               
        if (collapseEdge&&deleteVertex)
            result = CollapseMCEdgeToMCVertex(collapseEdge, deleteVertex);
        else
            result = false;

    } while (result);

    {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge22){
               if (_mcBody->G21()->GetArcs().size() <= 1)
                   break;
               if (CheckIfTopoExists(edge22) == false)
                   printf("MC Edge ID=%lu not OK\n", edge22->get_id());
           }
    }

    return counter;
}

bool MCAA::CollapseMCEdgeToMCVertex(MCEdge * collapseEdge, MCVertex * deleteVertex)
{
    bool result = true;

        result = true;  
        MCVertex * vertex[2];
        vertex[0] = (MCVertex * ) collapseEdge->get_cosommet1()->get_sommet();
        vertex[1] = (MCVertex * ) collapseEdge->get_cosommet2()->get_sommet();
        MCVertex * targetVertex = (vertex[0]!=deleteVertex)?vertex[0]:vertex[1];
        std::set <MCEdge*> setNewEdges, setOldEdges;
        std::set <MCEdge*>::iterator itEdge, itEdge2;
        _mcBody->ContractEdgeToVertex(collapseEdge, targetVertex, setNewEdges, setOldEdges);
        
        for (itEdge = setOldEdges.begin(); itEdge != setOldEdges.end(); itEdge++)
        {
            MCEdge * mcEdge = *itEdge;
            EC_Delete(mcEdge);

            std::vector <MG_SEGMENT*> vecSegs;
            TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = mcEdge->get_lien_maillage();
            TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
            MG_ELEMENT_MAILLAGE* element2;
            for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
            {
                if (_feMesh->contient(element2))
                        vecSegs.push_back((MG_SEGMENT*)element2);
            }

            for (std::vector<MG_SEGMENT *>::iterator it = vecSegs.begin();
                it != vecSegs.end();
                it++)
            {
                MG_SEGMENT * seg = *it;

                _feMesh->supprimer_mg_segmentid(seg->get_id());
                
                for (itEdge2 = setNewEdges.begin(); itEdge2 != setNewEdges.end(); itEdge2++)
                {
                    MCEdge * mcEdge2  = *itEdge2;
                    mcEdge2->get_lien_maillage()->supprimer(seg);
                }
            }
        }

        for (itEdge = setNewEdges.begin(); itEdge != setNewEdges.end(); itEdge++)
            DeleteFEMesh(*itEdge);
                
        DeleteFEMesh(deleteVertex);
        VC_Delete(deleteVertex);
        
        TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage2 = deleteVertex->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it2;
        MG_ELEMENT_MAILLAGE* element2;
        for (element2 = lien_maillage2->get_premier(it2); element2; element2 = lien_maillage2->get_suivant(it2) )
        {
            MCEdge * mcEdgeOfDeletedVertex=*(setNewEdges.rbegin());
            element2->change_lien_topologie2(mcEdgeOfDeletedVertex);
        }

        _mcBody->DeleteMCVertex(deleteVertex);
        for (itEdge = setOldEdges.begin(); itEdge != setOldEdges.end(); itEdge++)
        {
            MCEdge * mcEdge = *itEdge;
            _mcBody->DeleteMCEdge(mcEdge);
        }
                                                          
        for (itEdge = setNewEdges.begin(); itEdge != setNewEdges.end(); itEdge++)
            InitializeFEMesh(*itEdge);
                                                      
        for (itEdge = setNewEdges.begin(); itEdge != setNewEdges.end(); itEdge++)
            EC_Init(*itEdge);

        {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge22){
                   if (_mcBody->G21()->GetArcs().size() <= 1)
                       break;
                   if (CheckIfTopoExists(edge22) == false)
                           printf("MC Edge ID=%lu not OK\n", edge22->get_id());
               }
        }

                
/*    std::set<MG_SEGMENT*> setSegmentBadTopology;
    std::set<MG_NOEUD*> setNodeBadTopology;
    CheckMCMesh( this->GetMCTess(), setSegmentBadTopology, setNodeBadTopology );
    CheckMCMesh( this->GetFEMesh(), setSegmentBadTopology, setNodeBadTopology );*/

    return result;
}

//---------------------------------------------------------------------------
bool MCAA::Next2( MCTChanges * __mctChanges )
{
    bool bSimplifDone = false;
    if (!bSimplifDone && _minMCVertexCount < _mcBody->G10()->GetArcs().size() ) bSimplifDone = SuppressNextVertex(__mctChanges);
    if (bSimplifDone) _nbB++;
    if (!bSimplifDone && _minMCEdgeCount < _mcBody->G10()->GetNodes().size()) { bSimplifDone = SuppressNextEdge2(__mctChanges); if (bSimplifDone) _nbA++; }

    return bSimplifDone;
}
//---------------------------------------------------------------------------
bool MCAA::SuppressNextEdge2(MCTChanges * __mctChanges)
{
    GlobalEdgeCriteria * ec = GetHighestEdgeDeletionScore();

    if (ec == NULL)
        return false;

    bool bSimplifDone = (ec->DeletionScore() > 1E-3);

    if (bSimplifDone)
    {
        RTChanges rtChanges;
        DeleteEdge2(ec->GetEdge(), __mctChanges, &rtChanges);
    }

    return bSimplifDone;
}               
//---------------------------------------------------------------------------
void
MCAA::DeleteEdge2(MCEdge * __mcEdge, MCTChanges * __mctChanges, RTChanges * __rtChanges)
{
    int i;
    double splitScore = 0, splitPoint[3];
    RTChanges rtChanges;
    double deletionScore;

    GlobalEdgeCriteria * gec=_ecMap[__mcEdge];
    splitScore = gec->SplitScore(splitPoint);
    if (splitScore > 0)
    {
        MCTChanges mctChanges;
        RTChanges rtChanges;
        SplitEdge(__mcEdge, splitPoint, &mctChanges, &rtChanges);
        if (mctChanges.E.size() != 0) // Successful split !!!
        {
            for (std::multimap<MCEdge*,MCEdge*>::iterator itEdge = mctChanges.E.begin();
                itEdge !=  mctChanges.E.end();
                itEdge++)
            {
                if (itEdge->second != NULL) // edge created by SplitEdge()
                {
                    DeleteEdge2(itEdge->second, __mctChanges, __rtChanges);
                    __mctChanges->Merge(mctChanges);
                }
            }
        }
        else // failure : there is 0 edge after split
        {
            deletionScore=gec->DeletionScore();
            if (deletionScore>0)
                SuppressMCEdge(__mcEdge, &mctChanges);
            __mctChanges->Merge(mctChanges);
        }
    }
    else
    {
        MCTChanges mctChanges;
        deletionScore=gec->DeletionScore();
        if (deletionScore>0)
            SuppressMCEdge(__mcEdge, &mctChanges);
        __mctChanges->Merge(mctChanges);
    }
}
Revision:	569
Committed:	Thu Oct 16 14:36:31 2014 UTC (10 years, 6 months ago) by foucault
File size:	97734 byte(s)
Log Message:	Mise à jour pour CAD4FE (Gilles) : operation 1 (tentative)