/// \file contact.cpp
/// \brief Given a MSH file of Gmsh and a displacement file, look for patches of faces close to each other (one patch contains a set of faces close to each other). Then write a new MSH file adding these patches as physical surfaces.
///
/// Usage: ./contact <mesh file name> <displacement file name>
///
/// Requires: libgmsh-dev ligcgal-dev (ver. >= 5.0)
/// \author Leblanc Christophe.
/// \date 25/10/2020
/// \mail cleblancad@gmail.com

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include "gmsh.h"
#include "extract_surface_mesh_from_msh.hpp"
#include "extract_deformed_mesh.hpp"
#include "extract_close_faces.hpp"
#include "extract_physical_patches.hpp"

//
// Typedefs.
//
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef typename Kernel::Point_3  Point;
typedef CGAL::Surface_mesh<Point> Surface_mesh;
typedef ExtractSurfaceMeshFromMSH<Surface_mesh> MSHExtractor;
typedef typename MSHExtractor::Node_map Node_map;
typedef typename MSHExtractor::Face_map Face_map;
typedef ExtractDeformedMesh<Surface_mesh, Node_map> DeformExtractor;
typedef ExtractCloseFaces<Surface_mesh, Face_map> CloseFacesExtractor;
typedef ExtractPhysicalPatches<Surface_mesh, Face_map, 
  typename CloseFacesExtractor::Triangle_triangle_distance> PhysicalExtractor;

int main(int argc, char** argv)
{
  gmsh::initialize(argc, argv);

  //
  // Extract surface mesh from MSH file (order 1 tetrahedron).
  //
  if(argc < 3)
  {
    std::cout << "Usage: \"./contact <mesh file name> <displacement file name> <patch expand value (>= 0)> (optional)\"\n";
    return EXIT_FAILURE;
  }
  
  unsigned int patch_expand_value = 0;
  if(argc >= 4)
  {
    patch_expand_value = strtoul(argv[3], NULL, 0);
  }
  
  const std::string mesh_filename = argv[1];
  const std::string displacement_filename = argv[2];
  
  // Determine view step from displacement file name.
  int view_step = 0;
  {
    char str1[256], str2[256];
    const std::size_t pos = displacement_filename.find_last_of("/");
    const std::string filename = displacement_filename.substr(pos+1);
    
    if(sscanf(filename.c_str(), "%[a-z_]%d%s", str1, &view_step, str2) < 1)
    {
      std::cerr << "Error: could not find step in displacement step file." << std::endl;
      return EXIT_FAILURE;
    }
  }
  
  const int face_element_type = gmsh::model::mesh::getElementType("triangle", 1);
  const int element_type = gmsh::model::mesh::getElementType("tetrahedron", 1);
  
  MSHExtractor extractor;
  
  extractor.set_element_type(face_element_type); // warning: only implemented for triangles and quadrangles (order 1 and 2).
  extractor.set_mesh_file_name(mesh_filename);
  extractor.extract();
  
  // Get the surface mesh and the node map with the MSH file.
  Surface_mesh mesh = extractor.get_surface_mesh();
  Surface_mesh initial_mesh = mesh;
  
  if(!mesh.is_valid())
  {
    std::cerr << "Invalid mesh found, please check your mesh..." << std::endl;
    return EXIT_FAILURE;
  }
  
  if(!CGAL::is_closed(mesh))
  {
    std::cerr << "Warning: surface mesh is not closed. This may lead to unexpected behavior..." << std::endl;
  }
  
  std::ofstream file_extracted("extracted_mesh.off");
  file_extracted << mesh;
  
  //
  // Compute the deformed mesh.
  //
  DeformExtractor deform_extractor;
  
  deform_extractor.set_mesh_file_name(extractor.get_mesh_file_name());
  deform_extractor.set_displacement_file_name(displacement_filename);
  deform_extractor.set_surface_mesh(mesh);
  deform_extractor.set_node_map(extractor.get_node_map());
  deform_extractor.set_element_type(element_type);
  deform_extractor.set_view_step(view_step);
  deform_extractor.extract();
  
  //
  // Check which triangular faces are close to each other.
  //
  CloseFacesExtractor close_extractor;
  typename CloseFacesExtractor::Mapped_triangle_triangle_distance
    mapped_tri_tri_distances;
  
  close_extractor.set_surface_mesh(mesh);
  close_extractor.set_face_map(extractor.get_face_map());
  close_extractor.set_threshold_distance(0.1); // 10% of mean edge length of considered face.
  close_extractor.extract();
  
  mapped_tri_tri_distances = close_extractor.get_face_mapped_triangle_triangle_distances();
  std::cout << "Number of close faces found: " << mapped_tri_tri_distances.size() << "\n";
  
  //
  // Modify a copy of the MSH file so to add patches of physical faces in contact.
  //
  PhysicalExtractor physical_extractor;
  
  physical_extractor.set_mesh_file_name(deform_extractor.get_mesh_file_name());
  physical_extractor.set_triangle_triangle_distance(
    close_extractor.get_triangle_triangle_distances());
  physical_extractor.set_surface_mesh(initial_mesh); // ! Undeformed mesh for determining the boundaries.
  physical_extractor.set_element_type(extractor.get_element_type());
  physical_extractor.set_face_map(extractor.get_face_map());
  physical_extractor.exclude_boundaries();
  physical_extractor.set_add_surrounded_faces(true);
  physical_extractor.set_make_patches(true);
  physical_extractor.set_expand_patches(patch_expand_value);
  physical_extractor.extract();
  
  // Show the mapped entity patches.
  {
    std::size_t patch_count = 1;
    const typename PhysicalExtractor::Mapped_entity_patches mapped_entity_patches =
      physical_extractor.get_mapped_entity_patches();
    
    if(!mapped_entity_patches.empty())
    {
    std::cout << "\nPatches found:\n";
    for(const typename PhysicalExtractor::Mapped_entity_patch &patch :
        mapped_entity_patches)
    {
      std::cout << "Patch " << patch_count << ": [";
      
      for(const typename PhysicalExtractor::Mapped_entity_patch::value_type &entity : patch)
      {
        std::cout << entity << ", ";
      }
      
      std::cout << "]\n";
      ++patch_count;
    }
    }
  }
  
  // Compute the minimal distance between mesh elements.
  if(!mapped_tri_tri_distances.empty())
  {
    double average_distance = 0.0;
    std::size_t count = 0;
    
    for(const std::pair<typename  
          CloseFacesExtractor::Mapped_triangle_triangle_distance::key_type,
          typename CloseFacesExtractor::Mapped_triangle_triangle_distance::mapped_type> 
          it : mapped_tri_tri_distances)
    {
      const double distance = it.second.second;
      average_distance += distance;
      ++count;
    }
    
    average_distance /= count;
    std::cout << "Average distance between faces in patches: " << average_distance << "\n";
  }
  
  // Identified patches.
  {
    CGAL::Color fcolor_gray(125, 125, 125);
        
    typename Surface_mesh::Property_map<typename Surface_mesh::Face_index, CGAL::Color>
    face_color_map = mesh.add_property_map<
      typename Surface_mesh::Face_index, CGAL::Color>("f:color", fcolor_gray).first;
      
    typename Surface_mesh::Property_map<typename Surface_mesh::Face_index, CGAL::Color>
    face_color_map_init = initial_mesh.add_property_map<
      typename Surface_mesh::Face_index, CGAL::Color>("f:color", fcolor_gray).first;
      
    const typename PhysicalExtractor::Patches &patches = physical_extractor.get_patches();
    std::srand(std::time(nullptr));
    
    int r, g, b;
    for(const auto &it_patches : patches)
    {
      r = std::rand() / ((RAND_MAX + 1u) / 255);
      g = std::rand() / ((RAND_MAX + 1u) / 255);
      b = std::rand() / ((RAND_MAX + 1u) / 255);
      CGAL::Color fcolor(r, g, b);
    
      for(const auto &it_patch : it_patches)
      {
        face_color_map[it_patch] = fcolor;
        face_color_map_init[it_patch] = fcolor;
      }
    }
  
    std::ofstream file_deformed_patches("patches_deformed_mesh.off");
    file_deformed_patches << mesh;
    
    std::ofstream file_patches("patches_mesh.off");
    file_patches << initial_mesh;
    
    mesh.remove_property_map(face_color_map);
    initial_mesh.remove_property_map(face_color_map_init);
  }
  
  gmsh::finalize();
  return EXIT_SUCCESS;
}