/// \file extract_surface_mesh_from_msh.hpp
/// \brief Convert the surface of a MSH mesh to a CGAL Surface_mesh.
/// \author Leblanc Christophe.
/// \date 23/10/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 <CGAL/Surface_mesh.h>

#include "gmsh.h"

/// \brief Extract a CGAL Surface Mesh from a msh file.
/// \warning: only implemented for triangles and quadrangles (order 1 and 2).
template<typename Surface_meshT>
class ExtractSurfaceMeshFromMSH
{
public:
  typedef Surface_meshT Surface_mesh;
  typedef std::map<std::size_t, typename Surface_mesh::Vertex_index> Node_map; // Maps node tag of MSH to node index of surface mesh.
  typedef std::map<std::size_t, typename Surface_mesh::Face_index> Face_map; // Maps face tags of MSH to face index of surface mesh.
  
  /// \brief Constructor.
  ExtractSurfaceMeshFromMSH() {}
  
  /// \brief Destructor.
  ~ExtractSurfaceMeshFromMSH() {}
  
  /// \brief Copy constructor.
  ExtractSurfaceMeshFromMSH(const ExtractSurfaceMeshFromMSH &rhs):
    m_mesh(rhs.m_mesh), m_filename(rhs.m_filename), m_element_type(rhs.m_element_type),
    m_vhandle(rhs.m_vhandle), m_face_map(rhs.m_face_map) {}
  
  /// \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 Get the element type to extract from the MSH.
  int get_element_type() const
  { return m_element_type; }
  
  /// \brief Set the element type to extract from the MSH.
  void set_element_type(const int e)
  { m_element_type = e; }
  
  /// \brief Get the map between node tag of MSH and node index of surface mesh.
  const Node_map& get_node_map() const
  { return m_vhandle; }
  
  /// \brief Get the map between face tag of MSH and face index of surface mesh.
  const Face_map& get_face_map() const
  { return m_face_map; }
  
  /// \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;

  Surface_mesh m_mesh;
  std::string m_filename;
  int m_element_type;
  Node_map m_vhandle; // Maps node tag of MSH to node index of surface mesh.
  Face_map m_face_map; // Maps face tag of MSH to face index of surface mesh.
  
  // 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.
};

template<typename Surface_meshT>
std::map<unsigned int, unsigned int> ExtractSurfaceMeshFromMSH<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 Surface_meshT>
ExtractSurfaceMeshFromMSH<Surface_meshT>& ExtractSurfaceMeshFromMSH<Surface_meshT>::
operator=(const ExtractSurfaceMeshFromMSH &rhs)
{
  if(this != rhs)
  {
    m_mesh = rhs.m_mesh;
    m_filename = rhs.m_filename;
    m_element_type = rhs.m_element_type;
    m_vhandle = rhs.m_vhandle;
    m_face_map = rhs.m_face_map;
  }
  
  return this;
}


template<typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<Surface_meshT>::
extract()
{
  // Check element type.
  if(m_element_type == null_element_type())
  {
    std::cerr << "Error: element type not set." << std::endl;
    return;
  }

  // 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);
  }
  
  // Get the nodes.
  std::vector<std::size_t> node_tags;
  std::vector<double> coordinates, parametric_coordinates;
  gmsh::model::mesh::getNodesByElementType(m_element_type, node_tags, coordinates, parametric_coordinates);
  
  // Get the elements.
  std::vector<std::size_t> element_tags;
  std::vector<std::size_t> element_node_tags; 
  gmsh::model::mesh::getElementsByType(m_element_type, element_tags, element_node_tags);
  
  std::cout << "done.\n";
  
  //
  // Convert to surface mesh.
  //
  {
    std::cout << "Convert to surface mesh...\n";
  
    const unsigned int nb_nodes_per_element = m_element_type_nb_nodes_map.at(m_element_type);
    m_mesh.clear();
    
    // Copy vertices.
    for(std::size_t i = 0; i < node_tags.size(); ++i)
    {
      m_vhandle[node_tags[i]] = m_mesh.add_vertex(Point(
        coordinates[3*i+0],
        coordinates[3*i+1],
        coordinates[3*i+2]
      ));
    }
  
    // Copy face connectivity.
    std::vector<typename Surface_mesh::Vertex_index> face_vhandles;
    
    for(std::size_t i = 0; i < element_tags.size(); ++i)
    {
      face_vhandles.clear();
      
      for(unsigned int j = 0; j < nb_nodes_per_element; ++j)
      {
        face_vhandles.push_back(m_vhandle[element_node_tags[i*nb_nodes_per_element + j]]);
      }
      
      m_face_map[element_tags[i]] = m_mesh.add_face(face_vhandles);
    }
    
    std::cout << "done.\n";
  }
}




#endif // EXTRACT_SURFACE_MESH_FROM_MSH_HPP