/**============================================================================
 *
 * \file   gmsh_writer.hxx
 * \brief  Implementations for the class GmshWriter.
 *
 * \author Leblanc Christophe.
 * \date   23/02/2016
 * \email cleblancad@gmail.com
 *
 *===========================================================================*/

#ifndef GMSH_WRITER_HXX
#define GMSH_WRITER_HXX

#include "gmsh_writer.hpp"

namespace marinov {

template<typename Face_decimationT>
GmshWriter<Face_decimationT>::GmshWriter(const std::string &fileName, 
			      MeshType &globalMesh,
		              double scale, 
			      unsigned int printPrecision,
                              std::ostream &os):
  m_FileName(fileName),
  m_Mesh(&globalMesh),
  m_Scale(scale),
  m_PrintPrecision(printPrecision + 1),
  m_Verbose(false),
  m_Os(&os)
{}

template<typename Face_decimationT>
GmshWriter<Face_decimationT>::GmshWriter(const GmshWriter &rhs):
  m_FileName(rhs.m_FileName),
  m_Mesh(rhs.m_Mesh),
  m_Scale(rhs.m_Scale),
  m_PrintPrecision(rhs.m_PrintPrecision),
  m_Verbose(rhs.m_Verbose),
  m_Os(rhs.m_Os),
  m_VertexData(rhs.m_VertexData),
  m_EdgeData(rhs.m_EdgeData),
  m_FaceData(rhs.m_FaceData),
  m_MeshProperties(rhs.m_MeshProperties)
{}

template<typename Face_decimationT>
GmshWriter<Face_decimationT>& GmshWriter<Face_decimationT>::operator=(const GmshWriter &rhs)
{
  if(this != &rhs)
  {
    m_FileName   = rhs.m_FileName;
    m_Mesh       = rhs.m_Mesh;
    m_Scale      = rhs.m_Scale;
    m_PrintPrecision = rhs.m_PrintPrecision;
    m_Verbose = rhs.m_Verbose;
    m_Os = rhs.m_Os;
    m_VertexData = rhs.m_VertexData;
    m_EdgeData = rhs.m_EdgeData;
    m_FaceData = rhs.m_FaceData;
    m_MeshProperties = rhs.m_MeshProperties;
  }
  
  return *this;
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::Write()
{
  std::ofstream file;
  file.open(m_FileName.c_str());

  if(file.is_open())
    WriteHelper(file);
  else
    throw std::runtime_error("GmshWriter was unable to create file.");

  file.close();
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::WriteHelper(std::ostream &os)
{
  // Prints data to 'os'.
  PrintHeader(os);
  PrintVertices(os);
  PrintEdges(os);
  PrintFaces(os);
  PrintVolumes(os);
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::PrintHeader(std::ostream &os) const
{
  // Computes the current time.
  // Copied from: http://www.cplusplus.com/reference/ctime/asctime
  time_t rawtime;
  struct tm *timeinfo;
  time(&rawtime);
  timeinfo = localtime(&rawtime);
  
  os << "/*********************************************************************\n"
   << "*\n"
   << "* Geo file generated by Tomoprocess, part of CadXFem.\n"
   << "* Date: "
   << asctime(timeinfo)
   << "* Version: "
   << "1.0"
   << "\n*\n"
   << "*\n"
   << "* Leblanc Christophe\n"
   << "* cleblancad@gmail.com\n"
   << "* University of Liege\n"
   << "*\n"
   << "*********************************************************************/\n\n"
   << "\n// Global characteristic length.\n"
   << "scale " << " = "
   << std::fixed << std::setprecision(m_PrintPrecision)
   << m_Scale
   << "; // This is a default value. Change it according to your needs.\n";
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::PrintVertices(std::ostream &os)
{
  // Print a comment.
  os << "\n//\n// Vertex coordinates.\n//\n";
  std::size_t vertexCounter = 1;
  
  typename MeshType::Vertex_range::iterator v_it = m_Mesh->vertices().begin(),
      v_ite = m_Mesh->vertices().end();

  for(; v_it != v_ite; ++v_it) 
  {
    const typename MeshType::Point &p = m_Mesh->point(*v_it);

    os << "Point(" << vertexCounter << ") = {"
       << std::setprecision(m_PrintPrecision) << p[0] << ", " 
       << std::setprecision(m_PrintPrecision) << p[1] << ", " 
       << std::setprecision(m_PrintPrecision) << p[2] << ", "
       << "scale" << "};\n";

    m_VertexData[*v_it].num = vertexCounter;
    ++vertexCounter;
  }
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::PrintEdges(std::ostream &os)
{
  // Print a comment.
  os << "\n//\n// Line definitions.\n//\n";
  long int edgeCounter = 1;

  typename MeshType::Face_range::iterator f_it = m_Mesh->faces().begin(),
    f_ite = m_Mesh->faces().end();

  for(; f_it != f_ite; ++f_it) 
  {
    // Loop on the edges of the current face.
    CGAL::Edge_around_face_iterator<MeshType> fe_it, fe_ite;
    boost::tie(fe_it, fe_ite) = CGAL::edges_around_face(m_Mesh->halfedge(*f_it), *m_Mesh);
    
    for(; fe_it != fe_ite; ++fe_it)
    {
      if(!m_EdgeData[*fe_it].IsVisited()) // Do not add already visited edges.
      {
        // Set current edge as visited.
        m_EdgeData[*fe_it].SetVisited();
        m_EdgeData[*fe_it].num = edgeCounter;

        // Get the vertices handle connected to the current edge.
        const typename MeshType::Halfedge_index he = m_Mesh->halfedge(*fe_it, 0);
        typename MeshType::Vertex_index vh1 = m_Mesh->source(he),
          vh2 = m_Mesh->target(he);

        // Print edge/line.
        os << "Line(" << edgeCounter <<") = {"
           << (m_VertexData[vh1].num) << ", " << (m_VertexData[vh2].num) << "};\n";

        ++edgeCounter;
      }
    }
  }
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::PrintFaces(std::ostream &os)
{
  // Print a comment.
  os << "\n//\n// Surfaces definitions.\n//\n";
  std::size_t faceCounter = 1;

  //
  // Loop on faces of the global mesh.  
  //
  typename MeshType::Face_range::iterator f_it = m_Mesh->faces().begin(),
    f_ite = m_Mesh->faces().end();

  for(; f_it != f_ite; ++f_it)
  {
    // Check if the edges belonging to the current face
    // are correctly oriented for Gmsh.
    CheckEdgeOrientations(*m_Mesh, *f_it);

    // Loop on the edges of the current face.
    // Write line loops.
    os << "Line Loop(" << faceCounter << ") = {";

    CGAL::Edge_around_face_iterator<MeshType> fe_it, fe_ite;
    boost::tie(fe_it, fe_ite) = CGAL::edges_around_face(m_Mesh->halfedge(*f_it), *m_Mesh);
    
    os << (m_EdgeData[*fe_it].num * m_EdgeData[*fe_it].tag);  
    ++fe_it;

    std::size_t face_num = 1;
    for(; fe_it != fe_ite; ++fe_it, ++face_num)
    {
      if((face_num % 10) == 0) // Avoid too long line in the text file.
        os << ",\n";
      else
        os << ", ";

      os << (m_EdgeData[*fe_it].num * m_EdgeData[*fe_it].tag);
    }

    os << "};\n";
    ++faceCounter;

    // Write surfaces.
    m_FaceData[*f_it].num = faceCounter;

    // Plane surface.
    os << "Plane Surface(" << faceCounter << ") = {"
       << (faceCounter-1) << "};\n";

    ++faceCounter;
  }
}

template<typename Face_decimationT>
void GmshWriter<Face_decimationT>::PrintVolumes(std::ostream &os)
{
  // Print a comment.
  os << "\n//\n// Volumes definitions.\n//\n";

  //
  // Global mesh.
  //
  // Loop on faces of the current mesh.
  // Write surface loop.
  os << "Surface Loop(1) = {";

  typename MeshType::Face_range::iterator f_it = m_Mesh->faces().begin(),
    f_ite = m_Mesh->faces().end();

  if(f_it != f_ite) 
  {
    os << m_FaceData[*f_it].num;
    ++f_it;
  }

  int face_num = 1;
  for(; f_it != f_ite; ++f_it, ++face_num)
  {
    if((face_num % 10) == 0) // Avoid too long lines in the text file.
      os << ",\n";
    else
      os << ", ";

    os << m_FaceData[*f_it].num;
  }

  os << "};\n";

  // Write volume.
  os << "Volume(2) = {1};";


  // Write physical volume.
  os << "Physical Volume(\"Foam volume\") = {" 
     << "2" << "};\n";

  // Write physical faces.
  {
    std::map<BoundaryType, std::vector<long int>> physical_faces;

    typename MeshType::Face_range::iterator f_it = m_Mesh->faces().begin(),
      f_ite = m_Mesh->faces().end();

    for(; f_it != f_ite; ++f_it)
    {
      if(m_MeshProperties->at(*f_it).boundary != BoundaryType::BOUNDARY_UNKNOWN)
      {
        physical_faces[m_MeshProperties->at(*f_it).boundary].push_back(
          m_FaceData[*f_it].num);
      }
    }

    for(const std::pair<BoundaryType, std::vector<long int>> &it : physical_faces)
    {
      std::size_t face_counter = 0;
      typename std::vector<long int>::const_iterator it2 = it.second.begin(),
        ite2 = it.second.end();

      if(it2 == ite2) continue;

      os << "Physical Surface(";
         if(BoundaryType::BOUNDARY_NONE == it.first)
           os << "\"interior\"";
         else
           os << it.first;
       os << ") = {" << (*it2);
      ++it2;

      for(; it2 != ite2; ++it2)
      {
        if((face_counter % 10) == 0)
          os << ",\n";
        else
          os << ", ";

        os << (*it2);
        ++face_counter;
      }

      os << "};\n";
    }
  }
}

template<typename MeshT>
void GmshWriter<MeshT>::CheckEdgeOrientations(const MeshType &mesh,
  const typename MeshType::Face_index &f)
{
  // Check if the first two edges of the face are correctly oriented.
  CGAL::Edge_around_face_iterator<MeshType> fe_it, fe_ite;
  boost::tie(fe_it, fe_ite) = CGAL::edges_around_face(mesh.halfedge(f), mesh);

  CGAL::Edge_around_face_iterator<MeshType> fe_it2 = fe_it;

  /** Looking for the common vertice between the first
      two segments.

      4 possibilities: {a, b} and {c, a}
                       {a, b} and {a, c}
                       {b, a} and {c, a}
                       {b, a} and {a, c} (the correct one)
   */
  {
    long int a, b, c, d;

    // Get one halfedges of the first and second edge.
    typename MeshType::Halfedge_index h1 = mesh.halfedge(*fe_it, 0);
    ++fe_it2;
    typename MeshType::Halfedge_index h2 = mesh.halfedge(*fe_it2, 0);

    // Get the vertex numbers connected to the two halfedges.
    a = m_VertexData[mesh.source(h1)].num;
    b = m_VertexData[mesh.target(h1)].num;
    c = m_VertexData[mesh.source(h2)].num;
    d = m_VertexData[mesh.target(h2)].num;

    // First segment wrongly oriented,
    // second correctly oriented.
    if(a == c) 
    {
      m_EdgeData[*fe_it2].tag = 1;
      m_EdgeData[*fe_it].tag = -1;
    }
    // The two segments are wrongly oriented.
    else if(a == d) 
    {
      m_EdgeData[*fe_it2].tag = -1;
      m_EdgeData[*fe_it].tag = -1;
    }
    // The first segment is correctly oriented,
    // the second is wrongly oriented.
    else if(b == d) 
    {
      m_EdgeData[*fe_it2].tag = -1;
      m_EdgeData[*fe_it].tag = 1;
    }
    // Else: nothing to do: the segments are correctly oriented.
    else 
    {
      m_EdgeData[*fe_it2].tag = 1;
      m_EdgeData[*fe_it].tag = 1;
    }
  }

  /** Looking for the other segments.

      As the previous segment is correctly oriented (by above),
      it remains two possibilities:

        {b, a} and {c, a}
        {b, a} and {a, c} (the correct one)
  */
  long int a, b;
  CGAL::Edge_around_face_iterator<MeshType> fe_it3 = fe_it2; 
  ++fe_it3; // Third edge.

  while(fe_it3 != fe_ite) 
  {
    // Get halfedges of the preceeding edge and the current one.
    typename MeshType::Halfedge_index h1 = mesh.halfedge(*fe_it2, 0);
    int orientation = m_EdgeData[*fe_it2].tag; // Orientation of the preceeding edge.

    typename MeshType::Halfedge_index h2 = mesh.halfedge(*fe_it3, 0);

    // Get the last vertex number of the preceeding edge
    // and the first vertex number of the current edge.
    a = ( (orientation == 1) ? m_VertexData[mesh.target(h1)].num
                             : m_VertexData[mesh.source(h1)].num );

    b = m_VertexData[mesh.source(h2)].num;

    if(a != b) // Wrong orientation.
      m_EdgeData[*fe_it3].tag = -1;
    else // Correct orientation.
      m_EdgeData[*fe_it3].tag = 1;

    // Next edge.
    fe_it2 = fe_it3;
    ++fe_it3;
  }

  /* Example of general reordering:
       Given points number 2, 9, 8, 4.

       Line(4) =  {9, 2}; -> Line(-4) = {2, 9};
       Line(10) = {9, 8}; -> Line(10) = {9, 8};
       Line(11) = {8, 4}; -> Line(11) = {8, 4};
       Line(5)  = {2, 4}; -> Line(-5) = {4, 2};

       Exemple of reordering with the last element to be reordered:
       Given points numbers 1, 2, 3 and lines:

       Line(1) = {1, 2}; -> Line(1) = {1, 2};
       Line(2) = {2, 3}; -> Line(2) = {2, 3};
       Line(3) = {1, 3}; -> Line(-3) = {3, 1};

       Line(3) must be reordered by {3, 1} to match the
       orientation of Line(1) and Line(2).
   */
}


} // end namespace marinov

#endif // GMSH_WRITER_HXX