// 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 "surface.h"
#include "LVertex.h"



void LSurface::process_tangent_edges()
{
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    LVertex *v0 = ce->get_vertex ( 0 );
    LVertex *v1 = ce->get_vertex ( 1 );
    for ( int j = 0; j < v0->nbSurfZero(); j++ ) // should be efficient enough because of small size of array
      for ( int k = 0; k < v1->nbSurfZero(); k++ )
      {
        int surf = v0->surfZeroId ( j );
        if ( get_entity_manager()->is_virtual_face ( surf ) && ( surf == v1->surfZeroId ( k ) ) )
        {
          ce->set_tangent();
          break;
        }
      }
  }
}

void LSurface::extract_cutfaces ( Mesh *m )
{
  progress_bar pbar ( m->tetra_size() );
  m->set_all_vertex_flag ( -1 );
  pbar.start();
  for ( int i = 0; i < m->tetra_size(); i++ )
  {
    pbar.progress ( i );
    LTetra *t = m->get_tetra ( i );
    for ( int j = 0; j < 4; j++ )
    {
      std::vector<LVertex*> vx;
      vx.push_back ( t->getFaceVertex ( j, 0 ) );
      vx.push_back ( t->getFaceVertex ( j, 1 ) );
      vx.push_back ( t->getFaceVertex ( j, 2 ) );
      m->compute_common_zero_surf ( vx );

      if ( m->common_surf_size() >= 1 )
      {
        std::vector<int> lsurf;
        int insert = 0; // par defaut --> ls virtuelle
//                 int insert = 1; // par defaut --> all ls

        for ( int k = 0; k < m->common_surf_size(); k++ )
        {
          int surf = m->get_common_surf ( k );
          lsurf.push_back ( surf );
          if ( get_entity_manager()->is_real_face ( surf ) || get_entity_manager()->is_restore_face ( surf ) )
            insert = 1;
        }
        if ( insert )
        {
          LVertex *opp = t->getVertex ( j );
//        for (int k = 0; k < m->common_surf_size(); k++)
          if ( opp->checkSurf ( lsurf[0] ) )
          {
            int sign = 0;
            if ( opp->surfVal ( lsurf[0] ) > 0.0 )
              sign = 1;
            addFace ( t, j, lsurf, sign );
//      dumpVTK("testing");
//      break;
          }
        }
      }
    }
  }
  pbar.end();
}

void LSurface::build_surface ( Mesh *m )
{
  tprintf ( "Surface construction START\n" );

  m->set_all_vertex_flag ( -1 );
  m->set_all_tetra_flag ( -1 );

  extract_cutfaces ( m );

  if ( _opts->brep_topo )
    process_tangent_edges();

  debug_edge();
  dumpVTK ( "initial" );

  // suppression des bords libres
  if ( _opts->only_valid_topo )
    erode_surface();

#if 0
  debug_edge();
  dumpVTK ( "surface_inter" );
#endif

  // calcul des faces distinctes
  compute_face_topo();
  // calcul des aretes
  compute_edge_topo();

  debug_border_edge();

  if ( _nb_topo_face == 0 )
  {
    tprintf ( "Invalid topology --> open surface !\n" );
    exit ( 0 );
  }
  else
  {
    dumpVTK ( "final" );
    tprintf ( "Valid topology --> closed surface\n" );
  }

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

int LSurface::addFace ( LTetra *t, int face, std::vector<int> surf, int sign )
{
  assert ( t != NULL );
  assert ( face>=0 && face<4 );

  hashFaceKey *k1 = new hashFaceKey ( t, face );
  std::multimap<int, hashFaceKey*>::iterator it;
  std::pair<std::multimap<int, hashFaceKey*>::iterator,std::multimap<int, hashFaceKey*>::iterator> ret;
  ret = _face_hashtable.equal_range ( k1->key() );
  for ( it=ret.first; it!=ret.second; ++it )
  {
    hashFaceKey *k2 = it->second;
    if ( k2->equal ( k1 ) )
    {
      delete k1;
//      _face_hashtable.erase(it);
//      delete k2;
      return 0;
    }
  }
  // not found
  _face_hashtable.insert ( std::pair<int, hashFaceKey *> ( k1->key(), k1 ) );
  k1->setIndex ( _face_list.size() );
  CutFace *cf = new CutFace ( t, face, surf, sign );
  _face_list.push_back ( cf ); // adding face
  // edge hashing
  addEdges ( cf );
  return 1;
}

int LSurface::addEdges ( CutFace *f )
{
  assert ( f != NULL );
  int add = 0;
  LTetra *t = f->int_tetra();
  int face = f->int_face();
  assert ( t != NULL );
  for ( int i = 0; i < 3; i++ )
  {
    hashEdgeKey *k1 = new hashEdgeKey ( t, t->getFaceEdgeIndex ( face, i ) );
    std::multimap<int, hashEdgeKey*>::iterator it;
    std::pair<std::multimap<int, hashEdgeKey*>::iterator,std::multimap<int, hashEdgeKey*>::iterator> ret;
    ret = _edge_hashtable.equal_range ( k1->key() );
    int found = 0;
    for ( it=ret.first; it!=ret.second; ++it )
    {
      hashEdgeKey *k2 = it->second;
      if ( k2->equal ( k1 ) )
      {
        int index = k2->index();
        _edge_list[index]->addFace ( f, i );
        f->add_edge ( _edge_list[index] );
        delete k1;
        found = 1;
        break;
      }
    }
    if ( !found )
    {
      _edge_hashtable.insert ( std::pair<int, hashEdgeKey *> ( k1->key(), k1 ) );
      k1->setIndex ( _edge_list.size() );
      CutEdge *ce = new CutEdge ( f, i );
      f->add_edge ( ce );
      _edge_list.push_back ( ce );
      add = 1;
    }
  }
  return add;
}

void LSurface::addBorderEdge ( CutEdge *ce )
{
  assert ( ce != NULL );
  assert ( ce->get_vertex ( 0 ) != ce->get_vertex ( 1 ) );
  ////////////////////////////
  ce->set_flag ( 1 );
  LVertex *v0 = ce->get_vertex ( 0 );
  if ( v0->getFlag() == -1 )
  {
    _vertex_border_list.push_back ( v0 );
    v0->setFlag ( 0 );
  }
  _border_edge_hashtable.insert ( std::pair<LVertex *, CutEdge *> ( v0, ce ) );
  v0->incFlag();
  ////////////////////////////
  LVertex *v1 = ce->get_vertex ( 1 );
  if ( v1->getFlag() == -1 )
  {
    _vertex_border_list.push_back ( v1 );
    v1->setFlag ( 0 );
  }
  _border_edge_hashtable.insert ( std::pair<LVertex *, CutEdge *> ( v1, ce ) );
  v1->incFlag();
}

void LSurface::linkBorderEdges()
{
  std::multimap<LVertex *, CutEdge*>::iterator it;
  std::pair<std::multimap<LVertex *, CutEdge*>::iterator, std::multimap<LVertex *, CutEdge*>::iterator> ret;

  for ( int i = 0; i < _vertex_border_list.size(); i++ )
    assert ( _border_edge_hashtable.count ( _vertex_border_list[i] ) == _vertex_border_list[i]->getFlag() );

  debug_edge_vertices();

//     tprintf("--> remove orphan edges ...\n");

  // supprimer les aretes orphelines
  for ( int i = 0; i < _vertex_border_list.size(); i++ )
  {
//  printf("%d / %d\n", i, (int)_vertex_border_list.size());
    LVertex *node = _vertex_border_list[i];
//  printf("flag --> %d\n", node->getFlag());
    if ( node->getFlag() == 1 )
    {
      assert ( _border_edge_hashtable.count ( node ) == 1 );
      it = _border_edge_hashtable.find ( node );
      while ( it != _border_edge_hashtable.end() )
      {
        CutEdge *ce = it->second;
        _border_edge_hashtable.erase ( it );
        if ( ce->get_vertex ( 0 ) != node )
          ce->invert_dir();
        node->decFlag();
        assert ( node->getFlag() == 0 );
        ce->set_color ( -2 );
        assert ( ce->get_vertex ( 0 ) == node );
        node = ce->get_vertex ( 1 );
        if ( node->getFlag() >= 2 )
        {
          node->decFlag();
          ret = _border_edge_hashtable.equal_range ( node );
          for ( it = ret.first; it != ret.second; ++it )
            if ( it->second == ce )
            {
              _border_edge_hashtable.erase ( it );
              break;
            }
          if ( node->getFlag() >= 2 )
            break;
          it = _border_edge_hashtable.find ( node );
        }
        else
          break;
      }
    }
  }

//     tprintf("--> link open edges ...\n");

  // construire les aretes non fermees
  for ( int i = 0; i < _vertex_border_list.size(); i++ )
  {
    LVertex *node = _vertex_border_list[i];
    if ( node->getFlag() >= 3 )
    {
      it = _border_edge_hashtable.find ( node );
      while ( it != _border_edge_hashtable.end() )
      {
        CutEdge *ce = it->second;
        _border_edge_hashtable.erase ( it );
        if ( ce->get_color() == -1 )
        {
          if ( ce->get_vertex ( 0 ) != node )
            ce->invert_dir();
          CutEdge *n = NULL;
          LVertex *v = ce->get_vertex ( 1 );
          assert ( v != node );
          assert ( ce->get_color() != -2 );
          ce->set_color ( _nb_topo_edge );
          _border_edge_list.push_back ( ce );
          while ( v->getFlag() < 3 )
          {
            assert ( v->getFlag() == 2 );
            for ( int j = 0; j < 2; j++ )
            {
              it = _border_edge_hashtable.find ( v );
              assert ( it != _border_edge_hashtable.end() );
              if ( it->second != ce )
                n = it->second;
              _border_edge_hashtable.erase ( it );
            }
            if ( n->get_vertex ( 0 ) != v )
              n->invert_dir();
            assert ( n->get_vertex ( 0 ) == v );
            set_link ( ce, n );
            ce = n;
            assert ( ce->get_color() != -2 );
            ce->set_color ( _nb_topo_edge );
            _border_edge_list.push_back ( ce );
            v = ce->get_vertex ( 1 );
          }
          _nb_topo_edge++;
        }
        it = _border_edge_hashtable.find ( node );
      }
    }
  }

  if ( _border_edge_hashtable.size() )
  {
    tprintf ( "--> Warning : Edge loops found ...\n" );
  }

//     tprintf("--> link loops ...\n");

  // construire les aretes fermees
  for ( int i = 0; i < _vertex_border_list.size(); i++ )
  {
    LVertex *node = _vertex_border_list[i];
    if ( node->getFlag() == 2 )
    {
      it = _border_edge_hashtable.find ( node );
      while ( it != _border_edge_hashtable.end() )
      {
        CutEdge *ce = it->second;
        assert ( ce->get_color() != -2 );
        _border_edge_hashtable.erase ( it );
        if ( ce->get_color() == -1 )
        {
          if ( ce->get_vertex ( 0 ) != node )
            ce->invert_dir();
          CutEdge *n = NULL;
          LVertex *v = ce->get_vertex ( 1 );
          assert ( v != node );
          ce->set_color ( _nb_topo_edge );
          _border_edge_list.push_back ( ce );
          while ( v != node )
          {
            assert ( v->getFlag() == 2 );
            for ( int j = 0; j < 2; j++ )
            {
              it = _border_edge_hashtable.find ( v );
              assert ( it != _border_edge_hashtable.end() );
              if ( it->second != ce )
                n = it->second;
              _border_edge_hashtable.erase ( it );
            }
            if ( n->get_vertex ( 0 ) != v )
              n->invert_dir();
            assert ( n->get_vertex ( 0 ) == v );
            set_link ( ce, n );
            ce = n;
            assert ( ce->get_color() != -2 );
            ce->set_color ( _nb_topo_edge );
            _border_edge_list.push_back ( ce );
            v = ce->get_vertex ( 1 );
          }
          _nb_topo_edge++;
        }
        it = _border_edge_hashtable.find ( node );
      }
    }
  }
}

int LSurface::erode_surface ()
{
  tprintf ( "Erosion --> %d edges / %d faces to be processed\n", ( int ) _edge_list.size(), ( int ) _face_list.size() );

  int loop = 0;
  int max_loop = 10;
  int unfilled_color = -1;
  int erased_color = -2;

#if 1
  // building list of free faces
  std::vector<CutFace *> free_face_list;
  do
  {
    free_face_list.clear();

    if ( _opts->brep_topo || _opts->auto_topo )
    {
      for ( int i = 0; i < _face_list.size(); i++ )
      {
        CutFace *cf = _face_list[i];
        if ( cf->get_color() == unfilled_color )
        {
          for ( int j = 0; j < cf->edge_size(); j++ )
          {
            CutEdge *ce = cf->get_edge ( j );
            if ( ce->face_size() == 1 ) // free edge
            {
              cf->inc_free();
              free_face_list.push_back ( cf );
              break;
            }
          }
        }
      }
    }
    if ( _opts->recut_enable )
    {
      for ( int i = 0; i < _face_list.size(); i++ )
      {
        CutFace *cf = _face_list[i];
        if ( cf->get_color() == unfilled_color && get_entity_manager()->is_restore_face ( cf->surface ( 0 ) ) )
        {
          for ( int j = 0; j < cf->edge_size(); j++ )
          {
            CutEdge *ce = cf->get_edge ( j );
            if ( ce->face_size() == 1 ) // free edge
            {
              cf->inc_free();
              free_face_list.push_back ( cf );
              break;
            }
          }
        }
      }
    }

//         printf("face list size --> %d\n", (int)free_face_list.size());

    if ( free_face_list.size() )
    {
      int removed = 0;
      for ( int i = 0; i < free_face_list.size(); i++ )
      {
        CutFace *cf = free_face_list[i];
        if ( cf->get_color() == unfilled_color )
        {
          seed ( cf );
          removed += color ( erased_color, unfilled_color );
        }
      }

      // update
      for ( int i = 0; i < _edge_list.size(); i++ )
      {
        CutEdge *ce = _edge_list[i];
        for ( int j = 0; j < ce->face_size(); j++ )
        {
          CutFace *cf = ce->get_face ( j );
          if ( cf->get_color() == erased_color )
          {
            cf->set_free ( -1 );
            ce->erase ( j );
            j--;
          }
        }
      }
    }
  }
  while ( free_face_list.size() && ( loop++ < max_loop ) );

#else
  std::vector<CutFace *> free_face_list;
  max_loop = 200;
  do
  {
    free_face_list.clear();
    for ( int i = 0; i < _face_list.size(); i++ )
    {
      CutFace *cf = _face_list[i];
      if ( cf->get_color() == unfilled_color )
      {
        for ( int j = 0; j < cf->edge_size(); j++ )
        {
          CutEdge *ce = cf->get_edge ( j );
          if ( ce->face_size() == 1 ) // free edge
            cf->inc_free();
        }
      }
    }
    for ( int i = 0; i < _face_list.size(); i++ )
    {
      CutFace *cf = _face_list[i];
      if ( cf->get_color() == unfilled_color )
        if ( cf->free() != 0 )
          free_face_list.push_back ( cf );
    }

    if ( free_face_list.size() )
    {
      int removed = 0;
      for ( int i = 0; i < free_face_list.size(); i++ )
      {
        CutFace *cf = free_face_list[i];
        cf->set_color ( erased_color );
      }

      // update
      for ( int i = 0; i < _edge_list.size(); i++ )
      {
        CutEdge *ce = _edge_list[i];
        for ( int j = 0; j < ce->face_size(); j++ )
        {
          CutFace *cf = ce->get_face ( j );
          if ( cf->get_color() == erased_color )
          {
            cf->set_free ( -1 );
            ce->erase ( j );
            j--;
          }
        }
      }
    }
    char name[100];
    sprintf ( name ,"surface_inter_%d", loop );
    dumpVTK ( name );
  }
  while ( free_face_list.size() && ( loop++ < max_loop ) );
#endif

  // retour
  if ( loop == max_loop )
  {
    tprintf ( "Ending erosion at max iteration %d --> Something went wrong !!\n", loop );
    return 0;
  }

  return 1;
}

int LSurface::color ( int surf, int unfilled_surf )
{
  int count = 0;
  while ( _sc.size() )
  {
    CutFace *cf = _sc.front();
    if ( cf->get_color() == unfilled_surf )
    {
      cf->set_color ( surf );
      count++;
      for ( int i = 0; i < cf->edge_size(); i++ )
      {
        CutEdge *ce = cf->get_edge ( i );
        if ( ce->face_size() == 2 && !ce->is_tangent() )
        {
          if ( ( ce->get_face ( 0 ) != cf ) && ( ce->get_face ( 0 )->same_surface ( cf ) ) )
            _sc.push ( ce->get_face ( 0 ) );
          else if ( ce->get_face ( 1 )->same_surface ( cf ) )
            _sc.push ( ce->get_face ( 1 ) );
        }
      }
    }
    _sc.pop();
  }
  return count;
}

void LSurface::compute_face_topo()
{
  int unfilled_color = -1;
  std::vector<CutFace *> multi_surf;
  for ( int i = 0; i < _face_list.size(); i++ )
  {
    CutFace *cf = _face_list[i];
    if ( cf->surf_size() == 1 && cf->get_color() == unfilled_color ) // ne commence qu'à partir d'une face monosurface
    {
      seed ( cf );
      int nb = color ( _nb_topo_face, unfilled_color );
      if ( nb )
        _nb_topo_face++;
    }
    if ( cf->surf_size() > 1 ) // traitement ulterieur
      multi_surf.push_back ( cf );
  }
  // face multisurface
  for ( int i = 0; i < multi_surf.size(); i++ )
  {
    CutFace *cf = multi_surf[i];
    if ( cf->get_color() == unfilled_color )
    {
      seed ( cf );
      int nb = color ( _nb_topo_face, unfilled_color );
      if ( nb )
        _nb_topo_face++;
    }
  }
  tprintf ( "--> Found %d topologic faces\n", _nb_topo_face );

}

void LSurface::debug_edge()
{
  if ( !_opts->debug )
    return;

  FILE *fp = fopen ( "debug_edges.vtk", "w" );
  fprintf ( fp, "# vtk DataFile Version 2.0\n" );
  fprintf ( fp, "Really cool data\n" );
  fprintf ( fp, "ASCII\n" );
  fprintf ( fp, "DATASET POLYDATA\n" );
  int edge_size = 0;
  long seek;
  fprintf ( fp, "POINTS " );
  seek = ftell ( fp );
  fprintf ( fp, "%10d double\n", 0 );
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    if ( ce->face_size() )
    {
      LVertex *v = ce->get_vertex ( 0 );
      fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
      v = ce->get_vertex ( 1 );
      fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
      edge_size++;
    }
  }
  fseek ( fp, seek, SEEK_SET );
  fprintf ( fp, "%10d ", 2*edge_size );
  fseek ( fp, 0, SEEK_END );
  fprintf ( fp, "LINES %d %d\n", edge_size, 3*edge_size );
  int vertex_index = 0;
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    if ( ce->face_size() )
    {
      fprintf ( fp, "%d %d %d\n", 2, 2*vertex_index, 2*vertex_index+1 );
      vertex_index++;
    }
  }
  fprintf ( fp, "CELL_DATA %d\n",edge_size );
//     fprintf (fp, "CELL_DATA %d\n",(int) _edge_list.size());
  fprintf ( fp, "SCALARS edge_color int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    if ( ce->face_size() )
      fprintf ( fp, "%d\n", ce->face_size() );
  }
  fprintf ( fp, "SCALARS edge_tangent int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    if ( ce->face_size() )
      fprintf ( fp, "%d\n", ce->is_tangent() );
  }
  fclose ( fp );
}

void LSurface::debug_edge_vertices()
{
  if ( !_opts->debug )
    return;

  FILE *fp = fopen ( "debug_edge_vertices.vtk", "w" );
  fprintf ( fp, "# vtk DataFile Version 2.0\n" );
  fprintf ( fp, "Really cool data\n" );
  fprintf ( fp, "ASCII\n" );
  fprintf ( fp, "DATASET POLYDATA\n" );
  fprintf ( fp, "POINTS %d double\n", ( int ) _vertex_border_list.size() );
  for ( int i = 0; i < _vertex_border_list.size(); i++ )
  {
    LVertex *v = _vertex_border_list[i];
    fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
  }

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

//     fprintf (fp, "LINES %d %d\n",  (int)_border_edge_list.size(), 3*(int)_vertex_border_list.size());
//     for (int i = 0; i < _vertex_border_list.size(); i++)
//     {
//             fprintf (fp, "%d %d\n", 1, i);
//     }

  fprintf ( fp, "POINT_DATA %d\n", ( int ) _vertex_border_list.size() );
  fprintf ( fp, "SCALARS point_color int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < _vertex_border_list.size(); i++ )
  {
    LVertex *v = _vertex_border_list[i];
    fprintf ( fp, "%d\n", v->getFlag() );
  }

  fclose ( fp );
}

void LSurface::debug_border_edge()
{
  if ( !_opts->debug )
    return;

  FILE *fp = fopen ( "debug_border_edges.vtk", "w" );
  fprintf ( fp, "# vtk DataFile Version 2.0\n" );
  fprintf ( fp, "Really cool data\n" );
  fprintf ( fp, "ASCII\n" );
  fprintf ( fp, "DATASET POLYDATA\n" );

  fprintf ( fp, "POINTS %d double\n", 2* ( int ) _border_edge_list.size() );
  for ( int i = 0; i < _border_edge_list.size(); i++ )
  {
    CutEdge *ce = _border_edge_list[i];
    LVertex *v = ce->get_vertex ( 0 );
    fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
    v = ce->get_vertex ( 1 );
    fprintf ( fp, "%lf %lf %lf\n", v->x(), v->y(), v->z() );
  }
  fprintf ( fp, "LINES %d %d\n", ( int ) _border_edge_list.size(), 3* ( int ) _border_edge_list.size() );
  int c = 0;
  for ( int i = 0; i < _border_edge_list.size(); i++ )
  {
    CutEdge *ce = _border_edge_list[i];
    fprintf ( fp, "2 %d %d\n", c++, c++ );
  }
  fclose ( fp );
}



void LSurface::compute_edge_topo()
{
  int erased_color = -2;
  // selectionner les aretes
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    ce->set_flag ( 0 );
    // supprimer les faces invalides des aretes
    for ( int j = 0; j < ce->face_size(); j++ )
      if ( ce->get_face ( j )->get_color() == erased_color )
        ce->erase ( j-- );
    // marquer les aretes a chainer
    for ( int j = 0; j < ce->face_size(); j++ )
      if ( ce->get_face ( j )->get_color() != ce->get_face ( ( j+1 ) %ce->face_size() )->get_color() )
      {
        ce->set_flag ( 1 );
        break;
      }
    // cas des aretes de faces ouvertes
    if ( !_opts->only_valid_topo )
      if ( ce->face_size() == 1 )
        ce->set_flag ( 1 );
  }

  // ranger les aretes de bords
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    ce->reset_default();
    if ( ce->get_flag() )
      addBorderEdge ( ce );
  }

  debug_border_edge();

  // chainer les aretes
  linkBorderEdges();

  // resume
  tprintf ( "--> Found %d topologic edges\n", _nb_topo_edge );
  tprintf ( "    +--> %d linear edges\n", ( int ) _border_edge_list.size() );

}

void LSurface::clear()
{
  for ( int i = 0; i < _face_list.size(); i++ )
  {
    CutFace *cf = _face_list[i];
    if ( cf->free() )
      delete cf;
  }
  for ( int i = 0; i < _edge_list.size(); i++ )
  {
    CutEdge *ce = _edge_list[i];
    if ( !ce->face_size() )
      delete ce;
  }
  for ( std::multimap<int, hashFaceKey*>::iterator it = _face_hashtable.begin(); it != _face_hashtable.end(); it++ )
  {
    hashFaceKey *hek = it->second;
    delete hek;
  }
}

void LSurface::dumpVTK ( const char *ext, bool number )
{
  if ( !_opts->debug )
    return;

  static int count = 0;
  char filename[100];
  if ( number )
    sprintf ( filename, "%s_%s_surf_%d.vtk", _opts->basename, ext, count++ );
  else
    sprintf ( filename, "%s_%s_surf.vtk", _opts->basename, ext );
  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_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 < _face_list.size(); i++ )
  {
    CutFace *f = _face_list[i];
    if ( f->free() != -1 )
    {
      nb_f_v += 3;
      nb_f++;
      for ( int j = 0; j < 3; j++ )
      {
        LVertex *v = f->int_tetra()->getFaceVertex ( f->int_face(), j );
        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 < _face_list.size(); i++ )
  {
    CutFace *f = _face_list[i];
    if ( f->free() != -1 )
    {
      fprintf ( fp, "%d ", 3 );
      for ( int j = 0; j < 3; j++ )
      {
        LVertex *v = f->int_tetra()->getFaceVertex ( f->int_face(), j );
        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 color int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < _face_list.size(); i++ )
  {
    CutFace *f = _face_list[i];
    if ( f->free() != -1 )
      fprintf ( fp, "%d\n", f->get_color() );
  }
  fprintf ( fp, "SCALARS surf_number int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < _face_list.size(); i++ )
  {
    CutFace *f = _face_list[i];
    if ( f->free() != -1 )
      fprintf ( fp, "%d\n", f->surf_size() );
  }
  fprintf ( fp, "SCALARS surf int 1\n" );
  fprintf ( fp, "LOOKUP_TABLE default\n" );
  for ( int i = 0; i < _face_list.size(); i++ )
  {
    CutFace *f = _face_list[i];
    if ( f->free() != -1 )
      fprintf ( fp, "%d\n", f->surface ( 0 ) );
  }

  fclose ( fp );

  if ( _border_edge_list.size() )
  {
    sprintf ( filename, "%s_%s_edges.vtk", _opts->basename, ext );
    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" );
    ////////////////////////////
    // dump points
    fprintf ( fp, "POINTS %d double\n", ( int ) vx.size() );
    for ( int i = 0 ; i < vx.size(); i++ )
      fprintf ( fp, "%lf %lf %lf\n", vx[i]->x(), vx[i]->y(), vx[i]->z() );

    fprintf ( fp, "LINES %d %d\n", ( int ) _border_edge_list.size(), ( int ) _border_edge_list.size() *3 );

    for ( int i = 0; i < _border_edge_list.size(); i++ )
    {
      CutEdge *ce = _border_edge_list[i];
      fprintf ( fp, "%d ", 2 );
      LVertex *v0 = ce->get_vertex ( 0 );
      std::map<int, int >::iterator iter;
      iter = map_v.find ( v0->getIndex() );
      fprintf ( fp, "%d ", iter->second );
      LVertex *v1 = ce->get_vertex ( 1 );
      iter = map_v.find ( v1->getIndex() );
      fprintf ( fp, "%d ", iter->second );
      fprintf ( fp, "\n" );
    }

    fprintf ( fp, "CELL_DATA %d\n", ( int ) _border_edge_list.size() );
    fprintf ( fp, "SCALARS edge_color int 1\n" );
    fprintf ( fp, "LOOKUP_TABLE default\n" );
    for ( int i = 0; i < _border_edge_list.size(); i++ )
    {
      CutEdge *ce = _border_edge_list[i];
      fprintf ( fp, "%d\n", ce->get_color() );
    }
    fprintf ( fp, "SCALARS edge_tangent int 1\n" );
    fprintf ( fp, "LOOKUP_TABLE default\n" );
    for ( int i = 0; i < _border_edge_list.size(); i++ )
    {
      CutEdge *ce = _border_edge_list[i];
      fprintf ( fp, "%d\n", ce->is_tangent() );
    }

    fclose ( fp );
  }
}