/// \author Leblanc Christophe.
/// \date 18/11/2020
/// \mail cleblancad@gmail.com

// No guards.

template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
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
{
  extract_physicals_dim0(nodes_physical_groups);
  // Note: extract_physicals_dim1 not implemented.
  extract_physicals_dim2(faces_physical_groups);
  extract_physicals_dim3(volumes_physical_groups);
}


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
extract_physicals_dim0(
  std::map<std::size_t, std::vector<int> > &nodes_physical_groups // Node number, physical groups ids it belongs to.
) const
{
  // Get physical groups for entity points.
  std::map<int, std::vector<int> > entity_points_physical_groups;
  {
    gmsh::vectorpair physical_groups;
    gmsh::model::getPhysicalGroups(physical_groups, 0);
    
    // Get the entity point tags corresponding to the physical groups.
    for(std::size_t i = 0; i < physical_groups.size(); ++i)
    {
      std::vector<int> entity_point_tags;
      gmsh::model::getEntitiesForPhysicalGroup(0, physical_groups[i].second, entity_point_tags);
      
      for(std::size_t j = 0; j < entity_point_tags.size(); ++j)
        entity_points_physical_groups[entity_point_tags[j]].push_back(
          physical_groups[i].second);
    }
  }
  
  // Get the correspondance between entity (geometry) point tags and element (mesh) point tags.
  if(entity_points_physical_groups.empty())
    return; // Nothing to do.

  std::map<int, std::size_t> entity_to_element_map;
  {
    // Get the types of 0D entities.
    std::vector<int> element_types;
    gmsh::model::mesh::getElementTypes(element_types, 0);
    
    for(const int it_element_type : element_types)
    {
      for(const std::pair<int, std::vector<int> > &it_entity_points_tags :
          entity_points_physical_groups)
      {
        std::vector<std::size_t> element_tags, node_tags;
        
        gmsh::model::mesh::getElementsByType(it_element_type, element_tags,
          node_tags, it_entity_points_tags.first);
      
        if(node_tags.size() == 1) // Normaly there is only one element point for each entity point.
          entity_to_element_map[it_entity_points_tags.first] = node_tags[0];
      }
    }
  }

  // Get the physicals corresponding to the element faces.
  for(const std::pair<int, std::vector<int> > &it_entity_points :
      entity_points_physical_groups)
  {
    const std::size_t &corresponding_node = 
     entity_to_element_map.at(it_entity_points.first);
  
    nodes_physical_groups[corresponding_node] = it_entity_points.second;
  }
}


/* // Simpler version of the above function, but seem to have some issues with the function
   // "gmsh::model::mesh::getNodesForPhysicalGroup".
template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
extract_physicals_dim0(
  std::map<std::size_t, std::vector<int> > &nodes_physical_groups // Node number, physical groups ids it belongs to.
) const
{
  // Get physical groups for entity points.
  gmsh::vectorpair physical_groups;
  gmsh::model::getPhysicalGroups(physical_groups, 0);

  // Get the nodes corresponding to these physical groups.
  for(std::size_t i = 0; i < physical_groups.size(); ++i)
  {
    std::vector<double> coords;
    std::vector<std::size_t> node_tags;

    gmsh::model::mesh::getNodesForPhysicalGroup(0, physical_groups[i].second,
      node_tags, coords);
      
    for(std::size_t j = 0; j < node_tags.size(); ++j)
      nodes_physical_groups[node_tags[j]].push_back(physical_groups[i].second);
  }
}*/


template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
extract_physicals_dim2(
  std::map<std::size_t, std::pair<std::vector<int>, std::vector<std::size_t> > > &faces_physical_groups // Face number, physical groups ids it belongs to (first) and corresponding node ids (second).
) const
{
  // Get physical groups for entity surfaces.
  std::map<int, std::vector<int> > entity_faces_physical_groups; // Face entity number.
  {
    gmsh::vectorpair physical_groups;
    gmsh::model::getPhysicalGroups(physical_groups, 2);
  
    // Get the entity face tags corresponding to the physical groups.
    for(std::size_t i = 0; i < physical_groups.size(); ++i)
    {
      std::vector<int> entity_face_tags;
      gmsh::model::getEntitiesForPhysicalGroup(2, physical_groups[i].second, entity_face_tags);
    
      for(std::size_t j = 0; j < entity_face_tags.size(); ++j)
        entity_faces_physical_groups[entity_face_tags[j]].push_back(physical_groups[i].second);
    }
  }

  // Get the correspondance between entity (geometry) face tags and element (mesh) face tags.
  if(entity_faces_physical_groups.empty())
    return; // Nothing to do.
  
  std::map<int, std::pair<std::vector<std::size_t>, std::vector<std::size_t> > >
    entity_to_elements_map;
  {
    // Get the types of the 2D entities.
    std::vector<int> element_types;
    gmsh::model::mesh::getElementTypes(element_types, 2);
    
    for(const int it_element_type : element_types)
    {
      for(const std::pair<int, std::vector<int> > &it_entity_faces_tags : 
          entity_faces_physical_groups)
      {
        std::vector<std::size_t> element_tags, node_tags;
    
        gmsh::model::mesh::getElementsByType(it_element_type, element_tags,
          node_tags, it_entity_faces_tags.first);
          
        entity_to_elements_map[it_entity_faces_tags.first] = 
          std::make_pair(element_tags, node_tags);
      }
    }
  }
  
  // Get the physicals corresponding to the element faces.
  for(const std::pair<int, std::vector<int> > &it_entity_faces :
      entity_faces_physical_groups)
  {
    const std::vector<std::size_t> &corresponding_elements =
      entity_to_elements_map.at(it_entity_faces.first).first;
      
    const std::vector<std::size_t> &corresponding_elements_nodes =
      entity_to_elements_map.at(it_entity_faces.first).second;
      
    const unsigned int nb_nodes_per_element = corresponding_elements_nodes.size() /
      corresponding_elements.size();

    for(std::size_t i = 0; i < corresponding_elements.size(); ++i)
    {
      // Get the nodes of the current corresponding element.
      std::vector<std::size_t> nodes_of_element;
      nodes_of_element.reserve(nb_nodes_per_element);
      
      for(unsigned int j = 0; j < nb_nodes_per_element; ++j)
        nodes_of_element.push_back(corresponding_elements_nodes[ 
          i*nb_nodes_per_element + j]);
    
      faces_physical_groups[corresponding_elements[i]] = std::make_pair(it_entity_faces.second,
        nodes_of_element);
      
      // Code for debugging purposes.
      /*std::cout << "element face " << it_corr << " has physicals: [";
      for(std::size_t ii = 0; ii < it_entity_faces.second.size(); ++ii)
        std::cout << it_entity_faces.second[ii] << ", ";
      std::cout << "]\n";*/
    }
  }
}

template<typename KernelT, typename Surface_meshT>
void ExtractSurfaceMeshFromMSH<KernelT, Surface_meshT>::
extract_physicals_dim3(std::map<std::size_t, std::vector<int> > 
  &volumes_physical_groups) const
{
  // Get physical groups for entity volumes.
  std::map<int, std::vector<int> > entity_volumes_physical_groups;
  {
    gmsh::vectorpair physical_groups;
    gmsh::model::getPhysicalGroups(physical_groups, 3);
    
    // Get the entity volume tags corresponding to the physical groups.
    for(std::size_t i = 0; i < physical_groups.size(); ++i)
    {
      std::vector<int> entity_volume_tags;
      gmsh::model::getEntitiesForPhysicalGroup(3, physical_groups[i].second, 
        entity_volume_tags);
      
      for(std::size_t j = 0; j < entity_volume_tags.size(); ++j)
        entity_volumes_physical_groups[entity_volume_tags[j]].push_back(
          physical_groups[i].second);
    }
  }
  
  // Note: no notion of "mesh volume". "Entity volume tags = element volume tags".
  
  // Get the physicals corresponding to the volume.
  for(const std::pair<int, std::vector<int> > &it_entity_volumes :
      entity_volumes_physical_groups)
  {
    const int num_int = it_entity_volumes.first;
    const std::size_t num = (num_int < 0) ? static_cast<std::size_t>(-num_int) : 
      static_cast<std::size_t>(num_int);
      
    volumes_physical_groups[num] = it_entity_volumes.second;
  }
}