// file automatically generated by: ../../../rheolef/nfem/pbasis/make_basis_list_cxx.sh
///
/// 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
///
/// =========================================================================
#include "basis.h"
#include "equispaced.icc"
#include "basis_visu_gnuplot.icc"
#include "basis_on_pointset_evaluate.icc"
#include "rheolef/iorheo.h"
namespace rheolef {
using namespace std;

// --------------------------------------------------------------------------
// show polynomials
// --------------------------------------------------------------------------
template<class T>
void
basis_rep<T>::put (ostream& os, reference_element hat_K) const
{
  switch (valued_tag()) {
    case space_constant::scalar: put_scalar_valued (os, hat_K); break;
    case space_constant::vector: put_vector_valued (os, hat_K); break;
    default: error_macro ("unimplemented gnuplot visualization of "<<valued()<<"-valued polynomials");
  }
}
template<class T>
void
basis_rep<T>::put_scalar_valued (ostream& os, reference_element hat_K) const
{
   details::put(*this,os,hat_K); // shared with basis_raw
}
template<class T>
void
basis_rep<T>::put_vector_valued (ostream& os, reference_element hat_K) const
{
  typedef reference_element::size_type  size_type;
  bool verbose  = iorheo::getverbose(os);
  bool clean    = iorheo::getclean(os);
  bool execute  = iorheo::getexecute(os);
  size_type nsub = iorheo::getsubdivide(os);
  if (nsub <= 1) nsub = 20; // default value

  size_type loc_ndof = ndof(hat_K);
  string basename  = "basis-" + name() + "-" + hat_K.name();
  string filelist;
  // --------------------
  // .plot
  // --------------------
  string plot_name = basename + ".plot";
  string gdat_name = basename + ".gdat";
  ofstream plot (plot_name.c_str());
  cerr << "! file \"" << plot_name << "\" created" << endl;
  filelist += " \"" + plot_name + "\"";
  size_t d = hat_K.dimension();
  check_macro (d==1 || d==2, "unsupported dimension " << d);
  plot << "gdat = \"" << gdat_name << "\"" << endl
       << "set colors classic" << endl
       << "set size square" << endl
    ;
  // draw all edges
  point_basic<T> xmin, xmax;
  for (size_type mu = 0; mu < d; ++mu) {
    xmin[mu] =  1;
    xmax[mu] = -1;
  }
  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(1); loc_isid < loc_nsid; ++loc_isid) {
    size_type loc_ia = hat_K.subgeo_local_vertex (1, loc_isid, 0);
    size_type loc_ib = hat_K.subgeo_local_vertex (1, loc_isid, 1);
    const point& a = hat_K.vertex (loc_ia);
    const point& b = hat_K.vertex (loc_ib);
    plot << "set arrow from ";
    size_type d2 = max(d,size_type(2));
    for (size_type mu = 0; mu < d2; ++mu) {
      plot << a[mu];
      if (mu+1 != d2) plot << ", ";
      xmin[mu] = min(xmin[mu],T(a[mu]));
      xmax[mu] = max(xmax[mu],T(a[mu]));
    }
    plot << " to ";
    for (size_type mu = 0; mu < d2; ++mu) {
      plot << b[mu];
      if (mu+1 != d2) plot << ", ";
      xmin[mu] = min(xmin[mu],T(b[mu]));
      xmax[mu] = max(xmax[mu],T(b[mu]));
    }
    plot << " nohead lc 0 lw 1" << endl;
  }
  T delta = 0.1;
  plot << endl
       <<   "set xrange ["<<xmin[0]-delta<<":"<<xmax[0]+delta<<"]"<<endl;
  if (d >= 2) {
    plot << "set yrange ["<<xmin[1]-delta<<":"<<xmax[1]+delta<<"]"<<endl;
  } else {
    plot << "set yrange [-1:1]"<<endl;
  }
  if (d >= 3) {
    plot << "set zrange ["<<xmin[2]-delta<<":"<<xmax[2]+delta<<"]"<<endl
         << "set xyplane at "<< xmin[2]-delta << endl;
  }
  if (d == 1) {
    plot << "plot \\" << endl
      ;
    for (size_t loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
      plot << "  gdat u 1:"<<loc_idof+2 << " t \"L" << loc_idof+1 << "\" w l";
      if (loc_idof+1 != loc_ndof) { plot << ",\\"; }
      plot << endl;
    }
  } else if (d == 2) {
    plot << "pause_duration = 0.5" << endl
         << "d    = 2" << endl
         << "scale= 1./" << nsub << endl
         << "idx  = 0" << endl
         << "idxmax = " << loc_ndof << endl
         << "plot gdat i idx u 1:2:(scale*$3):(scale*$4) t sprintf(\"P%d\",idx+1) w vec lc 0 lw 2" << endl
         << "load \"" << basename << ".loop\"" << endl
      ;
    string loop_name = basename + ".loop";
    ofstream loop (loop_name.c_str());
    cerr << "! file \"" << loop_name << "\" created" << endl;
    filelist += " \"" + loop_name + "\"";
    loop << "idx = idx+1" << endl
         << "if (idx < idxmax) \\" << endl
         << "  pause -1 \"<return>\"; \\" << endl
         << "  replot; \\" << endl
         << "  reread" << endl
      ;
  }
  plot << "pause -1 \"<return>\"" << endl;
  plot.close();

  // --------------------
  // .gdat
  // --------------------
  ofstream gdat (gdat_name.c_str());
  cerr << "! file \"" << gdat_name << "\" created" << endl;
  filelist += " \"" + gdat_name + "\"";
  gdat << setprecision(std::numeric_limits<T>::digits10)
       << "# basis   " << name() << endl
       << "# element " << hat_K.name() << endl
       << "# degree  " << degree() << endl
    ;
  switch (hat_K.variant()) {
#ifdef TODO
    case reference_element::p: {
      break;
    } 
    case reference_element::e: {
      gdat << "# ndof " << loc_ndof << endl
           << "# x";
      for (size_t loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
        gdat << " L" << loc_idof+1;
      }
      gdat << endl;
      for (size_type i = 0; i <= nsub; i++) { 
        point_basic<T> hat_x (T(i)/T(nsub));
        eval (hat_K, hat_x, value);
	gdat << hat_x[0];
	for (size_t loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
	  for (size_t mu = 0; mu < d; ++mu) {
	    gdat << " " << value[loc_idof][mu];
          }
        }
        gdat << endl;
      }
      break;
    }
#endif // TODO
    case reference_element::t:
    case reference_element::q: {
      Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>              hat_node;
      Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic> value;
      pointset_lagrange_equispaced (hat_K, nsub, hat_node); 
      details::basis_on_pointset_evaluate (*this, hat_K, hat_node, value);
      size_t loc_nnod = value.rows();
      size_t loc_ndof = value.cols();
      for (size_t loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
        for (size_t loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
	  for (size_t mu = 0; mu < d; ++mu) {
	    gdat << hat_node[loc_inod][mu] << " ";
          }
	  for (size_t mu = 0; mu < d; ++mu) {
	    gdat << value(loc_inod,loc_idof)[mu] << " ";
          }
          gdat << endl;
        }
        gdat << endl << endl;
      }
      break;
    }
#ifdef TODO
    case reference_element::T: {
      for (size_type k = 0; k <= nsub; k++) {
        for (size_type j = 0; j+k <= nsub; j++) {
          for (size_type i = 0; i+j+k <= nsub; i++) { 
            point_basic<T> hat_x (T(i)/T(nsub), T(j)/T(nsub), T(k)/T(nsub));
            eval_lagrange (hat_x, hat_K, degree, sopt, inv_vdm, value);
            Lambda = max(Lambda, norm(value,1));
          }
        }
      }
      return Lambda;
    }
    case reference_element::P: {
      for (size_type k = 0; k <= degree; k++) { 
        for (size_type j = 0; j <= degree; j++) { 
          for (size_type i = 0; i+j <= degree; i++) { 
            point_basic<T> hat_x (T(i)/T(nsub), T(j)/T(nsub), 2*T(k)/T(nsub)-1);
            eval_lagrange (hat_x, hat_K, degree, sopt, inv_vdm, value);
            Lambda = max(Lambda, norm(value,1));
          }
        }
      }
      break;
    } 
    case reference_element::H: {
      for (size_type k = 0; k <= degree; k++) { 
        for (size_type j = 0; j <= degree; j++) { 
          for (size_type i = 0; i <= degree; i++) { 
            point_basic<T> hat_x (2*T(i)/T(nsub)-1, 2*T(j)/T(nsub), 2*T(k)/T(nsub)-1);
            eval_lagrange (hat_x, hat_K, degree, sopt, inv_vdm, value);
            Lambda = max(Lambda, norm(value,1));
          }
        }
      }
      break;
    } 
#endif // TODO
    default:
      error_macro ("unexpected element type `"<<hat_K.name()<<"'");
  } 
  // -----------
  if (execute) {
  // -----------
    string command = "gnuplot \"" + plot_name + "\"";
    if (verbose) clog << "! " << command << endl;
    int status = system (command.c_str());
  }
  // -----------
  if (clean) {
  // -----------
    string command = "/bin/rm -f " + filelist;
    if (verbose) clog << "! " << command << endl;
    int status = system (command.c_str());
  }
}
// --------------------------------------------------------------------------
// show nodes
// --------------------------------------------------------------------------
template<class T>
static
void
put_hat_node_header (std::ostream& plot, reference_element hat_K) 
{
  typedef reference_element::size_type size_type;
  size_type d = hat_K.dimension();
  plot << "set colors classic" << endl
       << "set nokey" << endl
       << "set noborder" << endl
       << "set noxtics" << endl
       << "set noytics" << endl
       << "set noztics" << endl
    ;
  if (d == 2){
    plot << "set size ratio -1 # 2d equal scales" << endl;
  } else if (d == 3) {
    plot << "set view 70, 140   # x to left, y to right, z to top" << endl;
    plot << "set view equal xyz # 3d equal scales" << endl;
  }
  plot << endl;

  // draw all edges
  point_basic<T> xmin, xmax;
  for (size_type mu = 0; mu < d; ++mu) {
    xmin[mu] =  1;
    xmax[mu] = -1;
  }
  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(1); loc_isid < loc_nsid; ++loc_isid) {
    size_type loc_ia = hat_K.subgeo_local_vertex (1, loc_isid, 0);
    size_type loc_ib = hat_K.subgeo_local_vertex (1, loc_isid, 1);
    const point& a = hat_K.vertex (loc_ia);
    const point& b = hat_K.vertex (loc_ib);
    plot << "set arrow from ";
    size_type d2 = max(d,size_type(2));
    for (size_type mu = 0; mu < d2; ++mu) {
      plot << a[mu];
      if (mu+1 != d2) plot << ", ";
      xmin[mu] = min(xmin[mu],T(a[mu]));
      xmax[mu] = max(xmax[mu],T(a[mu]));
    }
    plot << " to ";
    for (size_type mu = 0; mu < d2; ++mu) {
      plot << b[mu];
      if (mu+1 != d2) plot << ", ";
      xmin[mu] = min(xmin[mu],T(b[mu]));
      xmax[mu] = max(xmax[mu],T(b[mu]));
    }
    plot << " nohead lc 0 lw 1" << endl;
  }
  T delta = 0.1;
  plot << endl
       <<   "set xrange ["<<xmin[0]-delta<<":"<<xmax[0]+delta<<"]"<<endl;
  if (d >= 2) {
    plot << "set yrange ["<<xmin[1]-delta<<":"<<xmax[1]+delta<<"]"<<endl;
  } else {
    plot << "set yrange [-1:1]"<<endl;
  }
  if (d >= 3) {
    plot << "set zrange ["<<xmin[2]-delta<<":"<<xmax[2]+delta<<"]"<<endl
         << "set xyplane at "<< xmin[2]-delta << endl;
  }
}
// class-function for shrinkage of discontinuous element internal nodes
template<class T>
struct shrink_fun {
  point_basic<T> operator() (const point_basic<T>& x) const {
    return is_continuous ? x : x - a*(x-xc);
  }
  shrink_fun (bool ic, const point_basic<T>& xc1) 
  : is_continuous(ic), xc(xc1), a(0.15) {}
  bool is_continuous;
  point_basic<T> xc;
  T a;
};
template<class T>
void
basis_rep<T>::put_hat_node (std::ostream& os, reference_element hat_K) const
{
  using namespace std;
  bool verbose  = iorheo::getverbose(os);
  bool clean    = iorheo::getclean(os);
  bool execute  = iorheo::getexecute(os);

  _initialize_data_guard (hat_K);
  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_x = hat_node(hat_K);
  size_type d = hat_K.dimension();

  string basename  = "basis-" + name() + "-" + hat_K.name();
  string filelist;
  // --------------------
  // .gdat
  // --------------------
  string gdat_name = basename + ".gdat";
  ofstream gdat (gdat_name.c_str());
  cerr << "! file \"" << gdat_name << "\" created" << endl;
  filelist += " \"" + gdat_name + "\"";
  gdat << setprecision(std::numeric_limits<T>::digits10)
       << "# basis    " << name() << endl
       << "# element  " << hat_K.name() << endl
       << "# degree   " << degree() << endl
       << "# ndof     " << ndof(hat_K) << endl
       << "# raw_poly " << _sopt.get_raw_polynomial_name() << endl
       << "# nodeset  " << _sopt.get_node_name() << endl
       << "# nnode    " << hat_x.size() << endl
    ;
  point_basic<T> xc;
  reference_element_barycenter (hat_K, xc);
  shrink_fun<T> shrink (option().is_continuous(), xc);
  if (! option().is_trace_n()) {
    for (size_type dim = 0, d = hat_K.dimension(); dim <= d; ++dim) {
      for (size_type loc_inod = first_inod_by_dimension(hat_K,dim), loc_nnod = first_inod_by_dimension(hat_K,dim+1); loc_inod < loc_nnod; ++loc_inod) {
        point_basic<T> shrink_hat_x = shrink (hat_x [loc_inod]);
        shrink_hat_x.put (gdat, max(d,size_type(2)));
        gdat << " " << loc_inod << endl;
      }
      gdat << endl << endl;
    }
  } else {
    // TODO: loop on sides and shrink also inside each side
    size_type loc_inod = 0, loc_nnod = nnod(hat_K);
    size_type d = hat_K.dimension();
    for (size_type loc_isid = 0, loc_nsid = hat_K.n_side(); loc_isid < loc_nsid; ++loc_isid) {
      reference_element S = hat_K.side (loc_isid);
      // 1) compute the side barycenter
      point xs;
      size_type loc_nsidvert = hat_K.subgeo_size (d-1, loc_isid);
      for (size_type loc_jsidvert = 0; loc_jsidvert < loc_nsidvert; ++loc_jsidvert) {
        size_type iloc = hat_K.subgeo_local_vertex(d-1, loc_isid, loc_jsidvert);
        xs += hat_K.vertex (iloc);
      }
      xs /= T(loc_nsidvert);
      // 2) loop on side and shrink
      shrink_fun<T> shrink_s (option().is_continuous(), xs);
      size_type nns = nnod_on_subgeo_internal (d, S.variant());
      for (size_type ins = 0; ins < nns; ++ins, ++loc_inod) { 
        check_macro (loc_inod < loc_nnod, "invalid loc_inod");
        point_basic<T> shrink_hat_x = shrink_s (shrink (hat_x [loc_inod]));
        shrink_hat_x.put (gdat, max(d,size_type(2)));
        gdat << " " << loc_inod << endl;
      }
    }
  }
  // --------------------
  // .plot
  // --------------------
  string plot_name = basename + ".plot";
  ofstream plot (plot_name.c_str());
  cerr << "! file \"" << plot_name << "\" created" << endl;
  filelist += " \"" + plot_name + "\"";
  plot << "gdat = \"" << gdat_name << "\"" << endl
       << "set title \"" << name() << ": nodes\"" << endl
       ;
  put_hat_node_header<T> (plot, hat_K) ;
  string color[4] = {"#ff0000", "#008800", "#0000ff", "#ff00ff"};
  bool need_coma = false;
  size_type index = 0;
  plot << ((d==3) ? "s" : "") << "plot \\" << endl;
  if (! option().is_trace_n()) {
    for (size_type dim = 0, d = hat_K.dimension(); dim <= d; ++dim) {
      size_type nnod_dim = first_inod_by_dimension (hat_K, dim+1) - first_inod_by_dimension (hat_K, dim);
      if (nnod_dim == 0) continue;
      plot << " " << (need_coma ? ", " : " ") << "gdat \\" << endl
           << "  index " << index << " \\" << endl
           << "  w p ps 2 pt 7 lc rgb \"" << color[dim] << "\",\\" << endl
           << " gdat \\" << endl
           << "  index " << index << " \\" << endl
           << "  u 1:2:3" << ((d==3) ? ":4" : "") << " w labels offset char 2 lc 0 \\" << endl
        ;
      need_coma = true;
      index++;
    }
  } else {
    plot << " gdat \\" << endl
         << "  w p ps 2 pt 7 lc rgb \"" << color[0] << "\",\\" << endl
         << " gdat \\" << endl
         << "  u 1:2:3" << ((d==3) ? ":4" : "") << " w labels offset char 2 lc 0" << endl
        ;
  }
  plot << endl
       << "pause -1 \"<return>\"" << endl
    ;
  // -----------
  if (execute) {
  // -----------
    string command = "gnuplot \"" + plot_name + "\"";
    if (verbose) clog << "! " << command << endl;
    int status = system (command.c_str());
  }
  // -----------
  if (clean) {
  // -----------
    string command = "/bin/rm -f " + filelist;
    if (verbose) clog << "! " << command << endl;
    int status = system (command.c_str());
  }
}
// --------------------------------------------------------------------------
// show nodes
// --------------------------------------------------------------------------
template<class T>
void
basis_rep<T>::put_hat_node_on_side (
  std::ostream&                os,
  reference_element            hat_K,
  const side_information_type& sid) const
{
  using namespace std;
  check_macro (is_nodal(), "side-node visualization: only nodal basis supported");
  Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_inod;
  local_idof_on_side (hat_K, sid, loc_inod);

  bool verbose  = iorheo::getverbose(os);
  bool clean    = iorheo::getclean(os);
  bool execute  = iorheo::getexecute(os);

  _initialize_data_guard (hat_K);
  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_x = hat_node (hat_K);
  size_type d = hat_K.dimension();

  string basename  = "basis-" + name() + "-" + hat_K.name();
  string filelist;
  // --------------------
  // .gdat
  // --------------------
  string gdat_name = basename + ".gdat";
  ofstream gdat (gdat_name.c_str());
  cerr << "! file \"" << gdat_name << "\" created" << endl;
  filelist += " \"" + gdat_name + "\"";
  gdat << setprecision(std::numeric_limits<T>::digits10)
       << "# basis    " << name() << endl
       << "# element  " << hat_K.name() << endl
       << "# degree   " << degree() << endl
       << "# ndof     " << ndof(hat_K) << endl
       << "# raw_poly " << _sopt.get_raw_polynomial_name() << endl
       << "# nodeset  " << _sopt.get_node_name() << endl
       << "# nnode    " << hat_x.size() << endl
    ;
  for (size_type loc_sid_inod = 0, loc_sid_nnod = loc_inod.size(); loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
    hat_x [loc_inod[loc_sid_inod]].put (gdat, max(d,size_type(2)));
    gdat << " " << loc_inod[loc_sid_inod] << endl;
  }
  // --------------------
  // .plot
  // --------------------
  string plot_name = basename + ".plot";
  ofstream plot (plot_name.c_str());
  cerr << "! file \"" << plot_name << "\" created" << endl;
  filelist += " \"" + plot_name + "\"";
  plot << "gdat = \"" << gdat_name << "\"" << endl
       << "set title \"" << name() << ": nodes\"" << endl
       ;
  put_hat_node_header<T> (plot, hat_K) ;
  string color[4] = {"#ff0000", "#008800", "#0000ff", "#ff00ff"};
  bool need_coma = false;
  size_type index = 0;
  plot << ((d==3) ? "s" : "") << "plot \\" << endl
       << "  gdat \\" << endl
       << "    w p ps 2 pt 7 lc rgb \"" << color[d] << "\",\\" << endl
       << "  gdat \\" << endl
       << "    u 1:2:3" << ((d==3) ? ":4" : "") << " w labels offset char 2 lc 0" << endl
       << "pause -1 \"<return>\"" << endl
    ;
  // -----------
  if (execute) {
  // -----------
    string command = "gnuplot \"" + plot_name + "\"";
    if (verbose) clog << "! " << command << endl;
    int status = system (command.c_str());
  }
  // -----------
  if (clean) {
  // -----------
    string command = "/bin/rm -f " + filelist;
    if (verbose) clog << "! " << command << endl;
    int status = system (command.c_str());
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T)                                               \
template void basis_rep<Float>::put (ostream& os, reference_element hat_K) const; \
template void basis_rep<Float>::put_scalar_valued (ostream& os, reference_element hat_K) const; \
template void basis_rep<Float>::put_vector_valued (ostream& os, reference_element hat_K) const; \
template void basis_rep<Float>::put_hat_node (ostream& os, reference_element hat_K) const; \
template void basis_rep<Float>::put_hat_node_on_side (				\
  std::ostream&, reference_element, const side_information_type&) const;

_RHEOLEF_instanciation(Float)

} // namespace rheolef