/// \author Leblanc Christophe.
/// \date 18/11/2020
/// \mail cleblancad@gmail.com

// No guards.

template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_geo(Surface_mesh &mesh,
  const std::map<Vertex_index, std::vector<int> > &vertices_physical_groups,
  const std::map<Face_index, std::vector<int> > &faces_physical_groups,
  const std::map<Volume_index, std::vector<int> > &volumes_physical_groups) const
{
  // Checks.
  if(m_geo_filename.empty())
  {
    std::cerr << "Error: empty file name for the output geo file.\n";
    return;
  }
  
  if(mesh.number_of_vertices() == 0)
  {
    std::cerr << "Error: mesh is empty.\n";
    return;
  }

  // Open geo file.
  std::ofstream geo_file;
  geo_file.open(m_geo_filename);
  
  if(!geo_file.is_open())
  {
    std::cerr << "Error: could bot open file \'" << m_geo_filename << "\'\n";
    return;
  }

  // Write the geo file.
  write_header(geo_file);
  geo_file << "mesh_size = " << m_mesh_size << ";\n\n";
  
  // Vertices.
  Vertex_map_long vertex_num_map = mesh.template add_property_map<
    Vertex_index, long int>("v:num", 0).first;
  
  write_vertices(geo_file, mesh, vertex_num_map);
  
  // Edges.
  Edge_map_bool edge_visited_map = mesh.template add_property_map<
    Edge_index, bool>("e:is_visited", false).first;
  Edge_map_long edge_num_map = mesh.template add_property_map<
    Edge_index, long int>("e:num", 0).first;
  
  write_edges(geo_file, mesh, vertex_num_map, edge_visited_map, edge_num_map);
  
  // Plane surfaces.
  Face_map_long face_num_map = mesh.template add_property_map<
    Face_index, long>("f:num", 0).first;
  
  write_faces(geo_file, mesh, vertex_num_map, edge_num_map, face_num_map);
  
  // Volume.
  Volume_map_long volume_num_map = mesh.template add_property_map<
    Face_index, long>("vol:num", 0).first;
  
  write_volume(geo_file, mesh, face_num_map, volume_num_map);
  
  // Physical points.
  write_physical_points(geo_file, vertices_physical_groups, vertex_num_map);
  
  // Physical lines (warning: not implemented).
  write_physical_lines(geo_file);
  
  // Physical faces.
  write_physical_faces(geo_file, faces_physical_groups, face_num_map);
  
  // Physical volumes.
  write_physical_volumes(geo_file, volumes_physical_groups, volume_num_map);
  
  mesh.remove_property_map(vertex_num_map);
  mesh.remove_property_map(edge_visited_map);
  mesh.remove_property_map(edge_num_map);
  mesh.remove_property_map(face_num_map);
  mesh.remove_property_map(volume_num_map);
  geo_file.close();
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_header(std::ofstream &file) const
{
  time_t rawtime;
  struct tm* timeinfo;
  time(&rawtime);
  timeinfo = localtime(&rawtime);
  
  file << "/*********************************************************************\n"
     << "*\n"
     << "* Geo file generated from MSH file \"" << m_filename << "\".\n"
     << "* Date: "
     << asctime(timeinfo)
     << "*\n*\n"
     << "*\n"
     << "* Generator written by: Leblanc Christophe\n"
     << "* Email: cleblancad@gmail.com\n"
     << "* University of Liege\n"
     << "*\n"
     << "*********************************************************************/\n\n";
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_vertices(std::ofstream &file, Surface_mesh &mesh,
 Vertex_map_long &vertex_num_map) const
{
  file << "// Vertices.\n\n";
  
  long int vertex_counter = 1;
  typename Surface_mesh::Vertex_range::const_iterator v_it = mesh.vertices().begin(),
    v_ite = mesh.vertices().end();
    
  for(; v_it != v_ite; ++v_it, ++vertex_counter)
  {
    const Point &p = mesh.point(*v_it);
    vertex_num_map[*v_it] = vertex_counter;
    
    file << "Point(" << vertex_counter << ") = {"
         << std::setprecision(25) << p[0] << ", "
         << std::setprecision(25) << p[1] << ", "
         << std::setprecision(25) << p[2] << ", "
         << "mesh_size" << "};\n";
  }
  
  file << "\n";
}

template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_edges(std::ofstream &file, Surface_mesh &mesh,
  const Vertex_map_long &vertex_num_map,
  Edge_map_bool &edge_visited_map, Edge_map_long &edge_num_map) const
{
  file << "// Edges.\n\n";
  long int edge_counter = 1;
  
  typename Surface_mesh::Face_range::iterator f_it = mesh.faces().begin(),
    f_ite = mesh.faces().end();
  for(; f_it != f_ite; ++f_it)
  {
    // Loop around the edges of the current face.
    CGAL::Halfedge_around_face_iterator<Surface_mesh> hf_it, hf_ite;
    boost::tie(hf_it, hf_ite) = CGAL::halfedges_around_face(mesh.halfedge(*f_it), mesh);
  
    for(; hf_it != hf_ite; ++hf_it)
    {
      if(!edge_visited_map[mesh.edge(*hf_it)])
      {
        edge_visited_map[mesh.edge(*hf_it)] = true;
        edge_num_map[mesh.edge(*hf_it)] = edge_counter;
      
        const Vertex_index vh1 = mesh.source(*hf_it);
        const Vertex_index vh2 = mesh.target(*hf_it);
        
        file << "Line(" << edge_counter << ") = {"
             << vertex_num_map[vh1] << ", " << vertex_num_map[vh2] << "};\n";
        ++edge_counter;
      }
    }
  }
  
  file << "\n";
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_faces(std::ofstream &file, Surface_mesh &mesh,
  const Vertex_map_long &vertex_num_map,
  Edge_map_long &edge_num_map, Face_map_long &face_num_map) const
{
  file << "// Faces.\n\n";
  
  long int face_counter = 1;
  typename Surface_mesh::Face_range::iterator f_it = mesh.faces().begin(),
    f_ite = mesh.faces().end();
    
  for(; f_it != f_ite; ++f_it)
  {
    // Check if the edges belonging to the current face
    // are correctly oriented for Gmsh.
    check_edge_orientations(mesh, *f_it, vertex_num_map, edge_num_map);
    
    // Loop around the edges of the current face.
    file << "Curve Loop(" << face_counter << ") = {";
    CGAL::Halfedge_around_face_iterator<Surface_mesh> hf_it, hf_ite;
    boost::tie(hf_it, hf_ite) = CGAL::halfedges_around_face(mesh.halfedge(*f_it), mesh);
    
    file << edge_num_map[mesh.edge(*hf_it)];
    ++hf_it;
    
    int face_num = 1;
    for(; hf_it != hf_ite; ++hf_it)
    {
      if((face_num % 10) == 0) // Avoid too long lines in the geo file.
      {
        file << ",\n";
        face_num = 1;
      }
      else
      {
        file << ", ";
        ++face_num;
      }
    
      file << edge_num_map[mesh.edge(*hf_it)];
    }
  
    file << "};\n";
    ++face_counter;
    
    // Write surface.
    face_num_map[*f_it] = face_counter;
    file << "Plane Surface(" << face_counter << ") = {"
         << (face_counter-1) << "};\n";
    ++face_counter;
  }
  
  file << "\n";
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
check_edge_orientations(Surface_mesh &mesh, const Face_index &f,
  const Vertex_map_long &vertex_num_map, Edge_map_long &edge_num_map) const
{
  // Check if the first two edges of the face are correctly oriented.
  CGAL::Halfedge_around_face_iterator<Surface_mesh> hf_it, hf_ite, hf_it2, hf_it3;
  boost::tie(hf_it, hf_ite) = CGAL::halfedges_around_face(mesh.halfedge(f), mesh);
  hf_it2 = hf_it;
  
  /** Looking for the common vertice between the first
      two segments.

      4 possibilities: {a, b} and {c, a}
                       {a, b} and {a, c}
                       {b, a} and {c, a}
                       {b, a} and {a, c} (the correct one)
   */
   {
     long int a, b, c, d;
     
     // Get one halfedges of the first and second edge.
     // Warning: doing mesh.halfedge(mesh.edge(halfedge, 0))
     // is important here (erase halfedge-orientation which is non-existent in geo files)
     Halfedge_index he1 = mesh.halfedge(mesh.edge(*hf_it), 0);
     ++hf_it2;
     Halfedge_index he2 = mesh.halfedge(mesh.edge(*hf_it2), 0);
   
     // Get the vertex numbers connected to the two halfedges.
     a = vertex_num_map[mesh.source(he1)];
     b = vertex_num_map[mesh.target(he1)];
     c = vertex_num_map[mesh.source(he2)];
     d = vertex_num_map[mesh.target(he2)];
     
     // First segment wrongly oriented,
     // second correctly oriented.
     if(a == c)
     {
       edge_num_map[mesh.edge(*hf_it2)] = std::abs(edge_num_map[mesh.edge(*hf_it2)]);
       edge_num_map[mesh.edge(*hf_it)] = std::abs(edge_num_map[mesh.edge(*hf_it)]) * (-1);
     }
     // The two segments are wrongly oriented.
     else if(a == d)
     {
       edge_num_map[mesh.edge(*hf_it2)] = std::abs(edge_num_map[mesh.edge(*hf_it2)]) * (-1);
       edge_num_map[mesh.edge(*hf_it)] = std::abs(edge_num_map[mesh.edge(*hf_it)]) * (-1);
     }
     // The first segment is correctly oriented,
     // the second is wrongly oriented.
     else if(b == d)
     {
       edge_num_map[mesh.edge(*hf_it2)] = std::abs(edge_num_map[mesh.edge(*hf_it2)]) * (-1);
       edge_num_map[mesh.edge(*hf_it)] = std::abs(edge_num_map[mesh.edge(*hf_it)]);
     }
     // Else: nothing to do: the segments are correctly oriented.
     else
     {
       edge_num_map[mesh.edge(*hf_it2)] = std::abs(edge_num_map[mesh.edge(*hf_it2)]);
       edge_num_map[mesh.edge(*hf_it)] = std::abs(edge_num_map[mesh.edge(*hf_it)]);
     }
   }

   /** Looking for the other segments.

      As the previous segment is correctly oriented (by above),
      it remains two possibilities:

        {b, a} and {c, a}
        {b, a} and {a, c} (the correct one)
  */
  long int a, b;
  hf_it3 = hf_it2;
  ++hf_it3; // Third edge.
  
  while(hf_it3 != hf_ite)
  {
    // Get halfedges of the preceeding edge and the current one.
    // Warning: doing mesh.halfedge(mesh.edge(halfedge, 0))
    // is important here (erase halfedge-orientation which is non-existent in geo files).
    Halfedge_index he1 = mesh.halfedge(mesh.edge(*hf_it2), 0);
    Halfedge_index he2 = mesh.halfedge(mesh.edge(*hf_it3), 0);
    int orientation = (edge_num_map[mesh.edge(*hf_it2)] >= 0 ? 1 : -1);
  
    // Get the last vertex number of the preceeding edge
    // and the first vertex number of the current edge.
    a = ( (orientation == 1) ? vertex_num_map[mesh.target(he1)]
                             : vertex_num_map[mesh.source(he1)] );
    b = vertex_num_map[mesh.source(he2)];
  
    if(a != b) // Wrong orientation.
      edge_num_map[mesh.edge(*hf_it3)] = std::abs(edge_num_map[mesh.edge(*hf_it3)]) * (-1);
    // Correct orientation.
    else
      edge_num_map[mesh.edge(*hf_it3)] = std::abs(edge_num_map[mesh.edge(*hf_it3)]);
    
    // Next edge.
    hf_it2 = hf_it3;
    ++hf_it3;
  }
  
  /* Example of general reordering:
       Given points number 2, 9, 8, 4.

       Line(4) =  {9, 2}; -> Line(-4) = {2, 9};
       Line(10) = {9, 8}; -> Line(10) = {9, 8};
       Line(11) = {8, 4}; -> Line(11) = {8, 4};
       Line(5)  = {2, 4}; -> Line(-5) = {4, 2};

       Exemple of reordering with the last element to be reordered:
       Given points numbers 1, 2, 3 and lines:

       Line(1) = {1, 2}; -> Line(1) = {1, 2};
       Line(2) = {2, 3}; -> Line(2) = {2, 3};
       Line(3) = {1, 3}; -> Line(-3) = {3, 1};

       Line(3) must be reordered by {3, 1} to match the
       orientation of Line(1) and Line(2).
   */
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_volume(std::ofstream &file, Surface_mesh &mesh, 
  const Face_map_long &face_num_map, Volume_map_long &volume_num_map) const
{
  long int volume_counter = 1;
  
  file << "// Volume.\n\n";
  file << "Surface Loop(" << volume_counter << ") = {";
  
  typename Surface_mesh::Face_range::iterator f_it = mesh.faces().begin(),
    f_ite = mesh.faces().end();
    
  if(f_it != f_ite)
  {
    file << face_num_map[*f_it];
    ++f_it;
  }
  
  int face_num = 1;
  for(; f_it != f_ite; ++f_it)
  {
    if((face_num % 10) == 0) // Avoid too long lines in the geo file.
    {
      file << ",\n";
      face_num = 1;
    }
    else
    {
      file << ", ";
      ++face_num;
    }
    
    file << face_num_map[*f_it];
  }
  
  file << "};\n";
  ++volume_counter;
  
  // Write volume.
  file << "Volume(" << volume_counter << ") = {"
       << (volume_counter-1) << "};\n";
       
  // Warning: implementation for only one volume.
  volume_num_map[Face_index(0)] = volume_counter;
  
  file << "\n";
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_physical_points(std::ofstream &file,
  const std::map<Vertex_index, std::vector<int> > &vertices_physical_groups,
  const Vertex_map_long &vertex_num_map) const
{
  file << "// Physical points.\n\n";
  
  // Get physical <-> point list correspondance.
  std::map<int, std::vector<long int> > physical_groups;
  for(const std::pair<Vertex_index, std::vector<int> > &it : vertices_physical_groups)
  {
    for(std::size_t i = 0; i < it.second.size(); ++i)
      physical_groups[it.second[i]].push_back( vertex_num_map[it.first] );
  }
  
  // Write the physical points to the geo file.
  for(const std::pair<int, std::vector<long int> > &it : physical_groups)
  {
    if(it.second.empty()) continue;
  
    std::size_t i = 0;
    file << "Physical Point(" << it.first << ") = {";
    
    file << it.second[i];
    ++i;
    int vertex_num = 0;
    
    for(; i < it.second.size(); ++i)
    {
      file << ", " << it.second[i];
      
      // Avoid too long lines in the geo file.
      if(vertex_num < 10) ++vertex_num;
      else if((i + 1) < it.second.size())
      {
        vertex_num = 0;
        file << "\n";
      }
    }
    
    file << "};\n";
  }
  
  file << "\n";
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_physical_lines(std::ofstream &file) const
{
  // Get physical groups of entity lines.
  gmsh::vectorpair physical_groups;
  gmsh::model::getPhysicalGroups(physical_groups, 1);

  if(!physical_groups.empty())
  {
    std::cerr << "Warning: physical lines not implement and will be ignored.\n";
    file << "// Physical lines: warning not implemented and are ignored.\n";
    file << "\n";
  }
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_physical_faces(std::ofstream &file, 
  const std::map<Face_index, std::vector<int> > &faces_physical_groups,
  const Face_map_long &face_num_map) const
{
  file << "// Physical faces.\n\n";

  // Get physical <-> face list correspondance.
  std::map<int, std::vector<long int> > physical_groups;
  for(const std::pair<Face_index, std::vector<int> > &it : faces_physical_groups)
  {
    for(std::size_t i = 0; i < it.second.size(); ++i)
      physical_groups[it.second[i]].push_back( face_num_map[it.first] );
  }
  
  // Write the physical faces to the geo file.
  for(const std::pair<int, std::vector<long int> > &it : physical_groups)
  {
    if(it.second.empty()) continue;
  
    std::size_t i = 0;
    file << "Physical Surface(" << it.first << ") = {";
    
    file << it.second[i];
    ++i;
    int face_num = 0;
    
    for(; i < it.second.size(); ++i)
    {
      file << ", " << it.second[i];
      
      // Avoid too long lines in the geo file.
      if(face_num < 10) ++face_num;
      else if((i + 1) < it.second.size())
      {
        face_num = 0;
        file << "\n";
      }
    }
    
    file << "};\n";
  }
  
  file << "\n";
}


// Warning: implementation for only one volume.
template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
write_physical_volumes(std::ofstream &file,
  const std::map<Volume_index, std::vector<int> > &volumes_physical_groups,
  const Volume_map_long &volume_num_map) const
{
  file << "// Physical volumes.\n\n";
  
  if(volumes_physical_groups.size() > 1)
  {
    std::cerr << "Warning: more than one volume. Consider only first one, ignore the others.\n";
  }
  
  // Get physical <-> volume list correspondance.
  std::map<int, std::vector<long int> > physical_groups;
  {
    typename std::map<Volume_index, std::vector<int> >::const_iterator 
      it = volumes_physical_groups.begin();
      
    for(std::size_t i = 0; i < (it->second.size()); ++i)
      physical_groups[it->second[i]].push_back(volume_num_map[Face_index(0)]);
  }

  // Write the physical volumes to the geo file.
  for(const std::pair<int, std::vector<long int> > &it : physical_groups)
  {
    if(it.second.empty()) continue;
    
    std::size_t i = 0;
    file << "Physical Volume(" << it.first << ") = {";
    
    file << it.second[i];
    ++i;
    int volume_num = 0;
  
    for(; i < it.second.size(); ++i)
    {
      file << ", " << it.second[i];
      
      // Avoid too long lines in the geo file.
      if(volume_num < 10) ++volume_num;
      else if((i + 1) < it.second.size())
      {
        volume_num = 0;
        file << "\n";
      }
    }
  
    file << "};\n";
  }

  file << "\n";
}