// 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 "cadlevelset.h"

void cadLevelset::readLsnFile()
{
  int error = 0;
  char filename[127];
  sprintf ( filename, "%s.lsn", _opts->basename );
  tprintf ( "Reading Lsn file %s ....\n", filename );
  FILE *fp = fopen ( filename, "r" );

  if ( fp==NULL )
  {
    printf ( "Cannot open %s ...\n", filename );
    exit ( 0 );
  }

  int i,nb_ls=0;
  if ( fscanf ( fp, "%d", &nb_ls ) != 1 )
    error++;

  for ( i = 0; i < nb_ls; i++ )
  {
    DiscreteSurface *surf = new DiscreteSurface ( fp );
    _lsurf.push_back ( surf );
  }
  fclose ( fp );

  compute_bounds();

  int nb_nodes=0;
  tprintf ( "Discrete LS summary :\n" );
  tprintf ( "--> %d level-sets,\n", nb_ls );
  for ( i = 0; i < nb_ls; i++ )
  {
    tprintf ( "--> lsn[%d] contains %d nodes,\n", i, _lsurf[i]->size() );
    nb_nodes += _lsurf[i]->size();
  }
  tprintf ( "--> total : %d nodes.\n", nb_nodes );
}

int cadLevelset::intersected ( DiscreteSurface *surf, LTetra *t, int *out )
{
  int minus = 0, plus = 0;
  double projdist;
  double dist;
  for ( int i = 0; i < 4; i++ )
  {
    LVertex *v = t->getVertex ( i );
    if ( !surf->bbox_contain ( v ) )
      ( *out ) ++;
    if ( get_proj_dist ( v ) == INFINITY )
    {
      surf->distance ( v, dist, projdist );
      set_proj_dist ( v, projdist );
      set_dist ( v, dist );
    }

    if ( get_proj_dist ( v ) < 0.0 )
      minus++;
    else if ( get_proj_dist ( v ) > 0.0 )
      plus++;
    else
    {
      minus++;
      plus++;
    }
  }
//   printf("intersection --> %lf %lf %lf %lf\n", t->getVertex(0)->get_projdist(), t->getVertex(1)->get_projdist(), t->getVertex(2)->get_projdist(), t->getVertex(3)->get_projdist());
  if ( minus && plus )
    return 1;
  else return 0;
}

LTetra* cadLevelset::init_intersected ( DiscreteSurface *surf, LTetra *init )
{
  std::vector<LTetra *> flagged;
  double min_dist = INFINITY;
  LTetra *min_tetra = NULL;

  LVertex *spt;
  LVertex *sptn;
  double curv;

  surf->get_point ( 0, &spt, &sptn, &curv );

  if ( init->contain ( spt ) )
    return init;

  LTetra *t = init;

  srand ( 0 );

  while ( 1 )
  {
    int dir = rand() %4;
    LTetra *adj = t->getAdj ( dir );
    int opp = t->getOpp ( dir );
    if ( adj != NULL )
    {
      if ( adj->getFlag() == 0 )
      {
        if ( !init->contain ( spt ) )
        {
          t = adj;
        }
        else
        {
          for ( int j = 0; j < flagged.size(); j++ )
            flagged[j]->setFlag ( 0 );
          return adj;
        }
      }
      else
      {
        for ( int j = 0; j < flagged.size(); j++ )
          flagged[j]->setFlag ( 0 );
        return adj;
      }
    }
  }
}

LTetra* cadLevelset::search_intersected ( DiscreteSurface *surf, LTetra *init )
{
  std::vector<LTetra *> flagged;
  double min_dist = INFINITY;
  LTetra *min_tetra = NULL;

  int out = 0;
  if ( intersected ( surf, init, &out ) )
    return init;

  assert ( init );

  LTetra *t = init;

  while ( 1 )
  {
    min_dist = INFINITY;
    for ( int i = 0; i < 4; i++ )
    {
      assert ( t );
      LTetra *adj = t->getAdj ( i );
      int opp = t->getOpp ( i );
      if ( adj != NULL && adj->getFlag() == 0 )
      {
        out = 0;
        if ( !intersected ( surf, adj, &out ) )
        {
          double tpd = get_proj_dist ( adj->getVertex ( opp ) );
          if ( fabs ( tpd ) < min_dist )
          {
            min_dist = fabs ( tpd );
            min_tetra = adj;
          }
        }
        else
        {
          for ( int j = 0; j < flagged.size(); j++ )
            flagged[j]->resetFlag();
          return adj;
        }
      }
    }
    t->setFlag ( 1 );
//  printf("t --> %d\n", t->getIndex());
//  printf("min dist --> %lf\n", min_dist);
//  printf("min tetra --> %d\n", min_tetra->getIndex());
    flagged.push_back ( t );
    t = min_tetra;
  }
}

void cadLevelset::build_support ( DiscreteSurface *surf )
{
//     tprintf("--> building support ...\n");

  int local_idsurf =  _m->get_entity_manager()->add_face_id ( surf->get_id() );
  _m->get_entity_manager()->set_real_face ( local_idsurf );

  init_dist();

  for ( int i = 0; i < _m->tetra_size(); i++ )
    _m->get_tetra ( i )->setFlag ( 0 );

  std::queue<LTetra *> tq;
  std::vector<LTetra *> flagged;
  double min_dist = INFINITY;
  LTetra *min_tetra = NULL;

  LTetra *init = _m->get_vertex ( 0 )->getTetra();

//     printf("init tetra --> %d\n", init->getIndex());

  init = search_intersected ( surf, init );

  tq.push ( init );

  while ( !tq.empty() )
  {
    LTetra *t = tq.front();
    tq.pop();
    for ( int i = 0; i < 4; i++ )
    {
      LTetra *adj = t->getAdj ( i );
      if ( adj != NULL && adj->getFlag() == 0 )
      {
        int out = 0;
        if ( intersected ( surf, adj, &out ) )
        {
//        printf("out %d\n", out);
          if ( out < 4 )
          {
            double d = 0.0;
            for ( int j = 0; j < 4; j++ )
            {
              LVertex *v = t->getVertex ( j );
              if ( get_dist ( v ) > d )
                d = get_dist ( v );
            }
//      printf ( "longuest --> %lf / %lf\n", get_longuest(), d);

            if ( d <= 1.2*_m->longuest_edge() )
            {
              double pdt[4];
              for ( int j = 0; j < 4; j++ )
              {
                LVertex *v = adj->getVertex ( j );
                pdt[j] = get_proj_dist ( v );
              }
              _m->addLevelset ( local_idsurf, adj, pdt );
            }
          }

          tq.push ( adj );
          adj->setFlag ( 1 );
          flagged.push_back ( adj );
        }
      }
    }
  }
  // cleanup
  for ( int i = 0; i < flagged.size(); i++ )
    flagged[i]->setFlag ( 0 );
}


void cadLevelset::compute_edge_lenght ( GEdge *edge, double tol, std::list<LVertex>& vlist, std::list<double>& plist )
{
  std::list<double> llist;
  std::queue< std::list<LVertex>::iterator > itvlist;
  std::queue< std::list<double>::iterator > itplist;
  std::queue< std::list<double>::iterator > itllist;
  std::list<LVertex>::iterator itv;
  std::list<double>::iterator itp;
  std::list<double>::iterator itl;
  double t0 = /*dynamic_cast<OCCEdge *>(edge)*/edge->parBounds ( 1 ).low();
  double t1 = /*dynamic_cast<OCCEdge *>(edge)*/edge->parBounds ( 1 ).high();
  GPoint pt;
  LVertex d;
//     printf("t0 = %lf ---> t1 = %lf\n", t0, t1);
  pt = edge->point ( t0 );
  vlist.push_back ( LVertex ( pt.x(), pt.y(), pt.z() ) );
  plist.push_back ( t0 );
  pt = edge->point ( t1 );
  vlist.push_back ( LVertex ( pt.x(), pt.y(), pt.z() ) );
  plist.push_back ( t1 );
  d.setVector ( &vlist.front(), &vlist.back() );
  double tl = d.norm();
//     printf("Initial lenght --> %lf\n", tl);
  llist.push_back ( tl );
  itvlist.push ( vlist.begin() );
  itplist.push ( plist.begin() );
  itllist.push ( llist.begin() );
  while ( !itvlist.empty() )
  {
    itv = itvlist.front();
    LVertex v0 = *itv;
    std::advance ( itv, 1 );
    LVertex v1 = *itv;
    itp = itplist.front();
    t0 = *itp;
    std::advance ( itp, 1 );
    t1 = *itp;
//  printf("t0 = %lf ---> t1 = %lf\n", t0, t1);
    itl = itllist.front();
    double ol = *itl;
//  printf("old lenght --> %lf\n", ol);
    double tm = ( t0+t1 ) /2.0;
    pt = edge->point ( tm );
    LVertex vm ( pt.x(), pt.y(), pt.z() );
    d.setVector ( &v0, &vm );
    double nl1 = d.norm();
    d.setVector ( &vm, &v1 );
    double nl2 = d.norm();
//  printf("l1 = %lf -- l2 = %lf\n", nl1, nl2);
    assert ( nl1+nl2 - ol >= -tol );
    if ( nl1+nl2 - ol > tol*ol )
    {
      itv = vlist.insert ( itv, vm );
      itvlist.push ( itv );
      itp = plist.insert ( itp, tm );
      itplist.push ( itp );
      *itl = nl1;
      itl = llist.insert ( itl, nl2 );
      itllist.push ( itl );
    }
    else
    {
      itvlist.pop();
      itplist.pop();
      itllist.pop();
    }
//         tl = 0.0;
//         for (itl = llist.begin(); itl != llist.end(); itl++)
//             tl += *itl;
//  printf("Temp lenght --> %lf\n", tl);
  }
  tl = 0.0;
  for ( itl = llist.begin(); itl != llist.end(); itl++ )
    tl += *itl;
  edge->setLength ( tl );

//     for (itp = plist.begin(); itp != plist.end(); itp++)
//         printf("--> %lf\n", *itp);
//     printf("Final lenght --> %lf\n", edge->length());
  return;
}

DiscreteSurface *cadLevelset::add_levelset_if_tangetial_edge ( GEdge *edge, double curv_tol, double dotmin, int max_step )
{
  DiscreteSurface *surf = NULL;

  std::list<LVertex> vlist;
  std::list<double> plist;
  std::list<LVertex>::iterator itvlist;
  std::list<double>::iterator itplist;

  std::vector<GFace *> faces = edge->faces();
  std::vector<GFace *>::iterator itface = faces.begin();

  int surf0 = _gface_map.find ( *itface )->second;
  itface++;
  int surf1 = _gface_map.find ( *itface )->second;

  compute_edge_lenght ( edge, curv_tol, vlist, plist );
  assert ( vlist.size() >= 2 );

  int tangent = 0;
  int nb_step = 5;

  LVertex v0, v1;
  double p0, p1;
  double tstep;
  double dotmean = 0.0;

  LVertex uv;

  v0 = vlist.front();
  v1 = * ( ++vlist.begin() );
//     v1 = vlist.front();
  p0 = plist.front();
  p1 = * ( ++plist.begin() );
//     p1 = plist.front();
  itvlist = vlist.begin();
  itplist = plist.begin();
  tstep = ( plist.back() - plist.front() ) / ( double ) max_step;
//     printf("tstep = %lf\n", tstep);

  progress_bar pbar;

  for ( int i = 0; i < nb_step-1; i++ )
  {
    if ( tangent )
      pbar.progress ( i );

    double t = plist.front() + ( double ) i*tstep;
    while ( t >= p1 )
    {
      v0 = v1;
      v1 = * ( ++itvlist );
      p0 = p1;
      p1 = * ( ++itplist );
    }

    uv.setVector ( &v0, &v1 );
    double lt = p1-p0;
    double l = t-p0;
    double tt = l/lt;
    LVertex *pt = new LVertex ( v0.x() +tt* ( v1.x()-v0.x() ), v0.y() +tt* ( v1.y()-v0.y() ), v0.z() +tt* ( v1.z()-v0.z() ) );
    // on recupere les sommets les plus proches des surfaces formant l'arete

//         std::vector<Vertex *> results = _gb->get_closests(pt, surf_id);
    //         assert(results.size() == 2);
    LVertex *cl0 = new LVertex;
    LVertex *cl1 = new LVertex;
    LVertex *n0 = new LVertex;
    LVertex *n1 = new LVertex;
    _lsurf[surf0]->get_closest ( pt, &cl0, &n0 );

    _lsurf[surf1]->get_closest ( pt, &cl1, &n1 );

    // on analyse
    double dot = fabs ( n0->dotproduct ( n1 ) );
    dotmean += dot;
//  printf("dot --> %lf\n", dot);
    if ( tangent || dot > dotmin )
    {
      if ( !tangent )
      {
        v0 = vlist.front();
        v1 = vlist.front();
        p0 = plist.front();
        p1 = plist.front();
        i = -1;
        tangent = 1;
        surf = new DiscreteSurface ( NULL, _gface_map.size(), plist.front(), plist.back(), 0.0, 0.0, max_step );
        nb_step = max_step;
        pbar.start ( nb_step );;
      }
      else
      {
        LVertex *n = new LVertex;
        double curv = 0.0;
        n1->crossproduct ( &uv, n );
        surf->add_point ( pt, n, curv );
      }
    }
    else
      delete pt;
  }
  if ( surf )
    pbar.end();

  return surf;
}

void cadLevelset::compute_bounds()
{
  _min = new LVertex;
  _max = new LVertex;

  double global_bbmin[3] = {INFINITY, INFINITY, INFINITY};
  double global_bbmax[3] = {-INFINITY, -INFINITY, -INFINITY};

  for ( int i = 0; i < _lsurf.size(); i++ )
  {
    double surf_bbmin[3], surf_bbmax[3];
    _lsurf[i]->get_bbox ( surf_bbmin, surf_bbmax );

//     printf("bound --> %lf %lf %lf %lf %lf %lf\n", surf_bbmin[0], surf_bbmin[1], surf_bbmin[2], surf_bbmax[0], surf_bbmax[1], surf_bbmax[2]);

    if ( global_bbmin[0] > surf_bbmin[0] )
      global_bbmin[0] = surf_bbmin[0];
    if ( global_bbmin[1] > surf_bbmin[1] )
      global_bbmin[1] = surf_bbmin[1];
    if ( global_bbmin[2] > surf_bbmin[2] )
      global_bbmin[2] = surf_bbmin[2];
    if ( global_bbmax[0] < surf_bbmax[0] )
      global_bbmax[0] = surf_bbmax[0];
    if ( global_bbmax[1] < surf_bbmax[1] )
      global_bbmax[1] = surf_bbmax[1];
    if ( global_bbmax[2] < surf_bbmax[2] )
      global_bbmax[2] = surf_bbmax[2];
  }

  _min->setPosition ( global_bbmin[0], global_bbmin[1], global_bbmin[2] );
  _max->setPosition ( global_bbmax[0], global_bbmax[1], global_bbmax[2] );
}

void cadLevelset::build_surfaces()
{
  int surf_id = 0;
  progress_bar pbar ( 100 );

  tprintf ( "Surface building START\n" );

  int nb_face = _pModel->getNumFaces();
  assert(nb_face>0);
  pbar.start ( nb_face );
  for ( GModel::fiter it = _pModel->firstFace(); it != _pModel->lastFace(); ++it )
  {
    GFace *gf = *it;
//     printf("face %d --> %p\n", surf_id, gf);
    _gface_map.insert ( std::pair<GFace*, int> ( gf, surf_id ) );
    pbar.progress ( surf_id );
    double umin = gf->parBounds ( 0 ).low();
    double umax = gf->parBounds ( 0 ).high();
    double vmin = gf->parBounds ( 1 ).low();
    double vmax = gf->parBounds ( 1 ).high();
    DiscreteSurface *surf = new DiscreteSurface ( gf, surf_id, umin, umax, vmin, vmax, _opts->max_surf_sample );
    surf->generate();
    _lsurf.push_back ( surf );
    surf_id++;
  }
  pbar.end();

  compute_bounds();
//     dump_surfaces();

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

void cadLevelset::build_mesh_box()
{
  tprintf ( "Building mesh START\n" );

  double dl[3];
  dl[0] = _max->x() - _min->x();
  dl[0] *= _opts->enlarge_mesh_factor;
  dl[1] = _max->y() - _min->y();
  dl[1] *= _opts->enlarge_mesh_factor;
  dl[2] = _max->z() - _min->z();
  dl[2] *= _opts->enlarge_mesh_factor;
  _min->setPosition ( _min->x() - dl[0], _min->y() - dl[1], _min->z() - dl[2] );
  _max->setPosition ( _max->x() + dl[0], _max->y() + dl[1], _max->z() + dl[2] );

  tprintf ( "Bounds : %lf %lf %lf -- %lf %lf %lf\n", _min->x(), _min->y(), _min->z(), _max->x(), _max->y(), _max->z() );

  double step = _max->x() - _min->x();
  if ( _max->y() - _min->y() > step )
    step = _max->y() - _min->y();
  if ( _max->z() - _min->z() > step )
    step = _max->z() - _min->z();

  step /= ( double ) _opts->max_mesh_size;

  int nx = ( int ) ( ( _max->x() - _min->x() ) / step );
  int ny = ( int ) ( ( _max->y() - _min->y() ) / step );
  int nz = ( int ) ( ( _max->z() - _min->z() ) / step );

  double dx = ( _max->x() - _min->x() ) / ( double ) nx;
  double dy = ( _max->y() - _min->y() ) / ( double ) ny;
  double dz = ( _max->z() - _min->z() ) / ( double ) nz;

  tprintf ( "Mesh : %d (%lf) %d (%lf) %d (%lf)\n", nx, dx, ny, dy, nz, dz );

  _m->allocate_vertices ( ( nx+1 ) * ( ny+1 ) * ( nz+1 ) );

  for ( int k = 0; k < nz+1; k++ )
    for ( int j = 0; j < ny+1; j++ )
      for ( int i = 0; i < nx+1; i++ )
      {
        int index = i+j* ( nx+1 ) +k* ( nx+1 ) * ( ny+1 );
        LVertex *lv = _m->new_vertex ( _min->x() + ( double ) i*dx, _min->y() + ( double ) j*dy, _min->z() + ( double ) k*dz );
        _m->set_vertex ( index, lv );
//    printf("index = %d --> %p\n", index, v);
      }

  static const int tt[6][4] = { {0, 1, 2, 6}, {0, 5, 6, 4}, {0, 1, 6, 5}, {1, 7, 6, 5}, {1, 3, 2, 7}, {1, 2, 6, 7} };
  static const int pp[8][3] = { {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1} };

  progress_bar pbar;
  pbar.start ( nx*ny*nz );
  _m->allocate_tetra ( 6*nx*ny*nz );

  hashFaceStruct *hfs = new hashFaceStruct;

  for ( int k = 0; k < nz; k++ )
    for ( int j = 0; j < ny; j++ )
      for ( int i = 0; i < nx; i++ )
      {
        pbar.progress ( i+j*nx+k*nx*ny );
        for ( int tet = 0; tet < 6; tet++ )
        {
          int num;
          num = i+pp[tt[tet][0]][0] + ( j+pp[tt[tet][0]][1] ) * ( nx+1 ) + ( k+pp[tt[tet][0]][2] ) * ( nx+1 ) * ( ny+1 );
//        printf("%d %d %d %d\n", i, j, k, tet);
          LVertex *e0 = _m->get_vertex ( num );
          num = i+pp[tt[tet][1]][0] + ( j+pp[tt[tet][1]][1] ) * ( nx+1 ) + ( k+pp[tt[tet][1]][2] ) * ( nx+1 ) * ( ny+1 );
//        printf("%d %d %d %d\n", i, j, k, tet);
          LVertex *e1 = _m->get_vertex ( num );
          num = i+pp[tt[tet][2]][0] + ( j+pp[tt[tet][2]][1] ) * ( nx+1 ) + ( k+pp[tt[tet][2]][2] ) * ( nx+1 ) * ( ny+1 );
//        printf("%d %d %d %d\n", i, j, k, tet);
          LVertex *e2 = _m->get_vertex ( num );
          num = i+pp[tt[tet][3]][0] + ( j+pp[tt[tet][3]][1] ) * ( nx+1 ) + ( k+pp[tt[tet][3]][2] ) * ( nx+1 ) * ( ny+1 );
//        printf("%d %d %d %d\n", i, j, k, tet);
          LVertex *e3 = _m->get_vertex ( num );
          int index = tet+6*i+j*6*nx+k*6*nx*ny;
          LTetra *t = _m->new_tetra ( e0, e1, e2, e3, NULL );
          _m->set_tetra ( index, t );
          for ( int ii = 0; ii < 4; ii++ )
            hfs->addAdjFace ( t, ii );
        }
      }

  _m->compute_longuest_edge();

  pbar.end();
  tprintf ( "Mesh : %d cells, %d vertices\n", _m->tetra_size(), _m->vertex_size() );
  tprintf ( "Building mesh END\n" );
}


void cadLevelset::refine_mesh()
{
//      int max_refinement_level = 5;
//      double error_threshold = 1e-2;
  for ( int i = 0; i < _m->tetra_size(); i++ )
  {
    LTetra *t = _m->get_tetra ( i );
    t->addRefineField();
  }

  // refinement part
  _m->bucketize();
  for ( int i = 0; i < _lsurf.size(); i++ )
  {
    DiscreteSurface *surf = _lsurf[i];
    _m->add_included ( surf );
  }

  double max_curv_allowed = _m->evaluate_refinement_params();


  int loop_ref = 0;
  while ( _m->adapt ( max_curv_allowed ) )
  {
    tprintf ( "refinement loop --> %d\n", loop_ref );
    _m->refine_triangulation();
    loop_ref++;
  }

  _m->renumber();

//     char refinedname[100];
//     sprintf(refinedname, "%s_refined", _basename);
//     _m->write(refinedname);
}

void cadLevelset::process()
{
  // variables d'ajustement
  double tangent_dotproduct_threshold = 0.99;
  double curvature_threshold = 1e-3;
  // fin des variables d'ajustement

  std::map<GVertex*, int>::iterator itgv;
  std::map<GFace*, int>::iterator itgf;

  tprintf ( "Surface support building START\n" );
  tprintf ( "Model contains %d surfaces\n", ( int ) _lsurf.size() );

  int surf_id = 0;

  progress_bar pbar;
  pbar.start ( _lsurf.size() );
  for ( int i = 0; i < _lsurf.size(); i++ )
  {
    pbar.progress ( i );
    DiscreteSurface *surf = _lsurf[i];
    surf->bucketize ( _m->longuest_edge() );
//    surf->dump(i);
    build_support ( surf );
    surf_id = surf->get_id();
  }
  pbar.end();
  tprintf ( "Surface support building END\n" );

  if ( _opts->brep_topo )
  {
    tprintf ( "Tangent support building START\n" );
    int nb_edge = _pModel->getNumEdges();
    tprintf ( "Model contains %d edges\n", nb_edge );
    int tanget_edge_count = 0;
    int edge_id = 0;
    pbar.start ( nb_edge );
    for ( GModel::eiter it = _pModel->firstEdge(); it != _pModel->lastEdge(); ++it )
    {
      pbar.progress ( edge_id );
      // l'arete courante
      GEdge *edge = *it;
      std::list<LVertex> vlist;
      std::list<double> plist;
      std::list<LVertex>::iterator itvlist;
      std::list<double>::iterator itplist;
      DiscreteSurface *surf = add_levelset_if_tangetial_edge ( edge, curvature_threshold, tangent_dotproduct_threshold, _opts->max_surf_sample );
      if ( surf )
      {
        surf->bucketize ( _m->longuest_edge() );
//      surf->dump(-surf_id);
//      std::vector<LTetra*> tlist;
        build_support ( surf/*, tlist*/ );
        int s1 = _gface_map.find ( edge->faces().front() )->second;
        int s2 = _gface_map.find ( edge->faces().back() )->second;
        surf_id++;
        tanget_edge_count++;
      }
      edge_id++;
    }
    pbar.end();
    tprintf ( "Added %d tangent edges\n", tanget_edge_count );
    tprintf ( "Tangent support building END\n" );
  }

  _m->collect_edges();

#if 0
  std::vector<Vertex> vertices;
  std::vector<std::vector<int> > edges;
  std::vector<std::vector<int> > edge_faces;

  for ( GModel::viter it = _pModel->firstVertex(); it != _pModel->lastVertex(); ++it )
  {
    GVertex *gv = *it;
    _gvertex_map.insert ( std::pair<GVertex*, int> ( gv, vertex_id ) );
    vertices.push_back ( Vertex ( gv->x(), gv->y(), gv->z() ) );
    vertex_id++;
  }

  edges.resize ( _pModel->getNumEdges() );
  edge_faces.resize ( _pModel->getNumEdges() );

  for ( GModel::eiter it = _pModel->firstEdge(); it != _pModel->lastEdge(); ++it )
  {
    // l'arete courante
    GEdge *edge = ( *it );
    // le premier vertex
    GVertex *gv0 = edge->getBeginVertex();
    itgv = _gvertex_map.find ( gv0 );
    edges[edge_id].push_back ( itgv->second );
    // les suivants
    std::list<Vertex> vlist;
    std::list<double> plist;
//  compute_edge_lenght(edge, 1e-3, vlist, plist);
    cut_edge_with_mesh ( edge, _m, 1e-3, vlist, plist );
    vlist.erase ( vlist.begin() );
    vlist.erase ( --vlist.end() );
    for ( std::list<Vertex>::iterator itvx = vlist.begin(); itvx != vlist.end(); itvx++ )
    {
      vertices.push_back ( *itvx );
      edges[edge_id].push_back ( vertex_id );
      vertex_id++;
    }
    // le dernier vertex
    GVertex *gv1 = edge->getEndVertex();
    itgv = _gvertex_map.find ( gv1 );
    edges[edge_id].push_back ( itgv->second );

    std::list<GFace*> lf = edge->faces();
    for ( std::list<GFace*>::iterator itf = lf.begin(); itf != lf.end(); itf++ )
    {
      itgf = _gface_map.find ( *itf );
      edge_faces[edge_id].push_back ( itgf->second );
    }
    edge_id++;
  }

  _m->append_topo_vertex_size ( vertices.size() );
  for ( int i = 0; i < vertices.size(); i++ )
  {
    _m->append_topo_vertex ( i, vertices[i] );
  }

  _m->append_topo_edge_size ( edges.size() );
  for ( int i = 0; i < edges.size(); i++ )
  {
    _m->append_topo_edge ( i, edges[i], edge_faces[i] );
  }
#endif

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