CAD4FE/src/CAD4FE_MCBody.cpp

#include <vector>
#include <string>     
#include <sstream>   
#include <fstream>
#include <algorithm>
#include <math.h>
#pragma hdrstop

// Headers of Magic

#include "gestionversion.h"
#include "lc_point.h"
#include "step_point.h"

#include "CAD4FE_MCFace.h"
#include "CAD4FE_MCEdge.h"
#include "CAD4FE_PolyCurve.h"
#include "CAD4FE_PolySurface.h"
#include "CAD4FE_MCVertex.h"
#include "mg_maillage.h"       
#include "CAD4FE_mg_utils.h"     
#include "mg_gestionnaire.h"
#include "ot_mathematique.h"

// Headers of CAD4FE 
#include "CAD4FE_Geometric_Tools.h"
#include "CAD4FE_MakeLoops.h"
#include "CAD4FE_MCBody.h"
#include "CAD4FE_ColorMap.h"
#include "CAD4FE_Criteria.h"
#include "hypergraphlib_components.h"
#include "hypergraphlib_neato.h"
#include "hypergraphlib_findcycles.h"

using namespace CAD4FE;


void MCTChanges::Merge(MCTChanges& __mctChanges)
{
    /*int i;

    std::multimap<MCEdge*,MCEdge*>::iterator itEdge1, itEdge2;
    for ( itEdge1 = E.begin();
        itEdge1 != E.end();
        itEdge1++)
    {
        i=0;
        for ( itEdge2 = __mctChanges.E.find(itEdge1->second);
            itEdge2->first == itEdge1->second && itEdge2 != __mctChanges.E.end();
            itEdge2++)
        {
            if (itEdge2->second == itEdge2->first)
                continue;
            E.insert(std::make_pair(itEdge1->first, itEdge2->second));
            i++;
        }
        if (i>0){
            itEdge2 = itEdge1;
            itEdge2 --;
            E.erase(itEdge1);
            itEdge1 = itEdge2;
        }
    }
    
    for ( itEdge2 = __mctChanges.E.begin();
        itEdge2 != __mctChanges.E.end();
        itEdge2++)
    {
        itEdge1 = E.find(itEdge2->first);
        if (itEdge1 == E.end())
        {
            E.insert(std::make_pair(itEdge1->first,itEdge1->second));
        }
    }*/
}

//---------------------------------------------------------------------------
MCBody::MCBody(MG_GEOMETRIE* __refGeom, MG_VOLUME * __refBody)
{
    time = 0;
    _refBody = __refBody;
    _refGeom = __refGeom; 
    InitHyperGraphs();
}

//---------------------------------------------------------------------------
MCBody::~MCBody()
{
    delete _G10;
    delete _G20;
    delete _G21;
}

//---------------------------------------------------------------------------
void MCBody::Graph_SetUserData(MG_ELEMENT_TOPOLOGIQUE * __topo)
{
        Graph_SetUserData(_G10, __topo);
        Graph_SetUserData(_G20, __topo);
        Graph_SetUserData(_G21, __topo);
}

//---------------------------------------------------------------------------
void MCBody::Graph_SetUserData(Graph::Graph * __G, MG_ELEMENT_TOPOLOGIQUE * __topo)
{
        int id = __topo->get_id();

        switch (__topo->get_dimension())
        {
        case 0:// MCVertex  
        {
                // reference vertex in arc of g10 and g20
                if (__G == _G10 || __G == _G20 )
                        __G->GetArc(id)->SetUserData(__topo);
                                   
                break;
        }
        case 1:// MCEdge
        {
                // reference edge in node of g10 and arc of g21
                if (__G == _G10 )
                        __G->GetNode(id)->SetUserData(__topo);     
                else if (__G == _G21 )
                        __G->GetArc(id)->SetUserData(__topo);

                break;
        }
        case 2:// MCFace
        {
                // reference face in node of g20 and g21
                if (__G == _G20  || __G == _G21 )
                        __G->GetNode(id)->SetUserData(__topo);
                        
                break;
        }
        }
}

//---------------------------------------------------------------------------
MCFace * MCBody::AddMCFace(MG_FACE * __refFace)
{
       MCFace *newMCFace = new MCFace(__refFace);
       // Set the face ID
       _refGeom->ajouter_mg_surface(newMCFace->GetPolySurface());
       _refGeom->ajouter_mg_face(newMCFace);

       return newMCFace;
}

//---------------------------------------------------------------------------
MCFace * MCBody::AddMCFace( MCFace *__a,  MCFace *__b)
{
        std::ostringstream idOriginal;
        MCFace & mcFace1 = *__a;                                     
        MCFace & mcFace2 = *__b;
        MCFace *newMCFace = new MCFace(mcFace1, mcFace2);
       _refGeom->ajouter_mg_surface(newMCFace->GetPolySurface());
       _refGeom->ajouter_mg_face(newMCFace);
        return newMCFace;
}

//---------------------------------------------------------------------------
MCFace* MCBody::GetMCFace(const int __id) const
{
        Graph::Node * node = _G21->GetNode(__id);
        if (node)
                return GetMCFace(node);
        else
                printf("Error in GetMCFace ID = %d\n", __id);
          
        return NULL;
}

//---------------------------------------------------------------------------
bool MCBody::Contains(MCFace * __mcFace) const
{
    return (GetMCFace(__mcFace->get_id()) != NULL);
}

//---------------------------------------------------------------------------
MCFace* MCBody::GetMCFace(Graph::GraphObject * __graphObject)
{
        return (MCFace*) __graphObject->GetUserData();
}
//---------------------------------------------------------------------------
//
void MCBody::AddMCEdgeCovertices(MCEdge * __mcEdge, std::vector<MCVertex *> __mcVertex)
{
    int i;

    // Set the covertex ID and add it in the geometry manager
    MG_COSOMMET * covertex[2];
    for (i=0;i<2;i++)
        covertex[i] = _refGeom->ajouter_mg_cosommet(__mcEdge, __mcVertex[i]);

    if ( covertex[0]->get_t() > covertex[1]->get_t() && ! __mcEdge->GetPolyCurve()->est_periodique() )
        std::swap (covertex[0], covertex[1]);


    __mcEdge->changer_cosommet1(covertex[0]);
    __mcEdge->changer_cosommet2(covertex[1]);
}
//---------------------------------------------------------------------------
void MCBody::DeleteMCEdgeCovertices(MCEdge * __mcEdge)
{
    _refGeom->supprimer_mg_cosommet(__mcEdge->get_cosommet1());  
    _refGeom->supprimer_mg_cosommet(__mcEdge->get_cosommet2());
}
//---------------------------------------------------------------------------
void MCBody::DeleteMCEdge(MCEdge * __mcEdge)
{
        printf("Deleting MC Edge\t : ID \"%lu\" at adress %p\n",__mcEdge->get_id(), __mcEdge);
    _refGeom->supprimer_mg_courbeid(__mcEdge->GetPolyCurve()->get_id());
    _refGeom->supprimer_mg_areteid(__mcEdge->get_id());
}                              
//---------------------------------------------------------------------------
void MCBody::DeleteMCFace(MCFace * __mcFace)
{
        printf("Deleting MC Face\t : ID \"%lu\" at adress %p\n",__mcFace->get_id(), __mcFace);
        _refGeom->supprimer_mg_surfaceid(__mcFace->GetPolySurface()->get_id());
        _refGeom->supprimer_mg_faceid(__mcFace->get_id());
}                              
//---------------------------------------------------------------------------
void MCBody::DeleteMCVertex(MCVertex * __mcVertex)
{
//    printf("Deleting MC Vertex\t : ID \"%lu\" at adress %p\n",__mcVertex->get_id(), __mcVertex);

    _refGeom->supprimer_mg_sommetid(__mcVertex->get_id());
}
//---------------------------------------------------------------------------
MCEdge * MCBody::AddMCEdge(MG_ARETE * __refEdge)
{
        // Construct a mc edge
        MCEdge * newMCEdge = new MCEdge(__refEdge);

        // Set the newEdge ID                
    _refGeom->ajouter_mg_courbe(newMCEdge->GetPolyCurve());

    if (_refGeom->get_mg_courbeid( newMCEdge->GetPolyCurve()->get_id()) == NULL)
            printf("PB!\n");

        _refGeom->ajouter_mg_arete(newMCEdge);
        if (_refGeom->get_mg_areteid( newMCEdge->get_id()) == NULL)
            printf("PB!\n");

        return newMCEdge;
}                                    
//---------------------------------------------------------------------------                                        
MCEdge * MCBody::AddMCEdge(MCEdge *__A, MCEdge *__B)
{                               
    std::ostringstream idOriginal;
    idOriginal<< __A->get_idoriginal();
    idOriginal<<"+";
    idOriginal<< __B->get_idoriginal();

    // Create a polycurve
    MCEdge * newMCEdge = new MCEdge( idOriginal.str(), *__A, *__B); 

    // Set the newEdge ID
    _refGeom->ajouter_mg_courbe(newMCEdge->GetPolyCurve());

    if (_refGeom->get_mg_courbeid( newMCEdge->GetPolyCurve()->get_id()) == NULL)
            printf("PB!\n");

    _refGeom->ajouter_mg_arete(newMCEdge);
        if (_refGeom->get_mg_areteid( newMCEdge->get_id()) == NULL)
            printf("PB!\n");

    return newMCEdge;
}                                 
//---------------------------------------------------------------------------
MCEdge * MCBody::AddMCEdge(PolyCurve *__curve)
{
    std::ostringstream idOriginal;
    idOriginal<< "MCEdge" ;
    for (unsigned i=0; i<__curve->GetRefEdgeCount(); i++)
        idOriginal<<"-"<<__curve->GetRefEdge(i)->get_id();
    MCEdge * newMCEdge = new MCEdge(idOriginal.str(), __curve);
    
    // Set the newEdge ID                                          
    _refGeom->ajouter_mg_courbe(newMCEdge->GetPolyCurve());

    if (_refGeom->get_mg_courbeid( newMCEdge->GetPolyCurve()->get_id()) == NULL)
            printf("PB!\n");

    _refGeom->ajouter_mg_arete(newMCEdge);        
        if (_refGeom->get_mg_areteid( newMCEdge->get_id()) == NULL)
            printf("PB!\n");

    return newMCEdge;
}

//---------------------------------------------------------------------------
Graph::Graph * MCBody::G10() const
{
        return _G10;
}

//---------------------------------------------------------------------------
Graph::Graph * MCBody::G21() const
{
        return _G21;
}

//---------------------------------------------------------------------------
Graph::Graph * MCBody::G20() const
{
        return _G20;
}

//---------------------------------------------------------------------------
MCEdge* MCBody::GetMCEdge(int __id) const
{
        Graph::Node * node = _G10->GetNode(__id);

        if (node)
                return GetMCEdge(node);
        else
                cout << "GetMCEdge : Error: Did not find the MCEdge corresponding to the ID \""<<__id<<"\"\n";

        return NULL;
}

//---------------------------------------------------------------------------
MCEdge* MCBody::GetMCEdge(Graph::GraphObject * __graphObject)
{
        return (MCEdge *) __graphObject->GetUserData();
}

//---------------------------------------------------------------------------
bool MCBody::Contains(MCEdge * __mcEdge) const
{
    Graph::Node * n = _G10->GetNode(__mcEdge->get_id());

    if (!n)
        return false;

    return ( n->GetUserData() == __mcEdge );
}

//---------------------------------------------------------------------------
MCVertex * MCBody::AddMCVertex(MG_SOMMET * __refVertex)
{
        double pos[3]; __refVertex->get_point()->evaluer(pos);
        MCVertex *mcVertex = new MCVertex(__refVertex);

        // Add the vertex to the geometry manager in order to
        // set its ID
        _refGeom->ajouter_mg_sommet(mcVertex);

        return mcVertex;
}

//---------------------------------------------------------------------------
MCVertex* MCBody::GetMCVertex(int __id)const
{
        Graph::Arc * arc;
        arc = _G10->GetArc(__id);
        if (arc)
                return GetMCVertex(arc);
        else
        {
                        printf("Error in GetMCVertex (int id) method : can't find the arc ID %d in G10 Graph !\n", __id);
        }
        return NULL;
}

//---------------------------------------------------------------------------
bool MCBody::Contains(MCVertex * __mcVertex) const
{
    Graph::Arc * a = _G10->GetArc(__mcVertex->get_id());

    if (!a)
        return false;

    return (a->GetUserData() == __mcVertex );
}

//---------------------------------------------------------------------------
MCVertex * MCBody::GetMCVertex(Graph::GraphObject * __graphObject)
{
        return (MCVertex*) __graphObject->GetUserData();
}

//---------------------------------------------------------------------------
void MCBody::InitHyperGraphs()
{
        int j, k;
        _G21 = new Graph::Graph;
        _G10 = new Graph::Graph;
        _G20 = new Graph::Graph;

    std::set < MG_FACE * > lst_ref_faces;
    std::set < MG_ARETE * > lst_ref_edges;
    std::set < MG_SOMMET * > lst_ref_vertices;

    // initialize face, edge, and vertices lists
    unsigned nb_shells = _refBody->get_nb_mg_coquille ();
    for (unsigned it_shells = 0; it_shells < nb_shells; it_shells++)
    {
            MG_COQUILLE * shell = _refBody->get_mg_coquille(it_shells);
            unsigned nb_coface = shell->get_nb_mg_coface();
            for (unsigned it_coface = 0; it_coface < nb_coface; it_coface++)
            {
                    MG_FACE * face = shell->get_mg_coface(it_coface)->get_face();
                    lst_ref_faces.insert(face);

                    unsigned nb_loop = face->get_nb_mg_boucle();
                    for (unsigned it_loop = 0; it_loop < nb_loop; it_loop++)
                    {
                            MG_BOUCLE * loop = face->get_mg_boucle(it_loop);
                            unsigned nb_edge = loop->get_nb_mg_coarete();

                            for (unsigned it_edge = 0; it_edge < nb_edge; it_edge++)
                            {
                                    MG_ARETE * edge = loop->get_mg_coarete(it_edge)->get_arete();
                                    lst_ref_edges.insert(edge);

                                    lst_ref_vertices.insert (edge->get_cosommet1()->get_sommet());
                                    lst_ref_vertices.insert (edge->get_cosommet2()->get_sommet());
                            }
                    }
            }
    }
    

        // Correspondance map between Reference Entities and MC entities
        std::map <MG_FACE *, MCFace *> mapRefToMCFace;
        std::map <MG_ARETE *, MCEdge *> mapRefToMCEdge;
        std::map <MG_SOMMET *, MCVertex *> mapRefToMCVertex;

        // Create the MC Vertices
        for (std::set<MG_SOMMET*>::iterator it_vert = lst_ref_vertices.begin();
                it_vert != lst_ref_vertices.end(); it_vert++)
        {
                MG_SOMMET * refVertex = *it_vert;
                // Create the MC Vertex
                // and initialize the Indentifier using the Geometry manager
                MCVertex * mcVertex = AddMCVertex(refVertex);
                // Map the MC Vertex to its reference vertex in the correspondance map
                mapRefToMCVertex [ refVertex ] = mcVertex;
        }
                                     
        // Create the MC Edges
        for (std::set <MG_ARETE *>::iterator it_edge = lst_ref_edges.begin();
                it_edge != lst_ref_edges.end();
                it_edge ++)                        
        {
                MG_ARETE * refEdge = *it_edge;
                // Create the MC Edge
                // and initialize the Indentifier using the Geometry manager
                MCEdge * mcEdge = AddMCEdge(refEdge);
                // Map the MC Edge to its reference edge in the correspondance map
                mapRefToMCEdge [ refEdge ] = mcEdge;
        }
                                       
        // Create the MC Faces
        for (std::set<MG_FACE*>::iterator it_face = lst_ref_faces.begin();
                it_face != lst_ref_faces.end();
                it_face++)                     
        {
                MG_FACE * refFace = *it_face;
                // Create the MC Face
                // and initialize the Indentifier using the Geometry manager
                MCFace * mcFace = AddMCFace(refFace);
                // Map the MC Face to its reference Face in the correspondance map
                mapRefToMCFace [ refFace ] = mcFace;
        }

            // for each face of the reference topology
        for (std::set<MG_FACE*>::iterator it_face = lst_ref_faces.begin();
                it_face != lst_ref_faces.end();
                it_face++)
        {
                MG_FACE * refFace = *it_face;                    

                // Get the MCFace corresponding to the reference face
                MCFace * mcFace = mapRefToMCFace [ refFace ];

                        printf("Initializing\tMC Face\tID \"%lu\" \t(Ref Face ID \"%lu\")\n",mcFace->get_id(), refFace->get_id());

                        // Add a node in _G21 and _G20
                        _G21->AddNode (mcFace->get_id());
                        _G20->AddNode (mcFace->get_id());

                // initialize graph G21 and G20 nodes to face
                Graph_SetUserData(mcFace);  
        }
                                  
        // For each edge of refBody :
        // create an arc in G21 and connect-it to each nodes of adjacent faces
        // create a node in G1
        for (std::set <MG_ARETE *>::iterator it_edge = lst_ref_edges.begin();
                it_edge != lst_ref_edges.end();
                it_edge ++)
        {
                MG_ARETE * refEdge = *it_edge;

                // Get the MCEdge corresponding to the reference Edge
                MCEdge * mcEdge = mapRefToMCEdge [ refEdge ];

                        // printing the step number
                printf("Initializing\t MC Edge\t%d of %d : ID = \"%lu\" \t(Reference Edge ID \"%lu\")\n", std::distance(lst_ref_edges.begin(), it_edge)+1, lst_ref_edges.size(), mcEdge->get_id(), refEdge->get_id());

                        // create a node in G1

                        // add a node in _G10
                        _G10->AddNode (mcEdge->get_id());

                        // create an arc in G21 and connect-it to each nodes of adjacent faces
                        std::vector <int> nodes;
                for (int j=0; j < refEdge->get_nb_mg_coarete(); j++)
                {
                        MG_FACE * refFace = refEdge->get_mg_coarete(j)->get_boucle()->get_mg_face();
                        // mcFace = adjacent face
                        MCFace * mcFace = mapRefToMCFace [ refFace ];

                                    // connect arc to adjacent faces
                                    nodes.push_back(mcFace->get_id());
                }
                        _G21->AddArc(nodes, mcEdge->get_id());
                                   
                // initialize graph G21 arcs and G10 nodes to MC Edge
                Graph_SetUserData(mcEdge);
        }


    // For each ref vertex of reference Body :
    for (std::set<MG_SOMMET*>::iterator it_vertex = lst_ref_vertices.begin();
    it_vertex != lst_ref_vertices.end();
    it_vertex++)
    {
        MG_SOMMET * refVertex = *it_vertex;

        // Get the MCVertex corresponding to the reference Vertex
        MCVertex * mcVertex = mapRefToMCVertex [ refVertex ];

        // -------------- G20 Arc --------------------------------------------
            // create an arc in G20 and connect-it to each nodes of adjacent faces
        {
        std::vector <int> nodes;
        for (unsigned j = 0; j < refVertex->get_nb_mg_cosommet(); j++)
        {
            MG_ARETE * refEdge = refVertex->get_mg_cosommet(j)->get_arete();

            for (int k=0; k < refEdge->get_nb_mg_coarete(); k++)
            {
                    MG_FACE * refFace = refEdge->get_mg_coarete(k)->get_boucle()->get_mg_face();
                    MCFace * mcFace = mapRefToMCFace [ refFace ];

                    if (std::find (nodes.begin(), nodes.end(), mcFace->get_id()) == nodes.end()) 
                        nodes.push_back(mcFace->get_id());
            }
        }
         // create a graph-arc in G20
         _G20->AddArc(nodes, mcVertex->get_id());
         }

                           
        // -------------- G10 Arc --------------------------------------------  
        // create an arc in G10 and connect-it to each nodes of adjacent edges
        {
        std::vector <int> nodes;
        for (unsigned j = 0; j < refVertex->get_nb_mg_cosommet(); j++)
        {
            MG_ARETE * refEdge = refVertex->get_mg_cosommet(j)->get_arete();
            MCEdge * mcEdge = mapRefToMCEdge [ refEdge ];

            nodes.push_back(mcEdge->get_id());
        }
        // create a graph-arc in G10
        _G10->AddArc(nodes, mcVertex->get_id());
        }

        // initialize graph G10 arc and G20 arc to MC Vertex
        Graph_SetUserData(mcVertex);
    }

    // Add covertices for each edge
    {
    MCBODY_FOR_EACH_MCEDGE(this, mcEdge)
    {
        std::vector <MCVertex*> extremities;
        Edge_GetVertices(mcEdge, extremities);
        AddMCEdgeCovertices(mcEdge, extremities);
    }
    }
}

//---------------------------------------------------------------------------
int
MCBody::G2ToNeatoFile (Graph::Graph * __G2,
                                char * __filename,
                                char *__graphname)
{
    using std::cout;
    std::string graphNeato;
    std::ofstream neatoFile;
    neatoFile.open( __filename );
    graphNeato = NeatoGraph(*__G2);
    std::cout << graphNeato;
    neatoFile << graphNeato;
    neatoFile.flush();
    neatoFile.close(); //explicit

    return 0;
}
//---------------------------------------------------------------------------
// Create a MCVertex with a defined label
MCVertex * MCBody::CreateMCVertex (MG_SOMMET * __refVertex)
{
    // Create the MC Vertex in G20 and G10

    MCVertex * mcVertex = new MCVertex (__refVertex);

    _refGeom->ajouter_mg_sommet(mcVertex);
                                                                                   
    int id = mcVertex->get_id();
    Graph::Arc * G20Arc  = _G20->AddArc(id, 0 /*rank=0*/);
    Graph::Arc * G10Arc  = _G10->AddArc(id, 0 /*rank=0*/);

    Graph_SetUserData(mcVertex);

    return mcVertex;
}

//---------------------------------------------------------------------------
// Split an MC edge by inserting one MC vertex in its domain
void MCBody::SplitEdge(MCEdge * __mcEdge, double __xyz[3], MCTChanges * __mctChanges, RTChanges * __rtChanges)
{
    MG_ARETE *__origRefEdge;
    MG_SOMMET * __splitRefVertex;
    MG_ARETE *__splitRefEdges[2];
    MCVertex * __mcVertex;
    MCEdge *__result[2];
    int i;

    // store the extremities of the original edge in a vector 
    std::vector <MCVertex * > origEdgeExtremities;
    Edge_GetVertices ( __mcEdge, origEdgeExtremities );
                          
    // split curves
    PolyCurve * edgePolyCurve = __mcEdge->GetPolyCurve();
    PolyCurve * curves[2];
    SplitPolyCurve(edgePolyCurve, __xyz, _refBody, _refGeom, curves, &__origRefEdge, &__splitRefVertex, __splitRefEdges);
    if ( __splitRefVertex == NULL )
    {
        printf("Can't split MC Edge in MC body\n");
        return;
    }
    _refGeom->ajouter_mg_sommet(__splitRefVertex);

    __mcVertex = CreateMCVertex(__splitRefVertex);

    printf("Inserting MC Vertex\t : ID \"%lu\" \n", __mcVertex->get_id());
    
    // Create the two MC Edges
    for ( i=0; i<2; i++)
    {
        __result[i] = AddMCEdge(curves[i]);
        __result[i]->transfert_ccf(*__mcEdge);
    }

    // duplicate G21A to get 2 parallel arcs
    Graph::Arc * G21A[3];
    G21A[2] = _G21->GetArc(__mcEdge->get_id());
    int faceId[2];

    // save the G21 nodes adjacent to the mc Edges's arc 
    Graph::Node * G21N[2];
    i=0;
    for (Graph::Arc::MultimapNodesById::iterator itG21Node = G21A[2]->Nodes().begin();
        itG21Node != G21A[2]->Nodes().end();
        itG21Node++ )
    {
        faceId[i] = itG21Node->first;
        G21N[i++] = itG21Node->second;
    }


    // Set edges id
    int edgeId[3];
    edgeId[2] = __mcEdge->get_id();
    for (i = 0; i < 2; i++)
        edgeId[i] = __result[i]->get_id();

    // Duplicate the G21 arc of the split edge
    for (i = 0; i < 2; i++)
        G21A[i] = _G21->DuplicateArc(G21A[2], edgeId[i]);
    // Remove the G21 Arc of the split edge
    _G21->RemoveArc(G21A[2]);

    // split the G10 node of the split edge
    Graph::Arc * G10A[3];
    Graph::Node * G10N[3];
    G10A[2] = _G10->GetArc( __mcVertex->get_id() ); // the arc representing the split vertex
    G10N[2] = _G10->GetNode( __mcEdge->get_id() );
    // delete the G10 node of the old edge
    _G10->RemoveNode(G10N[2]);
    // Add two new nodes 
    for ( i=0; i<2; i++)
        G10N[i] = _G10->AddNode(edgeId[i]);
    // Add an arc linking the 2 nodes. This arc represents the split vertex
    // Link the G10 arc of the MC vertex to the 2 nodes representing the edges
    for (i=0; i<2; i++)
    {
        G10A[2]->Add(G10N[i]);
        G10N[i]->Add(G10A[2]);
    }
    
    for (i=0; i<2; i++)
    {
        MG_SOMMET * refVertex = NULL;
        if (curves[i]->get_sommet1()->get_id() != __mcVertex->GetRefVertex()->get_id())
            refVertex = curves[i]->get_sommet1();
        else
            refVertex = curves[i]->get_sommet2();
        MCVertex * V;
        
        for (int j=0; j<2; j++)
            if (origEdgeExtremities[j]->Contains( refVertex ))
                V = origEdgeExtremities[j];

        G10A[i] = _G10->GetArc(V->get_id());
    }


    // add new node[i] in arc[i]
    for ( i=0; i<2; i++)
    {
        G10A[i]->Add(G10N[i]);
        G10N[i]->Add(G10A[i]);
    }

    // Add an arc in G20
    Graph::Arc * G20Arc = _G20->GetArc(__mcVertex->get_id());
    Graph::Node * G20Node[2];
    for ( i=0; i<2; i++)
    {
        G20Node[i] = _G20->GetNode(faceId[i]);
        G20Arc->Add(G20Node[i]);
        G20Node[i]->Add(G20Arc);
    }

    printf("G10A[2]->Rank = %d\n", G10A[2]->Rank());
    printf("G20A[2]->Rank = %d\n", _G20->GetArc(__mcVertex->get_id())->Rank());

    /// Operations using the resulting Graph of this MC Edge
    /// split :
    /// Get the vertex extremities of both new MC Edge
    /// and Add 2 Covertices      
    for (i=0; i<2; i++)
    {
        MCEdge * newMCEdge = __result[i];

        std::vector<MCVertex*> mcVertexExtremities;
        Edge_GetVertices(newMCEdge, mcVertexExtremities);
        AddMCEdgeCovertices(newMCEdge, mcVertexExtremities);
        
        /// Set graph user Data of new MCEdge G21 Arc and G10 Node
        Graph_SetUserData(newMCEdge);
    }


    // Deleted and Created MC Topology
    if (__mctChanges)
    {
        // a new vertex has been created
        __mctChanges->V.insert( std::make_pair ((MCVertex*)NULL, __mcVertex));
        // __mcEdge has been replaces by 2 new edges                        
        __mctChanges->E.insert( std::make_pair (__mcEdge, (MCEdge*) NULL) );
        __mctChanges->E.insert( std::make_pair ((MCEdge*) NULL, __result[0]) );
        __mctChanges->E.insert( std::make_pair ((MCEdge*) NULL, __result[1]) );
    }

    // Deleted and created Reference Topology 
    if (__rtChanges)
    {
        // a new vertex __splitRefVertex has been created
        __rtChanges->V.insert( std::make_pair ((MCVertex*)NULL , __splitRefVertex) );
        // __origRefEdge has been replaces by 2 new edges __splitRefEdges
        __rtChanges->E.insert( std::make_pair (__origRefEdge , __splitRefEdges[0]) );
        __rtChanges->E.insert( std::make_pair (__origRefEdge , __splitRefEdges[1]) );
    }

}

//---------------------------------------------------------------------------
int
MCBody::G1ToNeatoFile (Graph::Graph * __G1,
                                char * __filename,
                                char *__graphname)
{
    using std::cout;
    std::ofstream neatoFile;
    neatoFile.open( __filename );
    neatoFile << "graph \""<< __graphname <<"\" {" << std::endl;
    neatoFile << "node [ shape=box , color = blue ]\n";

    //////// G1 graph
    GRAPH_FOR_EACH_ARC_CONST(__G1,itarcs)
        {
            neatoFile << "\""<<itarcs->first<<"\" [ shape=point, style=filled, fillcolor=red ]\n";

                ARC_FOR_EACH_NODE_CONST( itarcs->second, itnodes)
                {
                    neatoFile << "\""<< itnodes->first << "\"";
                    neatoFile << " -- ";
                    neatoFile << "\""<< itarcs->first << "\"";
                    neatoFile << " [ len=1  ] ";
                    neatoFile << std::endl;
                }
        }
    neatoFile << "}" << std::endl;
    neatoFile << std::endl;

    neatoFile.close(); //explicit

    return 0;
}

bool V2(MCBody * _mcBody, 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*) _mcBody->GetGeometry()->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 V(MCBody * _mcBody)
{
        int i=0;
    {MCBODY_FOR_EACH_MCEDGE(_mcBody, edge22){
               if (_mcBody->G21()->GetArcs().size() <= 1)
                   break;
               if (V2(_mcBody, edge22) == false)
                   printf("MC Edge ID=%lu not OK\n", edge22->get_id());
               i++;
           }
    }

}
//---------------------------------------------------------------------------
void
MCBody::SuppressMCEdge(MCEdge * __mcEdge, MCFace * __F[3])
{
        int i;
        int __id = __mcEdge->get_id();

        // Find the arc __id in G2V (G2 Virtual)
        Graph::Arc * arcG21 = _G21->GetArc(__id);
          //        V(this);

        ///// Process G2 Graph
        printf("Deleting MC Edge\t : ID \"%d\" at adress %p \n", __id, __mcEdge);

        if ( ! arcG21->IsLoop()  ) // The arc represents an boundary edge
        {
                // Get the two mc faces adjacent to the suppressed mc edge
                MCFace * face[3];
                int faceId[3];
                Graph::Node * faceNode[3];

                i = 0;
                Graph::Arc::MultimapNodesById::iterator itNode = arcG21->Nodes().begin();
                for (i=0; i<2; i++)
                {
                        faceId[i] = itNode->first;
                        faceNode[i] = itNode->second;
                        face[i] = GetMCFace(faceId[i]);

                        itNode++;
                }
                
                  //V(this);
                _G21->RemoveArc(__id);
                  //V(this);

                // Merge the nodes of these 2 mc face in G21
                // (contract the arc in G21)  
                MCFace* newMCFace = AddMCFace(face[0], face[1]);
                faceId[2] = newMCFace->get_id();
                faceNode[2] = _G21->MergeNodes(faceNode[0], faceNode[1], faceId[2]);

                // Merge the nodes of these 2 mc faces in G20    
                Graph::Node * mergedG20Node = _G20->MergeNodes( faceId[0], faceId[1], faceId[2]);
                  //V(this);

                // Remove the node in G10
                _G10->RemoveNode(__mcEdge->get_id());
                  //V(this);

                // Attach the MC Face to Graph nodes of G21 and G10
                Graph_SetUserData(newMCFace);
                  //V(this);

                printf ("Adding new merged face %d = { %d + %d }\n", faceId[2], faceId[0], faceId[1]);

                __F[0] =  face[0];
                __F[1] =  face[1];
                __F[2] =  newMCFace;
        }
        else // The arc represents an interior edge
        {
                Graph::Node *n = _G21->GetArc(__id)->First();
                
                // if the arc is a loop then delete it
                // in G2V : delete the g2v arc from the node and from the graph
                _G21->RemoveArc(__id);
                  //V(this);
                // in G1V : delete the g1v node from the graph
                _G10->RemoveNode(__id);
                  //V(this);
                                    
                __F[0] = 0;
                __F[1] = 0;
                __F[2] = GetMCFace(n);
                  //V(this);
        }
}
             
//---------------------------------------------------------------------------
void
MCBody::SuppressMCVertex(MCVertex * __mcVertex, MCTChanges * __mctChanges)
{                   
        printf("Deleting MC Vertex\t : ID \"%lu\" at adress %p\n",__mcVertex->get_id(), __mcVertex);

        // Find the arc __name in G1V (G1 Virtual)
        Graph::Arc *g1vArc = _G10->GetArc( __mcVertex->get_id() );
        MCVertex * mcVertex = __mcVertex;

        // Case of a B-Rep vertex separating only 2 edges
        if ( g1vArc->Rank() == 2 )
        {
                // Get the 2 MC Edges separated by the MC Vertex
                std::vector <MCEdge *> edge;
                Vertex_GetAdjacentEdges( mcVertex, edge );

                // Get the MC Faces Separated by the 2 MC Edges
                std::set <MCFace *> listMCFaces = Edge_GetAdjacentFaces(edge[0]);
                MCFace * mcFace[2];
                std::set<MCFace*>::iterator itr=listMCFaces.begin();
                for (unsigned i=0;i<2;i++)
                {
                        mcFace[i] = *itr;
                        itr++;
                }
                                         
                // Unify parallel G21 arcs

                // Create a unique arc in G2 Which connect n1 and n2
                Graph::Arc::MultimapNodesById newG21ArcNodes = _G21->GetArc( edge[0]->get_id() )->Nodes();      
                _G21->RemoveArc( edge[0]->get_id() );
                _G21->RemoveArc( edge[1]->get_id() );
                                                                  
                printf("the deleted vertex was adjacent to %d / %d faces\n",  _G20->GetArc( __mcVertex->get_id() )->Rank(), listMCFaces.size());
                printf("the deleted vertex was adjacent to %d edges\n",  _G10->GetArc( __mcVertex->get_id() )->Rank());
                // attach to the new node a edge attribute
                // constructed by the merge of the 2 edge attributes
                MCEdge * newMCEdge = AddMCEdge(edge[0], edge[1]);

                _G21->AddArc(newMCEdge->get_id(), newG21ArcNodes);

                // Merge nodes G10
                Graph::Node* n = _G10->MergeNodes(_G10->GetNode(edge[0]->get_id()), _G10->GetNode(edge[1]->get_id()), newMCEdge->get_id());
                // Remove G10 loop-arc to newNode
                _G10->RemoveArc( _G10->GetArc(__mcVertex->get_id()) );
                // Remove G20 arc
                _G20->RemoveArc( _G20->GetArc(__mcVertex->get_id()));

                /// Operations using the resulting Graph of this MC Vertex
                /// Suppression :
                /// Get the vertex extremities of this new MC Edge
                /// and Add 2 Covertices
                std::vector<MCVertex*> mcVertexExtremities;
                Edge_GetVertices(newMCEdge, mcVertexExtremities); 
                PolyCurve * newPolyCurve = newMCEdge->GetPolyCurve();

                if (newPolyCurve->est_periodique())
                    if (newPolyCurve->get_sommet1() != mcVertexExtremities[0]->GetRefVertex())
                        newPolyCurve->SetPeriodicPoleRefVertex(mcVertexExtremities[0]->GetRefVertex());

                AddMCEdgeCovertices(newMCEdge, mcVertexExtremities);
                /// Set graph user Data of new MCEdge G21 Arc and G10 Node
                Graph_SetUserData(newMCEdge);

                if (__mctChanges)
                {
                    // The mc vertex has been deleted
                    __mctChanges->V.insert( std::make_pair (__mcVertex, (MCVertex*) NULL) );
                    // The 2 adjacent edges have been merged together
                    __mctChanges->E.insert( std::make_pair (edge[0], newMCEdge) );
                    __mctChanges->E.insert( std::make_pair (edge[1], newMCEdge) );
                }
        }
        else if (g1vArc->Rank() == 0)
        {      // case of a vertex isolated in the domain of a face
               _G10->RemoveArc(_G10->GetArc(__mcVertex->get_id()));
               _G20->RemoveArc(_G20->GetArc(__mcVertex->get_id()));

                if (__mctChanges)
                {
                    // The mc vertex has been deleted
                    __mctChanges->V.insert( std::make_pair (__mcVertex, (MCVertex*) NULL) );
                }
        }
        else
        {
            printf("Warning : MCBody can't delete vertex %d because its rank equals %d\n", g1vArc->Id(), g1vArc->Rank());
        }
}

//---------------------------------------------------------------------------
void MCBody::MergeVertices(MCVertex * __deleteVertex, MCVertex * __targetVertex, std::set<MCEdge*> & edges)
{
    int V1_id, V2_id;
    MCVertex * V1, * V2;
    Graph::Arc * V1_G10Arc, *V2_G10Arc, *V1_G20Arc, *V2_G20Arc;
    V1=__deleteVertex; V1_id=V1->get_id();
    V2=__targetVertex; V2_id=V2->get_id();
    V1_G10Arc=_G10->GetArc(V1_id);
    V1_G20Arc=_G20->GetArc(V1_id);
    V2_G10Arc=_G10->GetArc(V2_id);
    V2_G20Arc=_G20->GetArc(V2_id);
    Graph::Arc::MultimapNodesById V1_G10Arc_MapNodes = V1_G10Arc->Nodes();
    Graph::Arc::MultimapNodesById V1_G20Arc_MapNodes = V1_G20Arc->Nodes();
    Graph::Arc::MultimapNodesById::iterator itNode;

    // Initialize edge list
    for (itNode = V1_G10Arc_MapNodes.begin(); itNode != V1_G10Arc_MapNodes.end(); itNode++)
    {
        Graph::Node * n = itNode->second;
        MCEdge * e = GetMCEdge(n);
        edges.insert(e);
    }
    // Replace V1 by V2 in G20 and G10
    for (itNode = V1_G10Arc_MapNodes.begin(); itNode != V1_G10Arc_MapNodes.end(); itNode++)
    {
        Graph::Node * n = itNode->second;
        V2_G10Arc->Add(n);
        n->Add(V2_G10Arc);
    }
    _G10->RemoveArc(V1_G10Arc);
    for (itNode = V1_G20Arc_MapNodes.begin(); itNode != V1_G20Arc_MapNodes.end(); itNode++)
    {
        Graph::Node * n = itNode->second;
        V2_G20Arc->Add(n);
        n->Add(V2_G20Arc);
    }
    _G20->RemoveArc(V1_G20Arc);

    // replace each covertices to V1 by covertex to V2
    int i,N;
    N=V1->get_nb_mg_cosommet();
    for (i=0; i<N; i++)
    {
        MG_COSOMMET * cov = V1->get_mg_cosommet(i);
        MCEdge * e = (MCEdge*) cov->get_arete();

        cov = _refGeom->ajouter_mg_cosommet(e,V2);

        if (e->get_cosommet1()->get_sommet()==V1)
            e->changer_cosommet1(cov);
        else if (e->get_cosommet2()->get_sommet()==V1)
            e->changer_cosommet2(cov);
    }

    V2->Merge(V1);
}
//---------------------------------------------------------------------------
/**
* Edge collapse operator
* arguments : name of the edge to be contracted, name of the vertex to contract the edge on
*
*/
void MCBody::ContractEdgeToVertex(MCEdge * __mcEdge, MCVertex * __targetVertex, std::set <MCEdge*> & newEdges, std::set<MCEdge*> & oldEdges)
{
    // E3 = mcedge V3 = target vertex V4 = delete vertex
        printf("Contracting\tMC Edge \"%lu\" to MC Vertex \"%lu\"\n", __mcEdge->get_id(), __targetVertex->get_id());

    MCVertex * V3 = __targetVertex, * V4;
    MCEdge * E3 = __mcEdge;
    int V4_id, V3_id, E3_id;
    Graph::Arc * V4_G10Arc, *V3_G10Arc, * V4_G20Arc, *V3_G20Arc, * E3_G21Arc;
    Graph::Node * E3_G10Node;

    V3_id = V3->get_id();
    V3_G10Arc = _G10->GetArc(V3_id);
    V3_G20Arc = _G20->GetArc(V3_id);

    E3_id = E3->get_id();
    E3_G10Node = _G10->GetNode(E3_id);
    E3_G21Arc = _G21->GetArc(E3_id);

    // Get G1 arc V4 the second incident arc to E3 different to V3
    for (Graph::Node::MultimapArcsById::const_iterator incidentArcIter = E3_G10Node->IncidentArcs().begin();
              incidentArcIter != E3_G10Node->IncidentArcs().end();
              incidentArcIter++)
              if (incidentArcIter->first != V3_id)
                        V4_G10Arc = incidentArcIter->second;
    V4 = GetMCVertex(V4_G10Arc);
    V4_id = V4->get_id();
    V4_G20Arc = _G20->GetArc(V4_id);

    // Merge E3 with edges that are adjacent to V4
    // get the nodes of the arc
        Graph::Arc::MultimapNodesById initialV4Nodes = V4_G10Arc->Nodes();
        for (Graph::Arc::MultimapNodesById::const_iterator itNode = initialV4Nodes.begin();
             itNode != initialV4Nodes.end();
             itNode++)
             {
                int id = itNode->first;

                                if (V4_G10Arc->IsLoopOfNode(id) || V3_G10Arc->IsLoopOfNode(E3_id)|| V4_G10Arc->IsLoopOfNode(E3_id) )
                                {
                                        printf("Warning while Contracting\tMC Edge \"%d\" to MC Vertex \"%d\" \nArc %d is a loop of node %d\n", 
                                                V4_id, E3_id, V3_id, id);
                        
                                        return;
                                }
                        }

            
        for (Graph::Arc::MultimapNodesById::const_iterator itNode = initialV4Nodes.begin();
             itNode != initialV4Nodes.end();
             itNode++)
        {
            oldEdges.insert(GetMCEdge(itNode->second));
        }

        for (Graph::Arc::MultimapNodesById::const_iterator itNode = initialV4Nodes.begin();
             itNode != initialV4Nodes.end();
             itNode++)
             {
                int id = itNode->first;
                Graph::Node * node = itNode->second;

                if (id != E3_id)
                {

                        // attach to the new node a MC edge attribute
                        // constructed by the merge of the 2 edge attributes
                        MCEdge * newEdge = AddMCEdge(GetMCEdge(id), E3);
                        int newId = newEdge->get_id();

                        Graph::Node * newG10Node = _G10->MergeNodes(node,  E3_G10Node, newId, true);
                                                printf("Edge Contraction : merging edges %d and %d in new edge %d\n", node->Id(), E3_id, newId);

                        // Create an edge in G2 named newVEdgeName and
                        // conecting nodes n1_g2nodes
                        Graph::Arc* newG21Arc = _G21->AddArc(newId, _G21->GetArc(id)->Nodes());

                        newEdges.insert(newEdge);
                }
             }
        for (Graph::Arc::MultimapNodesById::const_iterator itNode = V4_G20Arc->Nodes().begin();
             itNode != V4_G20Arc->Nodes().end();
             itNode++)
             {
                Graph::Node * n = itNode->second;
                int idN = itNode->first;
                if (E3_G21Arc->Nodes().find(idN) == E3_G21Arc->Nodes().end())
                {
                        V3_G20Arc->Add(n);
                        n->Add(V3_G20Arc);
                }
             }
        for (Graph::Arc::MultimapNodesById::const_iterator itNode = initialV4Nodes.begin();
             itNode != initialV4Nodes.end();
             itNode++)
        {
            // in G1 : remove the old edges adjacent to v4
            _G10->RemoveNode(itNode->first);

            // In G2 :
            // remove the old edges adjacent to v4
            _G21->RemoveArc(itNode->first);
        }

        _G20->RemoveArc(V4_id);
        _G10->RemoveArc(V4_id);

        for (std::set<MCEdge*>::iterator itEdge=newEdges.begin();itEdge!=newEdges.end();itEdge++)
        {
            MCEdge * newEdge = *itEdge;
            Graph_SetUserData(newEdge);
            std::vector<MCVertex*> mcVertexExtremities;
            Edge_GetVertices(newEdge, mcVertexExtremities);
            AddMCEdgeCovertices(newEdge, mcVertexExtremities);
        }

        V3->Merge(V4);
}

//---------------------------------------------------------------------------
void MCBody::Edge_GetVertices (MCEdge * __mcEdge, std::vector<MCVertex*> & __mcVertices)
{
        // Get the incident G10 arcs to G10 node
        // they compose the loops of the MC face

        Graph::Node::MultimapArcsById & arcsOfG10Node = _G10->GetNode( __mcEdge->get_id() )->IncidentArcs();

        // Get each MC edge from incident arcs
        for ( Graph::Node::MultimapArcsById::const_iterator it = arcsOfG10Node.begin();
                it != arcsOfG10Node.end(); it++)
        {
                __mcVertices.push_back(GetMCVertex (it->second));
        }
}

//---------------------------------------------------------------------------
std::set < MCFace * >
MCBody::Edge_GetAdjacentFaces(MCEdge * __mcEdge)
{
        std::set < MCFace * > result;

        // Get the G21 nodes to G21 arc
        // they compose the MC faces adjacent to the MC edge
        Graph::Arc *G21Arc = _G21->GetArc(__mcEdge->get_id());

        Graph::Arc::MultimapNodesById & nodesOfG21Arc = G21Arc->Nodes();


        for (Graph::Arc::MultimapNodesById::const_iterator it = nodesOfG21Arc.begin();
                it != nodesOfG21Arc.end();
                it++)
        {               

                MCFace *mcFace = GetMCFace(it->second);
                result.insert( mcFace );

        }

        return result;
}


//---------------------------------------------------------------------------
std::set < MCEdge * >
MCBody::Face_GetAdjacentEdges(MCFace * __mcFace)
{
        std::set < MCEdge * > result;

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

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


        for (Graph::Node::MultimapArcsById::const_iterator it = arcsOfG21Node.begin();
                it != arcsOfG21Node.end();
                it++)
        {               

                result.insert( GetMCEdge(it->second) );

        }

        return result;        
}

std::vector < std::vector < MCEdge * > >
MCBody::Face_GetUnsortedLoops(MCFace * __mcFace)
{
    std::set < std::set < Graph::Node * > > loops;
    std::vector < std::vector < MCEdge * > > result;

    // Extract G 10 subgraph of this face
    Graph::Graph * G10 = Face_G10SubGraph(__mcFace);

    // Extract Strongly Connected Components
    // The SCC represent the BRep loop edges !
    std::set < std::set < Graph::Node * > > scc;
    Graph::FindSCC(G10,scc);

    for (std::set < std::set < Graph::Node * > >::iterator it1 = scc.begin();
        it1 != scc.end();
        it1++
    )
    {
        std::vector <MCEdge * > loop;
        const std::set < Graph::Node * > & loopNodes = *it1;
        for ( std::set < Graph::Node * >::const_iterator it2 = loopNodes.begin();
            it2 != loopNodes.end();
            it2++)
        {
            Graph::Node * node = *it2;
            MCEdge * mcEdge = GetMCEdge( node );
            loop.push_back(mcEdge);
        }
        result.push_back(loop);
    }

    delete G10;
    return result;
}

std::vector < std::vector < MCEdge * > >
MCBody::Face_GetCycles(MCFace * __mcFace)
{
    std::set < std::set < Graph::Node * > > loops;
    std::vector < std::vector < MCEdge * > > result;

    // Extract G 10 subgraph of this face
    Graph::Graph * G10 = Face_G10SubGraph(__mcFace);

    // Extract Graph cycles
    // The graph cycles represent loops which toopology is a circle !
    std::vector < std::vector < Graph::Node * > > cycles;
    Graph::FindCycles(G10,cycles);

    for (std::vector < std::vector < Graph::Node * > >::iterator it1 = cycles.begin();
        it1 != cycles.end();
        it1++
    )
    {
        std::vector <MCEdge * > mcEdgeLoop;
        std::vector < Graph::Node * > & cycleNodes = *it1;
        for ( std::vector < Graph::Node * >::iterator it2 = cycleNodes.begin();
            it2 != cycleNodes.end();
            it2++)
        {
            Graph::Node * node = *it2;
            MCEdge * mcEdge = GetMCEdge( node );
            mcEdgeLoop.push_back(mcEdge);
        }
        result.push_back(mcEdgeLoop);
    }

    for (Graph::Graph::MapNodesById::const_iterator itNode = G10->GetNodes().begin();
        itNode != G10->GetNodes().end();
        itNode++)
    {
        Graph::Node * node = itNode->second;
        for (Graph::Node::MultimapArcsById::const_iterator itArc = node->IncidentArcs().begin();
            itArc != node->IncidentArcs().end();
            itArc ++)
        {
            Graph::Arc * arc = itArc->second;
            if ( arc->IsLoopOfNode(node->Id()) )
            {
                std::vector < MCEdge * > loop;
                MCEdge * mcEdge = GetMCEdge(node);
                loop.push_back(mcEdge);
                result.push_back(loop);
            }
        }
    }

    delete G10;
    return result;
}
//---------------------------------------------------------------------------
void MCBody::Face_GetMCVertices (MCFace * __mcFace, std::vector<MCVertex*> & __mcVertices) 
{
        // Get the incident G20 arcs to G20 node
        const Graph::Node::MultimapArcsById & arcsOfG20Node = _G20->GetNode(__mcFace->get_id())->IncidentArcs();

        // Get each MC edge from incident arcs
        for ( Graph::Node::MultimapArcsById::const_iterator it = arcsOfG20Node.begin();
                it != arcsOfG20Node.end(); it++)
        {
                __mcVertices.push_back(GetMCVertex (it->second));
        }
}

//---------------------------------------------------------------------------
void MCBody::Vertex_GetAdjacentEdges(MCVertex * __mcVertex, std::vector<MCEdge *> & __list)
{
        __list.clear();

        Graph::Arc * a = _G10->GetArc(__mcVertex->get_id());
        for (Graph::Arc::MultimapNodesById::const_iterator it = a->Nodes().begin();
                it != a->Nodes().end();
                it++)
        {
            MCEdge * mcEdge = GetMCEdge(it->second);
            __list.push_back(mcEdge);
        }
}

//---------------------------------------------------------------------------
void MCBody::Vertex_GetAdjacentEdges(MCVertex * __mcVertex, std::set<MCEdge *> & __list)
{
        __list.clear();

        Graph::Arc * a = _G10->GetArc( __mcVertex->get_id() );
        for (Graph::Arc::MultimapNodesById::const_iterator it = a->Nodes().begin();
                it != a->Nodes().end();
                it++)
        {
            MCEdge  * mcEdge = GetMCEdge(it->second);
            __list.insert(mcEdge);
        }
}


//---------------------------------------------------------------------------
void MCBody::Vertex_GetAdjacentFaces(MCVertex * __mcVertex, std::set<MCFace *> & __list)
{
        __list.clear();

        Graph::Arc * a = _G20->GetArc(__mcVertex->get_id());
        for (Graph::Arc::MultimapNodesById::const_iterator it = a->Nodes().begin();
                it != a->Nodes().end();
                it++)
        {
            MCFace * mcFace = GetMCFace(it->second);
            __list.insert(mcFace);
        }
}
//---------------------------------------------------------------------------
Graph::Graph * MCBody::Face_G10SubGraph(MCFace * __mcFace)
{
    std::set < Graph::Node * > nodes;
    Graph::Graph * result;
    
        // Get the G21 arcs to G21 node
        // they compose the MC edges adjacent to the MC face
        Graph::Node *G21Node = _G21->GetNode( __mcFace->get_id() );

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


        for (Graph::Node::MultimapArcsById::const_iterator it = arcsOfG21Node.begin();
                it != arcsOfG21Node.end();
                it++)
        {
                MCEdge * edge = GetMCEdge(it->second);
        Graph::Node * g10Node = _G10->GetNode(edge->get_id());
        nodes.insert(g10Node);
        }

    result = new Graph::Graph(*_G10,nodes);

    return result;
}
                                                                                           
//---------------------------------------------------------------------------
MG_COARETE * MCBody::CreateCoEdge(MCFace * mcFace, MG_BOUCLE * mcLoop, MCEdge * mcEdge, int sense)
{
    MG_COARETE * mcCoEdge = _refGeom->ajouter_mg_coarete (mcEdge, mcLoop, sense);
    mcLoop->ajouter_mg_coarete(mcCoEdge);
    return mcCoEdge;
}                               
//---------------------------------------------------------------------------
int MCBody::ExportBRep_MCEdgeSense(MCFace * mcFace, MCEdge * mcEdge)
{
                    
                // for each shell of the ref body
                unsigned nb_shells = _refBody->get_nb_mg_coquille ();
                for (unsigned it_shells = 0; it_shells < nb_shells; it_shells++)
                {
                        MG_COQUILLE * shell = _refBody->get_mg_coquille(it_shells);
                // for each coface->face F of the shell
                        unsigned nb_coface = shell->get_nb_mg_coface();
                        for (unsigned it_coface = 0; it_coface < nb_coface; it_coface++)
                        {
                                MG_FACE * face = shell->get_mg_coface(it_coface)->get_face();
                                             
                                if ( ! mcFace->GetPolySurface()->Contains(face) )
                                    continue;

                // for each loop L of this face
                                unsigned nb_loop = face->get_nb_mg_boucle();
                                for (unsigned it_loop = 0; it_loop < nb_loop; it_loop++)
                                {
                                        MG_BOUCLE * loop = face->get_mg_boucle(it_loop);
                // for each coedge coe of L
                                        unsigned nb_edge = loop->get_nb_mg_coarete();
                                        for (unsigned it_edge = 0; it_edge < nb_edge; it_edge++)
                                        {
                                                MG_COARETE * coedge = loop->get_mg_coarete(it_edge);
                                                MG_ARETE * edge = coedge->get_arete();

                                                int iRefEdge = -1;
                                                for (int itRefEdge = 0; itRefEdge != mcEdge->GetPolyCurve()->GetRefEdgeCount(); itRefEdge++)
                                                {

                                                   if (mcEdge->GetPolyCurve()->GetRefEdge(itRefEdge) == edge)
                                                   {
                                                        // Do not take reference edge if it is interior to the MC Face
                                                        // (e.g. boundary edge contracted along an interior edge)
                                                        int nbInteriorFaces = 0;
                                                        for (int it_face2 = 0; it_face2 < edge->get_nb_mg_coarete(); it_face2++)
                                                        {
                                                            MG_FACE * face2 = edge->get_mg_coarete(it_face2)->get_boucle()->get_mg_face();
                                                            if (mcFace->GetPolySurface()->Contains(face2))
                                                                nbInteriorFaces ++;
                                                        }
                                                        if (nbInteriorFaces  == 1)
                                                            iRefEdge = itRefEdge;
                                                    }
                                                }

                                                if ( iRefEdge != -1 )
                                                {
                                                    MG_SOMMET * polycurveStartRefVertex = mcEdge->GetPolyCurve()->GetRefVertex(iRefEdge);
                                                    MG_ARETE * polycurveStartRefEdge = mcEdge->GetPolyCurve()->GetRefEdge(iRefEdge);

                                                    // Get the start vertex of this coedge    
                                                    // coedgeStartVertex = start ref vertex of coe1
                                                    MG_SOMMET * coedgeStartVertex;
                                                    if ( coedge->get_orientation() == 1)
                                                        coedgeStartVertex = edge->get_cosommet1()->get_sommet();
                                                    else                                                          
                                                        coedgeStartVertex = edge->get_cosommet2()->get_sommet();

                                                    //            
                                                    // sense = ( coedgeStartVertex != v1 ) ? -1 : 1 ;
                                                    int sense;
                                                    if ( polycurveStartRefVertex != edge->get_cosommet1()->get_sommet() )
                                                        sense = -coedge->get_orientation();
                                                    else
                                                        sense = coedge->get_orientation();

                                                    /*OT_VECTEUR_3D tangentPc, tangentRefEdge;
                                                    PolyCurve * pc = mcEdge->GetPolyCurve();
                                                    MG_SOMMET * polycurveStartRefVertex = pc->GetRefVertex(iRefEdge);
                                                    MG_SOMMET * polycurveEndRefVertex = pc->GetRefVertex(iRefEdge+1);
                                                    double sMid = .5*(pc->RefVertex_GetS(polycurveStartRefVertex)+pc->RefVertex_GetS(polycurveEndRefVertex));
                                                    pc->deriver(sMid, tangentPc);
                                                    double tMidRefEdge, dtMidRefEdge; MG_ARETE * tmpRefEdge;
                                                    pc->Parameter_SToRefEdgeT(sMid, &tmpRefEdge, &tMidRefEdge, &dtMidRefEdge, false);
                                                    tmpRefEdge->deriver(tMidRefEdge, tangentRefEdge);

                                                    double fSense = tangentRefEdge * tangentPc;
                                                    int sense;
                                                    if ( fSense < 0 )
                                                        sense = -coedge->get_orientation();
                                                    else
                                                        sense = coedge->get_orientation(); */

                                                                            return sense;
                                                }
                                        } // end for coedge
                                } // end for face->loop
                        } // end for coface
                } // end shell

                return -10000;
}
MG_VOLUME * MCBody::ExportBRep()
{
    int i,j,k;

    std::ostringstream idBRepVolume;
    idBRepVolume << "MCBody of reference BRep "<< _refBody->get_idoriginal();
    _brep = new MG_VOLUME (idBRepVolume.str());
    _refGeom->ajouter_mg_volume(_brep);

    MG_COQUILLE * mcShell = new MG_COQUILLE(_brep);
    _refGeom->ajouter_mg_coquille(mcShell, _brep);

    {MCBODY_FOR_EACH_MCFACE(this, mcFace)
    {
        _refGeom->ajouter_mg_face(mcFace);        

        MG_COFACE * mcCoFace = _refGeom->ajouter_mg_coface(mcFace, mcShell, 1);          
        mcShell->ajouter_mg_coface(mcCoFace);  
        MG_SURFACE * surface = mcFace->get_surface();
        _refGeom->ajouter_mg_surface(surface);

        std::set <MCEdge*> mcEdgeArray = Face_GetAdjacentEdges(mcFace);
        std::vector <MG_FACE*> faces;std::vector <MG_ARETE*> edges;std::vector <int> edgesense;
        for (std::set <MCEdge*>::iterator itEdge = mcEdgeArray.begin();
            itEdge != mcEdgeArray.end(); itEdge++)
        {
            MCEdge * mcEdge = *itEdge;

                bool boundaryEdge = (!_G21->GetArc(mcEdge->get_id())->IsLoop());

            // If the MC Edge is a boundary, then create an oriented coedge
            if (boundaryEdge)
            {
                faces.push_back(mcFace);
                edges.push_back(mcEdge);
                int sense = ExportBRep_MCEdgeSense(mcFace,mcEdge);
                edgesense.push_back(sense);
            }
            // If the MC Edge is interior, then create two coedges
            // with same sense and inverse sense
            else
            {
                faces.push_back(mcFace); faces.push_back(mcFace);
                edges.push_back(mcEdge);edges.push_back(mcEdge);
                edgesense.push_back(+1);edgesense.push_back(-1);
            }
        }     
        std::vector < std::vector<CAD4FE::MakeLoops::CoEdge*> > loops;
        CAD4FE::MakeLoops makeLoops(faces,edges,edgesense);
        makeLoops.GetFaceLoops(mcFace,loops);
        for (unsigned k=0;k<loops.size();k++)
        {
            MG_BOUCLE * mcLoop = new MG_BOUCLE(mcFace);
            _refGeom->ajouter_mg_boucle(mcLoop, mcFace);
            for (unsigned l=0;l<loops[k].size();l++)
            {
                CAD4FE::MakeLoops::CoEdge* c = loops[k][l];
                CreateCoEdge(mcFace,mcLoop,(MCEdge*)c->Edge(),c->Sense());
            }
        }
    }}

    {MCBODY_FOR_EACH_MCVERTEX(this, mcVertex)
    {
        MCFace * mcFace = 0;
        int test = TestInteriorVertexIsolatedInFace(mcVertex, &mcFace);
        if (test)
        {
            ExportBRep_CreateInteriorVertexIsolatedInFace(mcVertex, mcFace);
        }
    }}

    return _brep;
}

int MCBody::TestInteriorVertexIsolatedInFace(MCVertex * __mcVertex, MCFace ** __mcFace)
{
    std::set<MCFace*> adjacentFaces;
    Vertex_GetAdjacentFaces(__mcVertex, adjacentFaces);
    std::set<MCEdge*> adjacentEdges;       
    Vertex_GetAdjacentEdges(__mcVertex, adjacentEdges);

    if ( adjacentFaces.size() == 1 && adjacentEdges.size() == 0)
    {
        *__mcFace = *(adjacentFaces.begin());
        return 1;
    }
    else
    {                    
        *__mcFace = 0;
        return 0;
    }
}

int MCBody::ExportBRep_CreateInteriorVertexIsolatedInFace(MCVertex * __mcVertex, MCFace * __mcFace)
{
/*
\item soit $RV$ le sommet r�f�rence du {\mcvertex} $V$,
\item g�n�rer une {\mcedge}, not�e $E$, n'ayant aucune ar�te r�f�rence,
  mais contenant sommet r�f�rence $RV$,
\item g�n�rer une coar�te $CE$ liant $E$ � $F$,
\item ajouter � $F$ une boucle $B$ contenant $CE$
*/
    PolyCurve * edgeCurve = new PolyCurve(__mcVertex->GetRefVertex());
    MCEdge * edge = AddMCEdge(edgeCurve);
    std::vector<MCVertex*> edgeVertices;
    for (int i=0;i<2;i++) edgeVertices.push_back(__mcVertex);
    AddMCEdgeCovertices(edge,edgeVertices);
    MG_BOUCLE * mcLoop = new MG_BOUCLE(__mcFace);
    _refGeom->ajouter_mg_boucle(mcLoop, __mcFace);
    CreateCoEdge(__mcFace, mcLoop, edge, +1);
    return 1;
}

void MCBody::GetFEVCount(int * __faceCount,int * __edgeCount,int * __vertexCount)
{
    * __faceCount = _G21->GetNodes().size();
    * __edgeCount = _G21->GetArcs().size();
    * __vertexCount = _G10->GetArcs().size();
}

void MCBody::GetRefBodyFEVCount(int * __faceCount,int * __edgeCount,int * __vertexCount)
{
    MG_UTILS::GetFEVCount(_refBody,__faceCount,__edgeCount,__vertexCount);
}

int MCBody::DeleteBRep(MG_GESTIONNAIRE * gest, MG_GEOMETRIE * mggeo)
{
    /// Lists of MCT entities
    std::set < CAD4FE::MCFace * > lst_mc_faces;
    std::set < CAD4FE::MCEdge * > lst_mc_edges;
    std::set < CAD4FE::MCVertex * > lst_mc_vertices;
    std::set < MG_COQUILLE * > lst_mc_shell;
    std::set < MG_BOUCLE * > lst_mc_loop;

    // initialize MC face, MC edge, MC loop, MC shell, and MC vertices lists
    unsigned nb_shells = _brep->get_nb_mg_coquille ();
    for (unsigned it_shells = 0; it_shells < nb_shells; it_shells++)
    {
            MG_COQUILLE * shell = _brep->get_mg_coquille(it_shells);
            lst_mc_shell.insert(shell);
            unsigned nb_coface = shell->get_nb_mg_coface();
            for (unsigned it_coface = 0; it_coface < nb_coface; it_coface++)
            {
                    CAD4FE::MCFace * face = (CAD4FE::MCFace *)shell->get_mg_coface(it_coface)->get_face();
                    lst_mc_faces.insert(face);

                    unsigned nb_loop = face->get_nb_mg_boucle();
                    for (unsigned it_loop = 0; it_loop < nb_loop; it_loop++)
                    {
                            MG_BOUCLE * loop = face->get_mg_boucle(it_loop);
                            lst_mc_loop.insert(loop);
                            unsigned nb_edge = loop->get_nb_mg_coarete();

                            for (unsigned it_edge = 0; it_edge < nb_edge; it_edge++)
                            {
                                    CAD4FE::MCEdge * edge = (CAD4FE::MCEdge*)loop->get_mg_coarete(it_edge)->get_arete();
                                    lst_mc_edges.insert(edge);

                                    lst_mc_vertices.insert ((CAD4FE::MCVertex*)edge->get_cosommet1()->get_sommet());
                                    lst_mc_vertices.insert ((CAD4FE::MCVertex*)edge->get_cosommet2()->get_sommet());
                            }
                    }
            }
    }

    for (unsigned i=0;i<gest->get_nb_mg_maillage();)
    {
        MG_MAILLAGE * tmpMesh = gest->get_mg_maillage(i);
        LISTE_MG_TRIANGLE::iterator itTri;
        for (MG_TRIANGLE * tri = tmpMesh->get_premier_triangle(itTri);
            tri; tri = tmpMesh->get_suivant_triangle(itTri))
            if (lst_mc_faces.end() != lst_mc_faces.find((CAD4FE::MCFace*)tri->get_lien_topologie()))
                tmpMesh->supprimer_mg_triangleid(tri->get_id());
        LISTE_MG_SEGMENT::iterator itSeg;
        for (MG_SEGMENT * seg = tmpMesh->get_premier_segment(itSeg);
            seg; seg=tmpMesh->get_suivant_segment(itSeg) )
            if (lst_mc_edges.end() != lst_mc_edges.find((CAD4FE::MCEdge*)seg->get_lien_topologie())
                || lst_mc_faces.end() != lst_mc_faces.find((CAD4FE::MCFace*)seg->get_lien_topologie())
            )
                tmpMesh->supprimer_mg_segmentid(seg->get_id());
        LISTE_MG_NOEUD::iterator itNo;
        for (MG_NOEUD * no = tmpMesh->get_premier_noeud(itNo);
            no; no=tmpMesh->get_suivant_noeud(itNo))
            if (
                lst_mc_vertices.end() != lst_mc_vertices.find((CAD4FE::MCVertex*)no->get_lien_topologie())
                || lst_mc_edges.end() != lst_mc_edges.find((CAD4FE::MCEdge*)no->get_lien_topologie())
                || lst_mc_faces.end() != lst_mc_faces.find((CAD4FE::MCFace*)no->get_lien_topologie())
            )
                tmpMesh->supprimer_mg_noeudid(no->get_id());

        if (tmpMesh->get_nb_mg_noeud() == 0)
            gest->supprimer_mg_maillage(i);
        else
            i++;
    }
    for (std::set < MG_COQUILLE * >::iterator itShell = lst_mc_shell.begin();
        itShell != lst_mc_shell.end();
        itShell++)
        mggeo->supprimer_mg_coquilleid((*itShell)->get_id());
    for (std::set<CAD4FE::MCFace *>::iterator itFace = lst_mc_faces.begin(); itFace != lst_mc_faces.end(); itFace++)
    {
        CAD4FE::MCFace* mcFace=*itFace;
        CAD4FE::PolySurface * ps = mcFace->GetPolySurface();
        mggeo->supprimer_mg_faceid(mcFace->get_id());
        mggeo->supprimer_mg_surfaceid(ps->get_id());
    }
    for (std::set < MG_BOUCLE * >::iterator itLoop = lst_mc_loop.begin();
        itLoop != lst_mc_loop.end();
        itLoop++)
        mggeo->supprimer_mg_boucleid((*itLoop)->get_id());

    for (std::set<CAD4FE::MCEdge *>::iterator itEdge = lst_mc_edges.begin(); itEdge != lst_mc_edges.end(); itEdge++)
    {
        CAD4FE::MCEdge* mcEdge=*itEdge;
        CAD4FE::PolyCurve * pc=mcEdge->GetPolyCurve();
        mggeo->supprimer_mg_areteid(mcEdge->get_id());
        mggeo->supprimer_mg_courbeid(pc->get_id());
    }
    for (std::set<CAD4FE::MCVertex *>::iterator itVertex = lst_mc_vertices.begin(); itVertex != lst_mc_vertices.end(); itVertex++)
    {
        CAD4FE::MCVertex * mcVertex = *itVertex;
        mggeo->supprimer_mg_sommetid(mcVertex->get_id());
    }     

    mggeo->supprimer_mg_volumeid(_brep->get_id());

    return 1;

}
Revision:	830
Committed:	Fri Sep 30 20:08:17 2016 UTC (8 years, 7 months ago) by francois
File size:	69182 byte(s)
Log Message:	Posttraitement de l'optimisation de topologie * transfert des conditions aux limites sur le maillage de peau via une geometrie virtuelle * ajout d'un fichier .volume qui suit l'évolution du volume au travers du processus de posttraitement