/// \file extract_surface_mesh_from_msh.hpp
/// \brief Convert the volume of a MSH mesh to a CGAL Surface_mesh.
/// \author Leblanc Christophe.
/// \date 18/11/2020
/// \mail cleblancad@gmail.com

#ifndef EXTRACT_SURFACE_MESH_FROM_MSH_HPP
#define EXTRACT_SURFACE_MESH_FROM_MSH_HPP

#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <map>
#include <set>
#include <array>
#include <ctime>

#include <CGAL/Surface_mesh.h>
#include <CGAL/box_intersection_d.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/repair.h>

#include "gmsh.h"

/// \brief Extract a CGAL Surface Mesh from a msh file.
/// \warning: only implemented for order 1 tetrahedra.
template<typename KernelT, typename Surface_meshT>
class ExtractSurfaceMeshFromMSH
{
public:
  typedef KernelT Kernel;
  typedef Surface_meshT Surface_mesh;
  
  /// \brief constructor.
  ExtractSurfaceMeshFromMSH(): m_paranoid(false) {}
  
  /// \brief Destructor.
  ~ExtractSurfaceMeshFromMSH() {}
  
  /// \brief Copy constructor.
  ExtractSurfaceMeshFromMSH(const ExtractSurfaceMeshFromMSH &rhs):
    m_mesh(rhs.m_mesh), m_filename(rhs.m_filename), m_paranoid(rhs.m_paranoid),
    m_geo_filename(rhs.m_geo_filename), m_mesh_size(rhs.m_mesh_size) {}
    
  /// \brief Copy operator.
  ExtractSurfaceMeshFromMSH& operator=(const ExtractSurfaceMeshFromMSH &rhs);
  
  /// \brief Return the extracted surface mesh.
  const Surface_mesh& get_surface_mesh() const
  { return m_mesh; }
  
  /// \brief Get the current file name.
  const std::string& get_mesh_file_name() const
  { return m_filename; }
  
  /// \brief Set the file name from which the MSH file is read.
  void set_mesh_file_name(const std::string &name)
  { m_filename = name; }
  
  /// \brief Set the paranoid mode on.
  void set_paranoid() { m_paranoid = true; }
  
  /// \brief Set the paranoid mode off.
  void unset_paranoid() { m_paranoid = false; }
  
  /// \brief Get the mode (paranoid or not).
  bool is_paranoid() const { return m_paranoid; }
  
  /// \brief Get the geo file name.
  const std::string& get_geo_file_name() const
  { return m_geo_filename; }
  
  /// \brief Set the geo file name.
  void set_geo_file_name(const std::string &name)
  { m_geo_filename = name; }
  
  /// \brief Get the mesh size use for the geo file.
  double get_mesh_size() const
  { return m_mesh_size; }
  
  /// \brief Set the mesh size for the geo file.
  void set_mesh_size(const double &s)
  { m_mesh_size = s; }
  
  /// \brief Extract surface mesh from MSH.
  void extract();
  
  /// \brief Get the null element type.
  static int null_element_type()
  { return std::numeric_limits<int>::min(); }
  
private:
  typedef typename Surface_mesh::Point Point;
  typedef typename Surface_mesh::Vertex_index Vertex_index;
  typedef typename Surface_mesh::Halfedge_index Halfedge_index;
  typedef typename Surface_mesh::Edge_index Edge_index;
  typedef typename Surface_mesh::Face_index Face_index;
  typedef          std::size_t              Volume_index;
  typedef typename Kernel::Triangle_3 Triangle;
  
  typedef typename Surface_mesh::template Property_map<Vertex_index, long int> Vertex_map_long;
  typedef typename Surface_mesh::template Property_map<Edge_index, bool> Edge_map_bool;
  typedef typename Surface_mesh::template Property_map<Edge_index, long int> Edge_map_long;
  typedef typename Surface_mesh::template Property_map<Face_index, long int> Face_map_long;
  typedef typename Surface_mesh::template Property_map<Face_index, long int> Volume_map_long;
  
  typedef std::pair<Triangle, std::array<std::size_t, 4> > Triangle_with_ids; // Triangle and its id followed by its vertex ids.
  typedef std::vector< Triangle_with_ids > Triangles_with_ids;
  typedef typename Triangles_with_ids::iterator Triangles_with_ids_iterator;
  typedef CGAL::Box_intersection_d::Box_with_handle_d<typename KernelT::FT, 3,
    Triangles_with_ids_iterator> Box;
  
  Surface_mesh m_mesh;
  std::string m_filename;
  bool m_paranoid;
  std::string m_geo_filename;
  double m_mesh_size;
  
  // Number of elements per element type.
  // https://gmsh.info/doc/texinfo/gmsh.html
  enum ElementType
  {
    TETRAHEDRON = 4,
    TETRAHEDRON10 = 11,
    HEXAHEDRON = 5,
    HEXAHEDRON20 = 17,
    HEXAHEDRON27 = 12,
    PRISM = 6,
    PRISM15 = 18,
    PRISM18 = 13,
    PYRAMID = 7,
    PYRAMID13 = 19,
    PYRAMID14 = 14,
    TRIANGLE = 2,
    TRIANGLE6 = 9,
    QUADRANGLE = 3,
    QUADRANGLE8 = 16
  };
  
  static std::map<unsigned int, unsigned int> m_element_type_nb_nodes_map; // Element type - number of nodes.
  
  /// \brief Write the geo file corresponding to the mesh.
  void 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;
  
  /// \brief Write a small header to the geo file.
  void write_header(std::ofstream &file) const;
  
  /// \brief Write vertices.
  void write_vertices(std::ofstream &file, Surface_mesh &mesh,
    Vertex_map_long &vertex_num_map) const;
  
  /// \brief Write edges.
  void 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;
    
  /// \brief Write faces.
  void 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;
  
  /// \brief Adapt edge orientation for faces in the geo file.
  void check_edge_orientations(Surface_mesh &mesh, const Face_index &f,
    const Vertex_map_long &vertex_num_map, Edge_map_long &edge_num_map) const;
    
  /// \brief Write the volume.
  void write_volume(std::ofstream &file, Surface_mesh &mesh, 
    const Face_map_long &face_num_map, Volume_map_long &volume_num_map) const;
    
  /// \brief Write physical points.
  void 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;
    
  /// \brief Write physical lines.
  /// \warning: not implemented.
  void write_physical_lines(std::ofstream &file) const;
  
  /// \brief Write physical faces.
  void 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;
    
  /// \brief Xrite physical volumes.
  void 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;
    
  /// \brief Extract physicals from the MSH file.
  void extract_physicals(std::map<std::size_t, std::vector<int> > &nodes_physical_groups,
    std::map<std::size_t, std::pair<std::vector<int>, std::vector<std::size_t> > > &faces_physical_groups,
    std::map<std::size_t, std::vector<int> > &volumes_physical_groups) const;

  /// \brief Extract physicals of dim 0 (a.k.a points) from the MSH file.
  void extract_physicals_dim0(std::map<std::size_t, std::vector<int> > &nodes_physical_groups) const;
  
  // Note: extraction of physicals of dim 1 not implemented.
  
  /// \brief Extract physicals of dim 2 (a.k.a faces) from the MSH file.
  void extract_physicals_dim2(std::map<std::size_t, std::pair<std::vector<int>, std::vector<std::size_t> > > &faces_physical_groups) const;
  
  /// \brief Extract physicals of dim 3 (a.k.a volumes) from the MSH file.
  void extract_physicals_dim3(std::map<std::size_t, std::vector<int> > 
    &volumes_physical_groups) const;
    
  /// \brief Try to repair a defective mesh (fill holes).
  void repair_mesh(Surface_mesh &mesh) const;
  
  class Report
  {
  public:
    std::vector<bool>* m_duplicated;
    bool m_paranoid; // Paranoid mode. Do not trust node connectivity.
  
    Report(std::vector<bool> &duplicated, bool paranoid = false): m_duplicated(&duplicated),
      m_paranoid(paranoid) {}
      
    void set_paranoid() { m_paranoid = true; }
    void unset_paranoid() { m_paranoid = false; }
    bool is_paranoid() const { return m_paranoid; }
      
    /// \brief Callback functor that computes closest triangles.
    void operator()(const Box &a, const Box &b)
    {
      // Note: the "box_self_intersection_d" function excludes self-intersections of boxes,
      // and records box intersections only once.
      const typename Triangle_with_ids::second_type &ids_a = a.handle()->second;
      const typename Triangle_with_ids::second_type &ids_b = b.handle()->second;
      
      if(!m_paranoid)
      {
        // Check if triangles have same last three ids.
        std::set<typename Triangle_with_ids::second_type::value_type> set_a, set_b;
        for(int i = 1; i < 4; ++i)
        {
          set_a.insert(ids_a[i]);
          set_b.insert(ids_b[i]);
        }
      
        if(set_a == set_b) // Yes, mark triangles as duplicated.
        {
          (*m_duplicated)[ids_a[0]] = (*m_duplicated)[ids_b[0]] = true;
        }
      }
      // Paranoid mode: do not relies on connectivity only, look at triangle's coordinates.
      else
      {
        constexpr typename Kernel::FT factor = 10.0;
        bool all_coordinates_equal = true;
        const typename Triangle_with_ids::first_type &tri_a = a.handle()->first;
        const typename Triangle_with_ids::first_type &tri_b = b.handle()->first;
        
        for(int i = 0; i < 3; ++i)
        {
          const Point &pa = tri_a[i];
          bool found_different_coordinates = true;
          
          for(int j = 0; j < 3; ++j)
          {
            const Point &pb = tri_b[j];
            
            if(nearly_equal(pa[0], pb[0], factor) &&
               nearly_equal(pa[1], pb[1], factor) &&
               nearly_equal(pa[2], pb[2], factor))
            {
              found_different_coordinates = false;
              break;
            }
          }
          
          if(found_different_coordinates)
          {
            all_coordinates_equal = false;
            break;
          }
        }
        
        if(all_coordinates_equal)
        {
          (*m_duplicated)[ids_a[0]] = (*m_duplicated)[ids_b[0]] = true;
        }
      }
    }
    
  private:
    // The machine epsilon has to be scaled to the magnitude of the values used
    // and multiplied by the desired precision in ULPs (units in the last place)
    // Source: https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon
    /// \brief Return true if x and y are "nearly" equal, false otherwise.
    template<class T>
    inline bool nearly_equal(const T &x, const T &y, const T factor = 1.0, 
      const int ulp = 2) const
    {
      return (std::abs(x-y) <= std::numeric_limits<T>::epsilon() * factor * std::abs(x+y) * ulp)
        // unless the result is subnormal
             || (std::abs(x-y) < std::numeric_limits<T>::min());
    }
  };
};

template<typename KernelT, typename Surface_meshT>
std::map<unsigned int, unsigned int> ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
  m_element_type_nb_nodes_map = {
  {TETRAHEDRON, 4},
  {TETRAHEDRON10, 10},
  {HEXAHEDRON, 8},
  {HEXAHEDRON20, 20},
  {HEXAHEDRON27, 27},
  {PRISM, 6},
  {PRISM15, 15},
  {PRISM18, 18},
  {PYRAMID, 5},
  {PYRAMID13, 13},
  {PYRAMID14, 14},
  {TRIANGLE, 3},
  {TRIANGLE6, 6},
  {QUADRANGLE, 4},
  {QUADRANGLE8, 8}
};


template<typename KernelT, typename Surface_meshT>
ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>& 
ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
operator=(const ExtractSurfaceMeshFromMSH &rhs)
{
  if(this != &rhs)
  {
    m_mesh = rhs.m_mesh;
    m_filename = rhs.m_filename;
    m_paranoid = rhs.m_paranoid;
    m_geo_filename = rhs.m_geo_filename;
    m_mesh_size = rhs.m_mesh_size;
  }
  
  return this;
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
extract()
{
  // Get the file stream on the file.
  if(m_filename.empty())
  {
    std::cerr << "Error: no file name set." << std::endl;
    return;
  }
  
  // Check if exists.
  std::ifstream file;
  file.open(m_filename.c_str(), std::ifstream::in);

  if(!file.is_open())
  {
    std::cerr << "Error: unable to open file: \"" << m_filename << "\"" << std::endl;
    return;
  }
  
  file.close();
  
  //
  // Get the mesh informations.
  //
  std::cout << "Read msh file...\n";
  // Open mesh file only if needed.
  {
    std::vector<std::string> models_open;
    gmsh::model::list(models_open);

    if(models_open.empty() || models_open[0].empty())
      gmsh::open(m_filename);
      
    gmsh::model::mesh::removeDuplicateNodes();
  }
  
  // Get the nodes.
  std::vector<std::size_t> node_tags;
  std::vector<double> coordinates, parametric_coordinates; // Coordinates[i] corresponds to node_tags[i].
  const int element_type = gmsh::model::mesh::getElementType("tetrahedron", 1);
  
  gmsh::model::mesh::getNodesByElementType(element_type, node_tags, coordinates,
    parametric_coordinates);
  
  // Get the elements (order 1 tetrahedra).
  /*  http://gmsh.info/doc/texinfo/gmsh.html: 
                 v
                .
              ,/
             /
           2                                     2
         ,/|`\                                 ,/|`\
       ,/  |  `\                             ,/  |  `\
     ,/    '.   `\                         ,6    '.   `5
   ,/       |     `\                     ,/       8     `\
 ,/         |       `\                 ,/         |       `\
0-----------'.--------1 --> u         0--------4--'.--------1
 `\.         |      ,/                 `\.         |      ,/
    `\.      |    ,/                      `\.      |    ,9
       `\.   '. ,/                           `7.   '. ,/
          `\. |/                                `\. |/
             `3                                    `3
                `\.
                   ` w
  */
  // Get element face nodes.
  std::vector<std::size_t> node_tags_on_faces; 
      // `node_tags_on_faces' contains the node tags of the faces
      // for all elements: [e1f1n1, ..., e1f1nFaceType, e1f2n1, ...].
  const int face_type = 3; // Triangular faces.
  
  gmsh::model::mesh::getElementFaceNodes(element_type, face_type, node_tags_on_faces);
  
  std::cout << "done.\n";
  
  //
  // Detect non-duplicated faces (surface faces).
  //
  std::vector<bool> duplicated_triangles(node_tags_on_faces.size() / 3, false);
  std::map<std::size_t, std::size_t> vhandle_tri;
  {
    std::cout << "Detect surface faces.\n";
    
    constexpr std::size_t number_of_vertices_per_face = 3;
    const std::size_t number_of_faces = node_tags_on_faces.size() /
      (number_of_vertices_per_face);
    
    Triangles_with_ids triangles_with_ids;
    triangles_with_ids.reserve(number_of_faces);
    
    for(std::size_t i = 0; i < node_tags.size(); ++i)
    {
      vhandle_tri[node_tags[i]] = i;
    }
    
    for(std::size_t i = 0; i < number_of_faces; ++i)
    {
      const std::size_t tags[] = {
        vhandle_tri[node_tags_on_faces[3*i+0]],
        vhandle_tri[node_tags_on_faces[3*i+1]],
        vhandle_tri[node_tags_on_faces[3*i+2]]
      };
    
      const Triangle triangle(
        Point(coordinates[3*tags[0]+0],
              coordinates[3*tags[0]+1],
              coordinates[3*tags[0]+2]
             ),
        Point(coordinates[3*tags[1]+0],
              coordinates[3*tags[1]+1],
              coordinates[3*tags[1]+2]
             ),
        Point(coordinates[3*tags[2]+0],
              coordinates[3*tags[2]+1],
              coordinates[3*tags[2]+2]
             )
      );
      
      const Triangle_with_ids triangle_with_ids = std::make_pair(
        triangle, std::array<std::size_t, 4>{i, node_tags_on_faces[3*i+0],    
          node_tags_on_faces[3*i+1], node_tags_on_faces[3*i+2]});
      
      triangles_with_ids.push_back(triangle_with_ids);
    }
    
    std::vector<Box> boxes;
    boxes.reserve(triangles_with_ids.size());
    {
      Triangles_with_ids_iterator it = triangles_with_ids.begin(),
        ite = triangles_with_ids.end();
        
      for(; it != ite; ++it)
      {
        CGAL::Bbox_3 bbox = it->first.bbox();
        
        // Compute bounding-box diagonal.
        typename Kernel::FT diagonal = 0.0;
        for(int l = 0; l < 3; ++l)
        {
          diagonal += (bbox.max(l) - bbox.min(l))
                    * (bbox.max(l) - bbox.min(l));
        }
        
        diagonal = sqrt(diagonal);
        
        if(diagonal <= std::numeric_limits<typename Kernel::FT>::epsilon())
          throw std::runtime_error("Error: msh2geo: found degenerated triangle.");
        
        if(std::abs(bbox.max(0) - bbox.min(0)) <= 
           10.0*std::numeric_limits<typename Kernel::FT>::epsilon())
        {
          bbox = CGAL::Bbox_3(bbox.min(0)-diagonal/100.0, bbox.min(1), bbox.min(2), 
                              bbox.max(0)+diagonal/100.0, bbox.max(1), bbox.max(2));
        }
      
        if(std::abs(bbox.max(1) - bbox.min(1)) <= 
           10.0*std::numeric_limits<typename Kernel::FT>::epsilon())
        {
          bbox = CGAL::Bbox_3(bbox.min(0), bbox.min(1)-diagonal/100.0, bbox.min(2), 
                              bbox.max(0), bbox.max(1)+diagonal/100.0, bbox.max(2));
        }
        
        if(std::abs(bbox.max(2) - bbox.min(2)) <= 
           10.0*std::numeric_limits<typename Kernel::FT>::epsilon())
        {
          bbox = CGAL::Bbox_3(bbox.min(0), bbox.min(1), bbox.min(2)-diagonal/100.0, 
                              bbox.max(0), bbox.max(1), bbox.max(2)+diagonal/100.0);
        }
        
        boxes.push_back( Box(bbox, it) );
      }
    }
    
    // Run the self intersection algorithm for boxes.
    Report report(duplicated_triangles);

    if(m_paranoid)
      report.set_paranoid();
      
    CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), report);
    
    std::size_t non_duplicated_count = 0;
    for(const bool b : duplicated_triangles)
      if(!b) ++non_duplicated_count;
    
    std::cout << "Found " << non_duplicated_count << " non duplicated triangles\n";
    std::cout << "done.\n";
  }
  
  //
  // Convert to surface mesh.
  //
  std::map<std::size_t, Vertex_index> node_map, node_map_base;
  std::map<Vertex_index, std::size_t> inverse_node_map, inverse_node_map_base;
  std::map<std::size_t, Face_index> face_map;
  std::map<Face_index, std::size_t> inverse_face_map;
  {
    std::cout << "Convert to surface mesh.\n";
    
    // Copy vertices.
    for(std::size_t i = 0; i < duplicated_triangles.size(); ++i)
    {
      if(!duplicated_triangles[i])
      {
        const std::size_t tags_base[] = {
          node_tags_on_faces[3*i+0],
          node_tags_on_faces[3*i+1],
          node_tags_on_faces[3*i+2]
        };
        
        const std::size_t tags[] = {
          vhandle_tri[tags_base[0]],
          vhandle_tri[tags_base[1]],
          vhandle_tri[tags_base[2]]
        };
        
        const Point ps[] = {
          Point(
            coordinates[3*tags[0]+0],
            coordinates[3*tags[0]+1],
            coordinates[3*tags[0]+2]
          ),
          Point(
            coordinates[3*tags[1]+0],
            coordinates[3*tags[1]+1],
            coordinates[3*tags[1]+2]
          ),
          Point(
            coordinates[3*tags[2]+0],
            coordinates[3*tags[2]+1],
            coordinates[3*tags[2]+2]
          )
        };
        
        // Add non-already added points.
        for(int j = 0; j < 3; ++j)
        {
          if(node_map.find(tags[j]) == node_map.end())
          {
            const Vertex_index vi = m_mesh.add_vertex(ps[j]);
            node_map[tags[j]] = vi;
            node_map_base[tags_base[j]] = vi;
            inverse_node_map[vi] = tags[j];
            inverse_node_map_base[vi] = tags_base[j];
          }
        }
        
      } // Triangle non-duplicated.
    } // Loop on triangles.
    
    // Copy face connectivity.
    for(std::size_t i = 0; i < duplicated_triangles.size(); ++i)
    {
      if(!duplicated_triangles[i])
      {
        const std::size_t tags[] = {
          vhandle_tri[node_tags_on_faces[3*i+0]],
          vhandle_tri[node_tags_on_faces[3*i+1]],
          vhandle_tri[node_tags_on_faces[3*i+2]]
        };
        
        std::vector<Vertex_index> face_vhandles = {
          node_map[tags[0]],
          node_map[tags[1]],
          node_map[tags[2]]
        };
      
        const Face_index face_index = m_mesh.add_face(face_vhandles);
        
        if(face_index == m_mesh.null_face())
        {
          std::cout << "Warning: failed to add face with node tags ["
                    << node_tags_on_faces[3*i+0] << ", " 
                    << node_tags_on_faces[3*i+1] << ", " 
                    << node_tags_on_faces[3*i+2]
                    << "]\n";
        }
        else
        {
          face_map[i] = face_index;
          inverse_face_map[face_index] = i;
        }
      } // Triangle non-duplicated.
    } // Loop on triangles.
    
    bool repair_tried = false;
    if(!CGAL::is_closed(m_mesh))
    {
      std::ofstream file_extracted("mesh_problems.off");
      file_extracted << m_mesh;
    
      // Try to repair mesh.
      std::cerr << "Mesh is not closed. Try to repair it. However, some physicals may be lost in the process...\n";
      repair_mesh(m_mesh);
      repair_tried = true;
    }
    
    if(!CGAL::is_closed(m_mesh))
    {
      std::cerr << "Error: repair failed. Check if the is any cavity in your mesh. Generated mesh may be invalid.\n";
    }
    else
    {
      // Orient mesh.
      if(repair_tried)
        std::cerr << "Repair successful!\n";
      
      CGAL::Polygon_mesh_processing::orient_to_bound_a_volume(m_mesh);
    }
    std::cout << "done.\n";
  }

  //
  // Extract physicals
  //
  std::map<Vertex_index, std::vector<int> > vertices_physical_groups;
  std::map<Face_index, std::vector<int> > faces_physical_groups;
  std::map<Volume_index, std::vector<int> > unique_volumes_physical_groups;
  {
    std::cout << "Extract physicals.\n";
  
    // Extract.
    std::map<std::size_t, std::vector<int> > nodes_physical_groups;
    std::map<std::size_t, std::pair<std::vector<int>, std::vector<std::size_t> > > 
      surfaces_physical_groups; // Face number, physical groups ids it belongs to (first) and corresponding node ids (second).
    std::map<std::size_t, std::vector<int> > volumes_physical_groups;
    
    extract_physicals(nodes_physical_groups, surfaces_physical_groups,
      volumes_physical_groups);
    
    // Map node tags to vertex indexes.
    //std::cout << "Vertices physical groups:\n";
    for(const std::pair<std::size_t, std::vector<int> > &it : nodes_physical_groups)
    {
      vertices_physical_groups[node_map_base.at(it.first)] = it.second;
      
      // Code for debugging purposes.
     /* std::cout << "Vertex " << node_map_base[it.first] << " has physical groups: {";
      for(std::size_t i = 0; i < it.second.size(); ++i)
        std::cout << it.second[i] << ", ";
      std::cout << "}\n";*/
    }
    
    // Map surface tags to face indexes.
    {
      // Construct node tags to surface tag map.
      std::map<std::set<Vertex_index>, std::size_t> nodes_to_face_map;
      for(const std::pair<std::size_t, std::pair<std::vector<int>, std::vector<std::size_t> > > 
          &it : surfaces_physical_groups)
      {
        std::set<Vertex_index> node_set;
        for(std::size_t i = 0; i < it.second.second.size(); ++i)
        {
          try
          {
            node_set.insert(node_map_base.at(it.second.second[i]));
          }
          catch(...)
          {
            std::cerr << "Warning: failed to add face " << it.first << " to physicals.\n";
          }
        }
        
        if(node_set.size() >= 3) // No face can have least than 3 nodes.
          nodes_to_face_map[node_set] = it.first;
      }
    
      // Retrieve the face indexes via their vertex indexes.
      typename Surface_mesh::Face_range::iterator f_it = m_mesh.faces().begin(),
        f_ite = m_mesh.faces().end();
      for(; f_it != f_ite; ++f_it)
      {
        // Get set of vertex indexes of the current face.
        std::set<Vertex_index> set;
        CGAL::Vertex_around_face_iterator<Surface_mesh> vf_it, vf_ite;
        boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(m_mesh.halfedge(*f_it), m_mesh);
        
        for(; vf_it != vf_ite; ++vf_it)
          set.insert(*vf_it);
        
        // Find the corresponding face tag and its associated physical groups.
        typename std::map<std::set<Vertex_index>, std::size_t>::const_iterator
          nodes_to_face_map_it = nodes_to_face_map.find(set);
          
        if(nodes_to_face_map_it != nodes_to_face_map.end())
        {
          // Found a corresponding face. Get the associated physical groups.
          const std::vector<int> &physical_groups = surfaces_physical_groups.at(
            nodes_to_face_map_it->second).first;
            
          faces_physical_groups[*f_it] = physical_groups;
          
          // Code for debugging purposes.
          /*std::cout << "face: " << (nodes_to_face_map_it->second) << " has physical groups: [";
          for(auto tt : physical_groups)
            std::cout << tt << ", ";
          std::cout << "]\n";*/
        } 
      }
    }
    
    // Map volume tags to volume indexes.
    // Warning: only one volume supported in this implementation.
    {
      if(volumes_physical_groups.size() > 1)
      {
        std::cerr << "More than one volume in the MSH file. This is not supported for the moment\n";
      }
      else
        unique_volumes_physical_groups = volumes_physical_groups;
    }
    
    std::cout << "done.\n";
  }
  
  //
  // Write geo file.
  //
  std::cout << "Export to geo file.\n";
  write_geo(m_mesh, vertices_physical_groups, faces_physical_groups,
    unique_volumes_physical_groups);
  std::cout << "done.\n";
}


// Code inspired from: https://doc.cgal.org/latest/Polygon_mesh_processing/Polygon_mesh_processing_2hole_filling_example_SM_8cpp-example.html#a12
template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
repair_mesh(Surface_mesh &mesh) const
{
  // Look for non-manifold vertices.
  std::vector<Vertex_index> non_manifold_vertices;
  {
    typename Surface_mesh::Vertex_range::const_iterator
      v_it = mesh.vertices().begin(), v_ite = mesh.vertices().end();
      
    for(; v_it != v_ite; ++v_it)
    {
      if(CGAL::Polygon_mesh_processing::is_non_manifold_vertex(*v_it, mesh))
      {
        non_manifold_vertices.push_back(*v_it);
      }
    }
  }
  
  // Duplicate non-manifold vertices.
  if(!non_manifold_vertices.empty())
  {
    std::cout << "Found " << non_manifold_vertices.size() << " non manifold " 
              << (non_manifold_vertices.size() == 1 ? "vertex" : "vertices") << ".\n";
    CGAL::Polygon_mesh_processing::duplicate_non_manifold_vertices(mesh);
  }

  // Collect one halfedge per boundary cycle.
  std::vector<Halfedge_index> border_cycles;
  CGAL::Polygon_mesh_processing::extract_boundary_cycles(mesh,
    std::back_inserter(border_cycles));

  for(Halfedge_index h : border_cycles)
  {
    std::vector<Face_index> patch_facets;
  
    CGAL::Polygon_mesh_processing::triangulate_hole(mesh,
      h,
      std::back_inserter(patch_facets),
      CGAL::Polygon_mesh_processing::parameters::
        use_delaunay_triangulation(true)
     );
  }
}

#include "write_geo.hpp"
#include "extract_physicals.hpp"

#endif // EXTRACT_SURFACE_MESH_FROM_MSH_HPP