~ubuntu-branches/ubuntu/trusty/rheolef/trusty

///
/// 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
///
/// =========================================================================
//
// visualization for vector fields
//  - with velocity arrow  (fluid)
//  - with deformed meshes (elasticity)
//
// author: Pierre.Saramito@imag.fr
//
// date: 12 february 2001
//
#include "rheolef/field.h"
#include "rheolef/geo-visu-aux.h"
using namespace std;
namespace rheolef {

Float maxi(Float a,Float b) { return (a>b)?a:b; }
int
field::gnuplot2d_velocity (const string& basename, plot_options& opt) const
{
  char filename [1000];
  sprintf (filename, "%s.plot", basename.c_str());
  ofstream plot (filename);
  space V0 = space(_V[0]);
  space V1 = space(_V[1]);	// TODO: can be better: avoid memory
  field u0 = get_comp(V0, 0);
  field u1 = get_comp(V1, 1);	// TODO: can be better: avoid memory 

  const geo& g = get_geo();
  if (opt.verbose) clog << "! file \"" << filename << "\" created.\n";
  plot << "#!gnuplot\n";
  if (!opt.bare)
  	plot << "set title \"" << basename << "\"" << endl;
  else
  	plot << "set nokey" << endl;
  if (opt.empty) 
   {
   	plot << "set noborder" << endl;
   	plot << "set noxtics" << endl;
   	plot << "set noytics" << endl;
   }
  if (opt.xmax==opt.xmin) 
   { 
   	opt.xmin=g.xmin(); 
	opt.xmax=g.xmax(); 
	point diag = opt.xmax-opt.xmin;
	Float w=opt.border*::sqrt(sqr(diag[0])+sqr(diag[1]));
	Float al=maxi(0,-u0.min());
	Float ar=maxi(0,u0.max());
	Float ab=maxi(0,-u1.min());
	Float at=maxi(0,u1.max());
	opt.xmax[0]+=w+ar*opt.vscale;
	opt.xmax[1]+=w+at*opt.vscale;
	opt.xmin[0]-=w+al*opt.vscale;
	opt.xmin[1]-=w+ab*opt.vscale;
	if (opt.xmirror) opt.xmin[0]=-opt.xmax[0];
	if (opt.ymirror) opt.xmin[1]=-opt.xmax[1];
   }
  if (opt.xmin[0]==opt.xmax[0])
   { 
	point diag = opt.xmax-opt.xmin;
	Float w=.05*::sqrt(sqr(diag[0])+sqr(diag[1]));
	Float al=maxi(0,-u0.min());
	Float ar=maxi(0,u0.max());
	opt.xmax[0]+=w+ar*opt.vscale;
	opt.xmin[0]-=w+al*opt.vscale;
   }
  if (opt.xmin[1]==opt.xmax[1])
   { 
	point diag = opt.xmax-opt.xmin;
	Float w=.05*::sqrt(sqr(diag[0])+sqr(diag[1]));
	Float al=maxi(0,-u1.min());
	Float ar=maxi(0,u1.max());
	opt.xmax[1]+=w+ar*opt.vscale;
	opt.xmin[1]-=w+al*opt.vscale;
   }
  if (opt.equal_scales)
   {
   	plot << "set size ratio -1" << endl;
   }

  if (opt.postscript && opt.noxplot)  
   {
	plot << "set output '" << basename << ".eps'" << endl;
	g.write_gnuplot_postscript_options(plot,opt);
   }
  ofstream gdat;
  string coma=",";
  plot << "vscale=" << opt.vscale << endl;
  if (opt.meshplot) plot << "load \"" << basename << ".mesh.plot\"" << endl << "replot "; 
  else 
   {
     if (!opt.keep_user_range) 
      {
     	plot << "set xrange ["<<opt.xmin[0]<<":"<<opt.xmax[0]<<"]" << endl;
     	plot << "set yrange ["<<opt.xmin[1]<<":"<<opt.xmax[1]<<"]" << endl;
      }
     plot << "plot ";
   }
  plot << "\"" << basename 
  		<< ".gdat\" using 1:2:(vscale*$4):(vscale*$5) title \"" 
		<< basename << "\" with vectors linetype 3 \\\n";
  if (opt.xmirror) plot << ", \"" << basename 
		<< ".gdat\" using (-$1):2:(-vscale*$4):(vscale*$5) notitle with vectors linetype 3 \\\n";
  if (opt.ymirror) plot << ", \"" << basename 
  	<< ".gdat\" using 1:(-$2):(vscale*$4):(-vscale*$5) notitle with vectors linetype 3 \\\n";
  if (opt.xmirror && opt.ymirror) plot << ", \"" << basename 
  	<< ".gdat\" using (-$1):(-$2):(-vscale*$4):(-vscale*$5) notitle with vectors linetype 3 \\\n";
  plot << endl;
   
  sprintf (filename, "%s.gdat", basename.c_str());
  gdat.open (filename);
  if (opt.verbose) clog << "! file \"" << filename << "\" created.\n";
  //
  // plot vectors
  //
  tiny_vector<size_type> dof0;
  tiny_vector<size_type> dof1;
  geo::const_iterator last_K = g.end();
  vector<bool> done(V0.size(), false);
  for (geo::const_iterator iter_K = g.begin(); iter_K != last_K; iter_K++) {
      const geo_element& K = *iter_K;
      V0.set_dof(K, dof0);
      V1.set_dof(K, dof1);
      for (geo::size_type i = 0; i < dof0.size(); i++) {
	  if (done[dof0(i)]) {
	     continue;
	  }
          point x_i = V0.x_dof(K, i);
          Float u0_i = u0.at(dof0(i));
          Float u1_i = u1.at(dof1(i));
          gdat << x_i << "\t" << u0_i << " " << u1_i << endl;
          done[dof0(i)] = true;
      }
  }
  //
  // end of plot
  //
  plot << endl; // because of the trailing "\"
  plot << endl; 
  if (opt.pause) plot << "pause -1 \"<return>\"\n";
  if (opt.postscript && (!opt.noxplot)) 
   {
	plot << "set output '" << basename << ".eps'" << endl;
   	g.write_gnuplot_postscript_options(plot,opt);
	plot << "replot";
   }
  plot.close();
  //
  // run gnuplot
  //
  if (opt.meshplot) 
   {
	plot_options optm(opt);
	optm.execute=false; 
	optm.postscript=false;
	optm.bare=true; 
	optm.clean=false;
	g.gnuplot2d (basename+".mesh", optm);
   }
  int status = 0;
  char command [1000];
  if (opt.execute) {
      sprintf (command, "gnuplot -persist %s.plot", basename.c_str());
      if (opt.verbose) clog << "! " << command << endl;
      status = system (command);
  }
  //
  // clear gnuplot data
  //
  if (opt.clean) {
      sprintf (command, "/bin/rm -f %s.plot %s.gdat", basename.c_str(), basename.c_str());
      if (opt.verbose) clog << "! " << command << endl;
      status |= system (command);
      if (opt.meshplot)
      {
      	sprintf (command, "/bin/rm -f %s.mesh.plot %s.mesh.gdat", basename.c_str(), basename.c_str());
      	if (opt.verbose) clog << "! " << command << endl;
      	status |= system (command);
      }
  }
  return status;
}

int
field::plotmtv_velocity (const string& basename, 
	bool execute, bool clean, bool verbose,
	const Float& vscale, bool put_domains) const
{
  string mtv_name = basename + ".mtv";
  ofstream mtv (mtv_name.c_str());
  if (verbose) clog << "! file \"" << mtv_name << "\" created.\n";
  int digits10 = numeric_limits<Float>::digits10;
  mtv << setprecision(digits10);
  //
  // header
  //
  mtv << "#!plotmtv\n"
      << "$ DATA = VECTOR\n"
      << "% comment    = \"rheolef\"\n"
      << "% toplabel   = \"" << basename << "\"\n"
      << "% subtitle   = \"\"\n"
      << "% equalscale = on\n"
      << "% vscale     = " << vscale << "\n"
      ;
  //
  // data
  //
  const geo& g = get_geo();
  space V0 = space(_V[0]);
  space V1 = space(_V[1]);	// TODO: can be better: avoid memory
  space V2;
  if (g.dimension() == 3) {
        V2 = _V[2];
  }
  field u0 = get_comp(V0, 0);
  field u1 = get_comp(V1, 1);	// TODO: can be better: avoid memory 
  field u2;
  if (g.dimension() == 3) {
        u2 = get_comp(V2, 2);
  }
  tiny_vector<size_type> dof0;
  tiny_vector<size_type> dof1;
  tiny_vector<size_type> dof2;
  geo::const_iterator last_K = g.end();
  vector<bool> done(V0.size(), false);
  for (geo::const_iterator iter_K = g.begin(); iter_K != last_K; iter_K++) {
      const geo_element& K = *iter_K;
      V0.set_dof(K, dof0);
      V1.set_dof(K, dof1);
      if (g.dimension() == 3) {
          V2.set_dof(K, dof2);
      }
      for (geo::size_type i = 0; i < dof0.size(); i++) {
	  if (done[dof0(i)]) {
	     continue;
	  }
          point x_i = V0.x_dof(K, i);
          Float u0_i = u0.at(dof0(i));
          Float u1_i = u1.at(dof1(i));
          Float u2_i = 0;
          if (g.dimension() == 3) {
                u2_i = u2.at(dof2(i));
	  }
          mtv << x_i << "\t" << u0_i << " " << u1_i << " " << u2_i << endl;
          done[dof0(i)] = true;
      }
  }
  //
  // plot domains if required
  //
  if (put_domains) 
   {
      mtv << "$ DATA = CURVE2D" << endl;
      g.put_mtv_domains(mtv,1);
   }
  mtv.close();
  //
  // run plotmtv
  //
  int status = 0;
  if (execute) {
      string option;
      if (g.dimension() == 3) {
 	 option = "-3d ";
      }
      if (put_domains) option += "-plotall ";
      string command;
      if (verbose) {
	  command = "plotmtv " + option + basename + ".mtv 1>&2";
          clog << "! " << command << endl;
      } else {
	  command = "plotmtv " + option + basename + ".mtv >/dev/null";
      }
      status = system (command.c_str());
  }
  //
  // clear plotmtv data
  //
  if (clean) {
      string command = "/bin/rm -f " + basename + ".mtv";
      if (verbose) clog << "! " << command << endl;
      status |= system (command.c_str());
  }
  return status;
}

ostream&
field::put_vtk_vector (ostream& vtk, const string& name, bool put_header) const
{
  const geo& g = get_geo();
  space V0 = space(_V[0]);
  space V1 = space(_V[1]);	// TODO: can be better: avoid memory
  space V2;
  if (g.dimension() == 3) {
        V2 = space(_V[2]);
  }
  field u0 = get_comp(V0, 0);
  field u1 = get_comp(V1, 1);	// TODO: can be better: avoid memory 
  field u2;
  //
  // put norm(u) as required scalar field for vector color
  //
  field norm_u;
  if (g.dimension() == 3) {
        u2 = get_comp(V2, 2);
	norm_u = sqrt(u0*u0+u1*u1+u2*u2);
  } else {
	norm_u = sqrt(u0*u0+u1*u1);
  }
  norm_u.put_vtk(vtk, name+"_norm", put_header);
  //
  // put vector data
  //
  vtk << "VECTORS " << name << " float" << endl;
  for (size_type i = 0; i < u0.size(); i++) {
      Float u0_i = u0.at(i);
      Float u1_i = u1.at(i);
      Float u2_i = 0;
      if (g.dimension() == 3) {
            u2_i = u2.at(i);
      }
      vtk << u0_i << " " << u1_i << " " << u2_i << endl;
  }
  vtk << endl;
  return vtk;
}
ostream&
field::put_vtk_tensor (ostream& vtk, const string& name, bool put_header) const
{
  const geo& g = get_geo();
  //
  // put norm(tau) as required scalar field for tensor color
  //
  field norm_tau;
  switch (g.dimension()) {
  case 1:
	norm_tau = abs(operator()(0,0)); 
	break;
  case 2:
	norm_tau = sqrt(
	            sqr(operator()(0,0)) 
		  + sqr(operator()(0,1)) 
		  + sqr(operator()(1,1)) 
		  + sqr(operator()(1,0)) 
	         );
	break;
  default:
	norm_tau = sqrt(
	            sqr(operator()(0,0)) 
		  + sqr(operator()(0,1)) 
		  + sqr(operator()(1,1)) 
		  + sqr(operator()(1,0)) 
		  + sqr(operator()(2,0)) 
		  + sqr(operator()(2,1)) 
		  + sqr(operator()(2,2)) 
		  + sqr(operator()(0,2)) 
		  + sqr(operator()(1,2)) 
	         );
	break;
  }
  norm_tau.put_vtk(vtk, name+"_norm", put_header);
  //
  // put tensors data
  //
  vtk << "TENSORS " << name << " float" << endl;
  switch (g.dimension()) {
  case 1: {
    field t00 = operator()(0,0);
    for (size_type i = 0; i < norm_tau.size(); i++) {
      vtk << t00.at(i) << " 0 0" << endl
          << "0 0 0" << endl
          << "0 0 0" << endl;
    }
    break;
  }
  case 2: {
    field t00 = operator()(0,0);
    field t01 = operator()(0,1);
    field t11 = operator()(1,1);
    field t10 = operator()(1,0);
    for (size_type i = 0; i < norm_tau.size(); i++) {
      vtk << t00.at(i) << " " << t01.at(i) << " 0" << endl
          << t10.at(i) << " " << t11.at(i) << " 0" << endl
          << "0 0 0" << endl;
    }
    break;
  }
  default: {
    field t00 = operator()(0,0);
    field t01 = operator()(0,1);
    field t11 = operator()(1,1);
    field t10 = operator()(1,0);
    field t02 = operator()(0,2);
    field t12 = operator()(1,2);
    field t22 = operator()(2,2);
    field t20 = operator()(2,0);
    field t21 = operator()(2,1);
    for (size_type i = 0; i < norm_tau.size(); i++) {
      vtk << t00.at(i) << " " << t01.at(i) << " " << t02.at(i) << endl
          << t10.at(i) << " " << t11.at(i) << " " << t12.at(i) << endl
          << t20.at(i) << " " << t21.at(i) << " " << t22.at(i) << endl;
    }
    break;
  }}
  vtk << endl;
  return vtk;
}

int
field::vtk_velocity (const string& basename, 
	bool execute, bool clean, bool verbose,
	const Float& vscale) const
{
  string vtk_name = basename + ".vtk";
  ofstream vtk (vtk_name.c_str());
  if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
  int digits10 = numeric_limits<Float>::digits10;
  vtk << setprecision(digits10);
  //
  // header
  //
  //
  // output geometry and data
  //
  const geo& g = get_geo();
  if (get_approx() == "P1") {
     vtk << vtkdata << get_geo();
  } else if (get_approx() == "P2") {
     vtk << lattice << vtkdata << get_geo();
  } else {
     fatal_macro ("vtk_vector: P1 or P2 approximation expected");
  }
  put_vtk_vector (vtk);
  vtk.close();
  //
  // output tcl script
  //
  string tcl_name = basename + ".tcl";
  ofstream tcl (tcl_name.c_str());
  if (verbose) clog << "! file \"" << tcl_name << "\" created.\n";
  tcl << setprecision(digits10);
 
  tcl << "# this is a tcl version of " << basename << ".vtk visualization\n"
      << "source \"" << _RHEOLEF_PKGDATADIR << "/vtk_interactor.tcl\"\n"
      << "#\n"
      << "# edit parameters:\n"
      << "#\n"
      << "set vscale          " << vscale << ";\n"
      << "set axis_color       \"0. 0. 0.\";\n"
      << "set outline_color    \"0. 0. 0.\";\n"
      << "set background_color \"1. 1. 1.\";\n"
      << "#\n"
      << "# Create rendering stuff\n"
      << "#\n"
      << "vtkRenderer ren1;\n"
      << "vtkRenderWindow render_window;\n"
      << "    render_window AddRenderer ren1;\n"
      << "vtkRenderWindowInteractor iren;\n"
      << "    iren SetRenderWindow render_window;\n"
      << "#\n"
      << "# load data\n"
      << "#\n"
      << "vtkUnstructuredGridReader reader;\n"
      << "    reader SetFileName \"" << basename << ".vtk\";\n"
      << "    reader Update;\n"
      << "#\n"
      << "# vector actor\n"
      << "#\n"
      << "vtkHedgeHog vector_filter;\n"
      << "    vector_filter SetInput [reader GetOutput];\n"
      << "    vector_filter SetScaleFactor $vscale;\n"
      << "\n"
      << "vtkLookupTable lookup_table;\n"
      << "    lookup_table Build;\n"
      << "\n"
      << "vtkPolyDataMapper vector_mapper;\n"
      << "    vector_mapper SetInput [vector_filter GetOutput];\n"
      << "    vector_mapper SetLookupTable lookup_table;\n"
      << "    vector_mapper ImmediateModeRenderingOn;\n"
      << "    eval vector_mapper SetScalarRange [[reader GetOutput] GetScalarRange];\n"
      << "\n"
      << "vtkActor vector_actor;\n"
      << "    vector_actor SetMapper vector_mapper;\n"
      << "    # eval [vector_actor GetProperty] SetColor 1 0.49 0.25;\n"
      << "\n"
      << "ren1 AddActor vector_actor;\n"
      << "#\n"
      << "# outline actor\n"
      << "#\n"
      << "vtkOutlineFilter outline_filter;\n"
      << "    outline_filter SetInput [reader GetOutput];\n"
      << "vtkPolyDataMapper outline_mapper;\n"
      << "    outline_mapper SetInput [outline_filter GetOutput];\n"
      << "vtkActor outline_actor;\n"
      << "    outline_actor SetMapper outline_mapper;\n"
      << "    eval [outline_actor GetProperty] SetColor $outline_color\n"
      << "set outline_prop [outline_actor GetProperty];\n"
      << "\n"
      << "ren1 AddActor outline_actor\n"
      << "#\n"
      << "# axis actor\n"
      << "#\n"
      << "vtkCubeAxesActor2D axis_actor;\n"
      << "    axis_actor SetInput [reader GetOutput];\n"
      << "    axis_actor SetCamera [ren1 GetActiveCamera];\n"
      << "    axis_actor SetLabelFormat \"%6.4g\";\n"
      << "    if {[get_vtk_major_version] >= 4 && [get_vtk_minor_version] > 0} {\n"
      << "        vtkTextProperty axis_actor_property\n"
      << "        axis_actor_property ShadowOn\n"
      << "        axis_actor SetAxisLabelTextProperty axis_actor_property\n"
      << "        axis_actor SetAxisTitleTextProperty axis_actor_property\n"
      << "    } else {\n"
      << "        axis_actor ShadowOn;\n"
      << "    }\n"
      << "    axis_actor SetFlyModeToOuterEdges;\n"
      << "    axis_actor SetFontFactor 0.8;\n"
      << "    eval [axis_actor GetProperty] SetColor $axis_color;\n"
      << "\n"
      << "ren1 AddProp  axis_actor\n"
      << "#\n"
      << "# color bar\n"
      << "#\n"
      << "vtkScalarBarActor bar_actor\n"
      << "    bar_actor SetLookupTable [vector_mapper GetLookupTable]\n"
      << "    bar_actor SetOrientationToHorizontal\n"
      << "    [bar_actor GetProperty] SetColor 0 0 0\n"
      << "    [bar_actor GetPositionCoordinate] SetCoordinateSystemToNormalizedViewport\n"
      << "    [bar_actor GetPositionCoordinate] SetValue 0.1 0.01\n"
      << "    bar_actor SetOrientationToHorizontal\n"
      << "    bar_actor SetWidth  0.8\n"
      << "    bar_actor SetHeight 0.10\n"
      << "    bar_actor SetTitle \"velocity module\"\n"
      << "    if {[get_vtk_major_version] >= 4 && [get_vtk_minor_version] > 0} {\n"
      << "        vtkTextProperty bar_text_property\n"
      << "        bar_text_property BoldOn\n"
      << "        bar_actor SetLabelTextProperty bar_text_property\n"
      << "        bar_actor SetTitleTextProperty bar_text_property\n"
      << "    } else {\n"
      << "        bar_actor BoldOn\n"
      << "    }\n"
      << "    eval [bar_actor GetProperty] SetColor $axis_color;\n"
      << "\n"
      << "ren1 AddActor bar_actor;\n"
      << "#\n"
      << "# render the image\n"
      << "#\n"
      << "render_window SetSize 500 500;\n"
      << "eval ren1 SetBackground $background_color;\n"
      << "if { [get_vtk_major_version] >= 4 && [get_vtk_minor_version] > 0} {\n"
      << "    iren AddObserver UserEvent {wm deiconify .vtkInteract};\n"
      << "} else {\n"
      << "    iren SetUserMethod {wm deiconify .vtkInteract};\n"
      << "}\n"
      << "render_window Render;\n"
      << "#\n"
      << "# prevent the tk window from showing up then start the event loop\n"
      << "#\n"
      << "wm withdraw .\n"
    ;
  tcl.close();
  //
  // run vtk
  //
  int status = 0;
  string command = "vtk " + tcl_name;
  if (execute) {
      if (verbose) clog << "! " << command << endl;
      status = system (command.c_str());
  }
  //
  // clear vtk data
  //
  command = "/bin/rm -f " + vtk_name + " " + tcl_name;
  if (clean) {
      if (verbose) clog << "! " << command << endl;
      status |= system (command.c_str());
  }
  return status;
}

int
field::plotmtv_deformation_P1 (const string& basename, 
	bool execute, bool clean, bool verbose,
	const Float& vscale) const
{
  if (get_approx() != "P1") {
      error_macro("field: " << get_approx()
                            << " plotmtv deformation output not implemented.");
  }
  string mtv_name = basename + ".mtv";
  ofstream mtv (mtv_name.c_str());
  if (verbose) clog << "! file \"" << mtv_name << "\" created.\n";
  int digits10 = numeric_limits<Float>::digits10;
  mtv << setprecision(digits10);
  //
  // header
  //
  const geo& g = get_geo();
  mtv << "#!plotmtv\n"
      << "$ DATA = CURVE" << g.dimension() << "D\n"
      << "% comment    = \"rheolef\"\n"
      << "% toplabel   = \"" << basename << "\"\n"
      << "% subtitle   = \"\"\n"
      << "% equalscale = on\n"
      ;
  //
  // plot elements
  //
  geo::const_iterator iter_e = g.begin();
  geo::const_iterator last_e = g.end();
  while (iter_e != last_e)
      plotmtv_element (mtv, *iter_e++, g.begin_node(), g.dimension());
  //
  // data
  //
  space V0 = space(_V[0]);
  space V1 = space(_V[1]);	// TODO: can be better: avoid memory
  space V2;
  if (g.dimension() == 3) {
        V2 = _V[2];
  }
  field u0 = get_comp(V0, 0);
  field u1 = get_comp(V1, 1);	// TODO: can be better: avoid memory 
  field u2;
  if (g.dimension() == 3) {
        u2 = get_comp(V2, 2);
  }
  tiny_vector<size_type> dof0;
  tiny_vector<size_type> dof1;
  tiny_vector<size_type> dof2;
  geo::const_iterator last_K = g.end();
  vector<bool>  done    (V0.size(), false);
  vector<point> new_node(V0.size());
  for (geo::const_iterator iter_K = g.begin(); iter_K != last_K; iter_K++) {
      const geo_element& K = *iter_K;
      V0.set_dof(K, dof0);
      V1.set_dof(K, dof1);
      if (g.dimension() == 3) {
          V2.set_dof(K, dof2);
      }
      for (geo::size_type i = 0; i < dof0.size(); i++) {
	  if (done[dof0(i)]) {
	     continue;
	  }
          point x_i = V0.x_dof(K, i);
          Float u0_i = u0.at(dof0(i));
          Float u1_i = u1.at(dof1(i));
          Float u2_i = 0;
          if (g.dimension() == 3) {
                u2_i = u2.at(dof2(i));
	  }
          point deplacement (u0_i, u1_i, u2_i);
          new_node [dof0(i)] = x_i + vscale*deplacement;
          done[dof0(i)] = true;
      }
  }
  //
  // plot deplaced mesh (P1 only)
  //
  iter_e = g.begin();
  last_e = g.end();
  while (iter_e != last_e)
      plotmtv_element (mtv, *iter_e++, new_node.begin(), g.dimension(), 2, 2);

  mtv.close();
  //
  // run plotmtv
  //
  int status = 0;
  if (execute) {
      string option;
      if (g.dimension() == 3) {
 	 option = "-3d ";
      }
      string command;
      if (verbose) {
	  command = "plotmtv " + option + basename + ".mtv 1>&2";
          clog << "! " << command << endl;
      } else {
	  command = "plotmtv " + option + basename + ".mtv >/dev/null";
      }
      status = system (command.c_str());
  }
  //
  // clear plotmtv data
  //
  if (clean) {
      string command = "/bin/rm -f " + basename + ".mtv";
      if (verbose) clog << "! " << command << endl;
      status |= system (command.c_str());
  }
  return status;
}

int
field::vtk_deformation (const string& basename, 
	bool execute, bool clean, bool verbose,
	const Float& vscale) const
{
  string vtk_name = basename + ".vtk";
  ofstream vtk (vtk_name.c_str());
  if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
  int digits10 = numeric_limits<Float>::digits10;
  vtk << setprecision(digits10);
  //
  // header
  //
  //
  // output geometry and data
  //
  const geo& g = get_geo();
  if (get_approx() == "P1") {
     vtk << vtkdata << get_geo();
  } else if (get_approx() == "P2") {
     vtk << lattice << vtkdata << get_geo();
  } else {
     fatal_macro ("vtk_vector: P1 or P2 approximation expected");
  }
  put_vtk_vector (vtk);
  vtk.close();
  //
  // output tcl script
  //
  string tcl_name = basename + ".tcl";
  ofstream tcl (tcl_name.c_str());
  if (verbose) clog << "! file \"" << tcl_name << "\" created.\n";
  tcl << setprecision(digits10);
  tcl << "# this is a tcl version of cube.vtk visualization\n"
      << "source \"" << _RHEOLEF_PKGDATADIR << "/vtk_interactor.tcl\"\n"
      << "#\n"
      << "# edit parameters:\n"
      << "#\n"
      << "set vscale           " << vscale << ";\n"
      << "set axis_color       \"0. 0. 0.\";\n"
      << "set outline_color    \"0. 0. 0.\";\n"
      << "set background_color \"1. 1. 1.\";\n"
      << "#\n"
      << "# Create rendering stuff\n"
      << "#\n"
      << "vtkRenderer ren1\n"
      << "vtkRenderWindow render_window\n"
      << "    render_window AddRenderer ren1\n"
      << "vtkRenderWindowInteractor iren\n"
      << "    iren SetRenderWindow render_window\n"
      << "#\n"
      << "#\n"
      << "# load data\n"
      << "#\n"
      << "vtkUnstructuredGridReader reader;\n"
      << "    reader SetFileName \"" << basename << ".vtk\";\n"
      << "    reader Update;\n"
      << "\n"
      << "vtkLookupTable lookup_table;\n"
      << "    lookup_table Build;\n"
      << "\n"
      << "\n"
      << "vtkWarpVector translate_filter\n"
      << "    translate_filter SetInput [reader GetOutput]\n"
      << "    translate_filter SetScaleFactor $vscale\n"
      << "vtkVectorNorm norm_filter\n"
      << "    norm_filter SetInput [translate_filter GetOutput]\n"
      << "vtkDataSetMapper deformation_mapper\n"
      << "    deformation_mapper SetInput [translate_filter GetOutput]\n"
      << "    deformation_mapper SetInput [norm_filter GetOutput]\n"
      << "    deformation_mapper SetLookupTable lookup_table;\n"
      << "    deformation_mapper ImmediateModeRenderingOn;\n"
      << "    eval deformation_mapper SetScalarRange [[reader GetOutput] GetScalarRange]\n"
      << "vtkActor deformation_actor\n"
      << "    deformation_actor SetMapper deformation_mapper\n"
      << "\n"
      << "ren1 AddActor deformation_actor\n"
      << "#\n"
      << "# create the outline\n"
      << "#\n"
      << "vtkOutlineFilter outline_filter\n"
      << "    outline_filter SetInput [reader GetOutput]\n"
      << "vtkPolyDataMapper outline_mapper\n"
      << "    outline_mapper SetInput [outline_filter GetOutput]\n"
      << "vtkActor outline_actor\n"
      << "    outline_actor SetMapper outline_mapper\n"
      << "    eval [outline_actor GetProperty] SetColor $outline_color\n"
      << "\n"
      << "ren1 AddActor outline_actor\n"
      << "#\n"
      << "# axis actor\n"
      << "#\n"
      << "vtkCubeAxesActor2D axis_actor;\n"
      << "    axis_actor SetInput [reader GetOutput];\n"
      << "    axis_actor SetCamera [ren1 GetActiveCamera];\n"
      << "    axis_actor SetLabelFormat \"%6.4g\";\n"
      << "    if {[get_vtk_major_version] >= 4 && [get_vtk_minor_version] > 0} {\n"
      << "        vtkTextProperty axis_actor_property\n"
      << "        axis_actor_property ShadowOn\n"
      << "        axis_actor SetAxisLabelTextProperty axis_actor_property\n"
      << "        axis_actor SetAxisTitleTextProperty axis_actor_property\n"
      << "    } else {\n"
      << "        axis_actor ShadowOn;\n"
      << "    }\n"
      << "    axis_actor SetFlyModeToOuterEdges;\n"
      << "    axis_actor SetFontFactor 0.8;\n"
      << "    eval [axis_actor GetProperty] SetColor $axis_color;\n"
      << "\n"
      << "ren1 AddProp  axis_actor\n"
      << "#\n"
      << "# color bar\n"
      << "#\n"
      << "vtkScalarBarActor bar_actor\n"
      << "    bar_actor SetLookupTable [deformation_mapper GetLookupTable]\n"
      << "    bar_actor SetOrientationToHorizontal\n"
      << "    [bar_actor GetProperty] SetColor 0 0 0\n"
      << "    [bar_actor GetPositionCoordinate] SetCoordinateSystemToNormalizedViewport\n"
      << "    [bar_actor GetPositionCoordinate] SetValue 0.1 0.01\n"
      << "    bar_actor SetOrientationToHorizontal\n"
      << "    bar_actor SetWidth  0.8\n"
      << "    bar_actor SetHeight 0.10\n"
      << "    bar_actor SetTitle \"deformation module\"\n"

      //
      // Version 4.2 of VTK introduces several modifications
      // on a couple of classes. Now for instance, fonts are
      // managed with a single class. In order to stay compatible
      // with older and newer versions, we use the code below
      //

      << "    if {[get_vtk_major_version] >= 4 && [get_vtk_minor_version] > 0} {\n"
      << "        vtkTextProperty bar_text_property\n"
      << "        bar_text_property BoldOn\n"
      << "        bar_actor SetLabelTextProperty bar_text_property\n"
      << "        bar_actor SetTitleTextProperty bar_text_property\n"
      << "    } else {\n"
      << "        bar_actor BoldOn\n"
      << "    }\n"
      << "    eval [bar_actor GetProperty] SetColor $axis_color;\n"
      << "\n"
      << "ren1 AddActor bar_actor;\n"
      << "#\n"
      << "# render the image\n"
      << "#\n"
      << "render_window SetSize 500 500;\n"
      << "eval ren1 SetBackground $background_color;\n"
      << "if { [get_vtk_major_version] >= 4 && [get_vtk_minor_version] > 0} {\n"
      << "    iren AddObserver UserEvent {wm deiconify .vtkInteract};\n"
      << "} else {\n"
      << "    iren SetUserMethod {wm deiconify .vtkInteract};\n"
      << "}\n"
      << "render_window Render;\n"
      << "#\n"
      << "# prevent the tk window from showing up then start the event loop\n"
      << "#\n"
      << "wm withdraw .\n"
     ;
  tcl.close();
  //
  // run vtk
  //
  int status = 0;
  string command = "vtk " + tcl_name;
  if (execute) {
      if (verbose) clog << "! " << command << endl;
      status = system (command.c_str());
  }
  //
  // clear vtk data
  //
  command = "/bin/rm -f " + vtk_name + " " + tcl_name;
  if (clean) {
      if (verbose) clog << "! " << command << endl;
      status |= system (command.c_str());
  }
  return status;
}
} // namespace rheolef