// 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 "mesh_cutter.h"
#ifdef USE_OLD_CUTTER
#include <MElementCut.h>
#endif
#include "MSubElement.h"
#include "discreteVertex.h"
#include <fstream>
#include <iostream>

void MeshCutter::build_GModel()
{
  tprintf ( "GModel building START\n" );

  std::map<LTetra*, MSubTetrahedron*> ltsbt;
  std::map<LTetra*, MSubTetrahedron*>::iterator itlt;
  std::map<MTetrahedron*, int> sbtv;
  std::map<MTetrahedron*, int>::iterator itmpv;

  std::map<TopoVertex*, GVertex*> tvgv;
  std::map<TopoVertex*, GVertex*>::iterator itv;
  std::map<TopoEdge*, GEdge*> tege;
  std::map<TopoEdge*, GEdge*>::iterator ite;
  std::map<TopoFace*, GFace*> tfgf;
  std::map<TopoFace*, GFace*>::iterator itf;
  std::map<int, GRegion*> vogr;
  std::map<int, GRegion*>::iterator itr;

  _mymesh->set_all_tetra_flag ( -1 );
  _mymesh->clear_all_tetra_child();

  _gm = new GModel ( "cut_mesh" );

  progress_bar pbar;

  EntityManager *em = _mymesh->get_entity_manager();

  tprintf ( "--> create regions ...\n" );

  int null_region = 0;
  GRegion *nullgr = new discreteRegion ( _gm, null_region );
//     nullgr->addPhysicalEntity(_mytopo->mpvols[null_region]);
  nullgr->addPhysicalEntity ( null_region );
  vogr.insert ( std::pair<int, GRegion*> ( null_region, nullgr ) );
  _gm->add ( nullgr );

  for ( int i = 0; i < _mytopo->topovolume_size(); i++ )
  {
    TopoVolume *tvol = _mytopo->get_topovolume ( i );
    GRegion *gr = new discreteRegion ( _gm, i+1 );
    if ( _opts->brep_topo )
    {
      assert ( tvol->getFlag() == 1 );
      assert ( tvol->potential_brep_index_size() == 1 );
      gr->addPhysicalEntity ( tvol->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      std::vector<int> list;
      if ( tvol->get_parent() != NULL )
      {
        TopoVolume *root = tvol->get_parent();
        int root_id_vol = root->get_physical ( 0 );
        assert ( em->is_restore_vol ( root_id_vol ) );
        list.push_back ( root_id_vol );
        dbgprintf ( 3, "--> root id vol = %d\n", root_id_vol );
        std::vector<int> assoc = em->get_associated_edges ( root_id_vol );
        list.insert ( list.begin(), assoc.begin(), assoc.end() );
      }

      assert ( tvol->physical_size() );
      int local_id_vol = tvol->get_physical ( 0 );
      list.push_back ( local_id_vol );
      dbgprintf ( 3, "--> local id vol = %d\n", local_id_vol );

      dbgprintf ( 3, "--> vol list size = %d\n", ( int ) list.size() );
      for ( int k = 0; k < list.size(); k++ )
      {
        local_id_vol = list[k];
        dbgprintf ( 3, "adding local_id_vol --> %d\n", local_id_vol );
        int idvol = em->get_vol_id ( local_id_vol );
        vogr.insert ( std::pair<int, GRegion*> ( local_id_vol, gr ) );
        dbgprintf ( 3, "vol --> %d\n", idvol );
        gr->addPhysicalEntity ( idvol );
      }
      _gm->add ( gr );
    }
  }

  tprintf ( "--> setting tetra children ...\n" );

  // cleaning the flags and association of children
  for ( int i = 0; i < _mymesh->tetra_size(); i++ )
  {
    LTetra *t = _mymesh->get_tetra ( i );
    t->setFlag ( -1 );
    for ( int j = 0; j < 4; j++ )
      t->getVertex ( j )->setFlag ( -1 );
    LTetra *at = t->getAncestor();
    if ( at )
    {
      at->setFlag ( -1 );
      for ( int j = 0; j < 4; j++ )
        at->getVertex ( j )->setFlag ( -1 );
      at->add_child ( t );
    }
  }

  // association of physicals on volumes : nullgr for parents, gr for others
  // creation of volume sub-elements for children and association of physicals
  for ( int i = 0; i < _mymesh->tetra_size(); i++ )
  {
    LTetra *t = _mymesh->get_tetra ( i );

    if ( t->getFlag() != -1 )
      continue;

    if ( t->getAncestor() )
      t = t->getAncestor();

    MTetrahedron *mt = NULL;
    if ( t->child_size() )
      mt = add_tetra ( t, nullgr );
    else
    {
      dbgprintf ( 4, "--> looking for volume %d\n", t->getVol() );
      itr = vogr.find ( t->getVol() );
      assert ( itr != vogr.end() );
      GRegion *gr = itr->second;
      mt = add_tetra ( t, gr );
    }

    for ( int j = 0; j < t->child_size(); j++ )
    {
      LTetra *tc = t->get_child ( j );
      assert ( tc != NULL );
      dbgprintf ( 4, "child tetra %p (%d) --> flag %d\n", tc, tc->getIndex(), tc->getFlag() );
      assert ( t != tc );
      itr = vogr.find ( tc->getVol() );
      assert ( itr != vogr.end() );
      GRegion *gr = itr->second;
      MSubTetrahedron* mtc = dynamic_cast<MSubTetrahedron*> ( add_tetra ( tc, gr ) );
      assert ( mtc );
      ltsbt.insert ( std::pair<LTetra*, MSubTetrahedron*> ( tc, mtc ) );
      mtc->setParent ( mt, !j );
    }
  }


  std::map< int, std::set<std::pair<int,int> > > map_edges; // map<physical, listofedges>
  std::map< LVertex*, std::map<int,int> > map_vertices; // map<vertex, map<idface,count> >
  std::map< int, std::vector<GVertex*> > map_endpoints; // map<physical, listofendpoints>

  if ( _opts->read_edges )
  {
    // identification of implicit edges
    for ( int i=0; i<_mymesh->tetra_size(); i++ )
    {
      LTetra *t = _mymesh->get_tetra ( i );
      for ( int j=0; j<6; j++ )
      {
        LVertex* v0 = t->getEdgeVertex(j,0);
        LVertex* v1 = t->getEdgeVertex(j,1);

        _mymesh->compute_common_zero_surf(v0,v1);
        if(_mymesh->common_surf_size())
        {
          if( v0->getIndex() > v1->getIndex() )
          {
            LVertex* tmp = v0;
            v0 = v1;
            v1 = tmp;
          }
        }
        for( int k=0; k<_mymesh->common_surf_size(); k++ )
        {
          int local_id_surf = _mymesh->get_common_surf(k);
          int idface = _mymesh->get_entity_manager()->get_face_id(local_id_surf);
          if (idface-idface%(1000000)==1000000)
          {
            std::pair<  std::map< int, std::set<std::pair<int,int> > >::iterator,bool  > ret_edge;
            std::pair< std::set<std::pair<int,int> >::iterator,bool> ret_pair;
            ret_edge = map_edges.insert( std::pair< int, std::set<std::pair<int,int> > >( idface, std::set<std::pair<int,int> >() ) ); // get the set if element 'idface' already existed
            ret_pair = ret_edge.first->second.insert(std::pair<int,int> (v0->getIndex(),v1->getIndex()));
//             printf("Insert in face %d the edge (%d,%d)\n",idface, v0->getIndex(),v1->getIndex());

            if (ret_pair.second)
            {
              std::pair< std::map< LVertex*, std::map<int,int> >::iterator,bool> ret_vertex;
              std::pair< std::map<int,int>::iterator,bool> ret_count;

              ret_vertex = map_vertices.insert(std::pair<LVertex*, std::map<int,int> >(v0, std::map<int,int>()) ); // get the set if element 'v0' already existed
              ret_count = ret_vertex.first->second.insert(std::pair<int,int>(idface,0) );
              ++ret_count.first->second;
//               printf("Vertex0 %d in face %d with counter=%d\n",v0->getIndex(),idface, ret_count.first->second);
              ret_vertex = map_vertices.insert(std::pair<LVertex*, std::map<int,int> >(v1, std::map<int,int>()) ); // get the set if element 'v1' already existed
              ret_count = ret_vertex.first->second.insert(std::pair<int,int>(idface,0) );
              ++ret_count.first->second;
//               printf("Vertex1 %d in face %d with counter=%d\n",v1->getIndex(),idface, ret_count.first->second);
            }
          }
        }
      }
    }

  }

  tprintf ( "--> mesh vertices referenced = %d\n", ( int ) _lmv.size() );
  tprintf ( "--> create vertices ...\n" );

  // mesh topovertices
  for ( int i = 0; i < _mytopo->topovertex_size(); i++ )
  {
    TopoVertex *tv = _mytopo->get_topovertex ( i );
    GVertex *gv = new discreteVertex ( _gm, i );
    if ( _opts->brep_topo )
    {
      assert ( tv->potential_brep_index_size() == 1 );
      gv->addPhysicalEntity ( tv->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      gv->addPhysicalEntity ( tv->get_physical ( 0 ) );
//             gv->addPhysicalEntity(i+1);
    }
    assert ( gv != NULL );
    tvgv.insert ( std::pair<TopoVertex*, GVertex*> ( tv, gv ) );
    MVertex *mv = get_vertex ( tv->get_vertex() );

    MPoint *mpt = NULL;

    LTetra *rt = tv->get_vertex()->getTetra();
    LTetra *at = rt->getAncestor();

    if ( at )
    {
      itlt = ltsbt.find ( rt );
      assert ( itlt != ltsbt.end() );
      MSubTetrahedron *st1 = itlt->second;
      mpt = new MSubPoint ( mv, 0, 0, false, get_tetra ( at ) );
    }
    else
    {
      mpt = new MPoint ( mv );
    }
    gv->points.push_back ( mpt );
    _gm->add ( gv );
  }

  if ( _opts->read_edges )
  {
    // select endpoints of implicite edges
    tprintf ( "--> select endpoints of implicite edges in %d vertices\n", (int)map_vertices.size() );
    std::map< LVertex*, std::map<int,int> >::iterator it_vertex;
    for (it_vertex=map_vertices.begin(); it_vertex!=map_vertices.end(); )
    {
      LVertex* v = it_vertex->first;
      std::map<int,int>::iterator it_physical;
      for (it_physical=it_vertex->second.begin(); it_physical!=it_vertex->second.end(); )
      {
        int idface = it_physical->first;
        int count = it_physical->second;
  //       printf("Vertex %d (%+lf,%+lf,%+lf) in face %d with counter=%d\n",v->getIndex(),v->x(),v->y(),v->z(),idface, count);

        if (count != 1)
          it_vertex->second.erase(it_physical++);
        else
          ++it_physical;
      }
      if (it_vertex->second.size()==0)
        map_vertices.erase(it_vertex++);
      else
        ++it_vertex;
    }

    // mesh endpoints of implicite edges and map endpoints
    if (map_vertices.size())
      tprintf ( "--> create endpoints of edges inserted = %d\n", (int)map_vertices.size() );
    int nbMPoints = _mytopo->topovertex_size();
    for (it_vertex=map_vertices.begin(); it_vertex!=map_vertices.end(); it_vertex++ )
    {
      LVertex* v = it_vertex->first;
      GVertex* gv = new discreteVertex ( _gm, ++nbMPoints );
      std::map<int,int>::iterator it_physical;
      for (it_physical=it_vertex->second.begin(); it_physical!=it_vertex->second.end(); it_physical++ )
      {
        int idface = it_physical->first;
  //       int count = it_physical->second;
  //       printf("Vertex %d (%+lf,%+lf,%+lf) in face %d with counter=%d\n",v->getIndex(),v->x(),v->y(),v->z(),idface, count);

        gv->addPhysicalEntity(idface+10000000*nbMPoints);

        std::pair< std::map< int, std::vector<GVertex*> >::iterator,bool> ret_endpoints;
        ret_endpoints = map_endpoints.insert(std::pair< int, std::vector<GVertex*> >(idface, std::vector<GVertex*>() ) );
        ret_endpoints.first->second.push_back(gv);
      }
      LTetra *rt = v->getTetra();
      LTetra *at = rt->getAncestor();
      MVertex *mv = get_vertex(v);
      MPoint *mpt = NULL;
      if (at)
      {
        mpt = new MSubPoint ( mv, 0, 0, false, get_tetra(at) );
      }
      else
      {
        mpt = new MPoint ( mv );
      }
      gv->points.push_back ( mpt );
      _gm->add ( gv );
    }

  }

  tprintf ( "--> create edges ...\n" );

  for ( int i = 0; i < _mytopo->topoedge_size(); i++ )
  {
    TopoEdge *te = _mytopo->get_topoedge ( i );
    itv = tvgv.find ( te->get_topovertex ( 0 ) );
    assert ( itv != tvgv.end() );
    GVertex *gv0 = itv->second;
    itv = tvgv.find ( te->get_topovertex ( 1 ) );
    assert ( itv != tvgv.end() );
    GVertex *gv1 = itv->second;
    GEdge *ge = new discreteEdge ( _gm, i, gv0, gv1 );
    if ( _opts->brep_topo )
    {
      assert ( te->potential_brep_index_size() == 1 );
      ge->addPhysicalEntity ( te->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      std::vector<int> list;
      if ( te->get_parent() != NULL )
      {
        TopoEdge *root = te->get_parent();
        int root_id_edge = root->get_physical ( 0 );
        assert ( em->is_restore_edge ( root_id_edge ) );
        list.push_back ( root_id_edge );
        dbgprintf ( 3, "--> root id edge = %d\n", root_id_edge );
        std::vector<int> assoc = em->get_associated_edges ( root_id_edge );
        dbgprintf ( 3, "--> assoc size = %d\n", ( int ) assoc.size() );
        list.insert ( list.begin(), assoc.begin(), assoc.end() );
      }

      int local_id_edge = te->get_physical ( 0 );
      list.push_back ( local_id_edge );
      dbgprintf ( 3, "--> local id edge = %d\n", local_id_edge );

      dbgprintf ( 3, "--> edge list size = %d\n", ( int ) list.size() );
      for ( int k = 0; k < list.size(); k++ )
      {
        local_id_edge = list[k];
        int edge = em->get_edge_id ( local_id_edge );
        dbgprintf ( 3, "adding local_id_edge --> %d (%d)\n", local_id_edge, edge );
        ge->addPhysicalEntity ( edge );
      }
    }
    assert ( ge != NULL );
    tege.insert ( std::pair<TopoEdge*, GEdge*> ( te, ge ) );
    for ( int j = 0; j < te->edge_size(); j++ )
    {
      CutEdge *ce = te->get_edge ( j, 1 );
      MVertex *mv0 = get_vertex ( ce->get_vertex ( 0 ) );
      MVertex *mv1 = get_vertex ( ce->get_vertex ( 1 ) );

      LTetra *rt = ce->get_rand_tetra();
      LTetra *at = rt->getAncestor();

      MLine *ml = NULL;

      if ( at )
      {
        itlt = ltsbt.find ( rt );
        assert ( itlt != ltsbt.end() );
        MSubTetrahedron *st1 = itlt->second;
        ml = new MSubLine ( mv0, mv1, 0, 0, false, get_tetra ( at ) );
      }
      else
      {
        ml = new MLine ( mv0, mv1 );
      }

      ge->addLine ( ml );
    }
    _gm->add ( ge );
  }

  if ( _opts->read_edges )
  {
    // mesh implicite edges
    if (map_edges.size())
      tprintf ( "--> create %d edges inserted\n", (int)map_edges.size() );
    int nbMLines = _mytopo->topoedge_size();
    std::map< int, std::set<std::pair<int,int> > >::iterator it_edge;
    std::map< int, std::vector<GVertex*> >::iterator it_endpoint=map_endpoints.begin();
    assert(map_edges.size()==map_endpoints.size());
    for (it_edge=map_edges.begin(); it_edge!=map_edges.end(); it_edge++, it_endpoint++ )
    {
      int idface = it_edge->first;
      assert(idface==it_endpoint->first);
      GVertex* gv0 = it_endpoint->second[0];
      GVertex* gv1 = it_endpoint->second[1];
      GEdge *ge = new discreteEdge ( _gm, ++nbMLines, gv0, gv1 );

      ge->addPhysicalEntity (idface);
      tprintf ( "     +--> edge inserted = %d (%d,%d) contains %d polylines\n", idface, gv0->getPhysicalEntities()[0], gv1->getPhysicalEntities()[0], (int)it_edge->second.size());

      std::set<std::pair<int,int> >::iterator it_pair;
      int i=0;
      for (it_pair=it_edge->second.begin(); it_pair!=it_edge->second.end(); it_pair++ )
      {
        LVertex *v0 = _mymesh->get_vertex(it_pair->first);
        LVertex *v1 = _mymesh->get_vertex(it_pair->second);

//         tprintf ( "     +--> %4d [ %8d -- %8d ] diff = %d\n", i++, v0->getIndex(), v1->getIndex(), v1->getIndex() - v0->getIndex());
//         if (v1->getIndex() - v0->getIndex()>1)
//           printf("    %d (%g %g %g)   %d (%g %g %g)\n", v0->getIndex(), v0->x(), v0->y(), v0->z(), v1->getIndex(), v1->x(), v1->y(), v1->z());

        // find common tetra
        LTetra *rt = NULL;
        Ball b0, b1;
        b0.set_vertex(v0); b1.set_vertex(v1);
        b0.build(); b1.build();
        for (int i=0; i<b0.size() && rt==NULL; ++i)
        {
          for (int j=0; j<b1.size() && rt==NULL; ++j)
          {
            LTetra *tb0, *tb1;
            int eb0, eb1;
            b0.get_tetra ( i, &tb0, &eb0 );
            b1.get_tetra ( j, &tb1, &eb1 );
            if (tb0==tb1)
            {
              rt = tb0;
              break;
            }
          }
        }
        assert(rt);
        LTetra *at = rt->getAncestor();
        MVertex *mv0 = get_vertex(v0);
        MVertex *mv1 = get_vertex(v1);
        MLine *ml = NULL;
        if (at)
        {
          ml = new MSubLine( mv0, mv1, 0, 0, false, get_tetra(at) );
        }
        else
        {
          ml = new MLine ( mv0, mv1 );
        }
        ge->addLine ( ml );
      }
      _gm->add ( ge );
    }
  }

  tprintf ( "--> create faces ...\n" );

  for ( int i = 0; i < _mytopo->topoface_size(); i++ )
  {
    TopoFace *tf = _mytopo->get_topoface ( i );
    GFace *gf = new discreteFace ( _gm, i );
    if ( _opts->brep_topo )
    {
      assert ( tf->getFlag() == 1 );
      assert ( tf->potential_brep_index_size() == 1 );
      gf->addPhysicalEntity ( tf->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      std::vector<int> list;
      if ( tf->get_parent() != NULL )
      {
        TopoFace *root = tf->get_parent();
        int root_id_surf = root->get_physical ( 0 );
//                 assert(em->is_restore_face(root_id_surf));
        list.push_back ( root_id_surf );
        dbgprintf ( 3, "--> root id face = %d\n", root_id_surf );
        std::vector<int> assoc = em->get_associated_faces ( root_id_surf );
        list.insert ( list.begin(), assoc.begin(), assoc.end() );
      }

      int local_id_surf = tf->get_physical ( 0 );
      list.push_back ( local_id_surf );
      dbgprintf ( 3, "--> initial id surf = %d\n", local_id_surf );

      dbgprintf ( 3, "--> surf list size = %d\n", ( int ) list.size() );
      for ( int k = 0; k < list.size(); k++ )
      {
        local_id_surf = list[k];
        int surf = em->get_face_id ( local_id_surf );
        dbgprintf ( 3, "adding local_id_surf --> %d (%d)\n", local_id_surf, surf );
        gf->addPhysicalEntity ( surf );
      }
    }
    assert ( gf != NULL );
    tfgf.insert ( std::pair<TopoFace*, GFace*> ( tf, gf ) );
    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      ite = tege.find ( te );
      assert ( ite != tege.end() );
      gf->addEmbeddedEdge ( ite->second );
    }
    std::set<MVertex*> mvl;
    for ( int j = 0; j < tf->face_size(); j++ )
    {
      CutFace *cf = tf->get_face ( j );
      MVertex *mv0 = get_vertex ( cf->get_face_vertex ( 0 ) );
      MVertex *mv1 = get_vertex ( cf->get_face_vertex ( 1 ) );
      MVertex *mv2 = get_vertex ( cf->get_face_vertex ( 2 ) );
      MTriangle *mt = NULL;
      LTetra* t_int = cf->int_tetra();
//             LTetra* t_ext = cf->ext_tetra();

      LTetra* ancestor_int = t_int->getAncestor();
//             LTetra* ancestor_ext = NULL;
//             if (t_ext)
//                 ancestor_ext = t_ext->getAncestor();

      if ( ancestor_int )
      {
        itlt = ltsbt.find ( t_int );
        assert ( itlt != ltsbt.end() );
        MSubTetrahedron *st1 = itlt->second;
//            mt = new MSubTriangle(mv0, mv1, mv2, 0, 0, false, st1);
        mt = new MSubTriangle ( mv0, mv1, mv2, 0, 0, false, get_tetra ( ancestor_int ) );
      }
      else
      {
        mt = new MTriangle ( mv0, mv1, mv2, 0, 0 );
      }

      gf->addTriangle ( mt );
    }
    _gm->add ( gf );
  }
  stat_GModel();
  tprintf ( "GModel building END\n" );
}




void MeshCutter::build_GModel_Old()
{
  tprintf ( "GModel building START\n" );
#ifdef DEFINE_OLD_CUTTER
  std::multimap<LTetra*, MPolyhedron*> ltmp;
  std::multimap<LTetra*, MPolyhedron*>::iterator itlt;
  std::pair<std::multimap<LTetra*, MPolyhedron*>::iterator, std::multimap<LTetra*, MPolyhedron*>::iterator> retlt;
  std::map<MPolyhedron*, int> mpv;
  std::map<MPolyhedron*, int>::iterator itmpv;

  std::map<TopoVertex*, GVertex*> tvgv;
  std::map<TopoVertex*, GVertex*>::iterator itv;
  std::map<TopoEdge*, GEdge*> tege;
  std::map<TopoEdge*, GEdge*>::iterator ite;
  std::map<TopoFace*, GFace*> tfgf;
  std::map<TopoFace*, GFace*>::iterator itf;
  std::map<int, GRegion*> vogr;
  std::map<int, GRegion*>::iterator itr;

  _mymesh->set_all_tetra_flag ( -1 );
  _mymesh->clear_all_tetra_child();

  _gm = new GModel ( "cut_mesh" );

  progress_bar pbar;

  tprintf ( "--> create regions ...\n" );
  for ( int i = 0; i < _mymesh->get_entity_manager()->vol_size(); i++ )
  {
    int vol = _mymesh->get_entity_manager()->get_vol_id ( i );
    GRegion *gr = new discreteRegion ( _gm, vol );
    gr->addPhysicalEntity ( _mytopo->mpvols[vol] );
    vogr.insert ( std::pair<int, GRegion*> ( vol, gr ) );
    _gm->add ( gr );
  }

  int null_region = _mymesh->get_entity_manager()->get_max_vol_id() +1;
  GRegion *nullgr = new discreteRegion ( _gm, null_region );
  nullgr->addPhysicalEntity ( _mytopo->mpvols[null_region] );

  vogr.insert ( std::pair<int, GRegion*> ( null_region, nullgr ) );
  _gm->add ( nullgr );

  std::set<LTetra*> ultlist;
  std::set<LTetra*> mltlist;
  std::vector<LTetra*> oltlist;

//     tprintf("--> checking uncutted tetras are lying on an interface ...\n");
//
//     for (int i = 0; i < _mytopo->topoface_size(); i++)
//     {
//         TopoFace *tf = _mytopo->get_topoface(i);
//
//  for (int j = 0; j < tf->face_size(); j++)
//         {
//             CutFace *cf = tf->get_face(j);
//             LTetra* t_int = cf->int_tetra();
//             LTetra* t_ext = cf->ext_tetra();
//             LTetra* ancestor_int = t_int->getAncestor();
//      if (!ancestor_int)
//      {
//        t_int->add_child(t_int);
//        mltlist.insert(t_int);
//      }
//      if (t_ext && !t_ext->getAncestor())
//      {
//        t_ext->add_child(t_ext);
//        mltlist.insert(t_ext);
//      }
//  }
//     }

  tprintf ( "--> setting tetra children ...\n" );

  for ( int i = 0; i < _mymesh->tetra_size(); i++ )
  {
    LTetra *t = _mymesh->get_tetra ( i );
    t->setFlag ( -1 );
    for ( int j = 0; j < 4; j++ )
    {
      t->getVertex ( j )->setFlag ( -1 );
    }
    LTetra *at = t->getAncestor();

    if ( at )
    {
      assert ( at != NULL );

      at->setFlag ( -1 );
      for ( int j = 0; j < 4; j++ )
        at->getVertex ( j )->setFlag ( -1 );

//     printf("tetra %p --> %d children\n", at, at->child_size());
      at->add_child ( t );
    }
  }

  for ( int i = 0; i < _mymesh->tetra_size(); i++ )
  {
    LTetra *t = _mymesh->get_tetra ( i );
    LTetra *at = t->getAncestor();
//  printf("tetra %p (%d) --> parent %p\n", t, t->getIndex(), at);
//  assert(at != NULL);

    if ( at )
    {
      if ( at->child_size() == 1 )
        mltlist.insert ( at );
      else
      {
        std::set<int> vollist;
        for ( int j = 0; j < at->child_size(); j++ )
          vollist.insert ( at->get_child ( j )->getVol() );

        if ( vollist.size() == 1 )
        {
          int vol = * ( vollist.begin() );
          int ult = 1;
          for ( int j = 0; j < at->child_size() && ult; j++ )
            for ( int k = 0; k < 4; k++ )
              if ( at->get_child ( j )->getAdj ( k ) != NULL && at->get_child ( j )->getAdj ( k )->getVol() != vol )
                ult = 0;
          if ( ult )
            ultlist.insert ( at );
          else
            mltlist.insert ( at );
        }
        else
          mltlist.insert ( at );
      }
    }
    else if ( !t->isParent() /*&& !t->child_size()*/ )
      oltlist.push_back ( t );
  }

  tprintf ( "--> create %d end tetras ...\n", ( int ) oltlist.size() );

  pbar.start ( oltlist.size() );
  for ( int i = 0; i < oltlist.size(); i++ )
  {
    pbar.progress ( i );
    LTetra *t = oltlist[i];
    itr = vogr.find ( t->getVol() );
    assert ( itr != vogr.end() );
    GRegion *gr = itr->second;
    add_tetra ( t, gr );
  }
  pbar.end();

  tprintf ( "--> create %d single volume tetra parents ...\n", ( int ) ultlist.size() );

  for ( std::set<LTetra*>::iterator it = ultlist.begin(); it != ultlist.end(); it++ )
  {
    LTetra *t = ( *it );
    assert ( t->child_size() );
    itr = vogr.find ( t->get_child ( 0 )->getVol() );
    assert ( itr != vogr.end() );
    GRegion *gr = itr->second;
    add_tetra ( t, gr );
  }

  tprintf ( "--> create %d multiple volume tetra parents + children ...\n", ( int ) mltlist.size() );
  int count = 0;
  pbar.start ( mltlist.size() );
  for ( std::set<LTetra*>::iterator it = mltlist.begin(); it != mltlist.end(); it++ )
  {
    pbar.progress ( count++ );
    LTetra *t = ( *it );
//  printf("tetra %p (%d) --> %d children\n", t, t->getIndex(), t->child_size());
//  for (int i = 0; i < t->child_size(); i++)
//    printf("child tetra %p --> index %d, flag %d\n", t->get_child(i), t->get_child(i)->getIndex(), t->get_child(i)->getFlag());

    assert ( t->child_size() );
    MTetrahedron *mt = add_tetra ( t, nullgr );
    std::set<int> vollist;
    for ( int i = 0; i < t->child_size(); i++ )
    {
      LTetra *tc = t->get_child ( i );
      assert ( tc != NULL );
//      printf("child tetra %p (%d) --> flag %d\n", tc, tc->getIndex(), tc->getFlag());
//      if (t != tc)
      MTetrahedron *mtc = add_tetra ( tc, nullgr );
      vollist.insert ( tc->getVol() );
    }
#if 1
    // polyhedra/polygon generation
    for ( std::set<int>::iterator it = vollist.begin(); it != vollist.end(); it++ )
    {
      int current_vol = ( *it );
      std::vector<MTetrahedron*> mtlist;
      for ( int i = 0; i < t->child_size(); i++ )
      {
        if ( t->get_child ( i )->getVol() == current_vol )
          mtlist.push_back ( get_tetra ( t->get_child ( i ) ) );
      }
      itr = vogr.find ( current_vol );
      assert ( itr != vogr.end() );
      GRegion *grc = itr->second;
      MPolyhedron *mp = new MPolyhedron ( mtlist, 0, 0, 0, mt );
      ltmp.insert ( std::pair<LTetra*, MPolyhedron*> ( t, mp ) );
      mpv.insert ( std::pair<MPolyhedron*, int> ( mp, current_vol ) );
      grc->polyhedra.push_back ( mp );
    }
#endif
  }
  pbar.end();

  tprintf ( "--> mesh vertices referenced = %d\n", ( int ) _lmv.size() );
  tprintf ( "--> create vertices ...\n" );

  for ( int i = 0; i < _mytopo->topovertex_size(); i++ )
  {
    TopoVertex *tv = _mytopo->get_topovertex ( i );
    GVertex *gv = NULL;
    if ( _opts->brep_topo )
    {
      assert ( tv->potential_brep_index_size() == 1 );
      gv = new discreteVertex ( _gm, tv->get_potential_brep_index ( 0 ) );
      gv->addPhysicalEntity ( tv->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      gv = new discreteVertex ( _gm, i );
      gv->addPhysicalEntity ( i+1 );
    }
    assert ( gv != NULL );
    tvgv.insert ( std::pair<TopoVertex*, GVertex*> ( tv, gv ) );
    MVertex *mv = get_vertex ( tv->get_vertex() );
    MPoint *mpt = new MPoint ( mv );
    gv->points.push_back ( mpt );
    _gm->add ( gv );
  }

  tprintf ( "--> create edges ...\n" );

  for ( int i = 0; i < _mytopo->topoedge_size(); i++ )
  {
    TopoEdge *te = _mytopo->get_topoedge ( i );
    itv = tvgv.find ( te->get_topovertex ( 0 ) );
    assert ( itv != tvgv.end() );
    GVertex *gv0 = itv->second;
    itv = tvgv.find ( te->get_topovertex ( 1 ) );
    assert ( itv != tvgv.end() );
    GVertex *gv1 = itv->second;
    GEdge *ge = NULL;
    if ( _opts->brep_topo )
    {
      assert ( te->potential_brep_index_size() == 1 );
      ge = new discreteEdge ( _gm, te->get_potential_brep_index ( 0 ), gv0, gv1 );
      ge->addPhysicalEntity ( te->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      ge = new discreteEdge ( _gm, i, gv0, gv1 );
//      printf("--> %d\n", te->potential_brep_index_size());
//      assert(te->potential_brep_index_size());
      int edge = i+1;
      if ( te->potential_brep_index_size() )
      {
        int local_id_edge = te->get_potential_brep_index ( 0 );
        edge = _mymesh->get_entity_manager()->get_edge_id ( local_id_edge );
      }
//      printf("--> local_id %d\n", local_id_edge);
      ge->addPhysicalEntity ( edge );
//      printf("--> edge %d\n", edge);
//          ge->addPhysicalEntity(i+1);
    }
    assert ( ge != NULL );
    tege.insert ( std::pair<TopoEdge*, GEdge*> ( te, ge ) );
    for ( int j = 0; j < te->edge_size(); j++ )
    {
      CutEdge *ce = te->get_edge ( j, 1 );
      MVertex *mv0 = get_vertex ( ce->get_vertex ( 0 ) );
      MVertex *mv1 = get_vertex ( ce->get_vertex ( 1 ) );
      MLine *ml = new MLine ( mv0, mv1 );
      ge->addLine ( ml );
    }
    _gm->add ( ge );
  }

  tprintf ( "--> create faces ...\n" );

  EntityManager *em = _mymesh->get_entity_manager();

  for ( int i = 0; i < _mytopo->topoface_size(); i++ )
  {
    TopoFace *tf = _mytopo->get_topoface ( i );
    GFace *gf = NULL;
    if ( _opts->brep_topo )
    {
      assert ( tf->getFlag() == 1 );
      assert ( tf->potential_brep_index_size() == 1 );
      gf = new discreteFace ( _gm, tf->get_potential_brep_index ( 0 ) );
      gf->addPhysicalEntity ( tf->get_potential_brep_index ( 0 ) /*+1*/ );
    }
    else
    {
      gf = new discreteFace ( _gm, i );

#if defined(NEW_PHYSICAL_NUMBERING)
      int local_id_surf = tf->get_surf ( 0 );
//      printf("--> initial id surf = %d\n", local_id_surf);
      std::vector<int> list;
      if ( em->is_restore_face ( local_id_surf ) )
      {
        list = em->get_associated_faces ( local_id_surf );
        int local_id_face = tf->get_potential_brep_index ( 0 );
        list.push_back ( local_id_face );
//    printf("--> initial id face = %d\n", local_id_face);
      }
      list.push_back ( local_id_surf );
//      printf("--> surf list size = %d\n", (int)list.size());
      for ( int k = 0; k < list.size(); k++ )
      {
        local_id_surf = list[k];
//                 printf("adding local_id_surf --> %d\n", local_id_surf);
        int surf = em->get_face_id ( local_id_surf );
//                 printf("surf --> %d\n", surf);
        gf->addPhysicalEntity ( surf );
      }
#else
      int local_id_surf = tf->get_surf ( 0 );
//                 printf("local_id_surf --> %d\n", local_id_surf);
      int surf = _mymesh->get_entity_manager()->get_face_id ( local_id_surf );
      gf->addPhysicalEntity ( surf );
      int maxid = _mymesh->get_entity_manager()->get_max_face_id();
      gf->addPhysicalEntity ( maxid+i+1 );
#endif

    }
    assert ( gf != NULL );
    tfgf.insert ( std::pair<TopoFace*, GFace*> ( tf, gf ) );
    for ( int j = 0; j < tf->topoedge_size(); j++ )
    {
      TopoEdge *te = tf->get_topoedge ( j );
      ite = tege.find ( te );
      assert ( ite != tege.end() );
      gf->addEmbeddedEdge ( ite->second );
    }
    std::set<MVertex*> mvl;
    for ( int j = 0; j < tf->face_size(); j++ )
    {
      CutFace *cf = tf->get_face ( j );
      MVertex *mv0 = get_vertex ( cf->get_face_vertex ( 0 ) );
      MVertex *mv1 = get_vertex ( cf->get_face_vertex ( 1 ) );
      MVertex *mv2 = get_vertex ( cf->get_face_vertex ( 2 ) );
      MTriangle *mt = NULL;
      LTetra* t_int = cf->int_tetra();
      LTetra* t_ext = cf->ext_tetra();
      LTetra* ancestor_int = t_int->getAncestor();
      LTetra* ancestor_ext = NULL;
      if ( t_ext )
        ancestor_ext = t_ext->getAncestor();

      if ( ancestor_int/* && (ancestor_ext == NULL || ancestor_int == ancestor_ext)*/ )
      {
//    assert(t_ext != ancestor_ext);
        int count_vol = ltmp.count ( ancestor_int );
        if ( count_vol == 1 ) // cas particulier
        {
          retlt = ltmp.equal_range ( ancestor_int );
          itlt = retlt.first;
          MPolyhedron *p1 = itlt->second;
          LTetra* t2 = cf->ext_tetra();
          LTetra* at2 = t_int->getAncestor();
          assert ( at2 != NULL );
          retlt = ltmp.equal_range ( at2 );
          itlt = retlt.first;
          MPolyhedron *p2 = itlt->second;
          mt = new MTriangleBorder ( mv0, mv1, mv2, 0, 0, p1, p2 );
        }
        else
        {
          retlt = ltmp.equal_range ( ancestor_int );
          itlt = retlt.first;
          MPolyhedron *p1 = NULL;
          MPolyhedron *p2 = NULL;

          for ( int k = 0; k < count_vol; k++ )
          {
            MPolyhedron *tmp = itlt->second;
            int vol = mpv.find ( tmp )->second;
            if ( t_int->getVol() == vol )
              p1 = tmp;
            if ( t_ext && t_ext->getVol() == vol )
              p2 = tmp;

            itlt++;
          }
          assert ( p1 != NULL );
//                     assert(p2 != NULL);
          assert ( p1 != p2 );
          mt = new MTriangleBorder ( mv0, mv1, mv2, 0, 0, p1, p2 );
        }
      }
      else
      {
//          // TODO : convert to MTriangleBorder
//                 mt = new MTriangle(mv0, mv1, mv2);

        printf ( "tetra int = %d\n", t_int->getIndex() );
        retlt = ltmp.equal_range ( t_int );
        itlt = retlt.first;
        MPolyhedron *p1 = NULL;
        MPolyhedron *p2 = NULL;

        p1 = itlt->second;
        if ( t_ext )
          p2 = itlt->second;
        assert ( p1 != NULL );
//                 assert(p2 != NULL);
        assert ( p1 != p2 );
        mt = new MTriangleBorder ( mv0, mv1, mv2, 0, 0, p1, p2 );

      }
      gf->addTriangle ( mt );
    }
    _gm->add ( gf );
  }

  stat_GModel();
#else
  tprintf ( "Warning - Old gmsh interface not compiled\n" );
#endif
  tprintf ( "GModel building END\n" );
}