// Cutmesh - Copyright (C) 2010-2018 T. Mouton, E. Bechet 
//
// See the LICENSE file for license information and contributions.
// bugs and problems to <thibaud.mouton@gmail.com>.

#include "topo.h"

void Topo::build_topo ( Mesh *m, LSurface *s )
{
  tprintf ( "Topology construction START\n" );
  m->set_all_vertex_flag ( -1 );
  _em = m->get_entity_manager();
  // table de flag des faces topologiques
  std::vector<int> flag_topo_face ( s->topo_face_size(), -1 );
  // remplissage de la topologie
  for ( int i = 0; i < s->border_edge_size(); i++ )
  {
    CutEdge *ce = s->get_border_edge ( i );
    int edge_color = ce->get_color();
    if ( ce->get_adj ( 0 ) == NULL )
    {
      // ajout d'une topoedge
      TopoEdge *te = new_topoedge ( -1, edge_color );
      // ajout des topofaces
      for ( int j = 0; j < ce->face_size(); j++ )
      {
        CutFace *cf = ce->get_face ( j );
        int face_color = cf->get_color();
        assert ( face_color != -1 );
        TopoFace *tf = NULL;
        if ( flag_topo_face[face_color] == -1 )
        {
          tf = new_topoface ( -1, face_color );
          flag_topo_face[face_color] = tf->index();
          for ( int k = 0; k < s->face_size(); k++ )
          {
            if ( s->get_face ( k )->get_color() == face_color )
              tf->add_face ( s->get_face ( k ) );
          }
        }
        else
        {
          int index = flag_topo_face[face_color];
          tf = get_topoface ( index );
        }
        tf->add_topoedge ( te );
        te->add_topoface ( tf );
      }
      // ajout des topovertex
      TopoVertex *tv = NULL;
      LVertex *v0 = ce->get_vertex ( 0 );
      if ( v0->getFlag() == -1 )
      {
        tv = new_topovertex ( v0, -1 );
        v0->setFlag ( tv->index() );
      }
      else
      {
        tv = get_topovertex ( v0->getFlag() );
      }
      te->set_topovertex ( 0, tv );
      tv->add_topoedge ( te );
      CutEdge *n = ce;
      CutEdge *p = n;
      while ( n != NULL ) // attention, ça boucle
      {
        te->add_edge ( n );
        p = n;
        n = n->get_adj ( 1 );
      }
      LVertex *v1 = p->get_vertex ( 1 );
      if ( v1->getFlag() == -1 )
      {
        tv = new_topovertex ( v1, -1 );
        v1->setFlag ( tv->index() );
      }
      else
      {
        tv = get_topovertex ( v1->getFlag() );
      }
      te->set_topovertex ( 1, tv );
      tv->add_topoedge ( te );
    }
  }

  tprintf ( "--> %4d topovertices\n", topovertex_size() );
  tprintf ( "--> %4d topoedges\n", topoedge_size() );
  tprintf ( "--> %4d topofaces\n", topoface_size() );

  // reoganisation des entites topologiques
  tprintf ( "Reorganizing faces ...\n" );
  for ( int i = 0; i < topoface_size(); i++ )
    get_topoface ( i )->build_internal ( get_entity_manager() );
  tprintf ( "Reorganizing edges ...\n" );
  for ( int i = 0; i < topoedge_size(); i++ )
    get_topoedge ( i )->build_internal ( get_entity_manager() );
  tprintf ( "Reorganizing vertices ...\n" );
  for ( int i = 0; i < topovertex_size(); i++ )
    get_topovertex ( i )->build_internal ( get_entity_manager() );

//             dump_bohren();
//        dump();


  tprintf ( "Topology construction END\n" );
}

void Topo::init_comparison()
{
  build_vertex_comparison_table();
  build_edge_comparison_table();
  build_face_comparison_table();
}

void TopoFace::build_internal ( EntityManager *em )
{
  double tolerance = 0.5; // hack tolerance
  // orientation des edges
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    int dir = te->get_direction ( color() );
    if ( !dir )
      set_reverse_dir ( i );
  }

  for ( int i = 0; i < topoedge_size()-2; i++ ) // pas optimal pour des grandes quantites, mais efficace pour des petites !
  {
    TopoVertex *tv = get_topovertex ( i, 1 );
    if ( tv != get_topovertex ( i+1, 0 ) )
      for ( int j = i+2; j < topoedge_size(); j++ )
        if ( tv == get_topovertex ( j, 0 ) )
          swap_topoedge ( i+1, j );
  }

  update_surf ( em );
}

void TopoFace::update_surf ( EntityManager *em )
{
  // calcul des surfaces presentes dans la topoface
  clear_surf();
  std::map<int, int> surf_map;
  std::map<int, int>::iterator it;
  std::multimap<int, int> reverse_surf_map;
  std::map<int, int>::reverse_iterator rit;
  for ( int i = 0; i < face_size(); i++ )
  {
    CutFace *cf = get_face ( i );
    for ( int j = 0; j < cf->surf_size(); j++ )
    {
      int surf = cf->surface ( j );
      it = surf_map.find ( surf );
      if ( it != surf_map.end() )
        it->second++;
      else
        surf_map.insert ( std::pair<int, int> ( surf, 1 ) );
    }
  }

  for ( it = surf_map.begin(); it != surf_map.end(); it++ )
    reverse_surf_map.insert ( std::pair<int, int> ( it->second, it->first ) ); // equivalent a sort

  int count_inserted = 0;
  int last_score = 0;
  for ( rit = reverse_surf_map.rbegin(); rit != reverse_surf_map.rend(); rit++ )
  {
    int surf = rit->second;
    if ( !count_inserted || rit->first == last_score ) // on insert seulement la plus significative
    {
      add_surf ( surf );
      last_score = rit->first;
      count_inserted++;
    }
  }
}

void TopoFace::update_topoedge()
{
  std::map<TopoEdge*, int> count_edge;
  std::map<TopoEdge*, int>::iterator itce;
  for ( int i = 0; i < child_size(); i++ )
  {
    TopoFace *c = get_child ( i );
    for ( int j = 0; j < c->topoedge_size(); j++ )
    {
      itce = count_edge.find ( c->get_topoedge ( j ) );
      if ( itce != count_edge.end() )
        itce->second++;
      else
        count_edge.insert ( std::pair<TopoEdge*, int> ( c->get_topoedge ( j ), 1 ) );
    }
  }
  for ( itce = count_edge.begin(); itce != count_edge.end(); itce++ )
  {
    if ( itce->second == 1 )
      add_topoedge ( itce->first );
  }
}

void TopoFace::add_common_surfaces()
{
  std::map<int, int> surf_map;
  std::map<int, int>::iterator it;
//     printf("child size --> %d\n", child_size());
  assert ( child_size() == 2 );
  clear_surf();

  for ( int i = 0; i < child_size(); i++ )
  {
    TopoFace *tf = get_child ( i );
    for ( int j = 0; j < tf->surf_size(); j++ )
    {
      int surf = tf->get_surf ( j );
//          printf("put surf --> %d\n",surf);
      it = surf_map.find ( surf );
      if ( it == surf_map.end() )
        surf_map.insert ( std::pair<int, int> ( surf, 1 ) );
      else
        it->second++;
    }
  }
  for ( it = surf_map.begin(); it != surf_map.end(); ++it )
  {
    if ( it->second == 2 )
    {
      int surf = it->first;
//          printf("add surf --> %d\n", surf);
      add_surf ( surf );
    }
  }
}

void TopoEdge::build_internal ( EntityManager *em )
{
  update_surf ( em );
}

void TopoEdge::update_surf ( EntityManager *em )
{
  double tolerance = 0.5;
  // calcul des surfaces presentes dans la topoedge
  clear_surf();
  std::map<int, int> surf_map;
  std::map<int, int>::iterator it;
  std::multimap<int, int> reverse_surf_map;
  std::map<int, int>::reverse_iterator rit;
  for ( int i = 0; i < vertex_size(); i++ )
  {
    LVertex *v = get_vertex ( i, 1 );
    std::set<int> surf_set;
    for ( int j = 0; j < v->nbSurfZero(); j++ )
    {
      int surf = v->surfZeroId ( j );
      if ( em->is_virtual_face ( surf ) )
      {
        int s1, s2;
        em->get_cutted_face_id ( surf, &s1, &s2 );
        if ( v->checkZero ( s1 ) || v->checkZero ( s2 ) ) // important --> assure que le vertex appartient bien a la surface virtuelle
        {
          surf_set.insert ( s1 );
          surf_set.insert ( s2 );
          surf_set.insert ( surf );
        }
      }
      else
        surf_set.insert ( surf );
    }
    for ( std::set<int>::iterator itset = surf_set.begin(); itset != surf_set.end(); itset++ )
    {
      int surf = *itset;
      it = surf_map.find ( surf );
      if ( it != surf_map.end() )
        it->second++;
      else
        surf_map.insert ( std::pair<int, int> ( surf, 1 ) );
    }
  }

//     printf("topoedge %d :\n", index());
  for ( it = surf_map.begin(); it != surf_map.end(); it++ )
  {
    reverse_surf_map.insert ( std::pair<int, int> ( it->second, it->first ) ); // equivalent a sort
//         printf("--> surf %d : score = %d\n", it->first, it->second);
  }

  int last_score = 0;
  for ( rit = reverse_surf_map.rbegin(); rit != reverse_surf_map.rend(); rit++ )
  {
    int surf = rit->second;
//         printf("trying surf %d with score %d\n", rit->second, rit->first);
    if ( rit->first >= tolerance*last_score )
    {
      if ( em->is_real_face ( surf ) || em->is_restore_face ( surf ) )
      {
        add_surf ( surf );
      }
      else
      {
//                 printf("set imposed topoedge %d\n", index());
        if ( rit->first == vertex_size() )
        {
          set_imposed();
          get_topovertex ( 0 )->set_imposed();
          get_topovertex ( 1 )->set_imposed();
        }
      }
      last_score = rit->first;
    }
  }
}

void TopoVertex::build_internal ( EntityManager *em )
{
  // calcul des surfaces presentes dans le topovertex
  LVertex *v = get_vertex();
  std::set<int> surf_map;
  std::set<int>::iterator it;
  int s1, s2;
  for ( int j = 0; j < v->nbSurfZero(); j++ )
  {
    int surf = v->surfZeroId ( j );
    if ( em->is_virtual_face ( surf ) )
    {
      em->get_cutted_face_id ( surf, &s1, &s2 );
      if ( v->checkZero ( s1 ) || v->checkZero ( s2 ) ) // important --> assure que le vertex appartient bien a la surface virtuelle
      {
        surf_map.insert ( s1 );
        surf_map.insert ( s2 );
      }
    }
    else
      surf_map.insert ( surf );
  }
  for ( it = surf_map.begin(); it != surf_map.end(); it++ )
    add_surf ( *it );
}

void Topo::build_face_tree()
{
  topoface_map_clear();

  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    for ( int j = 0; j < tf->surf_size(); j++ )
    {
      int surf = tf->get_surf ( j );
//             if (_em->is_restore_face(surf))
      topoface_map_insert ( tf->get_surf ( j ), tf );
    }
  }

  std::set<int> surf_list = get_topoface_key_list();
  std::set<int>::iterator itlist;

  for ( itlist = surf_list.begin(); itlist != surf_list.end(); ++itlist )
  {
    int surf = *itlist;
    std::set<TopoFace *> tfset = get_topoface_mapped_set ( surf );
    if ( tfset.size() > 1 )
    {
      TopoFace *root = new_root_topoface();
      root->add_surf ( surf );
      root->add_physical ( surf );
      std::set<TopoFace *>::iterator itset;
      for ( itset = tfset.begin(); itset!=tfset.end(); ++itset )
      {
        TopoFace *tf = *itset;
        root->add_child ( tf );
        tf->set_parent ( root );
      }
//             root->add_common_surfaces();
    }
  }
}


void Topo::build_edge_tree()
{
  topoedge_map_clear();

  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    // recuperation premiere arete
    CutEdge *ce = te->get_edge ( 0 );
    LVertex *v0 = ce->get_vertex ( 0 );
    LVertex *v1 = ce->get_vertex ( 1 );

    v0->computeRootParents();
    v1->computeRootParents();

    LVertex *rv0 = NULL;
    LVertex *rv1 = NULL;

    std::set<LVertex*> roots;

    if ( !v0->getRootParentSize() )
      roots.insert ( v0 );
    else
    {
      for ( int j = 0; j < v0->getRootParentSize(); j++ )
        roots.insert ( v0->getRootParent ( j ) );
    }
    if ( !v1->getRootParentSize() )
      roots.insert ( v1 );
    else
    {
      for ( int j = 0; j < v1->getRootParentSize(); j++ )
        roots.insert ( v1->getRootParent ( j ) );
    }

    if ( roots.size() == 2 )
    {
      std::set<LVertex*>::iterator itroots = roots.begin();
      rv0 = *itroots;
      itroots++;
      rv1 = *itroots;

      std::map<std::pair<int, int>, int>::iterator it;
//          printf("rv0 --> %d, rv1 --> %d\n", rv0->getIndex(), rv1->getIndex());

      if ( rv0->getIndex() < rv1->getIndex() )
        it = _previous_edge_map.find ( std::pair<int, int> ( rv0->getIndex(), rv1->getIndex() ) );
      else
        it = _previous_edge_map.find ( std::pair<int, int> ( rv1->getIndex(), rv0->getIndex() ) );

      if ( it != _previous_edge_map.end() )
      {
        int edge = it->second;
        topoedge_map_insert ( edge, te );
      }
    }
  }

  std::set<int> edge_list = get_topoedge_key_list();
  std::set<int>::iterator itlist;

  for ( itlist = edge_list.begin(); itlist != edge_list.end(); ++itlist )
  {
    int edge = *itlist;
    std::set<TopoEdge *> teset = get_topoedge_mapped_set ( edge );
    if ( teset.size() > 1 )
    {
//             printf("edge --> %d\n", edge);
      std::set<TopoFace*> faceset;
      TopoEdge *root = new_root_topoedge();
      root->add_physical ( edge );
      std::set<TopoEdge *>::iterator itset;
      for ( itset = teset.begin(); itset != teset.end(); ++itset )
      {
        TopoEdge *te = *itset;
//                 printf("topoedge --> %d\n", te->index());
        for ( int j = 0; j < te->topoface_size(); j++ )
        {
          TopoFace *tf = te->get_topoface ( j );
          if ( tf->get_parent() != NULL )
            faceset.insert ( tf->get_parent() );
        }
        root->add_child ( te );
        te->set_parent ( root );
      }

      for ( std::set<TopoFace*>::iterator it = faceset.begin(); it != faceset.end(); it++ )
      {
        TopoFace *tf = *it;
        root->add_topoface ( tf );
        tf->add_topoedge ( root );
      }
    }
    else if ( teset.size() == 1 )
    {
      std::set<TopoEdge *>::iterator itset;
      for ( itset = teset.begin(); itset!=teset.end(); ++itset )
      {
        TopoEdge *te = *itset;
        te->add_physical ( edge );
      }
    }
  }
}


void Topo::build_volume_tree()
{
  topovolume_map_clear();

  for ( int i = 0; i < topovolume_size(); i++ )
  {
    TopoVolume *tvol = get_topovolume ( i );
    LTetra *t = tvol->get_seed();
    dbgprintf ( 2, "topovolume %d (%d) --> old %d vol %d\n", i, tvol->physical_size(), t->getOldVol(), t->getVol() );
    topovolume_map_insert ( t->getOldVol(), tvol );
    topovolume_map_insert ( t->getVol(), tvol );
  }

  std::set<int> vol_list = get_topovolume_key_list();
  std::set<int>::iterator itlist;

  for ( itlist = vol_list.begin(); itlist != vol_list.end(); ++itlist )
  {
    int volume = *itlist;
    std::set<TopoVolume *> tfset = get_topovolume_mapped_set ( volume );
    dbgprintf ( 2, "volume tag %d --> %d subsets\n", volume, ( int ) tfset.size() );

    if ( tfset.size() > 1 )
    {
      TopoVolume *root = new_root_topovolume();
      root->add_physical ( volume );
      std::set<TopoVolume *>::iterator itset;
      for ( itset = tfset.begin(); itset!=tfset.end(); ++itset )
      {
        TopoVolume *tvol = *itset;
        root->add_child ( tvol );
        tvol->set_parent ( root );
      }
    }
  }
}


void Topo::dump_heritage()
{
  std::set<TopoFace *> facerootset;
  std::set<TopoFace *>::iterator itfr;

  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    TopoFace *tfp = tf->get_parent();
    if ( tfp != NULL )
      facerootset.insert ( tfp );
  }

  for ( itfr = facerootset.begin(); itfr != facerootset.end(); itfr++ )
  {
    TopoFace *tfp = *itfr;
    tprintf ( "--> Topoface %d\n", tfp->index() );
    for ( int j = 0; j < tfp->child_size(); j++ )
    {
      tprintf ( "     +--> Topoface %d face %d (%d)\n", tfp->get_child ( j )->index(), _em->get_face_id ( tfp->get_surf ( 0 ) ), tfp->get_surf ( 0 ) );
    }
  }

  std::set<TopoEdge *> edgerootset;
  std::set<TopoEdge *>::iterator iter;

  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    TopoEdge *tep = te->get_parent();
    if ( tep != NULL )
      edgerootset.insert ( tep );
  }

  for ( iter = edgerootset.begin(); iter != edgerootset.end(); iter++ )
  {
    TopoEdge *tep = *iter;
    if ( tep != NULL )
    {
      tprintf ( "--> Topoedge %d\n", tep->index() );
      for ( int j = 0; j < tep->child_size(); j++ )
      {
        tprintf ( "     +--> TopoEdge %d\n", tep->get_child ( j )->index() );
      }
    }
  }

#if 1
  dumpVTK ( "final" );
#endif
}

void Topo::build_hierarchy ( Mesh *m )
{
  build_face_tree();
  build_edge_tree();
  build_volume_tree();


  if ( !_opts->brep_topo )
  {
    // vertex physical assignment
    for ( int i = 0; i < topovertex_size(); i++ )
    {
      TopoVertex *tv = get_topovertex ( i );
      tv->add_physical ( i+1 );
    }
    // edge physical assignment
    for ( int i = 0; i < topoedge_size(); i++ )
    {
      TopoEdge *te = get_topoedge ( i );
      if ( te->get_parent() != NULL || !te->physical_size() )
      {
        int new_edge_id = get_entity_manager()->add_edge_id();
        te->add_physical ( new_edge_id );
      }
    }
    // face physical assignment
    for ( int i = 0; i < topoface_size(); i++ )
    {
      TopoFace *tf = get_topoface ( i );
      if ( tf->get_parent() != NULL )
      {
        int new_face_id = get_entity_manager()->add_face_id();
        tf->add_physical ( new_face_id );
      }
      else
        tf->add_physical ( tf->get_surf ( 0 ) );
    }
    // volume physical assignment
    for ( int i = 0; i < topovolume_size(); i++ )
    {
      TopoVolume *tvol = get_topovolume ( i );
      LTetra *t = tvol->get_seed();
      if ( tvol->get_parent() != NULL )
      {
        tvol->add_physical ( t->getVol() );
      }
      else
      {
        tvol->add_physical ( t->getOldVol() );
        m->seed ( t );
        m->color ( t->getOldVol(), t->getVol() );
      }
    }
  }

  dump_heritage();
}


int Topo::merge_edge ( TopoVertex *tv )
{
  assert ( tv->topoedge_size() == 2 );
  assert ( !tv->associated() );
  assert ( tv->getFlag() == -1 );


  TopoEdge *tep = tv->get_topoedge ( 0 );
  TopoEdge *tes = tv->get_topoedge ( 1 );

  assert ( tep->getFlag() != -1 && tes->getFlag() != -1 );

//     tep->dump();
//     tes->dump();

  if ( tep->getFlag() != tes->getFlag() )
    return 0;

//     if (tep->getFlag() != 1 || tes->getFlag() != 1)
//         return 0;

  assert ( tep->get_ref() == NULL );
  assert ( tes->get_ref() == NULL );

  if ( tep == tes )
    return 0;

//     printf("--> merging edge %d and edge %d :\n", tep->index(), tes->index());

  // trouve le topovertex commun
  int pos1 = -1;
  int pos2 = -1;
  for ( int i = 0; i < 2 && pos1 == -1; i++ )
    for ( int j = 0; j < 2 && pos1 == -1; j++ )
      if ( tep->get_topovertex ( i ) == tes->get_topovertex ( j ) )
      {
        pos1 = i;
        pos2 = j;
      }

  assert ( pos1 != -1 && pos2 != -1 );
  assert ( tep->get_topovertex ( pos1 ) == tes->get_topovertex ( pos2 ) );

  TopoVertex *tvc = tep->get_topovertex ( pos1 );
  TopoEdge *nte = new_topoedge ( tep->index(), tep->color() ); // cree une nouvelle edge
//      printf("--> create new edge %d\n", tep->index());

//     nte->setFlag(1);
  nte->setFlag ( tep->getFlag() );
  tep->setFlag ( -1 );
  tes->setFlag ( -1 );

  tep->set_ref ( nte );
  tes->set_ref ( nte );

  // copy attributes to new topoedge
  for ( int i = 0; i < tep->edge_size(); i++ )
    nte->add_edge ( tep->get_edge ( i, 1 ) );
  for ( int i = 0; i < tep->potential_brep_index_size(); i++ )
    nte->add_potential_brep_index ( tep->get_potential_brep_index ( i ) );
  for ( int i = 0; i < tes->potential_brep_index_size(); i++ )
    nte->add_potential_brep_index ( tes->get_potential_brep_index ( i ) );

  assert ( tes->get_topovertex ( 0 )->get_vertex() == tes->get_vertex ( 0, 1 ) );
  assert ( tes->get_topovertex ( 1 )->get_vertex() == tes->get_vertex ( tes->vertex_size()-1, 1 ) );
  assert ( tes->get_topovertex ( 1 )->get_vertex() == tes->get_vertex ( 0, 0 ) );
  assert ( tes->get_topovertex ( 0 )->get_vertex() == tes->get_vertex ( tes->vertex_size()-1, 0 ) );

  int xorpos = ( !pos1&&pos2 ) || ( pos1&&!pos2 ); // pos1 xor pos2
  TopoVertex *tvn0 = NULL;
  TopoVertex *tvn1 = NULL;
  if ( pos1 == 1 )
  {
    for ( int i = 0; i < tes->edge_size(); i++ )
    {
      CutEdge *ce = tes->get_edge ( i, xorpos );
      ce->set_color ( tep->color() );
      if ( pos2 )
        ce->invert_dir();
//          printf("add edge %d -- %d\n", ce->get_vertex(0)->getIndex(), ce->get_vertex(1)->getIndex());
      nte->add_edge ( ce );
    }
    tvn0 = tep->get_topovertex ( 0 );
    tvn1 = tes->get_topovertex ( ( pos2+1 ) %2 );
    nte->set_topovertex ( 0, tvn0 );
    nte->set_topovertex ( 1, tvn1 );
  }
  else
  {
    for ( int i = 0; i < tes->edge_size(); i++ )
    {
      CutEdge *ce = tes->get_edge ( i, xorpos );
      ce->set_color ( nte->color() );
      if ( !pos2 )
        ce->invert_dir();
//          printf("insert edge %d -- %d\n", ce->get_vertex(0)->getIndex(), ce->get_vertex(1)->getIndex());
      nte->insert_edge ( i, ce );
    }
    tvn0 = tep->get_topovertex ( 1 );
    tvn1 = tes->get_topovertex ( ( pos2+1 ) %2 );
    nte->set_topovertex ( 0, tvn1 );
    nte->set_topovertex ( 1, tvn0 );
  }

  assert ( tvn0 != NULL && tvn1 != NULL );

  to_be_updated ( tvn0 );
  to_be_updated ( tvn1 );

  // reassocie les aretes pour chaque topovertex
  tvn0->add_topoedge ( nte );
  tvn1->add_topoedge ( nte );
  for ( int i = 0; i < tvn0->topoedge_size(); i++ )
    if ( tvn0->get_topoedge ( i ) == tep )
    {
      tvn0->remove_topoedge ( i );
      break;
    }
  for ( int i = 0; i < tvn1->topoedge_size(); i++ )
    if ( tvn1->get_topoedge ( i ) == tes )
    {
      tvn1->remove_topoedge ( i );
      break;
    }

  nte->update_surf ( _em );
  topoedge_map_remove_surf ( tep );
  topoedge_map_remove_surf ( tes );
  topoedge_map_insert_surf ( nte );

  assert ( nte->get_topovertex ( 0 )->get_vertex() == nte->get_vertex ( 0, 1 ) );
  assert ( nte->get_topovertex ( 1 )->get_vertex() == nte->get_vertex ( nte->vertex_size()-1, 1 ) );

//     nte->dump();

//     _lock = 1;

  return 1;
}

int Topo::merge_face ( TopoEdge *te )
{
  assert ( te->topoface_size() == 2 );
  assert ( te->getFlag() == -1 );
  assert ( !te->associated() );

  TopoFace *tfp = te->get_topoface ( 0 );
  TopoFace *tfs = te->get_topoface ( 1 );

  assert ( tfp != NULL );
  assert ( tfs != NULL );
  assert ( tfp->getFlag() == 1 );
  assert ( tfs->getFlag() == 1 );

//     printf("Edge %d --> merging topoface %d with topoface %d\n", te->index(), tfp->index(), tfs->index());

//     tfp->dump();
//     tfs->dump();

  if ( tfp == tfs ) // automerge --> pas de probleme de sens
  {
//      printf("--> automerge topoface %d\n", tfp->index());
    int eid0 = -1, eid1 = -1;
    std::vector<TopoEdge *> list;
    for ( int i = 0; i < tfp->topoedge_size() && eid0 == -1; i++ )
    {
      if ( tfp->get_topoedge ( i ) == te )
        eid0 = i;
    }
    for ( int i = 0; i < tfp->topoedge_size() && eid1 == -1; i++ )
    {
      if ( tfp->get_topoedge ( eid0+i+1 ) == te )
      {
        eid1 = ( eid0+i+1 ) %tfp->topoedge_size();
      }
      else
      {
        list.push_back ( tfp->get_topoedge ( eid0+i ) );
        tfp->remove_topoedge ( eid0+i-- );
      }
    }

    assert ( eid0 != -1 && eid1 != -1 );

    for ( int i = 0; i < list.size(); i++ )
      tfp->add_topoedge ( list[i] );

    for ( int i = 0; i < tfp->topoedge_size(); i++ )
    {
//          assert(tfp->get_topoedge(i)->get_ref() == NULL);
      if ( tfp->get_topoedge ( i ) == te )
        tfp->remove_topoedge ( i-- );
    }
//         printf("merging result :\n");
//         tfp->dump();
  }
  else
  {
    TopoFace *ntf = new_topoface ( tfp->index(), tfp->color() ); // cree une nouvelle face
//         printf("--> create new face %d\n", ntf->index());

    ntf->setFlag ( tfp->getFlag() );
    tfp->setFlag ( -1 );
    tfs->setFlag ( -1 );

    tfp->set_ref ( ntf );
    tfs->set_ref ( ntf );

    // copy attributes to new topoface
    for ( int i = 0; i < tfp->potential_brep_index_size(); i++ )
      ntf->add_potential_brep_index ( tfp->get_potential_brep_index ( i ) );
    for ( int i = 0; i < tfs->potential_brep_index_size(); i++ )
      ntf->add_potential_brep_index ( tfs->get_potential_brep_index ( i ) );

    int eid0 = -1, eid1 = -1;
    for ( int i = 0; i < tfp->topoedge_size() && eid0 == -1; i++ )
      if ( tfp->get_topoedge ( i ) == te )
        eid0 = i;
    for ( int i = 0; i < tfs->topoedge_size() && eid1 == -1; i++ )
      if ( tfs->get_topoedge ( i ) == te )
        eid1 = i;

    assert ( eid0 != -1 && eid1 != -1 );

    for ( int i = 1; i < tfp->topoedge_size(); i++ )
    {
      TopoEdge *te = tfp->get_topoedge ( ( eid0+i ) %tfp->topoedge_size() );
      ntf->add_topoedge ( te );
      for ( int j = 0; j < te->topoface_size(); j++ )
        if ( te->get_topoface ( j ) == tfp )
          te->remove_topoface ( j-- );
      te->add_topoface ( ntf );
    }
    for ( int i = 1; i < tfs->topoedge_size(); i++ )
    {
      TopoEdge *te = tfs->get_topoedge ( ( eid1+i ) %tfs->topoedge_size() );
      ntf->add_topoedge ( te );
      for ( int j = 0; j < te->topoface_size(); j++ )
        if ( te->get_topoface ( j ) == tfs )
          te->remove_topoface ( j-- );
      te->add_topoface ( ntf );
    }
    for ( int i = 0; i < tfp->face_size(); i++ )
    {
      CutFace *cf = tfp->get_face ( i );
      cf->set_color ( ntf->color() );
      ntf->add_face ( cf );
    }
    for ( int i = 0; i < tfs->face_size(); i++ )
    {
      CutFace *cf = tfs->get_face ( i );
      cf->set_color ( ntf->color() );
      ntf->add_face ( cf );
    }

    ntf->update_surf ( _em );

    topoface_map_remove_surf ( tfp );
    topoface_map_remove_surf ( tfs );
    topoface_map_insert_surf ( ntf );
    tfp->clear_all();
    tfs->clear_all();
//      getchar();
//         printf("merging result :\n");
//         ntf->dump();
  }
  te->clear_all();
  assert ( te->topoface_size() == 0 );
  return 1;
}

void Topo::build_vertex_comparison_table()
{
  for ( int i = 0; i < topovertex_size(); i++ )
  {
    TopoVertex *tv = get_topovertex ( i );
    tv->setFlag ( -1 );
    topovertex_map_insert_surf ( tv );
  }
}

void Topo::build_edge_comparison_table()
{
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    te->setFlag ( -1 );
    topoedge_map_insert_surf ( te );
  }
}

void Topo::build_face_comparison_table()
{
//     printf("topoface struct size --> %d\n", topoface_size());
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    tf->setFlag ( -1 );
    topoface_map_insert_surf ( tf );
  }
}

void Topo::topovertex_query ( std::vector<int> vx_surf, std::vector<TopoVertex*> &result )
{
  std::multimap<int, TopoVertex*>::iterator itmap;
  std::pair<std::multimap<int, TopoVertex*>::iterator, std::multimap<int, TopoVertex*>::iterator> ret;
  std::map<TopoVertex*, int> count_map;
  std::map<TopoVertex*, int>::iterator itcount;

//     printf("surf -->");
  for ( int i = 0; i < vx_surf.size(); i++ )
  {
//      printf(" %d", vx_surf[i]);
    ret = _key_vx_map.equal_range ( vx_surf[i] );
    for ( itmap=ret.first; itmap!=ret.second; ++itmap )
    {
      itcount = count_map.find ( itmap->second );
      if ( itcount == count_map.end() )
        count_map.insert ( std::pair<TopoVertex*, int> ( itmap->second, 1 ) );
      else
        itcount->second++;
      if ( itcount->second == vx_surf.size() )
        result.push_back ( itmap->second );
    }
  }
//     printf("\n");
//     assert(result.size());
}

void Topo::topoedge_query ( std::vector<int> edge_surf, std::vector<TopoEdge*> &result )
{
  std::multimap<int, TopoEdge*>::iterator itmap;
  std::pair<std::multimap<int, TopoEdge*>::iterator, std::multimap<int, TopoEdge*>::iterator> ret;
  std::map<TopoEdge*, int> count_map;
  std::map<TopoEdge*, int>::iterator itcount;

//     printf("surf -->");
  for ( int i = 0; i < edge_surf.size(); i++ )
  {
//         printf(" %d", edge_surf[i]);
    ret = _key_edge_map.equal_range ( edge_surf[i] );
    for ( itmap=ret.first; itmap!=ret.second; ++itmap )
    {
      itcount = count_map.find ( itmap->second );
      if ( itcount == count_map.end() )
        count_map.insert ( std::pair<TopoEdge*, int> ( itmap->second, 1 ) );
      else
        itcount->second++;
      if ( itcount->second == edge_surf.size() )
      {
        assert ( itmap->second->get_ref() == NULL );
        result.push_back ( itmap->second );
      }
    }
  }
//     printf("\n");
//     assert(result.size());
}

void Topo::topoface_query ( int face_surf, std::vector<TopoFace*> &tf_result )
{
//     printf("face surf --> %d\n", face_surf);
  std::multimap<int, TopoFace*>::iterator itmap;
  std::pair<std::multimap<int, TopoFace*>::iterator, std::multimap<int, TopoFace*>::iterator> ret;

  ret = _key_face_map.equal_range ( face_surf );
  for ( itmap=ret.first; itmap!=ret.second; ++itmap )
    tf_result.push_back ( itmap->second );
}

void Topo::read_topo_file ( char *name )
{
  int error = 0;
  char filename[100];
  sprintf ( filename, "%s.topo", name );
  tprintf ( "Reading topo file %s ....\n", filename );
  FILE *fp = fopen ( filename, "r" );
  _brm = new BRepModel;
  // lecture des sommets
  int nbvx;
  if ( fscanf ( fp, "%d", &nbvx ) != 1 )
    error++;
  tprintf ( "--> %4d topovertex\n", nbvx );
  for ( int i = 0; i < nbvx; i++ )
  {
    std::vector<int> surf;
    int idvx;
    int idsurf;
    double x, y, z;
    if ( fscanf ( fp, "%d", &idvx ) != 1 )
      error++;
    if ( fscanf ( fp, "%lf %lf %lf", &x, &y, &z ) != 3 )
      error++;
    for ( int j = 0; j < 3; j++ )
    {
      if ( fscanf ( fp, "%d", &idsurf ) != 1 )
        error++;
      idsurf = get_entity_manager()->get_face_id ( idsurf );
      surf.push_back ( idsurf );
    }
    std::sort ( surf.begin(), surf.end() );
    BRepEntity *be = new BRepEntity ( idvx );
    for ( int j = 0; j < 3; j++ )
      be->add_surf ( surf[j] );
    _brm->add_vertex ( be );
    LVertex *v = new LVertex ( x, y, z );
    _brep_vertex_coord.insert ( std::pair<BRepEntity*, LVertex*> ( be, v ) );
  }

  // lecture des edges
  int nbedge;
  if ( fscanf ( fp, "%d", &nbedge ) != 1 )
    error++;
  tprintf ( "--> %4d topoedges\n", nbedge );
  for ( int i = 0; i < nbedge; i++ )
  {
    int idedge;
    std::pair<int, int> idsurf;
    if ( fscanf ( fp, "%d %d %d", &idedge, &idsurf.first, &idsurf.second ) != 3 ) // WARNING --> parfois plus que 3 !!
      error++;
    idsurf.first = get_entity_manager()->get_face_id ( idsurf.first );
    idsurf.second = get_entity_manager()->get_face_id ( idsurf.second );
    if ( idsurf.first > idsurf.second )
    {
      int tmp = idsurf.first;
      idsurf.first = idsurf.second;
      idsurf.second = tmp;
    }
    BRepEntity *be = new BRepEntity ( idedge );
    be->add_surf ( idsurf.first );
    be->add_surf ( idsurf.second );
    _brm->add_edge ( be );
    std::pair<int, int> idendv;
    if ( fscanf ( fp, "%d %d", &idendv.first, &idendv.second ) != 2 )
      error++;
    be->add_adj ( idendv.first );
    be->add_adj ( idendv.second );
  }

  // lecture des faces
  int nbface;
  if ( fscanf ( fp, "%d", &nbface ) != 1 )
    error++;
  tprintf ( "--> %4d topofaces\n", nbface );
  for ( int i = 0; i < nbface; i++ )
  {
    int idface;
    int nbedgeface;
    int idedgeface;
    if ( fscanf ( fp, "%d %d", &idface, &nbedgeface ) != 2 )
      error++;
    idface = get_entity_manager()->get_face_id ( idface );
    BRepEntity *be = new BRepEntity ( idface );
    for ( int j = 0; j < nbedgeface; j++ )
    {
      if ( fscanf ( fp, "%d", &idedgeface ) != 1 )
        error++;
      be->add_adj ( idedgeface );
    }
    be->add_surf ( idface );
    _brm->add_face ( be );
  }
  fclose ( fp );

  assert ( !error );
}

void Topo::get_topo_from_gmodel ( GModel *gm )
{
  _brm = new BRepModel;

  std::map<GFace*, int> gface_map;
  std::map<GVertex*, int> gvertex_map;
  std::map<GEdge*, int> gedge_map;

  int idvx = 0;
  for ( GModel::viter it = gm->firstVertex(); it != gm->lastVertex(); it++ )
  {
    GVertex *gv = *it;
    gvertex_map.insert ( std::pair<GVertex*, int> ( *it, idvx ) );
    LVertex *v = new LVertex ( gv->x(), gv->y(), gv->z() );
    std::set<int> flist;
    std::vector<GEdge*> gle = gv->edges();
    for ( std::vector<GEdge*>::iterator itle = gle.begin(); itle != gle.end(); itle++ )
    {
      GEdge* ge = *itle;
      std::vector<GFace*> glf = ge->faces();
      for ( std::vector<GFace*>::iterator itlf = glf.begin(); itlf != glf.end(); itlf++ )
      {
        std::map<GFace*, int>::iterator itmf = gface_map.find ( *itlf );
        assert ( itmf != gface_map.end() );
        flist.insert ( itmf->second );
      }
    }
    std::vector<int> surf;
    for ( int i = 0; i < 3; i++ )
    {
      int idsurf = *flist.begin();
      flist.erase ( flist.begin() );
      idsurf = get_entity_manager()->get_face_id ( idsurf );
      surf.push_back ( idsurf );
    }
    std::sort ( surf.begin(), surf.end() );
    BRepEntity *be = new BRepEntity ( idvx );
    for ( int j = 0; j < 3; j++ )
      be->add_surf ( surf[j] );
    _brm->add_vertex ( be );
    _brep_vertex_coord.insert ( std::pair<BRepEntity*, LVertex*> ( be, v ) );
    idvx++;
  }

  int idedge = 0;
  for ( GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++ )
  {
    GEdge *ge = *it;
    std::vector<int> surf;
    gedge_map.insert ( std::pair<GEdge*, int> ( *it, idedge ) );
    std::vector<GFace*> glf = ge->faces();
    std::vector<GFace*>::iterator itlf;
    for ( itlf = glf.begin(); itlf != glf.end(); itlf++ )
    {
      std::map<GFace*, int>::iterator itmf = gface_map.find ( *itlf );
      assert ( itmf != gface_map.end() );
//             fprintf(fp, "%d ", itmf->second);
      int idsurf = itmf->second;
      idsurf = get_entity_manager()->get_face_id ( idsurf );
      surf.push_back ( idsurf );
    }
    std::sort ( surf.begin(), surf.end() );
    BRepEntity *be = new BRepEntity ( idedge );
    for ( int j = 0; j < 2; j++ )
      be->add_surf ( surf[j] );

    GVertex *gv = ge->getBeginVertex();
    std::map<GVertex*, int>::iterator itmv = gvertex_map.find ( gv );
    assert ( itmv != gvertex_map.end() );
//         fprintf(fp, "%d ", itmv->second);
    be->add_adj ( itmv->second );
    gv = ge->getEndVertex();
    itmv = gvertex_map.find ( gv );
    assert ( itmv != gvertex_map.end() );
//         fprintf(fp, "%d\n", itmv->second);
    be->add_adj ( itmv->second );
    idedge++;
  }

  int idface = 0;
  for ( GModel::fiter it = gm->firstFace(); it != gm->lastFace(); it++ )
  {
    GFace *gf = *it;
    int local_idface;
    gface_map.insert ( std::pair<GFace*, int> ( *it, idface ) );
    std::vector<GEdge*> gle = gf->edges();
//         fprintf(fp, "%d %d\n", idface, (int)gle.size());
    local_idface = get_entity_manager()->get_face_id ( idface );
    BRepEntity *be = new BRepEntity ( local_idface );
    std::vector<GEdge*>::iterator itlfe;
    for ( itlfe = gle.begin(); itlfe != gle.end(); itlfe++ )
    {
      std::map<GEdge*, int>::iterator itme = gedge_map.find ( *itlfe );
      assert ( itme != gedge_map.end() );
      be->add_adj ( itme->second );
//             fprintf(fp, "%d ", itme->second);
    }
    be->add_surf ( idface );
    _brm->add_face ( be );
//         fprintf(fp, "\n");
    idface++;
  }
}

int Topo::find_edge_occurrence ( TopoEdge *te, BRepEntity *bv0, BRepEntity *bv1, TopoVertex **tv0, TopoVertex **tv1, int &couple, int &ambiguous )
{
  int count = 0;
  int v0 = bv0->index();
  int v1 = bv1->index();
  int om[2][2] = {{0, 0}, {0, 0}};
  couple = 0;
  ambiguous = 0;

  for ( int i = 0; i < 2; i++ )
  {
    TopoVertex *tv = te->get_topovertex ( i );
    for ( int j = 0; j < tv->potential_brep_index_size(); j++ )
    {
      if ( tv->get_potential_brep_index ( j ) == v0 )
      {
        om[i][0]++;
        count++;
      }
      else if ( tv->get_potential_brep_index ( j ) == v1 )
      {
        om[i][1]++;
        count++;
      }
    }
    if ( count == 0 ) // pas la peine de continuer
      return 0;
  }
//     if (count == 2)
//     {
  if ( om[0][0] && om[1][1] )
  {
    *tv0 = te->get_topovertex ( 0 );
    *tv1 = te->get_topovertex ( 1 );
    couple++;
  }
  if ( om[0][1] && om[1][0] )
  {
    *tv0 = te->get_topovertex ( 1 );
    *tv1 = te->get_topovertex ( 0 );
    couple++;
  }
//     }

  if ( count > 2 )
    ambiguous = 1;

  if ( couple )
    return 1;
  else
    return 0;
}

void Topo::update_all_vertex()
{
  while ( !_vertex_inspection_queue.empty() )
  {
    TopoVertex *tv = _vertex_inspection_queue.front();
    _vertex_inspection_queue.pop();
    update_topovertex ( tv );
  }
}

void Topo::update_all_edge()
{
  while ( !_edge_inspection_queue.empty() )
  {
    TopoEdge *te = _edge_inspection_queue.front();
    _edge_inspection_queue.pop();
    update_topoedge ( te );
  }
}

void Topo::update_topoedge ( TopoEdge *te )
{
//     printf("--> update topoedge %d", te->index());

  for ( int i = 0; i < te->topoface_size(); i++ )
  {
    if ( te->get_topoface ( i )->potential_brep_index_size() == 0 )
      te->get_topoface ( i )->setFlag ( -1 );
    if ( te->get_topoface ( i )->getFlag() == -1 )
    {
      topoface_map_remove_surf ( te->get_topoface ( i ) );
      te->remove_topoface ( i-- );
    }
  }

//     printf(" : %d faces --> ", te->topoface_size());
//     for (int i = 0; i < te->topoface_size(); i++)
//         printf(" %d", te->get_topoface(i)->index());
//     printf("\n");

  TopoEdge *ter = te;
  while ( ter->get_ref() != NULL )
    ter = ter->get_ref();

  switch ( te->topoface_size() )
  {
    case 0 :
      assert ( te->getFlag() == -1 );
      return;
    case 1 :
//      printf("--> removing topoface %d\n", te->get_topoface(0)->index());
      assert ( !te->get_topoface ( 0 )->associated() );
      te->get_topoface ( 0 )->setFlag ( -1 );
      te->remove_topoface ( 0 );
      assert ( te->getFlag() == -1 );
      return;
    case 2:
      if ( ter->getFlag() == -1 )
        merge_face ( te );
      return;
    default :
      return;
  }
}

void Topo::update_topovertex ( TopoVertex *tv )
{
//     printf("--> update topovertex %d", tv->index());
  for ( int i = 0; i < tv->topoedge_size(); i++ )
  {
    if ( tv->get_topoedge ( i )->potential_brep_index_size() == 0 )
      tv->get_topoedge ( i )->setFlag ( -1 );
    if ( tv->get_topoedge ( i )->getFlag() == -1 )
    {
      topoedge_map_remove_surf ( tv->get_topoedge ( i ) );
      tv->remove_topoedge ( i-- );
    }
  }
//     printf(" : %d edges --> ", tv->topoedge_size());
//         for (int i = 0; i < tv->topoedge_size(); i++)
//     printf(" %d", tv->get_topoedge(i)->index());
//     printf("\n");

  switch ( tv->topoedge_size() )
  {
    case 0 :
      tv->setFlag ( -1 );
      return;
    case 1 :
//      printf("--> removing topoedge %d\n", tv->get_topoedge(0)->index());
      assert ( !tv->get_topoedge ( 0 )->associated() );
      tv->get_topoedge ( 0 )->setFlag ( -1 );
      tv->remove_topoedge ( 0 );
      tv->setFlag ( -1 );
      return;
    case 2:
      tv->setFlag ( -1 );
      merge_edge ( tv );
      return;
    default :
      return;
  }
}

void Topo::init_clean_edge()
{
  // pretraitement, on supprime les entites ne correspondant a aucune specification topologique
  for ( int i = 0; i < _brm->vertex_size(); i++ )
  {
    std::vector<int> surf;
    std::vector<TopoVertex*> result;
    BRepEntity * bv = _brm->get_vertex ( i );
    for ( int j = 0; j < bv->surf_size(); j++ )
      surf.push_back ( bv->get_surf ( j ) );
    topovertex_query ( surf, result );
    for ( int j = 0; j < result.size(); j++ )
    {
      result[j]->add_potential_brep_index ( bv->index() );
      _entity_to_vertex.insert ( std::pair<BRepEntity *, TopoVertex*> ( bv, result[j] ) );
    }
  }

  for ( int i = 0; i < _brm->edge_size(); i++ )
  {
    std::vector<int> surf;
    std::vector<TopoEdge*> result;
    BRepEntity * be = _brm->get_edge ( i );
    for ( int j = 0; j < be->surf_size(); j++ )
      surf.push_back ( be->get_surf ( j ) );
    topoedge_query ( surf, result );
    for ( int j = 0; j < result.size(); j++ )
    {
      result[j]->add_potential_brep_index ( be->index() );
      _entity_to_edge.insert ( std::pair<BRepEntity *, TopoEdge*> ( be, result[j] ) );
    }
  }

  // traitement inverse --> verifie les incompatibilites
  for ( int i = 0; i < topoedge_size() ; i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->surf_size() > 2 )
      check_compatibility ( te );
  }

  for ( int i = 0; i < topoedge_size() ; i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( !te->potential_brep_index_size() )
      te->setFlag ( -1 );
    else if ( te->imposed() )
      te->setFlag ( -2 );
    else
      te->setFlag ( 1 );
  }

  for ( int i = 0; i < topovertex_size() ; i++ )
  {
    TopoVertex *tv = get_topovertex ( i );
    if ( !tv->potential_brep_index_size() )
    {
      tv->setFlag ( -1 );
      to_be_updated ( tv );
    }
    else if ( tv->imposed() )
      tv->setFlag ( -2 );
    else
      tv->setFlag ( 1 );
  }

//     dumpVTK("init");

  update_all_vertex();
}

void Topo::init_clean_face()
{
  for ( int i = 0; i < topoface_size() ; i++ )
    assert ( get_topoface ( i )->getFlag() == -1 );

  // pretraitement, on supprime les entites ne correspondant a aucune specification topologique
  for ( int i = 0; i < _brm->face_size(); i++ )
  {
    int surf;
    std::vector<TopoFace*> result;

    BRepEntity * be = _brm->get_face ( i );
    surf = be->get_surf ( 0 );
    topoface_query ( surf, result );
//      printf("result size --> %d\n", (int)result.size());
    for ( int j = 0; j < result.size(); j++ )
    {
      result[j]->add_potential_brep_index ( be->index() );
      _entity_to_face.insert ( std::pair<BRepEntity *, TopoFace*> ( be, result[j] ) );
    }
  }

  for ( int i = 0; i < topoface_size() ; i++ )
  {
    TopoFace *tf = get_topoface ( i );
    assert ( tf->potential_brep_index_size() );
    tf->setFlag ( 1 );
    for ( int i = 0; i < tf->topoedge_size(); i++ )
    {
      TopoEdge *te = tf->get_topoedge ( i );
      to_be_updated ( te );
    }
  }

  update_all_edge();
}

void Topo::associate_face ( TopoFace *tf, BRepEntity *bef )
{
//     printf("--> associate topoface %d with brep face %d\n", tf->index(), bef->index());
  assert ( bef != NULL );
  assert ( !tf->associated() );
  assert ( tf->getFlag() == 1 );
  assert ( tf->get_ref() == NULL );

  std::set<TopoFace*> rlist;

  _retef = _entity_to_face.equal_range ( bef );

  std::vector<std::multimap<BRepEntity *, TopoFace*>::iterator> itlist;

  for ( _itef = _retef.first; _itef != _retef.second; ++_itef )
  {
    TopoFace *tef = _itef->second;
    std::vector<TopoFace*> tmplist;
    while ( tef->get_ref() != NULL )
    {
      tmplist.push_back ( tef );
      tef = tef->get_ref();
    }
    tmplist.push_back ( tef );
    if ( tef != tf )
    {
      for ( int i = 0; i < tmplist.size(); i++ )
        rlist.insert ( tmplist[i] );
    }
    if ( _itef->second != tf )
      itlist.push_back ( _itef );
  }

  for ( int i = 0; i < itlist.size(); i++ )
    _entity_to_face.erase ( itlist[i] );

  for ( std::set<TopoFace*>::iterator itset = rlist.begin(); itset != rlist.end(); itset++ )
  {
    TopoFace *tef = *itset;
    tef->remove_potential_brep_index ( bef->index() );
    for ( int i = 0; i < tef->topoedge_size(); i++ )
      if ( tef->get_topoedge ( i )->getFlag() == -1 )
        to_be_updated ( tef->get_topoedge ( i ) );
  }

  update_all_edge();

  tf->set_associated();
  tf->clear_potential_brep_index();
  tf->add_potential_brep_index ( bef->index() );

  assert ( tf->getFlag() == 1 );

  _nb_associated_face++;
//     _lock = 1;
}

void Topo::associate_edge ( TopoEdge *te, BRepEntity *bre, TopoVertex *tv0, BRepEntity *brv0, TopoVertex *tv1, BRepEntity *brv1 )
{
//     printf("--> associate topoedge %d with brep edge %d\n", te->index(), bre->index());
  assert ( bre != NULL );
  assert ( !te->associated() );
//     assert(te->getFlag() != -1);
  assert ( te->getFlag() == 1 );
  assert ( tv0->getFlag() == 1 );
  assert ( tv1->getFlag() == 1 );
  assert ( te->get_ref() == NULL );

  std::set<TopoEdge*> rlist;

  _retee = _entity_to_edge.equal_range ( bre );

  std::vector<std::multimap<BRepEntity *, TopoEdge*>::iterator> itlist;

  for ( _itee = _retee.first; _itee != _retee.second; ++_itee )
  {
    TopoEdge *tec = _itee->second;
    std::vector<TopoEdge*> tmplist;
    while ( tec->get_ref() != NULL )
    {
      tmplist.push_back ( tec );
      tec = tec->get_ref();
    }
    tmplist.push_back ( tec );
    if ( tec != te )
    {
      for ( int i = 0; i < tmplist.size(); i++ )
        rlist.insert ( tmplist[i] );
    }
    if ( _itee->second != te )
      itlist.push_back ( _itee );
  }

  for ( int i = 0; i < itlist.size(); i++ )
    _entity_to_edge.erase ( itlist[i] );

  for ( std::set<TopoEdge*>::iterator itset = rlist.begin(); itset != rlist.end(); itset++ )
  {
    TopoEdge *tec = *itset;
    tec->remove_potential_brep_index ( bre->index() );
    to_be_updated ( tec->get_topovertex ( 0 ) );
    to_be_updated ( tec->get_topovertex ( 1 ) );
  }

  update_all_vertex();

  bre->setFlag ( 1 );
  te->setFlag ( 1 );
  te->set_associated();
  te->clear_potential_brep_index();
  te->add_potential_brep_index ( bre->index() );

//     printf("tv0 : brepindex = %d\n", tv0->get_potential_brep_index(0));
  if ( !tv0->associated() )
    associate_vertex ( tv0, brv0 );
  else
  {
    assert ( tv0->potential_brep_index_size() == 1 );
    assert ( tv0->get_potential_brep_index ( 0 ) == brv0->index() );
  }

//      printf("tv1 : brepindex = %d\n", tv1->get_potential_brep_index(0));
  if ( !tv1->associated() )
    associate_vertex ( tv1, brv1 );
  else
  {
    assert ( tv1->potential_brep_index_size() == 1 );
    assert ( tv1->get_potential_brep_index ( 0 ) == brv1->index() );
  }

  _nb_associated_edge++;
//     _lock = 1;
}

void Topo::associate_vertex ( TopoVertex *tv, BRepEntity *brv )
{
//     printf("--> asssociate topovertex %d with brep vertex %d\n", tv->index(), brv->index());
  assert ( brv != NULL );
  assert ( !tv->associated() );
  assert ( tv->getFlag() == 1 );

  std::vector<std::multimap<BRepEntity *, TopoVertex*>::iterator> itlist;
  // supprime les sommets referencant ce sommet brep
  _retev = _entity_to_vertex.equal_range ( brv );
  for ( _itev=_retev.first; _itev!=_retev.second; ++_itev )
  {
    TopoVertex *tvc = _itev->second;
//      printf("trying to remove index %d to topovertex %d\n", brv->index(), _itev->second->index());
    if ( tvc != tv )
    {
//          printf("remove index %d to vertex %d\n", brv->index(), tvc->index());
      tvc->remove_potential_brep_index ( brv->index() );
      itlist.push_back ( _itev );
    }
  }

  for ( int i = 0; i < itlist.size(); i++ )
    _entity_to_vertex.erase ( itlist[i] );

  brv->setFlag ( 1 );
  tv->set_associated();
  tv->clear_potential_brep_index();
  tv->add_potential_brep_index ( brv->index() );

  _nb_associated_vertex++;
//     _lock = 1;
}

BRepEntity *Topo::get_closest_vertex ( TopoVertex *tv )
{
  double min = INFINITY;
  BRepEntity *cl = NULL;
  for ( int i = 0; i < tv->potential_brep_index_size(); i++ )
  {
    BRepEntity *be = _brm->get_vertex_by_index ( tv->get_potential_brep_index ( i ) );
    std::map<BRepEntity*, LVertex*>::iterator it = _brep_vertex_coord.find ( be );
    assert ( it != _brep_vertex_coord.end() );
    LVertex vv;
    vv.setVector ( tv->get_vertex(), it->second );
    double l = vv.sq_norm();
    if ( l < min )
    {
      min = l;
      cl = be;
    }
  }
  assert ( cl != NULL );
  return cl;
}

TopoVertex *Topo::get_closest_vertex ( BRepEntity *brv )
{
//     printf("looking for vertex closer than brep vertex %d\n", brv->index());
  double min = INFINITY;
  TopoVertex *tv = NULL;
  std::map<BRepEntity*, LVertex*>::iterator it = _brep_vertex_coord.find ( brv );
  assert ( it != _brep_vertex_coord.end() );
  _retev = _entity_to_vertex.equal_range ( brv );
  for ( _itev=_retev.first; _itev!=_retev.second; ++_itev )
  {
    TopoVertex *tvc = _itev->second;
//      printf("testing vertex %d (flag %d)\n", tvc->index(), tvc->getFlag());
    if ( tvc->getFlag() == 1 && !tvc->associated() )
    {
      LVertex vv;
      vv.setVector ( tvc->get_vertex(), it->second );
      double l = vv.sq_norm();
      if ( l < min )
      {
        min = l;
        tv = tvc;
      }
    }
  }
//     assert(tv != NULL);
  return tv;
}

void Topo::loop_detect ( TopoFace *tf )
{
  std::vector<std::vector<TopoEdge*> > edge_loops;
  int offset = -1;
  for ( int i = 0; i < tf->topoedge_size(); i++ )
  {
    TopoEdge *te0 = tf->get_topoedge ( i );
    TopoEdge *te1 = tf->get_topoedge ( ( i+1 ) %tf->topoedge_size() );

    if ( te0->get_topovertex ( 1 ) != te1->get_topovertex ( 0 ) &&
         te0->get_topovertex ( 1 ) != te1->get_topovertex ( 1 ) &&
         te0->get_topovertex ( 0 ) != te1->get_topovertex ( 0 ) &&
         te0->get_topovertex ( 0 ) != te1->get_topovertex ( 1 )
       )
    {
      offset = i+1;
      break;
    }
  }

  if ( offset == -1 ) // no loop
    return;

  edge_loops.resize ( 1 );
  for ( int i = 0; i < tf->topoedge_size(); i++ )
  {
    TopoEdge *te0 = tf->get_topoedge ( ( offset+i ) %tf->topoedge_size() );
    TopoEdge *te1 = tf->get_topoedge ( ( offset+i+1 ) %tf->topoedge_size() );

    edge_loops.back().push_back ( te0 );

    if ( te0->get_topovertex ( 1 ) != te1->get_topovertex ( 0 ) &&
         te0->get_topovertex ( 1 ) != te1->get_topovertex ( 1 ) &&
         te0->get_topovertex ( 0 ) != te1->get_topovertex ( 0 ) &&
         te0->get_topovertex ( 0 ) != te1->get_topovertex ( 1 )
       )
    {
      edge_loops.resize ( edge_loops.size() +1 );
    }
  }

  for ( int i = 0; i < edge_loops.size(); i++ )
  {
    for ( int j = 0; j < edge_loops[i].size(); j++ )
    {
      if ( edge_loops[i][j] == edge_loops[i][ ( j+1 ) %edge_loops[i].size() ] )
        edge_loops[i].erase ( edge_loops[i].begin() +j-- );
    }
  }
//     tf->dump();
  tf->clear_topoedge();
  for ( int i = 0; i < edge_loops.size(); i++ )
  {
    std::vector<TopoEdge*> tmp = edge_loops[i];
    for ( int j = 0; j < tmp.size(); j++ )
    {
      tf->add_topoedge ( tmp[j] );
    }
  }
//     tf->dump();
}

void Topo::final_clean_face()
{
  for ( int i = 0; i < topoface_size(); i++ ) // suppression des faces inutiles
  {
    if ( !get_topoface ( i )->associated() )
      remove_topoface ( i-- );
  }

  for ( int i = 0; i < topoface_size(); i++ ) // on remplace les aretes detruites par fusion par l'arete resultante
  {
    TopoFace *tf = get_topoface ( i );
    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      if ( te->get_ref() != NULL )
      {
        while ( te->get_ref() != NULL )
        {
          assert ( te->getFlag() == -1 );
          te = te->get_ref();
        }

        tf->remove_topoedge ( j );
        tf->insert_topoedge ( j, te, 1 );
      }
    }

    loop_detect ( tf );

    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      if ( tf->get_topoedge ( j ) == tf->get_topoedge ( ( j+1 ) %tf->topoedge_size() ) )
        tf->remove_topoedge ( j-- );
    }
  }

  for ( int i = 0; i < topoface_size(); i++ ) // supprime a priori les faces donc une arete est non valide
  {
    TopoFace *tf = get_topoface ( i );
//      tf->dump();
    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      if ( te->getFlag() == -1 )
      {
        remove_topoface ( i-- );
        break;
      }
    }
  }

  for ( int i = 0; i < topoedge_size(); i++ )
    if ( get_topoedge ( i )->getFlag() == -1 )
      remove_topoedge ( i-- );
}

int Topo::find_edge_association ( BRepEntity *bee, int nb_clone, std::vector<TopoEdge *> te_result, int &valid )
{
  // tri des resultats
  int v0 = bee->get_adj ( 0 );
  int v1 = bee->get_adj ( 1 );
  BRepEntity *bev0 = _brm->get_vertex_by_index ( v0 );
  BRepEntity *bev1 = _brm->get_vertex_by_index ( v1 );
  std::set<TopoEdge *> count_edge;

  TopoEdge *tefound = NULL;
  TopoVertex *tv0 = NULL;
  TopoVertex *tv1 = NULL;
  int couple;
  int ambiguous;
  int univocal;

  valid = 0;

  for ( int i = 0; i < te_result.size(); i++ )
  {
    TopoEdge *te = te_result[i];

    if ( te->getFlag() == 1 )
    {
//             printf("--> result %d = topoedge %d\n", i, te->index());
      assert ( te->get_ref() == NULL );
//             te->dump();

      if ( find_edge_occurrence ( te, bev0, bev1, &tv0, &tv1, couple, ambiguous ) )
      {
//                 printf("  --> OK\n");
        valid += couple;
        tefound = te;
        if ( tv0->potential_brep_index_size() == 1 || tv1->potential_brep_index_size() == 1 )
          univocal = 1;
        else
          univocal = 0;
      }
    }
  }
//     printf("count = %d, nb_clone = %d, univocal = %d\n", valid, nb_clone, univocal);
  if ( valid == 1 && ( ( nb_clone == 1 ) || ( nb_clone > 1 && univocal ) ) )  // association, sinon indetermination
  {
    if ( tv0->getFlag() == 1 && tv1->getFlag() == 1 )
    {
//         printf("found valid edge --> %d +--> associating ...\n", tefound->index());
      associate_edge ( tefound, bee, tv0, bev0, tv1, bev1 );
      return 1;
    }
  }
  return 0;
}

void Topo::find_best_candidate ( BRepEntity *bee, std::vector<TopoEdge *> te_result, TopoEdge **te, TopoVertex **tv0, TopoVertex **tv1 )
{
//     printf("--> find candidates\n");

  std::map<BRepEntity *, LVertex*>::iterator it;
  int ibev0 = bee->get_adj ( 0 );
  int ibev1 = bee->get_adj ( 1 );

  BRepEntity *bev0 = _brm->get_vertex_by_index ( ibev0 );
  BRepEntity *bev1 = _brm->get_vertex_by_index ( ibev1 );

  it = _brep_vertex_coord.find ( bev0 );
  assert ( it != _brep_vertex_coord.end() );
  LVertex *v0 = it->second;

  it = _brep_vertex_coord.find ( bev1 );
  assert ( it != _brep_vertex_coord.end() );
  LVertex *v1 = it->second;

//     printf("edge %d --> bev0 = %d -- bev1 = %d\n", bee->index(), ibev0, ibev1);

  *te = NULL;
  double min = INFINITY;

  for ( int i = 0; i < te_result.size(); i++ )
  {
    TopoEdge *temp = te_result[i];
    if ( temp->getFlag() != -1 )
    {
//             printf("--> testing topoedge %d\n", temp->index());
//          temp->dump();
      int om[2][2] = {{0, 0}, {0, 0}};
      int count = 0;
      for ( int j = 0; j < 2; j++ )
      {
        TopoVertex *tv = temp->get_topovertex ( j );
        for ( int k = 0; k < tv->potential_brep_index_size(); k++ )
        {
          if ( tv->get_potential_brep_index ( k ) == ibev0 )
          {
            om[j][0]++;
            count++;
          }
          else if ( tv->get_potential_brep_index ( k ) == ibev1 )
          {
            om[j][1]++;
            count++;
          }
        }
      }

      if ( count == 2 )
      {
        TopoVertex *t0 = temp->get_topovertex ( 0 );
        TopoVertex *t1 = temp->get_topovertex ( 1 );
        LVertex vv;

        if ( om[0][0] && om[1][1] )
        {
          double min1 = INFINITY;
          double min2 = INFINITY;
          vv.setVector ( t0->get_vertex(), v0 );
          min1 = vv.sq_norm();
          vv.setVector ( t1->get_vertex(), v1 );
          min2 = vv.sq_norm();
//                     printf("testing min (%f) = %f + %f\n", min, min1, min2);
          if ( min1+min2 < min )
          {
            min = min1+min2;
            *te = temp;
            *tv0 = temp->get_topovertex ( 0 );
            *tv1 = temp->get_topovertex ( 1 );
          }
        }
        if ( om[0][1] && om[1][0] )
        {
          double min1 = INFINITY;
          double min2 = INFINITY;
          vv.setVector ( t0->get_vertex(), v1 );
          min1 = vv.sq_norm();
          vv.setVector ( t1->get_vertex(), v0 );
          min2 = vv.sq_norm();
//                     printf("testing min (%f) = %f + %f\n", min, min1, min2);
          if ( min1+min2 < min )
          {
            min = min1+min2;
            *te = temp;
            *tv0 = temp->get_topovertex ( 1 );
            *tv1 = temp->get_topovertex ( 0 );
          }
        }
      }
    }
  }
}

int Topo::choose_edge ( std::vector<BRepEntity*> be_list, std::vector<TopoEdge *> te_result )
{
//      printf("--> choosing edges\n");

//     for (int i = 0; i < be_list.size(); i++)
//     {
//         BRepEntity *be0 = be_list[i];
//      printf("brep edge %d : v %d -- v %d\n", be0->index(), be0->get_adj(0), be0->get_adj(1));
//     }

  TopoEdge *te = NULL;
  TopoVertex *tv0 = NULL;
  TopoVertex *tv1 = NULL;

  int success = 0;
  std::vector<TopoEdge *> te_find;
  std::vector<TopoVertex *> tv0_find;
  std::vector<TopoVertex *> tv1_find;

  for ( int i = 0; i < be_list.size(); i++ )
  {
    BRepEntity *bre = be_list[i];

    int ibev0 = bre->get_adj ( 0 );
    int ibev1 = bre->get_adj ( 1 );
    BRepEntity *bev0 = _brm->get_vertex_by_index ( ibev0 );
    BRepEntity *bev1 = _brm->get_vertex_by_index ( ibev1 );

    find_best_candidate ( bre, te_result, &te, &tv0, &tv1 );

    if ( te != NULL )
    {
      te->setFlag ( 1 );
      tv0->setFlag ( 1 );
      tv1->setFlag ( 1 );
//             printf("found topoedge --> %d\n", te->index());
      associate_edge ( te, bre, tv0, bev0, tv1, bev1 );
      success = 1;
    }
  }

  return success;
}

void Topo::check_compatibility ( TopoEdge * te )
{
  if ( te->surf_size() > 2 && te->get_topovertex ( 0 )->getFlag() == -1 && te->get_topovertex ( 1 )->getFlag() == -1 )
  {
//      printf("--> checking te %d, surf size %d\n", te->index(), te->surf_size());
    for ( int i = 0; i < te->surf_size()-1; i++ )
      for ( int j = i+1; j < te->surf_size(); j++ )
      {
        std::vector<int> surf;
        surf.push_back ( te->get_surf ( i ) );
        surf.push_back ( te->get_surf ( j ) );
        std::sort ( surf.begin(), surf.end() );
//              printf("surf %d -- %d\n", te->get_surf(i), te->get_surf(j));
        std::vector<BRepEntity*> result;
        _brm->query ( surf, result );
        if ( !result.size() )
        {
          te->setFlag ( -1 );
          te->clear_potential_brep_index();
          to_be_updated ( te->get_topovertex ( 0 ) );
          to_be_updated ( te->get_topovertex ( 1 ) );
          return;
        }
      }
  }
}

void Topo::delete_imposed ( TopoEdge *te )
{
//     printf("--> unlock topoedge %d\n", te->index());
//     te->dump();
  assert ( te->getFlag() == -2 );
  TopoVertex *tv0 = te->get_topovertex ( 0 );
  TopoVertex *tv1 = te->get_topovertex ( 1 );
  te->setFlag ( -1 );
  if ( !tv0->associated() )
  {
    assert ( tv0->getFlag() == -2 );
    tv0->setFlag ( 1 );
    to_be_updated ( tv0 );
  }
  if ( !tv1->associated() )
  {
    assert ( tv1->getFlag() == -2 );
    tv1->setFlag ( 1 );
    to_be_updated ( tv1 );
  }
}

void Topo::clean_query_results ( std::vector<TopoEdge *> &te_result )
{
  for ( int i = 0; i < te_result.size(); i++ )
  {
    TopoEdge *te = te_result[i];
    TopoVertex *tv0 = te->get_topovertex ( 0 );
    TopoVertex *tv1 = te->get_topovertex ( 1 );
    if ( tv0->getFlag() != -1 )
      to_be_updated ( tv0 );
    if ( tv1->getFlag() != -1 )
      to_be_updated ( tv1 );
  }

  update_all_vertex();

  for ( int i = 0; i < te_result.size(); i++ )
  {
    if ( te_result[i]->getFlag() == -1 || te_result[i]->associated() )
      te_result.erase ( te_result.begin() +i-- );
  }
}

int Topo::try_edge_association ( BRepEntity *bee )
{
  std::vector<int> surf;
  std::vector<TopoEdge *> te_result;
  std::vector<BRepEntity *> be_result;

  printf ( "## BRep topoedge %d (v %d -- v %d) --> surf", bee->index(), bee->get_adj ( 0 ), bee->get_adj ( 1 ) );
  for ( int i = 0; i < bee->surf_size(); i++ )
  {
    printf ( " %d", bee->get_surf ( i ) );
    surf.push_back ( bee->get_surf ( i ) );
  }
  printf ( "\n" );

  topoedge_query ( surf, te_result );
  _brm->query ( bee, be_result );

  printf ( "te_query --> %d result\n", ( int ) te_result.size() );
  printf ( "be_query --> %d result\n", ( int ) be_result.size() );

  assert ( te_result.size() >= be_result.size() );

  clean_query_results ( te_result );

  int assoc_found;
  int success = find_edge_association ( bee, be_result.size(), te_result, assoc_found );

//     printf("assoc found = %d -- te result size = %d\n", assoc_found, (int)te_result.size());

  if ( success )
  {
    return 1;
  }
  else if ( assoc_found >/*=*/ be_result.size() ) // les arete imposees sont inutiles --> on debloque les topovertex pour fusion
  {
//      printf("--> unlock topoedge\n");
    for ( int i = 0; i < te_result.size(); i++ )
    {
      TopoEdge *te = te_result[i];
      if ( te->getFlag() != -1 )
        if ( te->imposed() && !te->associated() ) // unlock vertices
          delete_imposed ( te );
    }
    update_all_vertex();
  }
  else if ( ( assoc_found < be_result.size() ) || ( assoc_found > 1 && assoc_found == be_result.size() ) ) // on ne trouve pas assez d'assoc potentielle --> on debloque les arete imposees
  {
//      printf("--> insert topoedge\n");
    success = choose_edge ( be_result, te_result );
  }
  return success;
}

void Topo::process_edges()
{
  // construire une liste d'edge a associer
  std::queue<BRepEntity *> tequeue;
  for ( int i = 0; i < _brm->edge_size() ; i++ )
  {
    _brm->get_edge ( i )->setFlag ( -1 );
    tequeue.push ( _brm->get_edge ( i ) );
  }

  int loop = 0;
  int max_loop = 20;

  while ( !tequeue.empty() )
  {
    int queue_size = tequeue.size();
//      _lock = 0;

    for ( int i = 0; i < queue_size; i++ )
    {
      BRepEntity *bee = tequeue.front();
      tequeue.pop();
//             bee->setFlag(-1);
      if ( bee->getFlag() == -1 )
      {
        try_edge_association ( bee );
//                 printf("--> flag = %d -- %d associated edges -- queue size = %d\n", bee->getFlag(), _nb_associated_edge, (int)tequeue.size());
        if ( bee->getFlag() == -1 )
          tequeue.push ( bee );
      }
    }

    tprintf ( "--> %d associated edges\n", _nb_associated_edge );
    tprintf ( "--> %d not associated edges\n", ( int ) tequeue.size() );

    dumpVTK ( "inter" );

    if ( loop++ > max_loop )
    {
      printf ( "--> /!\\ <-- Cannot find edges !!\n" );
      printf ( "Exiting ...\n" );
      exit ( 1 );
    }

//         getchar();
  }

  dumpVTK ( "inter" );

  // suppression des topovertex inutiles
  for ( int i = 0; i < topovertex_size(); i++ )
    if ( get_topovertex ( i )->potential_brep_index_size() == 0 )
      remove_topovertex ( i-- );
}

int Topo::find_face_association ( BRepEntity *bef, std::vector<TopoFace *> tf_result )
{
//     for (int i = 0; i < tf_result.size(); i++)
//     {
//         TopoFace *tf = tf_result[i];
//         tf->dump();
//     }
  for ( int i = 0; i < tf_result.size(); i++ )
  {
    if ( tf_result[i]->getFlag() == -1 )
      tf_result.erase ( tf_result.begin() +i-- );
  }
  for ( int i = 0; i < tf_result.size(); i++ )
  {
    TopoFace *tf = tf_result[i];
    std::map<int, int> index_map;
    std::map<int, int>::iterator it;
    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      while ( te->get_ref() != NULL )
        te = te->get_ref();
      if ( te->associated() )
      {
        assert ( te->potential_brep_index_size() == 1 );
        int brep_edge_index = te->get_potential_brep_index ( 0 );
        it = index_map.find ( brep_edge_index );
        if ( it == index_map.end() )
          index_map.insert ( std::pair<int, int> ( brep_edge_index, 1 ) );
      }
      else
      {
        index_map.clear();
      }
    }
    if ( index_map.size() == bef->adj_size() )
    {
      for ( int j = 0; j < bef->adj_size(); j++ )
      {
        it = index_map.find ( bef->get_adj ( j ) );
        if ( it == index_map.end() )
          break;
        else
          it->second++;
      }
      int valid = 1;
      for ( it = index_map.begin(); it != index_map.end() && valid; it++ )
      {
        if ( it->second != 2 )
          valid = 0;
      }
      if ( valid )
      {
        associate_face ( tf, bef );
        return 1;
      }

    }
  }
  return 0;
}

int Topo::try_face_association ( BRepEntity *bef )
{
  // compute key
//     printf("brep face %d surf %d --> %d edges :", bef->index(), bef->get_surf(0), bef->adj_size());

  std::vector<TopoFace *> tf_result;

  topoface_query ( bef->get_surf ( 0 ), tf_result );

//     printf("--> found %d results\n", (int)tf_result.size());

  int success = 0;
  success = find_face_association ( bef, tf_result );
  return success;
}


void Topo::process_faces()
{
  // construire une liste de faces a associer
  std::queue<BRepEntity *> tfqueue;
  for ( int i = 0; i < _brm->face_size() ; i++ )
    tfqueue.push ( _brm->get_face ( i ) );

  int loop = 0;
  int max_loop = 20;

  while ( !tfqueue.empty() )
  {
    int queue_size = tfqueue.size();
//         _lock = 0;

    for ( int i = 0; i < queue_size; i++ )
    {
      BRepEntity *bef = tfqueue.front();
      tfqueue.pop();
      bef->setFlag ( -1 );

      if ( try_face_association ( bef ) )
        bef->setFlag ( 1 );
      else
        tfqueue.push ( bef );
//          tprintf("--> %d associated faces\n", _nb_associated_face);
//          dumpVTK("inter");
    }

    tprintf ( "--> %d associated faces\n", _nb_associated_face );
    tprintf ( "--> %d not associated faces\n", ( int ) tfqueue.size() );
//         tprintf("--> difference = %d\n", queue_size-(int)tfqueue.size());

    if ( loop++ > max_loop )
    {
      printf ( "--> /!\\ <-- Cannot find faces !!\n" );
      printf ( "Exiting ...\n" );
      exit ( 1 );
    }
#if 0
    dumpVTK ( "inter" );
#endif
//      getchar();
  }

  final_clean_face();

#if 1
  dumpVTK ( "inter" );
#endif
}


void Topo::clean_topo()
{
  tprintf ( "Topology cleaning START\n" );
  init_comparison();
  init_clean_edge();
//     dumpVTK("init");
  // deuxieme etape --> association/suppression/fusion edges
  process_edges();
  // troisieme etape --> association/suppression/fusion faces
  init_clean_face();

#if 0
  dumpVTK ( "inter" );
#endif

  process_faces();
  dump();
#if 1
  dumpVTK ( "final" );
#endif

  tprintf ( "Topology cleaning END\n" );
}

void Topo::build_volume ( Mesh *m )
{
  std::vector<LTetra *> bt;
  int unfilled_color = -1;
//     int ext_color = -2;



#if 0
  int count_vol;
  //  id of the unique volume (from original brep info)
  // if there is no volume or more than one volume, the id is set to 1
  //        tprintf("volsize %d\n",_brm->volume_size());
  int idvol=1;

  if ( _brm )
    idvol=_brm->volume_size() ==1?_brm->get_volume ( 0 )->index() :1;

//     tprintf("idvol %d\n",idvol);

  int nextidvol=idvol+1;
  mpvols[count_vol]=idvol;
#endif

  int current_vol_id = -1;

  tprintf ( "Volume building START\n" );

  for ( int i = 0; i < m->tetra_size(); i++ )
    m->get_tetra ( i )->setVol ( unfilled_color );

  // suppression des adjacences aux interfaces
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    for ( int j = 0; j < tf->face_size(); j++ )
    {
      CutFace *cf = tf->get_face ( j );
      if ( cf->ext_tetra() )
        cf->set_no_adj();
      if ( _opts->brep_topo || _opts->auto_topo )
      {
        bt.push_back ( cf->int_tetra() );
      }
      else if ( _opts->recut_enable )
      {
        if ( _em->is_restore_face ( tf->get_surf ( 0 ) ) )
          bt.push_back ( cf->int_tetra() );
      }
    }
  }

  tprintf ( "--> %d border tetra\n", ( int ) bt.size() );

  int count_vol = 0;

  for ( int i = 0; i < bt.size(); i++ )
  {
    LTetra *t = bt[i];
    assert ( !t->isParent() );
    if ( ( t->getVol() == unfilled_color ) )
    {
      m->seed ( t );
      current_vol_id = _em->add_vol_id();
      dbgprintf ( 2, "local vol id --> %d\n", current_vol_id );
      int nb = m->color ( current_vol_id, unfilled_color );
      if ( nb )
      {
        TopoVolume *tvol = new_topovolume ( count_vol, current_vol_id );
        tvol->set_seed ( t );
        if ( _opts->brep_topo )
          tvol->add_potential_brep_index ( current_vol_id );
        count_vol++;
//                 mpvols[count_vol]=nextidvol++;
      }
    }
  }

  tprintf ( "--> %d internal volume\n", count_vol );

  if ( count_vol == 0 )
  {
    tprintf ( "--> nothing to build !\n" );
    exit ( 0 );
  }

  // coloration externe
  int nb_ext_tetra = 0;
  for ( int i = 0; i < m->tetra_size(); i++ )
  {
    LTetra *t = m->get_tetra ( i );
    assert ( !t->isParent() );
    if ( t->getVol() == unfilled_color )
    {
      m->seed ( t );
      current_vol_id = _em->add_vol_id();
      dbgprintf ( 2, "local vol id --> %d\n", current_vol_id );
      int nb = m->color ( current_vol_id, unfilled_color );
      if ( nb )
      {
        TopoVolume *tvol = new_topovolume ( count_vol, current_vol_id );
        tvol->set_seed ( t );
        if ( _opts->brep_topo )
          tvol->add_potential_brep_index ( current_vol_id );
        nb_ext_tetra += nb;
        count_vol++;
//                 mpvols[count_vol]=nextidvol++;

      }
    }
  }

  tprintf ( "--> %d external tetrahedra\n", nb_ext_tetra );

  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    CutFace *cf = tf->get_face ( 0 );
    if ( ( cf->ext_tetra() && cf->int_tetra()->getVol() != cf->ext_tetra()->getVol() ) || !cf->ext_tetra() )
    {
      TopoVolume *tvol = NULL;
      int int_vol_id = cf->int_tetra()->getVol();
      tvol = get_topovolume_by_color ( int_vol_id );
      tvol->add_topoface ( tf );
      tf->set_inside_topovolume ( tvol );

      if ( tf->get_parent() )
      {
        TopoVolume *root_tvol = tf->get_parent()->get_inside_topovolume();
        root_tvol->add_child ( tvol );
        tvol->set_parent ( root_tvol );
      }

      if ( cf->ext_tetra() )
      {
        int ext_vol_id = cf->ext_tetra()->getVol();
        tvol = get_topovolume_by_color ( ext_vol_id );
        tvol->add_topoface ( tf );
        tf->set_outside_topovolume ( tvol );
      }
    }
  }

  // retablissement des adjacences
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    for ( int j = 0; j < tf->face_size(); j++ )
    {
      CutFace *cf = tf->get_face ( j );
      if ( cf->ext_tetra() )
        cf->set_adj();
    }
  }

//     for (int i = 0; i < topovolume_size(); i++)
//     {
//         TopoVolume *tvol = get_topovolume(i);
//         printf("topovolume %d parent %p --> %d --> %d faces\n", i, tvol->get_parent(), tvol->index(), tvol->face_size());
//         for (int j = 0; j < tvol->face_size(); j++)
//         {
//             printf("%s of topoface %d\n", tvol->get_face_orientation(j) ? "ext" : "int", tvol->get_face(j)->index());
//         }
//     }

  m->writeMeshVTK ( "final" );

  tprintf ( "Volume building END\n" );
}

void Topo::dump()
{
  printf ( "Topology :\n" );
  printf ( "%d faces \n", topoface_size() );
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    printf ( "Face %d (color %d) --> %d edges\n", i, tf->color(), tf->topoedge_size() );
    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      printf ( "---> Edge %d (color %d) (v0 %d -- v1 %d) --> %d vertex\n", j, te->color(), te->get_topovertex ( 0 )->index(), te->get_topovertex ( 1 )->index(), te->edge_size() );
    }
  }
//    getchar();
}

void Topo::export_dot_graph ( char *name )
{
  char filename[100];
  sprintf ( filename, "%s_edge_graph.dot", name );
  tprintf ( "Writing entities file %s ....\n", filename );
  FILE *fp = fopen ( filename, "w" );

  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    fprintf ( fp, "V_%d -> V_%d [ label = \"E_%d\"  ];\n", te->get_topovertex ( 0 )->index(), te->get_topovertex ( 1 )->index(), te->index() );
  }
  fclose ( fp );
}

void Topo::export_existing_entities ( char *name )
{
  char filename[100];
  sprintf ( filename, "%s_cut.ent", name );
  tprintf ( "Writing entities file %s ....\n", filename );
  FILE *fp = fopen ( filename, "w" );
  fprintf ( fp, "%d\n", topoface_size() ); // nb faces
  std::set<LVertex*> vset;
  for ( int i = 0; i < topoface_size(); i++ )
  {
    vset.clear();
    TopoFace *tf = get_topoface ( i );
    int real_id = -1;
    int surf = -1;
    if ( _opts->brep_topo )
    {
      real_id = tf->get_potential_brep_index ( 0 );
      surf = _em->get_lookup_face_id ( real_id );
    }
    else if ( _opts->auto_topo )
      surf = real_id = tf->get_surf ( 0 );

    for ( int j = 0; j < tf->face_size(); j++ )
    {
      CutFace *cf = tf->get_face ( j );
      vset.insert ( cf->get_face_vertex ( 0 ) );
      vset.insert ( cf->get_face_vertex ( 1 ) );
      vset.insert ( cf->get_face_vertex ( 2 ) );
      vset.insert ( cf->int_vertex() );
      if ( cf->ext_vertex() )
        vset.insert ( cf->ext_vertex() );
    }
    fprintf ( fp, "%d %d\n", real_id, ( int ) vset.size() );
    for ( std::set<LVertex*>::iterator it = vset.begin(); it != vset.end(); it++ )
    {
      fprintf ( fp, " %d %lf", ( *it )->getIndex(), ( *it )->surfVal ( surf ) );
    }
    fprintf ( fp, "\n" );
  }

  fprintf ( fp, "%d\n", topoedge_size() ); // nb edges
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );

    int real_idedge = -1;
    if ( _opts->brep_topo )
      te->get_potential_brep_index ( 0 );
    else if ( _opts->auto_topo )
      te->color();
//      int real_id0 = te->get_topoface(0)->get_potential_brep_index(0);
//         int real_id1 = te->get_topoface(1)->get_potential_brep_index(0);

    fprintf ( fp, "%d %d\n", real_idedge, te->vertex_size() );
    for ( int j = 0; j < te->vertex_size(); j++ )
      fprintf ( fp, " %d", te->get_vertex ( j, 1 )->getIndex() );
    fprintf ( fp, "\n" );
  }
//
//     for (int i = 0; i < topovertex_size(); i++)
//     {
//        std::set<int> sset;
//         TopoVertex *tv = get_topovertex(i);
//      sset.insert(tv->get_topoedge(0)->get_topoface(0)->get_potential_brep_index(0));
//      sset.insert(tv->get_topoedge(0)->get_topoface(1)->get_potential_brep_index(0));
//      sset.insert(tv->get_topoedge(1)->get_topoface(0)->get_potential_brep_index(0));
//      sset.insert(tv->get_topoedge(1)->get_topoface(1)->get_potential_brep_index(0));
//      sset.insert(tv->get_topoedge(2)->get_topoface(0)->get_potential_brep_index(0));
//      sset.insert(tv->get_topoedge(2)->get_topoface(1)->get_potential_brep_index(0));
//      std::set<int>::iterator it = sset.begin();
//         fprintf (fp, "%d %d %d %d %d\n", tv->get_potential_brep_index(0), tv->get_vertex()->getIndex(), (*it++), (*it++), (*it));
//     }
  fclose ( fp );
}

void Topo::dumpVTK ( const char *name )
{
  if ( !_opts->debug )
    return;

  char filename[100];
  sprintf ( filename, "%s_%s_topoface.vtk", _opts->basename, name );
  tprintf ( "Writing %s ....\n", filename );
  FILE *fp = fopen ( filename, "w" );
  fprintf ( fp, "# vtk DataFile Version 2.0\n" );
  fprintf ( fp, "Really cool data\n" );
  fprintf ( fp, "ASCII\n" );
  fprintf ( fp, "DATASET POLYDATA\n" );
  ////////////////////////////
  // map vertices
  int nb_v = 0;
  int nb_e = 0;
  int nb_f = 0;
  int nb_f_v = 0;
  long seek;
  // dump points
  fprintf ( fp, "POINTS " );
  seek = ftell ( fp );
  fprintf ( fp, "%10d double\n", 0 );
  std::map<int, int > map_v;
  std::map<int, int >::iterator it;
  std::vector<LVertex *> vx;
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    if ( tf->getFlag() == 1 )
    {
      for ( int j = 0; j < tf->face_size(); j++ )
      {
        CutFace *cf = tf->get_face ( j );
        nb_f_v += 3;
        nb_f++;
        for ( int k = 0; k < 3; k++ )
        {
          LVertex *v = cf->int_tetra()->getFaceVertex ( cf->int_face(), k );
          it = map_v.find ( v->getIndex() );
          if ( it == map_v.end() )
          {
            map_v.insert ( std::pair<int, int> ( v->getIndex(), nb_v++ ) );
            fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
            vx.push_back ( v );
          }
        }
      }
    }
  }

  fseek ( fp, seek, SEEK_SET );
  fprintf ( fp, "%10d ", nb_v );
  fseek ( fp, 0, SEEK_END );

  //dump polygons
  fprintf ( fp, "POLYGONS %d %d\n", nb_f, nb_f_v+nb_f );

  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    if ( tf->getFlag() == 1 )
    {
      for ( int j = 0; j < tf->face_size(); j++ )
      {
        CutFace *cf = tf->get_face ( j );
        fprintf ( fp, "%d ", 3 );
        for ( int k = 0; k < 3; k++ )
        {
          LVertex *v = cf->int_tetra()->getFaceVertex ( cf->int_face(), k );
          std::map<int, int >::iterator iter = map_v.find ( v->getIndex() );
          fprintf ( fp, "%d ", iter->second );
        }
        fprintf ( fp, "\n" );
      }
    }
  }

  fprintf ( fp, "CELL_DATA %d\n", nb_f );
  fprintf ( fp, "SCALARS face_color int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    if ( tf->getFlag() == 1 )
    {
      for ( int j = 0; j < tf->face_size(); j++ )
      {
        CutFace *cf = tf->get_face ( j );
        fprintf ( fp, "%d\n", cf->get_color() );
      }
    }
  }

  fprintf ( fp, "SCALARS face_index int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    if ( tf->getFlag() == 1 )
    {
      for ( int j = 0; j < tf->face_size(); j++ )
      {
        CutFace *cf = tf->get_face ( j );
        fprintf ( fp, "%d\n", tf->index() );
      }
    }
  }

  fprintf ( fp, "SCALARS face_associate int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );
    if ( tf->getFlag() == 1 )
    {
      for ( int j = 0; j < tf->face_size(); j++ )
      {
        CutFace *cf = tf->get_face ( j );
        fprintf ( fp, "%d\n", tf->associated() );
      }
    }
  }

  fclose ( fp );

  sprintf ( filename, "%s_%s_topoedge.vtk", _opts->basename, name );
  tprintf ( "Writing %s ....\n", filename );
  fp = fopen ( filename, "w" );
  fprintf ( fp, "# vtk DataFile Version 2.0\n" );
  fprintf ( fp, "Really cool data\n" );
  fprintf ( fp, "ASCII\n" );
  fprintf ( fp, "DATASET POLYDATA\n" );

  // map vertices
  nb_v = 0;
  nb_e = 0;
  nb_f_v = 0;
  // dump points
  fprintf ( fp, "POINTS " );
  seek = ftell ( fp );
  fprintf ( fp, "%10d double\n", 0 );
  map_v.clear();;
  vx.clear();
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->getFlag() != -1 )
    {
      nb_e += te->edge_size();
      for ( int j = 0; j < te->vertex_size(); j++ )
      {
        LVertex *v = te->get_vertex ( j, 1 );
        it = map_v.find ( v->getIndex() );
        if ( it == map_v.end() )
        {
          map_v.insert ( std::pair<int, int> ( v->getIndex(), nb_v++ ) );
          fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
          vx.push_back ( v );
        }
      }
    }
  }

  fseek ( fp, seek, SEEK_SET );
  fprintf ( fp, "%10d ", nb_v );
  fseek ( fp, 0, SEEK_END );

//     fprintf (fp, "VERTICES %d %d\n", (int)vx.size(), (int)vx.size()*2);
//
//     for (int i = 0; i < vx.size(); i++)
//         fprintf (fp, "%d %d\n", 1, i);

  fprintf ( fp, "LINES %d %d\n", nb_e, nb_e*3 );

  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->getFlag() != -1 )
    {
      for ( int j = 0; j < te->edge_size(); j++ )
      {
        CutEdge* ce = te->get_edge ( j );
        fprintf ( fp, "%d ", 2 );
        LVertex *v0 = ce->get_vertex ( 0 );
        std::map<int, int >::iterator iter;
        iter = map_v.find ( v0->getIndex() );
        assert ( iter != map_v.end() );
        fprintf ( fp, "%d ", iter->second );
        LVertex *v1 = ce->get_vertex ( 1 );
        iter = map_v.find ( v1->getIndex() );
        assert ( iter != map_v.end() );
        fprintf ( fp, "%d ", iter->second );
        fprintf ( fp, "\n" );
      }
    }
  }

  fprintf ( fp, "CELL_DATA %d\n", nb_e );
  fprintf ( fp, "SCALARS edge_color int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->getFlag() != -1 )
    {
      for ( int j = 0; j < te->edge_size(); j++ )
      {
        CutEdge* ce = te->get_edge ( j );
//                 fprintf (fp, "%d\n", ce->get_color());
        fprintf ( fp, "%d\n", te->index() );
      }
    }
  }
  fprintf ( fp, "SCALARS edge_index int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->getFlag() != -1 )
    {
      for ( int j = 0; j < te->edge_size(); j++ )
      {
        CutEdge* ce = te->get_edge ( j );
        fprintf ( fp, "%d\n", te->index() );
      }
    }
  }
  fprintf ( fp, "SCALARS edge_associate int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->getFlag() != -1 )
    {
      for ( int j = 0; j < te->edge_size(); j++ )
      {
        CutEdge* ce = te->get_edge ( j );
        fprintf ( fp, "%d\n", te->associated() );
      }
    }
  }
  fprintf ( fp, "SCALARS edge_brep_index int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < topoedge_size(); i++ )
  {
    TopoEdge *te = get_topoedge ( i );
    if ( te->getFlag() != -1 )
    {
      for ( int j = 0; j < te->edge_size(); j++ )
      {
        CutEdge* ce = te->get_edge ( j );
        fprintf ( fp, "%d\n", te->get_potential_brep_index ( 0 ) );
      }
    }
  }
  for ( int i = 0; i < _em->face_size(); i++ )
  {
    int surf = _em->get_face_id ( i );
    fprintf ( fp, "SCALARS edge_surf_%d int 1\n", surf );
    fprintf ( fp, "LOOKUP_TABLE default\n" );
    for ( int j = 0; j < topoedge_size(); j++ )
    {
      TopoEdge *te = get_topoedge ( j );
      if ( te->getFlag() != -1 )
      {
        int found = 0;
        for ( int k = 0; k < te->surf_size() && !found; k++ )
          if ( te->get_surf ( k ) == surf )
            found = 1;
        for ( int k = 0; k < te->edge_size(); k++ )
        {
          fprintf ( fp, "%d\n", found );
        }
      }
    }
  }
  fclose ( fp );
}


void Topo::dump_bohren()
{
  char filename[100];
  sprintf ( filename, "surface_tri.surf" );
  tprintf ( "Writing %s ....\n", filename );
  FILE *fp = fopen ( filename, "w" );

  std::vector<LVertex *> vx;
  int nb_v = 0;
  std::map<int, int > map_v;
  std::map<int, int >::iterator it;

  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );

    for ( int j = 0; j < tf->face_size() ; j++ )
    {
      CutFace *cf = tf->get_face ( j );
      for ( int k = 0; k < 3; k++ )
      {
        LVertex *v = cf->get_face_vertex ( k );
        it = map_v.find ( v->getIndex() );
        if ( it == map_v.end() )
        {
          map_v.insert ( std::pair<int, int> ( v->getIndex(), nb_v++ ) );
          vx.push_back ( v );
        }
      }
    }
  }

  fprintf ( fp, "%d\n", ( int ) vx.size() );
  for ( int j = 0 ; j < vx.size(); j++ )
    fprintf ( fp, "%.15lf %.15lf %.15lf\n", vx[j]->x(), vx[j]->y(), vx[j]->z() );

  fprintf ( fp, "%d\n", topoedge_size() );
  for ( int i = 0; i < topoedge_size() ; i++ )
  {
    TopoEdge *te = get_topoedge ( i );
//         if (te->color() != -1)
//         {
    fprintf ( fp, "%d %d\n", te->color(), te->vertex_size() );
//         int dir = te->get_direction(tf->color());
    for ( int j = 0; j < te->vertex_size(); j++ )
    {
      LVertex *v = te->get_vertex ( j, 0/*dir*/ );
      it = map_v.find ( v->getIndex() );
      fprintf ( fp, "%d ", it->second );
    }
//         }
    fprintf ( fp, "\n" );
  }

  fprintf ( fp, "%d\n", topoface_size() );
  for ( int i = 0; i < topoface_size(); i++ )
  {
    TopoFace *tf = get_topoface ( i );

    fprintf ( fp, "%d %d\n", tf->color(), tf->face_size() );
    for ( int j = 0; j < tf->face_size() ; j++ )
    {
      CutFace *cf = tf->get_face ( j );

      LVertex *v0 = cf->get_face_vertex ( 0 );
      it = map_v.find ( v0->getIndex() );
      fprintf ( fp, "%d ", it->second );
      LVertex *v1 = cf->get_face_vertex ( 1 );
      it = map_v.find ( v1->getIndex() );
      fprintf ( fp, "%d ", it->second );
      LVertex *v2 = cf->get_face_vertex ( 2 );
      it = map_v.find ( v2->getIndex() );
      fprintf ( fp, "%d\n", it->second );
    }

    fprintf ( fp, "%d\n", tf->topoedge_size() );
    for ( int j = 0; j < tf->topoedge_size() ; j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      fprintf ( fp, "%d ", te->color() );
    }
    fprintf ( fp, "\n" );
  }
  fclose ( fp );
}