/// \file extract_deformed_mesh.hpp
/// \brief Given a CGAL Surface_mesh converted from the surface of a MSH file (see "extract_surface_mesh_form_msh.hpp") and a deformation MSH file, compute the corresponding deformed CGAL Surface_mesh.
/// \author Leblanc Christophe.
/// \date 23/10/2020
/// \mail cleblancad@gmail.com

#ifndef EXTRACT_DEFORMED_MESH_HPP
#define EXTRACT_DEFORMED_MESH_HPP

#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <map>
#include <set>
#include <array>
#include <CGAL/Surface_mesh.h>

#include "gmsh.h"

/// \brief Compute the deformed CGAL surface mesh from a gmsh displacement file.
/// Note: the surface mesh must have been extract using the class "ExtractSurfaceMeshFromMSH"
///   with the correct corresponding msh file to the displacement file.
template<typename Surface_meshT, typename Node_mapT>
class ExtractDeformedMesh
{
public:
  typedef Surface_meshT Surface_mesh;
  typedef Node_mapT Node_map;
  
  /// \brief Default constructor.
  ExtractDeformedMesh(): m_mesh(nullptr), m_node_map(nullptr),
    m_view_step(null_view_step()), m_element_type(null_element_type()) {}
  
  /// \brief Copy constructor.
  ExtractDeformedMesh(const ExtractDeformedMesh &rhs);
  
  /// \brief Destructor.
  ~ExtractDeformedMesh() {}
  
  /// \brief Copy operator.
  ExtractDeformedMesh& operator=(const ExtractDeformedMesh &rhs);
 
  /// \brief Get the file name for the displacement mesh.
  const std::string& get_displacement_file_name() const
  { return m_displacement_file_name; }
  
  /// \brief Set the file name for the displacement mesh.
  void set_displacement_file_name(const std::string &name)
  { m_displacement_file_name = name; }
  
  /// \brief Get the deformed mesh.
  const Surface_mesh& get_surface_mesh() const
  { return *m_mesh; }

  /// \brief Set the mesh to deform.
  void set_surface_mesh(Surface_mesh &mesh)
  { m_mesh = &mesh; }
  
  /// \brief Get the node map giving the correspondance between the
  ///  CGAL surface mesh and the GMSH MSH mesh. (See class "ExtractSurfaceMeshFromMSH").
  const Node_map& get_node_map() const
  { return *m_node_map; }

  /// \brief Set the node map giving the correspondance between the
  ///  CGAL surface mesh and the GMSH MSH mesh. (See class "ExtractSurfaceMeshFromMSH").  
  void set_node_map(const Node_map &map)
  { m_node_map = &map; }
  
  /// \brief Get the mesh filename.
  const std::string& get_mesh_file_name() const
  { return m_mesh_file_name; }
  
  /// \brief Set the mesh filename.
  void set_mesh_file_name(const std::string name)
  { m_mesh_file_name = name; }
  
  /// \brief Get the view step.
  int get_view_step() const
  { return m_view_step; }
  
  /// \brief Set the view step.
  void set_view_step(const int s)
  { m_view_step = s; }
  
  /// \brief Get the element type.
  int get_element_type() const
  { return m_element_type; }
  
  /// \brief Set the element type.
  void set_element_type(const int e)
  { m_element_type = e; }
 
  /// \brief Compute (extract) the deformed mesh.
  void extract();
  
  /// \brief Get the null view step.
  static int null_view_step()
  { return -1; }
  
  /// \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_displacement_file_name;
  std::string m_mesh_file_name;
  const Node_map* m_node_map;
  int m_view_step;
  int m_element_type;
  
  // 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, typename Node_mapT>
std::map<unsigned int, unsigned int> ExtractDeformedMesh<Surface_meshT, Node_mapT>::
  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, typename Node_mapT>
ExtractDeformedMesh<Surface_meshT, Node_mapT>::
ExtractDeformedMesh(const ExtractDeformedMesh &rhs)
{
  *m_mesh = *rhs.m_mesh; // Full copy as the mesh can be modified.
  m_displacement_file_name = rhs.m_displacement_file_name;
  m_mesh_file_name = rhs.m_mesh_file_name;
  m_node_map = rhs.m_node_map; // Shallow copy as the node map is not modified.
  m_view_step = rhs.m_view_step;
  m_element_type = rhs.m_element_type;
}


template<typename Surface_meshT, typename Node_mapT>
ExtractDeformedMesh<Surface_meshT, Node_mapT>&
ExtractDeformedMesh<Surface_meshT, Node_mapT>::
operator=(const ExtractDeformedMesh &rhs)
{
  if(this != &rhs)
  {
    *m_mesh = *rhs.m_mesh; // Full copy as the mesh can be modified.
    m_displacement_file_name = rhs.m_displacement_file_name;
    m_mesh_file_name = rhs.m_mesh_file_name;
    m_node_map = rhs.m_node_map; // Shallow copy as the node map is not modified.
    m_view_step = rhs.m_view_step;
    m_element_type = rhs.m_element_type;
  }
  
  return *this;
}


template<typename Surface_meshT, typename Node_mapT>
void ExtractDeformedMesh<Surface_meshT, Node_mapT>::
extract()
{
  // Some checks...
  if(m_view_step == null_view_step() || m_view_step < 0)
  {
    std::cerr << "Error: invalid view step." << std::endl;
    return;
  }
  
  if(m_element_type == null_element_type())
  {
    std::cerr << "Error: invalid element type." << std::endl;
    return;
  }
  
  if(m_mesh == nullptr)
  {
    std::cerr << "Error: mesh to deform not set." << std::endl;
    return;
  }

  if(m_node_map == nullptr)
  {
    std::cerr << "Error: node map not set." << std::endl;
    return;
  }
  
  // Get the file stream on the file.
  if(m_displacement_file_name.empty())
  {
    std::cerr << "Error: no displacement file name set." << std::endl;
    return;
  }
  
  if(m_mesh_file_name.empty())
  {
    std::cerr << "Error: no mesh file name set." << std::endl;
    return;
  }
  
  // Check if exists.
  std::ifstream file;
  file.open(m_displacement_file_name.c_str(), std::ifstream::in);

  if(!file.is_open())
  {
    std::cerr << "Error: unable to open file: \"" << m_displacement_file_name << "\"" << std::endl;
    return;
  }
  
  file.close();

  // Read the deformed mesh.
  std::cout << "Read displacement 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_mesh_file_name);
  }
  
  gmsh::merge(m_displacement_file_name);
  
  // Get the informations on the displacement.
  std::vector<int> view_tags;
  std::string view_data_type;
  double view_time;
  int view_num_components;
  std::vector<std::size_t> view_element_tags; // Tags of, e.g. tetrahedra.
  std::vector< std::vector<double> > view_node_data; // Contains displacements.
  
  gmsh::view::getTags(view_tags);
  gmsh::view::getModelData(view_tags[0], m_view_step, view_data_type,
    view_element_tags, view_node_data, view_time, view_num_components);
  
  // Get node tags corresponding to each face of the elements.
  std::vector<std::size_t> node_tags;
  std::vector<std::size_t> element_tags;
  gmsh::model::mesh::getElementsByType(m_element_type, element_tags, node_tags);
  
  std::cout << "done.\n";
  
  // Apply displacement on the surface mesh.
  // view_element_tags = [e1, ..., eM] 
  // view_node_data[e1]: displacements of nodes [n1, ..., nN] (4 nodes for order 1 tetrahedra).
  // node_tags = [e1n1, e1n2, ..., e1nN, e2n1, ...] (e: element, n: node).
  std::cout << "Apply displacement...\n";
  const unsigned int nb_nodes_per_element = m_element_type_nb_nodes_map.at(m_element_type);
  std::set<std::size_t> node_moved; // Set of moved nodes (displacement applied).
  
  for(std::size_t e = 0; e < view_element_tags.size(); ++e)
  {
    // Node displacements for the current element.
    const std::vector<double> &element_nodes_displacements = view_node_data[e];
  
    for(std::size_t n = 0; n < nb_nodes_per_element; ++n)
    {
      // Get the tag of node n for element e.
      const std::size_t current_node_tag = node_tags[e*nb_nodes_per_element + n];
      
      // Check if the node has been already moved.
      if(node_moved.find(current_node_tag) != node_moved.end())
        continue;
      else
        node_moved.insert(current_node_tag);
        
      // Check if the current node tag corresponds to a surface node.
      if(m_node_map->find(current_node_tag) == m_node_map->end())
        continue;
    
      // Get the displacement of node n for element e.
      std::vector<double> current_node_displacement(view_num_components);
      
      for(int i = 0; i < view_num_components; ++i)
        current_node_displacement[i] = element_nodes_displacements[n*view_num_components + i];
    
      // Apply the displacement on the (not already moved) node.
      Point &p = m_mesh->point( m_node_map->at(current_node_tag) );
      const Point p_moved(p[0] + current_node_displacement[0],
                          p[1] + current_node_displacement[1],
                          p[2] + current_node_displacement[2]);
                          
      p = p_moved;
    } // Loop on nodes of the current element.
  
  } // Loop on elements.
  
  std::cout << "done.\n";
}

#endif // EXTRACT_DEFORMED_MESH_HPP