// 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"
#include <fstream>
#include <iostream>

int MeshCutter::restore()
{
  char filename[100];
  sprintf ( filename, "%s.ent", _opts->basename );
  tprintf ( "Reading entities file %s ....\n", filename );
  FILE *fp = fopen ( filename, "r" );

  int error = 0;

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

  tprintf ( "--> %d face to recover ...\n", nb_face );

  for ( int i = 0; i < nb_face; i++ )
  {
    int idsurf;
    int local_id;
    if ( fscanf ( fp, "%d", &idsurf ) != 1 )
      error++;

    local_id = _mymesh->get_entity_manager()->add_face_id ( idsurf );
    _mymesh->get_entity_manager()->set_restore_face ( local_id );

    int nbv;
    if ( fscanf ( fp, "%d", &nbv ) != 1 )
      error++;
    for ( int j = 0; j<nbv; j++ )
    {
      int id;
      double dist;
      if ( fscanf ( fp, "%d %lf", &id, &dist ) != 2 )
        error++;
      LVertex *v = _mymesh->get_vertex ( id );
      dist *= -1.0; // WARNING --> hack
      if ( !v->checkSurf ( local_id ) )
      {
        v->addSurf ( local_id, dist );
//                 if (fabs(dist) < TOL)
        if ( dist == 0.0 )
          v->addSurfZero ( local_id );
      }
    }
  }

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

  tprintf ( "--> %d edges to recover ...\n", nb_edge );

  for ( int i = 0; i < nb_edge; i++ )
  {
    int idedge;
    int local_id;
    int nbv;
    if ( fscanf ( fp, "%d %d", &idedge, &nbv ) != 2 )
      error++;

    local_id = _mytopo->get_entity_manager()->add_edge_id ( idedge );

    int id;
    if ( fscanf ( fp, "%d", &id ) != 1 )
      error++;
    LVertex *v0 = _mymesh->get_vertex ( id );
    LVertex *v1 = NULL;
    for ( int j = 1; j<nbv; j++ )
    {
      if ( fscanf ( fp, "%d", &id ) != 1 )
        error++;
      v1 = _mymesh->get_vertex ( id );
      _mytopo->add_existing_edge ( v0, v1, local_id );
      v0 = v1;
    }
  }
  return 0;
}

MeshCutter::MeshCutter ( options *opt )
{
  _opts = opt;
  _gm = NULL;
  _mymesh = new Mesh ( _opts );
  if ( _opts->generate_mesh )
    _mycadls = new cadLevelset ( _opts, _mymesh );
  else
    _mycadls = NULL;
  _mysurface = new LSurface ( _opts, _mymesh->get_entity_manager() );
  _mytopo = new Topo ( _opts, _mymesh->get_entity_manager() );
  _myquery = NULL;
}

void MeshCutter::load_GModel_geometry ( GModel *gm, int face_max_sample )
{
  int count_face = 0;
  for ( GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it )
  {
    GFace *gf = *it;
    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, count_face++, umin, umax, vmin, vmax, face_max_sample );
    surf->generate();
  }
}

int MeshCutter::process()
{
  if ( _opts->read_fiber_surface )
  {
    _mycadls->readLsnFile();
  }

  if ( _opts->load_gmodel_geometry )
  {
    GModel *gm = _mycadls->get_gmodel();
    _mycadls->read_step();
    _mycadls->build_surfaces();
    _mytopo->get_topo_from_gmodel ( gm );
  }

  if ( _opts->generate_mesh )
  {
    assert ( !_opts->read_meshfile );
    assert ( !_opts->read_lsfile );
    assert ( !_opts->load_ibrep );

    _mycadls->build_mesh_box();
    if ( _opts->refine_mesh_curvature )
      _mycadls->refine_mesh();
    _mycadls->process();
  }

#if defined(CUTMESH_ENABLE_IBREP)
  IBrep ib;
  if ( _opts->load_ibrep )
  {
    assert ( !_opts->read_meshfile );
    assert ( !_opts->read_lsfile );

    char filename[100];
    sprintf ( filename, "%s.ibrep", _opts->basename );
    std::ifstream inputf ( filename );
    inputf >> ib;
    _mymesh->ImportIbrep ( &ib );
    _mytopo->ImportIbrep ( &ib );
  }
#endif

  if ( _opts->load_gmodel_mesh && !_opts->library_mode )
  {
    char filename[100];
    sprintf ( filename, "%s.msh", _opts->basename );
    GModel *gm = new GModel;
    gm->readMSH ( filename );
    load_GModel_mesh ( gm );
  }

  if ( _opts->read_meshfile )
    _mymesh->readMeshFile ( _opts->basename );

//     if (_opts->recut_enable)
//       restore();

  if ( _opts->read_lsfile && !_opts->library_mode )
    _mymesh->readLevelsetFile ( _opts->basename );

  if ( _opts->read_topo_file && !_opts->library_mode )
    _mytopo->read_topo_file ( _opts->basename );

  _mymesh->cut_mesh();
//     _mymesh->optimize();

//     _mymesh->writeMesh();

  _mysurface->build_surface ( _mymesh );

  _mytopo->build_topo ( _mymesh, _mysurface );
  _mysurface->clear();

  if ( _opts->brep_topo )
    _mytopo->clean_topo();

//     if (_opts->brep_topo || _opts->auto_topo)
//         _mytopo->export_existing_entities(_opts->basename);

  _mytopo->build_volume ( _mymesh );

  if ( _opts->recut_enable )
    _mytopo->build_hierarchy ( _mymesh );

  _myquery = new Topology_query ( _mytopo );

#ifdef DEFINE_OLD_CUTTER
  build_GModel()_Old;

  if ( _opts->write_gmodel && !_opts->library_mode )
  {
    GModel *gm = get_gmodel();
    char filename[100];
    sprintf ( filename, "%s_cut.msh", _opts->basename );
    tprintf ( "Writing %s\n", filename );
    gm->writeMSH ( filename );
  }
#else
  build_GModel();
  if ( _opts->write_gmodel && !_opts->library_mode )
  {
    GModel *gm = get_gmodel();
    char filename[100];
    sprintf ( filename, "%s_cut.msh", _opts->basename );
    tprintf ( "Writing %s\n", filename );
    gm->writeMSH ( filename, 3 );
  }
#endif
//     _myquery->minishell();
  return 0;
}

void MeshCutter::physical_query ( const char *query, std::vector<int> phys[4] )
{
  std::vector<char*> exp;
  _myquery->split_line ( query, exp );
  _myquery->set_expression ( exp );
  if ( exp.size() )
  {
    _myquery->set_expression ( exp );
    _myquery->compute_expression();
    _myquery->dump ( _myquery->get_result() );
  }
  phys[3] = _myquery->get_result_volumes_physicals();
  phys[2] = _myquery->get_result_faces_physicals();
  phys[1] = _myquery->get_result_edges_physicals();
  phys[0] = _myquery->get_result_vertices_physicals();
}


void MeshCutter::load_GModel_mesh ( GModel *gm )
{
  _gm = gm;

  stat_GModel();

  //  get vertices
  std::vector<GEntity*> gm_ent;
  _gm->getEntities ( gm_ent );
  int num = 0;
  for ( unsigned int i = 0; i < gm_ent.size(); i++ )
  {
    for ( unsigned int j = 0; j < gm_ent[i]->getNumMeshVertices(); j++ )
    {
      MVertex *mv = gm_ent[i]->getMeshVertex ( j );

      if ( !_mymesh->get_vertex_from_lookup ( mv->getNum() ) )
      {
        LVertex *lv = _mymesh->new_vertex ( mv->x(), mv->y(), mv->z() );
        _mymesh->add_vertex ( lv, mv->getNum() );
      }
    }
  }

  tprintf ( "mesh --> %d vertices\n", _mymesh->vertex_size() );

  std::map<int , LTetra*> parents;
  std::map<int , LTetra*>::iterator itp;

  std::map<std::set<LVertex *> , MTetrahedron*> vertex_tetra_map;
  std::map<std::set<LVertex *> , MTetrahedron*>::iterator itvt;

  //  get tetras

  std::multimap<int, GRegion*> gregionphysmap;
  std::set<int> gregionphysset;
  for ( GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it )
  {
    GRegion *gr = *it;
    assert ( gr->getPhysicalEntities().size() );
    for ( int i = 0; i < gr->getPhysicalEntities().size(); i++ )
    {
      int idvol = gr->getPhysicalEntities() [i];
      dbgprintf ( 3, "vol tag --> %d\n", idvol );
      int local_id = _mymesh->get_entity_manager()->add_vol_id ( idvol );
      dbgprintf ( 3, "local id vol --> %d\n", local_id );
      _mymesh->get_entity_manager()->set_restore_vol ( local_id );
      gregionphysmap.insert ( std::pair<int, GRegion*> ( local_id, gr ) );
      gregionphysset.insert ( local_id );
    }
  }

  std::set<int>::iterator itgrps;
  for ( itgrps = gregionphysset.begin(); itgrps != gregionphysset.end(); itgrps++ )
  {
    int id = ( *itgrps );
//         dbgprintf(3, "vol id --> %d\n", id);
    if ( gregionphysmap.count ( id ) == 1 )
    {
      dbgprintf ( 3, "adding volume %d\n", id);
      std::multimap<int, GRegion*>::iterator it = gregionphysmap.equal_range ( id ).first;
      GRegion *gr = it->second;
      dbgprintf ( 3, "--> %d physicals\n", (int)gr->getPhysicalEntities().size());
      for ( int i = 0; i < gr->getPhysicalEntities().size(); i++ )
      {
        int idvol = gr->getPhysicalEntities() [i];
        int local_id = _mymesh->get_entity_manager()->add_vol_id ( idvol );
        if ( local_id != id )
          _mymesh->get_entity_manager()->associate_vol ( id, local_id );
      }

      for ( int i = 0; i < gr->tetrahedra.size(); i++ )
      {
        MTetrahedron *mt = gr->tetrahedra[i];
        LTetra *ltp = NULL;

        if ( mt->getTypeForMSH() == MSH_TET_SUB )
        {
          MSubTetrahedron *mtc = dynamic_cast<MSubTetrahedron*> ( mt );
          assert ( mtc->getParent() );
          if ( mtc->ownsParent() )
          {
            MTetrahedron *mtp = ( MTetrahedron * ) mtc->getParent();
            LVertex *vp0 = _mymesh->get_vertex_from_lookup ( mtp->getVertex ( 0 )->getNum() );
            LVertex *vp1 = _mymesh->get_vertex_from_lookup ( mtp->getVertex ( 1 )->getNum() );
            LVertex *vp2 = _mymesh->get_vertex_from_lookup ( mtp->getVertex ( 2 )->getNum() );
            LVertex *vp3 = _mymesh->get_vertex_from_lookup ( mtp->getVertex ( 3 )->getNum() );

            LTetra *ltp = _mymesh->new_tetra ( vp0, vp1, vp2, vp3, NULL );
            ltp->setOldVol ( -1 );
          }
        }

        LVertex *v0 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 0 )->getNum() );
        LVertex *v1 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 1 )->getNum() );
        LVertex *v2 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 2 )->getNum() );
        LVertex *v3 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 3 )->getNum() );

        LTetra *lt = _mymesh->new_tetra ( v0, v1, v2, v3, ltp );
        lt->setOldVol ( id );
        _mymesh->add_tetra ( lt, mt->getNum() );
      }
    }
  }

  tprintf ( "tetra size --> %d\n", _mymesh->tetra_size() );

  hashFaceStruct *hfs = new hashFaceStruct;

  for ( int i = 0; i < _mymesh->tetra_size(); i++ )
  {
    LTetra *t = _mymesh->get_tetra ( i );
    assert ( !t->isParent() );
    for ( int j = 0; j < 4; j++ )
      hfs->addFace ( t, j );
  }

  std::multimap<int, GFace*> gfacephysmap;
  std::set<int> gfacephysset;
  for ( GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it )
  {
    GFace *gf = *it;
    assert ( gf->getPhysicalEntities().size() );
    for ( int i = 0; i < gf->getPhysicalEntities().size(); i++ )
    {
      int idsurf = gf->getPhysicalEntities() [i];
      dbgprintf ( 3, "face tag --> %d\n", idsurf );
      int local_id = _mymesh->get_entity_manager()->add_face_id ( idsurf );
      dbgprintf ( 3, "local id face --> %d\n", local_id );
      _mymesh->get_entity_manager()->set_restore_face ( local_id );
      gfacephysmap.insert ( std::pair<int, GFace*> ( local_id, gf ) );
      gfacephysset.insert ( local_id );
    }
  }

  std::set<int>::iterator itgfps;
  for ( itgfps = gfacephysset.begin(); itgfps != gfacephysset.end(); itgfps++ )
  {
    int id = ( *itgfps );
    if ( gfacephysmap.count ( id ) == 1 )
    {
      dbgprintf ( 3, "adding face %d\n", id);
      std::multimap<int, GFace*>::iterator it = gfacephysmap.equal_range ( id ).first;
      GFace *gf = it->second;
      dbgprintf ( 3, "--> %d physicals\n", (int)gf->getPhysicalEntities().size());
      for ( int i = 0; i < gf->getPhysicalEntities().size(); i++ )
      {
        int idsurf = gf->getPhysicalEntities() [i];
        int local_id = _mymesh->get_entity_manager()->add_face_id ( idsurf );

        if ( local_id != id )
          _mymesh->get_entity_manager()->associate_face ( id, local_id );
      }
      for ( int i = 0; i < gf->triangles.size(); i++ )
      {
        MTriangle *mt = gf->triangles[i];

        LVertex *v0 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 0 )->getNum() );
        LVertex *v1 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 1 )->getNum() );
        LVertex *v2 = _mymesh->get_vertex_from_lookup ( mt->getVertex ( 2 )->getNum() );

        hashFaceKey *hk = hfs->getFace ( v0->getIndex(), v1->getIndex(), v2->getIndex() );

        assert ( hk );

        LTetra *t = hk->tetra();
        int face = hk->face();

        assert ( t != NULL );

        if ( t->getAdj ( face ) )
        {
          // check if face is correctly oriented
          for ( int j = 0; j < 3; j++ )
          {
            if ( t->getFaceVertex ( face, j ) == v0 )
            {
              if ( t->getFaceVertex ( face, ( j+1 ) %3 ) != v1 )
              {
                LTetra *tmp_t = t->getAdj ( face );
                int tmp_face = t->getOpp ( face );
                t = tmp_t;
                face = tmp_face;
                break;
              }
            }
          }
        }

        for ( int j = 0; j < 3; j++ )
        {
          LVertex *v = t->getFaceVertex ( face, j );
          if ( !v->checkSurf ( id ) )
          {
            v->addSurf ( id, 0.0 );
            v->addSurfZero ( id );
          }
        }

        LVertex *vf = t->getVertex ( face );
        if ( !vf->checkSurf ( id ) )
          vf->addSurf ( id, -1.0 ); // < 0 --> inside

        LTetra *adj = t->getAdj ( face );
        int opp = t->getOpp ( face );

        if ( adj )
        {
          LVertex *vfo = adj->getVertex ( opp );
          if ( !vfo->checkSurf ( id ) )
          {
            vfo->addSurf ( id, 1.0 ); // > 0 --> outside
          }
        }
      }
    }
  }

  _mymesh->compute_adjacencies();

  std::multimap<int, GEdge*> gedgephysmap;
  std::set<int> gedgephysset;
  for ( GModel::eiter it = _gm->firstEdge(); it != _gm->lastEdge(); ++it )
  {
    GEdge *ge = *it;
    assert ( ge->getPhysicalEntities().size() );
    for ( int i = 0; i < ge->getPhysicalEntities().size(); i++ )
    {
      int idedge = ge->getPhysicalEntities() [i];
      dbgprintf ( 3, "edge tag --> %d\n", idedge );
      int local_id = _mymesh->get_entity_manager()->add_edge_id ( idedge );
      dbgprintf ( 3, "local id edge --> %d\n", local_id );
      _mymesh->get_entity_manager()->set_restore_edge ( local_id );
      gedgephysmap.insert ( std::pair<int, GEdge*> ( local_id, ge ) );
      gedgephysset.insert ( local_id );
    }
  }

  std::set<int>::iterator itgeps;
  for ( itgeps = gedgephysset.begin(); itgeps != gedgephysset.end(); itgeps++ )
  {
    int id = ( *itgeps );
    if ( gedgephysmap.count ( id ) == 1 )
    {
      dbgprintf ( 3, "adding edge %d\n", id );
      std::multimap<int, GEdge*>::iterator it = gedgephysmap.equal_range ( id ).first;
      GEdge *ge = it->second;
      dbgprintf ( 3, "--> %d physicals\n", ( int ) ge->getPhysicalEntities().size() );
      for ( int i = 0; i < ge->getPhysicalEntities().size(); i++ )
      {
        int idedge = ge->getPhysicalEntities() [i];
        int local_id = _mymesh->get_entity_manager()->add_edge_id ( idedge );
        if ( local_id != id )
          _mymesh->get_entity_manager()->associate_edge ( id, local_id );
      }
      for ( int i = 0; i < ge->lines.size(); i++ )
      {
        MLine *ml = ge->lines[i];

        LVertex *v0 = _mymesh->get_vertex_from_lookup ( ml->getVertex ( 0 )->getNum() );
        LVertex *v1 = _mymesh->get_vertex_from_lookup ( ml->getVertex ( 1 )->getNum() );
        _mytopo->add_existing_edge ( v0, v1, id );
      }
    }
  }
}

MVertex *MeshCutter::get_vertex ( LVertex *v )
{
  int index = v->getFlag();
//     printf("index --> %d\n", index);
  assert ( index != -1 );
  return _lmv[index];
}

MTetrahedron *MeshCutter::get_tetra ( LTetra *t )
{
  int index = t->getFlag();
  assert ( index != -1 );
  return _lmt[index];
}

MVertex *MeshCutter::add_vertex ( LVertex *v, GRegion *gr )
{
  MVertex *mv = NULL;
  int index = v->getFlag();
  if ( index == -1 )
  {
    mv = new MVertex ( v->x(), v->y(), v->z() /*, NULL, v->getIndex()+1*/ );
//         printf("index = %d\n", mv->getIndex());
    v->setFlag ( _lmv.size() );
    _lmv.push_back ( mv );
    gr->addMeshVertex ( mv );
  }
  else
  {
    mv = _lmv[index];
  }
  return mv;
}

MTetrahedron *MeshCutter::add_tetra ( LTetra *t, GRegion *gr )
{
//     printf("add tetra %d\n", t->getIndex());
  assert ( t->getFlag() == -1 );
  std::vector<MVertex*> mvl;
  for ( int j = 0; j < 4; j++ )
  {
    LVertex *v = t->getVertex ( j );
    MVertex *mv = add_vertex ( v, gr );
    mvl.push_back ( mv );
  }
  MTetrahedron *mt = NULL;
  if ( t->getParent() )
    mt = new MSubTetrahedron ( mvl,false,mt );
  else
    mt = new MTetrahedron ( mvl, 0, 0 );
//     MTetrahedron *mt = new MTetrahedron(mvl, 0, 0);
  t->setFlag ( _lmt.size() );
//     printf("tetra %p --> flag %d\n", t, t->getFlag());
  _lmt.push_back ( mt );
  gr->tetrahedra.push_back ( mt );
  return mt;
}

void MeshCutter::stat_GModel()
{
  // summary
  tprintf ( "GModel contains :\n" );
  tprintf ( "--> %lu vertices\n", _gm->getNumVertices() );
  tprintf ( "--> %lu edges\n", _gm->getNumEdges() );
  tprintf ( "--> %lu faces\n", _gm->getNumFaces() );
  for ( GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it )
  {
    tprintf ( "     +--> face --> %lu triangles\n", ( *it )->triangles.size() );
  }
  tprintf ( "--> %lu regions\n", _gm->getNumRegions() );
  for ( GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); it++ )
  {
    tprintf ( "     +--> region --> %lu tetrahedra\n", ( *it )->tetrahedra.size() );
  }
  tprintf ( "--> %lu mesh elements\n", _gm->getNumMeshElements() );
  tprintf ( "--> %lu mesh vertices\n", _gm->getNumMeshVertices() );
}