/// \file extract_physical_patches.hpp
/// \brief Given a set of mesh triangles, a map of triangle ids and the corresponding MSH file, 
///   construct physical patches and return a copy of the MSH file with these 
///   physical patches added.
/// \author Leblanc Christophe.
/// \date 25/10/2020
/// \mail cleblancad@gmail.com

#ifndef EXTRACT_PHYSICAL_PATCHES_HPP
#define EXTRACT_PHYSICAL_PATCHES_HPP

#include <iostream>
#include <string.h>
#include <fstream>
#include <queue>
#include <vector>
#include <map>
#include <limits>
#include <set>
#include <algorithm>

#include <CGAL/boost/graph/selection.h>

#include "gmsh.h"

template<typename Surface_meshT, typename Face_mapT,
  typename Triangle_triangle_distanceT>
class ExtractPhysicalPatches
{
public:
  typedef Surface_meshT Surface_mesh;
  typedef Face_mapT Face_map;
  typedef Triangle_triangle_distanceT Triangle_triangle_distance;
  
  typedef typename Surface_mesh::Face_index Face_index;
  typedef typename Surface_mesh::Vertex_index Vertex_index;
  
  typedef std::vector<Face_index> Patch;
  typedef std::vector<Patch> Patches;
  
  typedef std::vector<std::size_t> Mapped_patch;
  typedef std::vector<Mapped_patch> Mapped_patches;
  
  typedef std::vector<int> Mapped_entity_patch;
  typedef std::vector<Mapped_entity_patch> Mapped_entity_patches;

  /// \brief Default constructor.
  ExtractPhysicalPatches(): m_tri_tri_map(nullptr),
    m_surface_mesh(nullptr), m_face_map(nullptr), m_add_boundaries(true),
    m_element_type(null_element_type()), m_add_surrounded_faces(false),
    m_make_patches(true),
    m_expand_face_selection(0) {}
  
  /// \brief Copy constructor.
  ExtractPhysicalPatches(const ExtractPhysicalPatches &rhs):
    m_tri_tri_map(rhs.m_tri_tri_map), // Shallow copy.
    m_mesh_file_name(rhs.m_mesh_file_name),
    m_surface_mesh(rhs.m_surface_mesh), // Shallow copy.
    m_face_map(rhs.m_face_map), // Shallow copy.
    m_add_boundaries(rhs.m_add_boundaries),
    m_patches(rhs.m_patches),
    m_element_type(rhs.m_element_type),
    m_add_surrounded_faces(rhs.m_add_surrounded_faces),
    m_make_patches(rhs.m_make_patches),
    m_expand_face_selection(rhs.m_expand_face_selection)
  {}

  /// \brief Destructor.
  ~ExtractPhysicalPatches() {}
  
  /// \brief copy operator.
  ExtractPhysicalPatches& operator=(const ExtractPhysicalPatches &rhs);

  /// \brief Get the mesh file name.
  const std::string& get_mesh_file_name() const
  { return m_mesh_file_name; }
  
  /// \brief Set the mesh file name.
  void set_mesh_file_name(const std::string &name)
  { m_mesh_file_name = name; }
  
  /// brief Get the map of triangles in contact.
  const Triangle_triangle_distance& get_triangle_triangle_distance() const
  { return *m_tri_tri_map; }
  
  /// \brief Set the map of triangles in contact.
  void set_triangle_triangle_distance(const Triangle_triangle_distance& ttd)
  { m_tri_tri_map = &ttd; }
  
  /// \brief Get if surface boundaries are added or not in the physical surface patches.
  bool get_boundaries() const
  { return m_add_boundaries; }
  
  /// \brief Surface boundaries in contact are added to the physical patches.
  void include_boundaries()
  { m_add_boundaries = true; }
  
  /// \brief Surface boundaries in contact are excluded from the physical patches.
  void exclude_boundaries()
  { m_add_boundaries = false; }
  
  /// \brief Get the surface mesh.
  const Surface_mesh& get_surface_mesh() const
  { return *m_surface_mesh; }
  
  /// \brief Set the surface mesh.
  void set_surface_mesh(const Surface_mesh &mesh)
  { m_surface_mesh = &mesh; }
  
  /// \brief Get the node map.
  const Face_map& get_face_map() const
  { return *m_face_map; }
  
  /// \brief Set the node map.
  void set_face_map(const Face_map& map)
  { m_face_map = &map; }
  
  /// \brief Get found patches of surfaces in contact.
  const Patches& get_patches() const
  { return m_patches; }
  
  /// \brief Get the number of patches.
  std::size_t number_of_patches() const
  { return m_patches.size(); }
  
  /// \brief Get the mapped patches.
  Mapped_patches get_mapped_patches() const;
  
  /// \brief Retrieve the entities patches corresponding to the element patches.
  Mapped_entity_patches get_mapped_entity_patches(const Mapped_patches &mapped_patches) const;
  
  /// \brief Get the mapped entity patches.
  Mapped_entity_patches get_mapped_entity_patches() const
  {
    Mapped_entity_patches ret;
    const Mapped_patches mapped_patches = get_mapped_patches();
    
    if(mapped_patches.empty()) return ret;
    ret = get_mapped_entity_patches(mapped_patches);
    
    return ret;
  }
  
  /// \brief Get the element type.
  int get_element_type() const
  { return m_element_type; }
  
  /// \brief Set the element type.
  void set_element_type(int e)
  { m_element_type = e; }
  
  /// \brief Get if faces surrounded by patch faces are added to this patch.
  bool get_add_surrounded_faces() const
  { return m_add_surrounded_faces; }
  
  /// \brief Set if faces surrounded by patch faces are added to this patch.
  void set_add_surrounded_faces(bool add)
  { m_add_surrounded_faces = add; }
  
  /// \brief Get if patches are made.
  bool get_make_patches() const
  { return m_make_patches; }
  
  /// \brief Set if sets of connected patches should be made.
  void set_make_patches(const bool p)
  { m_make_patches = p; }
  
  /// \brief Get the expand patches parameter.
  unsigned int get_expand_patches() const
  { return m_expand_face_selection; }
  
  /// \brief Set the expand patches parameter.
  void set_expand_patches(const unsigned int e)
  { m_expand_face_selection = e; }
  
  /// \brief Extract patches of physical given in the triangles map.
  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;
  
  enum Boundary_type
  {
    BOUNDARY_UNKNOWN = 0, ///< Do not known if belongs to a boundary or not.
    BOUNDARY_LEFT,
    BOUNDARY_RIGHT,
    BOUNDARY_BOTTOM,
    BOUNDARY_TOP,
    BOUNDARY_FRONT,
    BOUNDARY_REAR,
    BOUNDARY_NONE ///< Do not belong to a boundary.
  };

  const Triangle_triangle_distance* m_tri_tri_map;
  std::string m_mesh_file_name;
  const Surface_mesh* m_surface_mesh;
  const Face_map* m_face_map;
  bool m_add_boundaries;
  Patches m_patches;
  int m_element_type;
  bool m_add_surrounded_faces;
  bool m_make_patches;
  unsigned int m_expand_face_selection;
  
  /// \brief Extract patches of connected surfaces in contact.
  void extract_patches(const Surface_mesh &mesh,
    Patches &patches,
    const Triangle_triangle_distance &tri_tri_map) const;
    
  // Add faces completely surrouned by faces in the current patch.
  void add_surrounded_faces_by_patch_faces(const Surface_mesh &mesh,
    Patch &current_patch) const;
    
  /// \brief Insert into the queue the list of face in contact (in tri_tri_map)
  /// and that share at least common vertex with the current triangle.
  void get_neighbourhood_contact(const Surface_mesh &mesh,
    std::queue<Face_index> &queue,
    const Triangle_triangle_distance &tri_tri_map, 
    std::set<Face_index> &visited,
    Patch &patch,
    const std::array<double, 3> &min_box_point, 
    const std::array<double, 3> &max_box_point) const;
    
  /// \brief Compute the mesh box.
  void get_surface_mesh_box(std::array<double, 3> &min, std::array<double, 3> &max,
    const Surface_mesh &mesh) const;
    
  /// \brief Tell if a face index is in the triangle-triangle map.
  bool is_in_tri_tri_map(const Face_index &f, const Triangle_triangle_distance &tri_tri_map) const;
  
  /// \brief Tell if a face is located on a mesh boundary.
  bool is_on_boundary(const Surface_mesh &mesh, const Face_index &f,
    const std::array<double, 3> &min, const std::array<double, 3> &max) const;
    
  /// \brief Expand patches.
  void expand_face_selection(const Surface_mesh &mesh, 
    const unsigned int expand, Patches &patches) const;
  
  /// \brief Get the maximum id of the 2D physical groups.
  int get_max_physical_group_id(const gmsh::vectorpair &existing_physical_groups) const;
  
  /// \brief Remove faces from boundaries.
  void remove_faces_from_boundaries(const Surface_mesh &mesh, Patches &patches) const;
  
  // 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());
  }
  
  /// \class Vector_pmap_wrapper
  /// \brief Wrapper for the cgal expand_face_selection function.
  /// Source: https://doc.cgal.org/latest/Polygon_mesh_processing/Polygon_mesh_processing_2corefinement_difference_remeshed_8cpp-example.html#a7
  template<typename MeshT>
  struct Vector_pmap_wrapper
  {
    std::vector<bool> &vect;
    typedef typename boost::graph_traits<MeshT>::face_descriptor
      face_descriptor;
      
    Vector_pmap_wrapper(std::vector<bool> &v): vect(v) {}
    
    friend bool get(const Vector_pmap_wrapper &m, face_descriptor f)
    { return m.vect[f]; }
    
    friend void put(const Vector_pmap_wrapper &m, face_descriptor f, bool b)
    { m.vect[f] = b; }
  };
};


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>& 
ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
operator=(const ExtractPhysicalPatches &rhs)
{
  if(this != &rhs)
  {
    m_tri_tri_map = rhs.m_tri_tri_map; // Shallow copy.
    m_mesh_file_name = rhs.m_mesh_file_name;
    m_surface_mesh = rhs.m_surface_mesh; // Shallow copy.
    m_face_map = rhs.m_face_map; // Shallow copy.
    m_add_boundaries = rhs.m_add_boundaries;
    m_patches = rhs.m_patches;
    m_element_type = rhs.m_element_type;
    m_add_surrounded_faces = rhs.m_add_surrounded_faces;
    m_make_patches = rhs.m_make_patches;
    m_expand_face_selection = rhs.m_expand_face_selection;
  }
  
  return *this;
}

template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
extract()
{
  // Checks.
  if(m_tri_tri_map == nullptr)
  {
    std::cerr << "Error: triangle-triangle map not set." << std::endl;
    return;
  }
  
  if(m_surface_mesh == nullptr)
  {
    std::cerr << "Error: surface mesh not set." << std::endl;
    return;
  }
  
  if(m_face_map == nullptr)
  {
    std::cerr << "Error: face map not set." << std::endl;
    return;
  }
  
  // Get the file stream on the file.
  if(m_mesh_file_name.empty())
  {
    std::cerr << "Error: no file name set." << std::endl;
    return;
  }
  
  // Check if exists.
  std::ifstream file;
  file.open(m_mesh_file_name.c_str(), std::ifstream::in);

  if(!file.is_open())
  {
    std::cerr << "Error: unable to open file: \"" << m_mesh_file_name << "\"" << std::endl;
    return;
  }
  
  file.close();
  
  //
  // Perform operations.
  //
  // 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);
  
  std::cout << "Extract patches...\n";
  
  // Get the patches of adjacent surfaces in contact.
  if(m_make_patches)
  {
    extract_patches(*m_surface_mesh, m_patches, *m_tri_tri_map);
  }
  else
  {
    // Compute mesh box.
    std::array<double, 3> min_box_point, max_box_point;
    get_surface_mesh_box(min_box_point, max_box_point, *m_surface_mesh);
  
    for(const std::pair<Face_index, std::pair<Face_index, double> > &it : 
      *m_tri_tri_map)
    {
      // Skip boundary faces if asked by the user.
      if(!m_add_boundaries && 
         (is_on_boundary(*m_surface_mesh, it.first, min_box_point, max_box_point) ||
          is_on_boundary(*m_surface_mesh, it.second.first, min_box_point, max_box_point)) )
      {
        continue;
      }
    
      Patch patch(1);
      patch[0] = it.first;
      m_patches.push_back(patch);
      
      patch[0] = it.second.first;
      m_patches.push_back(patch);
    }
  }
  
  // Should have an even number of patches (always as they are associated two-by-two).
  if(m_patches.size() % 2 != 0)
    throw std::runtime_error("Number of patches not even.");
  
  // Expand face selection.
  if(m_expand_face_selection > 0)
    expand_face_selection(*m_surface_mesh, m_expand_face_selection, m_patches);
 
  // Remove faces from boundaries if asked by the user.
  if(!m_add_boundaries)
    remove_faces_from_boundaries(*m_surface_mesh, m_patches);
  
  std::cout << "Found " << m_patches.size() << " patch" 
            << (m_patches.size() > 1 ? "es" : "") << ". ("
            << (m_add_boundaries ? "including" : "excluding")
            << " boundaries)\n";
  std::cout << "done.\n";
  
  // Get the corresponding mesh element patches in the corresponding MSH file.
  const Mapped_patches mapped_patches = get_mapped_patches();
  
  if(mapped_patches.empty())
  {
    std::cout << "No patch found. Nothing to do and no MSH file written.\n";
    return;
  }
  
  // Retrieve the entities patches corresponding to the element patches
  // (needed as the physicals are defined at the entity (geometric) level).
  // element <-> mesh
  // entity  <-> geometry
  Mapped_entity_patches mapped_entity_patches = get_mapped_entity_patches(mapped_patches);
  
  // Write physicals (if any).
  std::cout << "Modifying \"" << models_open[0] << "\" by adding patches found...\n";
  
  // Get already existing surface physical groups.
  gmsh::vectorpair existing_physical_groups;
  gmsh::model::getPhysicalGroups(existing_physical_groups, 2);
  
  // Get the biggest physical group id for 2D entities.
  int max_physical_group_id = get_max_physical_group_id(existing_physical_groups);
  const int start_max_physical_group_id = max_physical_group_id;
  
  // Add the patches as new 2D physical groups.
  for(const Mapped_entity_patch mapped_entity_patch : mapped_entity_patches)
  {
    if(!mapped_entity_patch.empty()) // Empty mapped entity patch should not happen normaly...
    {
      gmsh::model::addPhysicalGroup(2, mapped_entity_patch, max_physical_group_id);
    
      std::cout << "Adding to model physical group id: " << max_physical_group_id << "\n";
      ++max_physical_group_id;
    }
  }
  
  const std::string filename = models_open[0] + "_modified.msh";
  gmsh::write(filename);
  
  std::cout << "Modified MSH file written to: " << filename << "\n";  
  std::cout << "done.\n";
  
  // Tell what python code to write.
  {
    std::cout << "\n\nAdd the following code to your python file:\n\n";
    std::cout << "\ncontact_list = []\n\n";
    std::cout << "for i in range(" 
              << start_max_physical_group_id
              << ", "
              << (max_physical_group_id-1)
              << "+1, 2):\n";
    std::cout << "\tcontact_list.append(dG3DNodeOnSurfaceContactDomain(1000, 2, i, 2, i+1, penalty, thick))\n";
    std::cout << "\tcontact_list.append(dG3DNodeOnSurfaceContactDomain(1000, 2, i+1, 2, i, penalty, thick))\n\n";
    std::cout << "for contact in contact_list:\n";
    std::cout << "\tmysolver.defoDefoContactInteraction(contact)\n";
    std::cout << "\n\n";
  }
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
extract_patches(const Surface_mesh &mesh,
  Patches &patches,
  const Triangle_triangle_distance &tri_tri_map) const
{
  // Compute mesh box.
  std::array<double, 3> min_box_point, max_box_point;
  get_surface_mesh_box(min_box_point, max_box_point, *m_surface_mesh);
  
  // Determine patches.
  std::set<Face_index> visited;
  for(const std::pair<Face_index, std::pair<Face_index, double> > &it : 
      tri_tri_map)
  {
    // Skip already visited triangles.
    if(visited.find(it.first) != visited.end())
      continue;
      
    // Skip boundary faces if asked by the user.
    // If face alone has closest faces on boundaries, ignore both groups of faces.
    if(!m_add_boundaries && 
       (is_on_boundary(mesh, it.first, min_box_point, max_box_point) ||
        is_on_boundary(mesh, it.second.first, min_box_point, max_box_point)) )
      continue;
    
    // First patch.
    Patch current_patch;
  
    std::queue<Face_index> queue;
    queue.push(it.first);
    visited.insert(it.first);
    
    current_patch.push_back(it.first);
    
    while(!queue.empty())
    {
      // Insert into the queue the list of non-visited faces in contact (in tri_tri_map)
      // and that share at least a common edge with the current triangle.
      get_neighbourhood_contact(mesh, queue, tri_tri_map, visited, current_patch,
        min_box_point, max_box_point);    
    } // While queue is not empty.
  
    // Add faces completely surrouned by faces in the current patch.
    if(m_add_surrounded_faces)
      add_surrounded_faces_by_patch_faces(mesh, current_patch);
  
    patches.push_back(current_patch);
    
    // Second patch.
    current_patch.clear();
    while(!queue.empty()) queue.pop();
    
    queue.push(it.second.first);
    visited.insert(it.second.first);
    current_patch.push_back(it.second.first);
    
    while(!queue.empty())
    {
      // Insert into the queue the list of non-visited faces in contact (in tri_tri_map)
      // and that share at least a common edge with the current triangle.
      get_neighbourhood_contact(mesh, queue, tri_tri_map, visited, current_patch,
        min_box_point, max_box_point);    
    } // While queue is not empty.
  
    // Add faces completely surrouned by faces in the current patch.
    if(m_add_surrounded_faces)
      add_surrounded_faces_by_patch_faces(mesh, current_patch);
  
    patches.push_back(current_patch);
  } // Loop on triangles in contact not already visited.
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
get_neighbourhood_contact(const Surface_mesh &mesh,
   std::queue<Face_index> &queue,
   const Triangle_triangle_distance &tri_tri_map, 
   std::set<Face_index> &visited,
   Patch &patch,
   const std::array<double, 3> &min_box_point, 
   const std::array<double, 3> &max_box_point) const
{
  const Face_index current_tri_id = queue.front();
  queue.pop();
  
  // Loop on the neighbouring faces of the current triangular face.
  {
    if(current_tri_id == Surface_mesh::null_face())
      return;
  
    CGAL::Face_around_face_iterator<Surface_mesh> ff_it, ff_ite;
    const typename Surface_mesh::Halfedge_index he_tri = mesh.halfedge(current_tri_id);

    if(mesh.is_border(he_tri))
      return;
    
    boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(he_tri, mesh);
    
    for(; ff_it != ff_ite; ++ff_it)
    {
      // Exclude visited faces.
      if(visited.find(*ff_it) != visited.end())
        continue;
        
      // Exclude faces on boundaries if asked by the user.
      if(!m_add_boundaries && is_on_boundary(mesh, *ff_it, min_box_point, max_box_point))
        continue;
        
      // Exclude faces that are not in contact (i.e. not if the tri-tri map).
      if(!is_in_tri_tri_map(*ff_it, tri_tri_map))
        continue;
      
      // Add the current face to the queue.
      queue.push(*ff_it);
      visited.insert(*ff_it);
      patch.push_back(*ff_it);
    } // Loop on neighbouring faces.
  }
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
get_surface_mesh_box(std::array<double, 3> &min, std::array<double, 3> &max,
  const Surface_mesh &mesh) const
{
  // Loop on vertices.
  if(mesh.number_of_vertices() == 0) // Nothing to do.
  {
    min[0] = min[1] = min[2] = max[0] = max[1] = max[2] =
      std::numeric_limits<double>::quiet_NaN();
  }

  typename Surface_mesh::Vertex_range::const_iterator v_it = mesh.vertices().begin(),
    v_ite = mesh.vertices().end();
    
  Point p = mesh.point(*v_it);
  min[0] = max[0] = p[0];
  min[1] = max[1] = p[1];
  min[2] = max[2] = p[2];

  for(; v_it != v_ite; ++v_it)
  {
    p = mesh.point(*v_it);
  
    for(int i = 0; i < 3; ++i)
    {
      if(min[i] > p[i]) min[i] = p[i];
      if(max[i] < p[i]) max[i] = p[i];
    }
  }
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
add_surrounded_faces_by_patch_faces(const Surface_mesh &mesh,
  Patch &patch) const
{
  // Find faces adjacent to the current patch.
  std::set<Face_index> adjacent_faces;
  Patch patch_copy = patch;
  bool patch_modified = false;
  
  for(const Face_index &f_patch : patch)
  {
    if(f_patch == Surface_mesh::null_face())
      continue;
  
    CGAL::Face_around_face_iterator<Surface_mesh> ff_it, ff_ite;
    const typename Surface_mesh::Halfedge_index he_patch = mesh.halfedge(f_patch);
    
    if(mesh.is_border(he_patch))
      continue;
    
    boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(he_patch, mesh);
    
    for(; ff_it != ff_ite; ++ff_it)
      adjacent_faces.insert(*ff_it);
  }
  
  // Find surrounded faces.
  for(const Face_index &f_adj : adjacent_faces)
  {
    // Loop on connected faces of f_patch.
    if(f_adj == Surface_mesh::null_face())
      continue;
    
    CGAL::Face_around_face_iterator<Surface_mesh> ff_it, ff_ite;
    const typename Surface_mesh::Halfedge_index he_adj = mesh.halfedge(f_adj);
    
    if(mesh.is_border(he_adj))
      continue;
    
    boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(he_adj, mesh);
    bool surrounded = true;
  
    for(; ff_it != ff_ite; ++ff_it)
    {
      // Adjacent face not completely surrounded, exclude it.
      if(std::find(patch.begin(), patch.end(), *ff_it) == patch.end())
      {
        surrounded = false;
        break;
      }
    }
    
    // Add new face to patch, avoiding duplicates.
    if(surrounded && 
       std::find(patch_copy.begin(), patch_copy.end(), f_adj) == patch_copy.end()) 
    {
      patch_modified = true;
      patch_copy.push_back(f_adj);
    }
  }
  
  if(patch_modified)
    patch = patch_copy;
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
bool ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
is_in_tri_tri_map(const Face_index &f, const Triangle_triangle_distance &tri_tri_map) const
{
  for(const std::pair<std::size_t, std::pair<std::size_t, double> > it :
      tri_tri_map)
  {
    if(f == it.first || f == it.second.first)
      return true;
  }

  return false;
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
bool ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
is_on_boundary(const Surface_mesh &mesh, const Face_index &f,
  const std::array<double, 3> &min_c, const std::array<double, 3> &max_c) const
{
  const double factor = 10.0;

  // Loop on the vertices of the face.
  if(f == Surface_mesh::null_face())
    return true;
  
  CGAL::Vertex_around_face_iterator<Surface_mesh> vf_it, vf_ite;
  const typename Surface_mesh::Halfedge_index he = mesh.halfedge(f);
  
  if(mesh.is_border(he))
    return true;
  
  boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(he, mesh);
  std::map< Vertex_index, std::set<Boundary_type> > v_boundaries;
  unsigned int count = 0;
  
  for(; vf_it != vf_ite; ++vf_it, ++count)
  {
        const Point &p = mesh.point(*vf_it);

        if(nearly_equal(p[0], min_c[0], factor))
        {
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_LEFT);
        }

        if(nearly_equal(p[0], max_c[0], factor))
        {
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_RIGHT);
        }

        if(nearly_equal(p[1], min_c[1], factor))
        {
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_FRONT);
        }

        if(nearly_equal(p[1], max_c[1], factor))
        {
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_REAR);
        }

        if(nearly_equal(p[2], min_c[2], factor))
        {
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_BOTTOM);
        }

        if(nearly_equal(p[2], max_c[2], factor))
        {
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_TOP);
        }
  }
  
  if(v_boundaries.size() == count) // All face's vertices on a boundary.
    return true;
  else
    return false;
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
typename ExtractPhysicalPatches<Surface_meshT, Face_mapT,
 Tri_tri_distanceT>::Mapped_patches 
ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
get_mapped_patches() const
{
  Mapped_patches mapped_patches;

  if(m_face_map == nullptr)
  {
    std::cerr << "Error: no node map set." << std::endl;
    return mapped_patches;
  }

  // Inverse face map.
  std::map<Face_index, std::size_t> inverse_face_map;
  
  for(const std::pair<std::size_t, Face_index> &it : *m_face_map)
    inverse_face_map[it.second] = it.first;
  
  // Construct the mapped patches.
  mapped_patches.reserve(m_patches.size());
  for(const Patch& patch : m_patches)
  {
    Mapped_patch mapped_patch;
    mapped_patch.reserve(patch.size());
  
    for(const Face_index &f : patch)
    {
      mapped_patch.push_back(inverse_face_map[f]);
    }
    
    mapped_patches.push_back(mapped_patch);
  }
  
  return mapped_patches;
}


// Retrieve the entities patches corresponding to the element patches.
template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
typename ExtractPhysicalPatches<Surface_meshT, Face_mapT,
 Tri_tri_distanceT>::Mapped_entity_patches 
ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
get_mapped_entity_patches(const Mapped_patches &mapped_patches) const
{
  Mapped_entity_patches mapped_entity_patches;

  // Check.
  if(m_element_type == null_element_type())
  {
    std::cerr << "Error: no element type set." << std::endl;
    return mapped_entity_patches;
  }
  
  // Shortcut.
  if(mapped_patches.empty())
  {
    return mapped_entity_patches;
  }

  // Get the tags of the 2D entities.
  gmsh::vectorpair entity_tags;
  gmsh::model::getEntities(entity_tags, 2);
  
  // For each 2D entity, get the corresponding elements inside it.
  std::map<std::size_t, int> element_to_entity_map; // Stores for each element (mesh) the corresponding entity (geometry).
  
  for(const typename gmsh::vectorpair::value_type &it_entity : entity_tags)
  {
    std::vector<std::size_t> element_tags;
    std::vector<std::size_t> node_tags;
  
    gmsh::model::mesh::getElementsByType(m_element_type,
      element_tags, node_tags, it_entity.second);
      
    for(std::size_t i = 0; i < element_tags.size(); ++i)
      element_to_entity_map[element_tags[i]] = it_entity.second;
    
    // Code the debugging purposes.
    /*std::cout << "entity tag: " << it_entity.second << " has elements: [";
    for(std::size_t i = 0; i < element_tags.size(); ++i)
      std::cout << element_tags[i] << ", ";
    std::cout << "]\n";*/
  }
  
  // Construct the mapped entity patches.
  for(const Mapped_patch &mapped_patch : mapped_patches)
  {
    Mapped_entity_patch mapped_entity_patch;
    
    for(const typename Mapped_patch::value_type &it : mapped_patch)
    {
      mapped_entity_patch.push_back( element_to_entity_map.at(it) );
    }
    
    mapped_entity_patches.push_back(mapped_entity_patch);
  }
  
  return mapped_entity_patches;
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
expand_face_selection(const Surface_mesh &mesh, 
  const unsigned int expand, Patches &patches) const
{
  for(Patch &patch : patches)
  {
    Patch new_patch = patch;
    
    std::vector<bool> is_selected(mesh.number_of_faces(), false);
    
    for(const typename Patch::value_type &f : new_patch)
    {
      is_selected[f] = true;
    }
    
    CGAL::expand_face_selection(new_patch, mesh, expand,
      Vector_pmap_wrapper<Surface_mesh>(is_selected),
      std::back_inserter(new_patch));
    
    patch = new_patch;
  }
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
int ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
get_max_physical_group_id(const gmsh::vectorpair &existing_physical_groups) const
{
  int max_physical_group_id = std::numeric_limits<int>::min();
  for(std::size_t i = 0; i < existing_physical_groups.size(); ++i)
  {
    if(max_physical_group_id < existing_physical_groups[i].second)
    {
      max_physical_group_id = existing_physical_groups[i].second;
    }
  }
  
  ++max_physical_group_id;
  
  return max_physical_group_id;
}


template<typename Surface_meshT, typename Face_mapT, typename Tri_tri_distanceT>
void ExtractPhysicalPatches<Surface_meshT, Face_mapT, Tri_tri_distanceT>::
remove_faces_from_boundaries(const Surface_mesh &mesh, Patches &patches) const
{
  // Compute mesh box.
  std::array<double, 3> min_box_point, max_box_point;
  get_surface_mesh_box(min_box_point, max_box_point, *m_surface_mesh);
  
  typename Patches::iterator it_patches = patches.begin();
    
  for(; it_patches != patches.end(); )
  {
    Patch &patch = *it_patches;
    
    typename Patch::iterator it = patch.begin();
    for(; it != patch.end(); )
    {
      if(is_on_boundary(mesh, *it, min_box_point, max_box_point))
        it = patch.erase(it);
      else
        ++it;
    }
    
    if(patch.empty())
      it_patches = patches.erase(it_patches);
    else
      ++it_patches;
  }
}

#endif // EXTRACT_PHYSICAL_PATCHES_HPP