///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
///
/// =========================================================================
//
// gnuplot geo visualisation
//
// author: Pierre.Saramito@imag.fr
//
// date: 16 sept 2011
//
# include "rheolef/field.h"
# include "rheolef/field_wdof_indirect.h"
# include "rheolef/field_evaluate.h"
# include "rheolef/piola_util.h"
# include "rheolef/rheostream.h"
# include "field_seq_visu_gnuplot_internal.h"
# include "rheolef/rounder.h"
# include "rheolef/interpolate.h"

namespace rheolef {

template <class T>
odiststream& visu_gnuplot (odiststream& ops, const geo_basic<T,sequential>& omega);

// ----------------------------------------------------------------------------
// puts for one element
// ----------------------------------------------------------------------------
template<class T>
static
void
put_edge (
  std::ostream&                    gdat,
  const geo_basic<T,sequential>&   omega,
  const geo_element&               K,
  const field_basic<T,sequential>& uh,
  const fem_on_pointset<T>&        fops,
  size_t                           my_order,
  bound_type<T>&                   bbox)
{
  using namespace std;
  typedef typename geo_basic<T,sequential>::size_type size_type;
  typedef point_basic<size_type>                      ilat;
  size_type dim = omega.dimension();
  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = fops.get_piola_on_pointset().get_piola (omega, K);
  Eigen::Matrix<T,Eigen::Dynamic,1> value;
  field_evaluate (uh, fops, omega, K, value);
  for (size_type loc_idof = 0, loc_ndof = piola.size(); loc_idof < loc_ndof; loc_idof++) {
    bbox.update (piola[loc_idof].F, value[loc_idof]);
  }
  for (size_type i = 0; i <= my_order; i++) {
    size_type loc_idof = reference_element_e::ilat2loc_inod (my_order, ilat(i));
    piola[loc_idof].F.put (gdat, dim); gdat << " " << value[loc_idof] << endl;
  }
  gdat << endl;
  gdat << endl;
}
template<class T>
static
void
put_triangle (
  std::ostream&                    gdat,
  const geo_basic<T,sequential>&   omega,
  const geo_element&               K,
  const field_basic<T,sequential>& uh,
  const fem_on_pointset<T>&        fops,
  size_t                           my_order,
  bound_type<T>&                   bbox)
{
  using namespace std;
  typedef typename geo_basic<T,sequential>::size_type size_type;
  typedef point_basic<size_type>                      ilat;
  size_type dim = omega.dimension();
  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = fops.get_piola_on_pointset().get_piola (omega, K);
  Eigen::Matrix<T,Eigen::Dynamic,1> value;
  field_evaluate (uh, fops, omega, K, value);
  for (size_type loc_idof = 0, loc_ndof = piola.size(); loc_idof < loc_ndof; loc_idof++) {
    bbox.update (piola[loc_idof].F, value[loc_idof]);
  }
  for (size_type j = 0; j <= my_order; j++) {
    for (size_type i1 = 0; i1 <= my_order; i1++) {
      size_type i = std::min(i1, my_order-j);
      size_type iloc = reference_element_t::ilat2loc_inod (my_order, ilat(i, j));
      piola[iloc].F.put (gdat, dim); gdat << " " << value[iloc] << endl;
    }
    gdat << endl;
  }
  gdat << endl << endl;
}
template<class T>
static
void
put_quadrangle (
  std::ostream&                    gdat,
  const geo_basic<T,sequential>&   omega,
  const geo_element&               K,
  const field_basic<T,sequential>& uh,
  const fem_on_pointset<T>&        fops,
  size_t                           my_order,
  bound_type<T>&                   bbox)
{
  using namespace std;
  typedef typename geo_basic<T,sequential>::size_type size_type;
  typedef point_basic<size_type>                      ilat;
  size_type dim = omega.dimension();
  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = fops.get_piola_on_pointset().get_piola (omega, K);
  Eigen::Matrix<T,Eigen::Dynamic,1> value;
  field_evaluate (uh, fops, omega, K, value);
  for (size_type loc_idof = 0, loc_ndof = piola.size(); loc_idof < loc_ndof; loc_idof++) {
    bbox.update (piola[loc_idof].F, value[loc_idof]);
  }
  for (size_type j = 0; j < my_order+1; j++) {
    for (size_type i = 0; i < my_order+1; i++) {
      size_type loc_idof00 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j));
      piola[loc_idof00].F.put (gdat, dim); gdat << " " << value[loc_idof00] << endl;
    }
    gdat << endl;
  }
  gdat << endl << endl;
}
template<class T>
void
put (
  std::ostream&                    gdat,
  const geo_basic<T,sequential>&   omega,
  const geo_element&               K,
  const field_basic<T,sequential>& uh,
  const fem_on_pointset<T>&        fops,
  size_t                           my_order,
  bound_type<T>&                   bbox)
{
  switch (K.variant()) {
   case reference_element::e: put_edge       (gdat, omega, K, uh, fops, my_order, bbox); break;
   case reference_element::t: put_triangle   (gdat, omega, K, uh, fops, my_order, bbox); break;
   case reference_element::q: put_quadrangle (gdat, omega, K, uh, fops, my_order, bbox); break;
   default: error_macro ("unsupported element variant `"<<K.variant()<<"'");
  }
}
// ----------------------------------------------------------------------------
// scalar field puts
// ----------------------------------------------------------------------------
template <class T>
odiststream&
visu_gnuplot_scalar (odiststream& ods, const field_basic<T,sequential>& uh)
{
  using namespace std;
  typedef typename field_basic<T,sequential>::float_type float_type;
  typedef typename geo_basic<float_type,sequential>::size_type size_type;
  typedef point_basic<size_type>                      ilat;
  ostream& os = ods.os();
  bool verbose  = iorheo::getverbose(os);
  bool clean    = iorheo::getclean(os);
  bool execute  = iorheo::getexecute(os);
  bool fill     = iorheo::getfill(os);    // show grid or fill elements
  bool elevation = iorheo::getelevation(os);
  bool color   = iorheo::getcolor(os);
  bool gray    = iorheo::getgray(os);
  bool black_and_white = iorheo::getblack_and_white(os);
  bool reader_on_stdin = iorheo::getreader_on_stdin(os);
  string format   = iorheo::getimage_format(os);
  string basename = iorheo::getbasename(os);
  size_type subdivide = iorheo::getsubdivide(os);
  size_type n_isovalue = iorheo::getn_isovalue(os);
  size_type n_isovalue_negative = iorheo::getn_isovalue_negative(os);
  string outfile_fmt = "";
  string tmp = get_tmpdir() + "/";
  if (!clean) tmp = "";

  const geo_basic<float_type,sequential>& omega = uh.get_geo();
  size_type dim     = omega.dimension();
  size_type map_dim = omega.map_dimension();
  size_type nv      = omega.sizes().ownership_by_dimension[0].size();
  size_type nedg    = omega.sizes().ownership_by_dimension[1].size();
  size_type nfac    = omega.sizes().ownership_by_dimension[2].size();
  size_type nvol    = omega.sizes().ownership_by_dimension[3].size();
  size_type ne      = omega.sizes().ownership_by_dimension[map_dim].size();

  const basis_basic<float_type>& b = uh.get_space().get_basis();
  if (subdivide == 0) { // subdivide is unset: use default
    subdivide = std::max(omega.order(), subdivide);
    subdivide = std::max(b.degree (), subdivide);
  }
  basis_basic<T> subdivide_pointset ("P"+std::to_string(subdivide));
  piola_on_pointset<T> pops; pops.initialize (omega.get_piola_basis(), subdivide_pointset, integrate_option());
  fem_on_pointset<T>   fops; fops.initialize (b, pops);

  bound_type<T> bbox;
  bbox.xmin = omega.xmin();
  bbox.xmax = omega.xmax();
  bbox.umin = uh.min();
  bbox.umax = uh.max();

  std::vector<T> values;
  if (n_isovalue_negative != 0) {
    for (size_t i = 0; i <= n_isovalue_negative; i++) {
      values.push_back (bbox.umin*(n_isovalue_negative - i)/n_isovalue_negative);
    }
    for (size_t i = 1; i <= n_isovalue - n_isovalue_negative; i++) {
      values.push_back (bbox.umax*i/(n_isovalue - n_isovalue_negative));
    }
  }
  string filelist;
  //
  // output .gdat
  //
  string filename = tmp+basename + ".gdat";
  string gdatname = filename;
  filelist = filelist + " " + filename;
  ofstream gdat (filename.c_str());
  if (verbose) clog << "! file \"" << filename << "\" created.\n";
  gdat << setprecision(numeric_limits<float_type>::digits10);
  size_type used_dim = (fill ? map_dim : 1);
  for (size_type ie = 0, ne = omega.size(used_dim); ie < ne; ie++) {
    const geo_element& K = omega.get_geo_element(used_dim,ie);
    put (gdat, omega, K, uh, fops, subdivide, bbox);
  }
  gdat.close();
  //
  // rounding bounds
  //
  T eps = 1e-7;
  bbox.umin = floorer(eps) (bbox.umin);
  bbox.umax = ceiler(eps)  (bbox.umax);
  if (fabs(bbox.umax - bbox.umin) < eps) { // empty range: constant field
    bbox.umax = 1.1*bbox.umax;
    bbox.umin = 0.9*bbox.umin;
  }
  //
  // output .plot
  //
  filename = tmp+basename + ".plot";
  filelist = filelist + " " + filename;
  ofstream plot (filename.c_str());
  if (verbose) clog << "! file \"" << filename << "\" created.\n";

  plot << "#!gnuplot" << endl
       << setprecision(numeric_limits<float_type>::digits10);
  if (format != "") {
    outfile_fmt = basename + "." + format;
    string terminal = format;
    if (terminal == "ps")  {
      terminal = "postscript eps";
      if (color) terminal += " color";
    }
    if (terminal == "jpg") terminal = "jpeg";
    if (terminal == "jpeg" || terminal == "png" || terminal == "gif") {
      terminal += " crop";
    }
    plot << "set terminal " << terminal    << endl
         << "set output \"" << outfile_fmt << "\"" << endl;
  }
  plot << "umin = " << bbox.umin << endl
       << "umax = " << bbox.umax << endl
       << "n_isovalue  = " << n_isovalue << endl
       << "uincr = 1.0*(umax-umin)/n_isovalue" << endl;
  if (dim == 2) {
    if (bbox.xmin[0] >= bbox.xmax[0]) plot << "#";
    plot << "set xrange [" << bbox.xmin[0] << ":" << bbox.xmax[0] << "]" << endl;
    if (bbox.xmin[1] >= bbox.xmax[1]) plot << "#";
    plot << "set yrange [" << bbox.xmin[1] << ":" << bbox.xmax[1] << "]" << endl;
    if (bbox.umin >= bbox.umax) plot << "#";
    plot << "set zrange [umin:umax]" << endl;
    if (bbox.xmin[0] >= bbox.xmax[0] || bbox.xmin[1] >= bbox.xmax[1]) {
      plot << "set size square" << endl;
    } else {
      plot << "set size ratio -1 # equal scales" << endl
           << "set noxtics" << endl
           << "set noytics" << endl;
    }
    if (elevation) {
      plot << "set xyplane at umin-0.1*(umax-umin)" << endl;
    } else {
      plot << "set view map" << endl;
    }
    if (map_dim == 2) {
      plot << "unset surface" << endl
           << "set contour base" << endl;
      if (values.size() != 0) {
        plot << "set cntrparam levels discrete ";
        for (size_t i = 0, n = values.size(); i < n; i++) {
          plot << values[i];
          if (i+1 != n) plot << ", ";
        }
        plot << endl;
      } else {
        plot << "eps = (umax-umin)*1e-4"<< endl;
        plot << "uincr_eps = 1.0*(umax-umin-2*eps)/n_isovalue"<< endl;
        plot << "set cntrparam levels incremental umin+eps,uincr_eps,umax-eps"<< endl;
      }
    }
    if (black_and_white) {
      plot << "set cbtics out scale 0.5" << endl;
    }
    plot << "set cbtics uincr" << endl
         << "set cbrange [umin:umax]" << endl;
    if (gray) {
      plot << "set palette gray" << endl;
    } else if (color) {
      plot << "set palette rgbformulae 33,13,-4" << endl;
    } else { // bw
      plot << "set palette rgbformulae 0,0,0" << endl;
    }
    plot << "set palette maxcolors n_isovalue+1" << endl
         << "set nokey" << endl;
  } else if (dim == 3 && map_dim == 2) {
    // field defined on a 3D surface
    point_basic<T> dx = 0.1*(omega.xmax() - omega.xmin());
    T dx_max = max(dx[0],max(dx[1],dx[2]));
    if (dx_max == 0) dx_max = 0.1;
    dx[0] = max(dx[0],dx_max);
    if (omega.dimension() >= 2) dx[1] = max(dx[1],dx_max);
    if (omega.dimension() == 3) dx[2] = max(dx[2],dx_max);
    point_basic<T> xmin = omega.xmin() - dx;
    point_basic<T> xmax = omega.xmax() + dx;
    // TODO: visu contours discrets comme en 2d classique
    plot << "set xrange [" << xmin[0] << ":" << xmax[0] << "]" << endl
         << "set yrange [" << xmin[1] << ":" << xmax[1] << "]" << endl
         << "set zrange [" << xmin[2] << ":" << xmax[2] << "]" << endl
         << "set xyplane at " << xmin[2] << endl
         << "set view equal xyz # equal scales" << endl
         << "set view 70,120" << endl
         << "n_isovalue  = " << n_isovalue << endl
	 << "n_subdivide = 40 " << endl
         << "uincr = 1.0*(umax-umin)/n_isovalue" << endl
         << "set cbtics uincr" << endl
         << "set pm3d interpolate n_subdivide,n_subdivide corners2color mean" << endl
         << "set palette rgbformulae 33,13,-4 maxcolors n_isovalue" << endl;
    if (format != "") {
      plot << "set noxlabel" << endl
           << "set noylabel" << endl
           << "set nozlabel" << endl;
    } else {
      plot << "set xlabel \"x\"" << endl
           << "set ylabel \"y\"" << endl
           << "set zlabel \"z\"" << endl;
    }
  }
  if (dim == 1) {
    if (color) {
      plot << "set colors classic" << endl;
    }
    plot << "plot \"" << gdatname << "\" notitle with lines lw 2";
    if (gray || black_and_white) {
      plot << " lc 0" << endl;
    }
    plot << endl;
  } else {
    if (!fill && dim == 2) {
      plot << "plot";
    } else {
      plot << "splot";
    } 
    plot << " \"" << gdatname << "\" notitle";
    if (map_dim == 2) {
      if (black_and_white) {
        plot << " with lines palette lw 2" << endl;
      } else {
        plot << " with lines palette" << endl;
      }
    } else { // a 2d line
        plot << " with lines palette lw 2" << endl;
    }
  }
  //
  // end of plot
  //
  if (format == "" && !reader_on_stdin) {
    plot << "pause -1 \"<return>\"\n";
  }
  plot.close();
  //
  // run gnuplot
  //
  int status = 0;
  string command;
  if (execute) {
      command = "gnuplot ";
      if (reader_on_stdin) command += "-persist ";
      command += tmp + basename + ".plot";
      if (verbose) clog << "! " << command << endl;
      cin.sync();
      status = system (command.c_str());
      if (format != "") {
        check_macro (file_exists (outfile_fmt), "! file \"" << outfile_fmt << "\" creation failed");
        if (verbose) clog << "! file \"" << outfile_fmt << "\" created" << endl;
      }
  }
  //
  // clear gnuplot data
  //
  if (clean) {
      command = "/bin/rm -f " + filelist;
      if (verbose) clog << "! " << command << endl;
      status = system (command.c_str());
  }
  return ods;
}
// ----------------------------------------------------------------------------
// vector field puts
// ----------------------------------------------------------------------------
template <class T>
odiststream&
visu_gnuplot_vector (odiststream& ods, const field_basic<T,sequential>& uh)
{
  typedef typename geo_basic<T,sequential>::size_type size_type;
  Float vscale = iorheo::getvectorscale(ods.os());
  const geo_basic<T,sequential>& omega = uh.get_geo();
  // TODO: non-isoparam => interpolate
  check_macro (omega.get_piola_basis().degree() == uh.get_space().get_basis().degree(),
	"gnuplot vector: unsupported non-isoparametric approx " << uh.get_space().get_basis().name());
  size_type n_elt = omega.size();
  size_type d     = omega.dimension();
  T diam_omega = norm (omega.xmax() - omega.xmin()); // characteristic length in omega
  T h_moy = diam_omega/pow(n_elt,1./d);
  space_basic<T,sequential> Xh (uh.get_geo(), "P"+std::to_string(uh.get_space().degree()));
  field_basic<T,sequential> norm_uh = interpolate(Xh, norm(uh));
  T norm_max_uh = norm_uh.max_abs();      // get max vector length
  if (norm_max_uh + 1 == 1) norm_max_uh = 1;
#ifdef TODO
  T scale = vscale*(h_moy/norm_max_uh);
#endif // TODO
  size_type n_comp = uh.get_space().get_basis().size();
  disarray<point_basic<T>, sequential> x = omega.get_nodes();
  for (size_type inod = 0, nnod = x.size(); inod < nnod; inod++) {
    point_basic<T> vi;
    for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
      size_type idof = n_comp*inod + i_comp;
      vi[i_comp] = uh.dof(idof);
    }
    x[inod] += vscale*vi;
  }
  geo_basic<T,sequential> deformed_omega = omega;
  deformed_omega.set_nodes(x);
  space_basic<T,sequential> deformed_Vh (deformed_omega, norm_uh.get_space().get_basis().name());
  field_basic<T,sequential> deformed_norm_uh (deformed_Vh);
  std::copy (norm_uh.begin_dof(), norm_uh.end_dof(), deformed_norm_uh.begin_dof());
  visu_gnuplot (ods, deformed_norm_uh);
  return ods;
}
// ----------------------------------------------------------------------------
// switch
// ----------------------------------------------------------------------------
template <class T>
odiststream&
visu_gnuplot (odiststream& ods, const field_basic<T,sequential>& uh)
{
  switch (uh.get_space().valued_tag()) {
    case space_constant::scalar: visu_gnuplot_scalar (ods, uh); break;
    case space_constant::vector: visu_gnuplot_vector (ods, uh); break;
    default: error_macro ("do not known how to print " << uh.valued() << "-valued field");
  }
  return ods;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciate(T) \
template odiststream& visu_gnuplot<T> (odiststream&, const field_basic<T,sequential>&);

_RHEOLEF_instanciate(Float)
#undef _RHEOLEF_instanciate
} // rheolef namespace