field_seq_visu_vtk_paraview.cc
Go to the documentation of this file.
//
// paraview vtk visualization
//
// author: Pierre.Saramito@imag.fr
//
// date: 12 may 1997 update: 23 oct 2011
//
#include "rheolef/field_expr.h"
#include "rheolef/piola_util.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/iofem.h"
#include "rheolef/interpolate.h"
#include "rheolef/render_option.h"
using namespace std;
namespace rheolef {
// ----------------------------------------------------------------------------
// field puts
// ----------------------------------------------------------------------------
// extern:
template <class T> odiststream& field_put_vtk (odiststream&, const field_basic<T,sequential>&);
template <class T>
odiststream&
visu_vtk_paraview (odiststream& ops, const field_basic<T,sequential>& uh)
{
//
// 0) prerequises
//
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 = ops.os();
//
// 1) set render options
//
render_option popt;
popt.valued = uh.get_space().valued();
popt.fill = iorheo::getfill(os); // isocontours or color fill
popt.elevation = iorheo::getelevation(os);
popt.color = iorheo::getcolor(os);
popt.gray = iorheo::getgray(os);
popt.black_and_white = iorheo::getblack_and_white(os);
popt.showlabel = iorheo::getshowlabel(os);
popt.stereo = iorheo::getstereo(os);
popt.volume = iorheo::getvolume(os);
popt.iso = true;
popt.cut = iorheo::getcut(os);
popt.grid = iorheo::getgrid(os);
popt.format = iorheo::getimage_format(os);
popt.mark = iorheo::getmark(os);
bool is_scalar = (popt.valued == "scalar");
if (popt.mark == "") popt.mark = "scalar";
if (popt.mark != "" && !is_scalar) popt.mark = "|"+popt.mark +"|";
popt.style = (popt.valued == "vector") ? (iorheo::getvelocity(os) ? "velocity" : "deformation") : "none";
popt.scale = iorheo::getvectorscale(os);
popt.origin = iofem::getorigin(os);
popt.normal = iofem::getnormal(os);
popt.resolution = iofem::getresolution(os);
popt.n_isovalue = iorheo::getn_isovalue(os);
popt.n_isovalue_negative = iorheo::getn_isovalue_negative(os);
popt.isovalue = iorheo::getisovalue(os); // isovalue is always a Float
popt.label = iorheo::getlabel(os);
const geo_basic<float_type,sequential>& omega = uh.get_geo();
size_type dim = omega.dimension();
size_type map_dim = omega.map_dimension();
popt.view_2d = (dim == 2);
popt.view_map = (dim > map_dim);
#if (_RHEOLEF_PARAVIEW_VERSION_MAJOR >= 5) && (_RHEOLEF_PARAVIEW_VERSION_MINOR >= 5)
// paraview version >= 5.5 has high order elements
popt.high_order = (omega.order() > 1 || uh.get_space().degree() > 1);
#else
popt.high_order = false;
#endif
#if (_RHEOLEF_PARAVIEW_VERSION_MAJOR == 5) && (_RHEOLEF_PARAVIEW_VERSION_MINOR == 7)
// paraview version == 5.7 has opacity bug
popt.have_opacity_bug = true;
#else
popt.have_opacity_bug = false;
#endif
popt.xmin = uh.get_geo().xmin();
popt.xmax = uh.get_geo().xmax();
field_basic<T,sequential> uh_scalar;
if (popt.valued == "scalar") {
uh_scalar = uh;
} else {
size_type k = uh.get_space().degree();
space_basic<T,sequential> T0h (uh.get_geo(), "P"+itos(k));
uh_scalar = interpolate(T0h, norm(uh));
}
// TODO: u_min=min(uh.u.min,uh.b.min) instead of interpolate norm ?
popt.f_min = uh_scalar.min();
popt.f_max = uh_scalar.max();
//
// 2) output data
//
bool verbose = iorheo::getverbose(os);
bool clean = iorheo::getclean(os);
bool execute = iorheo::getexecute(os);
string basename = iorheo::getbasename(os);
string outfile_fmt = "";
string tmp = get_tmpdir() + "/";
if (!clean) tmp = "";
string filelist;
string filename = tmp+basename + ".vtk";
filelist = filelist + " " + filename;
ofstream vtk_os (filename.c_str());
odiststream vtk (vtk_os);
if (verbose) clog << "! file \"" << filename << "\" created.\n";
field_put_vtk (vtk, uh);
vtk.close();
//
// 3) create python data file
//
std::string py_name = filename = tmp+basename + ".py";
filelist = filelist + " " + filename;
ofstream py (filename.c_str());
if (verbose) clog << "! file \"" << filename << "\" created.\n";
point xmin = uh.get_geo().xmin(),
xmax = uh.get_geo().xmax();
py << popt
<< endl
<< "paraview_field_" << popt.valued << "(paraview, \"" << tmp+basename << "\", opt)" << endl
<< endl
;
py.close();
//
// 4) run pyton
//
int status = 0;
string command;
if (execute) {
string prog = (popt.format == "") ? "paraview --script=" : "pvbatch --use-offscreen-rendering ";
command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " " + prog + py_name;
if (popt.format != "") command = "DISPLAY=:0.0 " + command;
if (popt.stereo && popt.format == "") command = command + " --stereo";
if (verbose) clog << "! " << command << endl;
status = system (command.c_str());
}
//
// 4) clear vtk data
//
if (clean) {
command = "/bin/rm -f " + filelist;
if (verbose) clog << "! " << command << endl;
status = system (command.c_str());
}
return ops;
}
// --------------------------------------------------------------------
// plane cut : via vtk-paraview script
// --------------------------------------------------------------------
template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
template <class T>
field_basic<T,sequential>
paraview_plane_cut (
const field_basic<T,sequential>& uh,
const point_basic<T>& origin,
const point_basic<T>& normal)
{
//
// 1) prerequises
//
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 typename geo_basic<float_type,sequential>::node_type node_type;
typedef point_basic<size_type> ilat;
ostream& os = std::cout;
bool verbose = iorheo::getverbose(os);
bool clean = iorheo::getclean(os);
bool execute = iorheo::getexecute(os);
string basename = iorheo::getbasename(os);
string tmp = get_tmpdir() + "/";
if (!clean) tmp = "";
//
// 2) output data
//
string filelist;
string vtk_name = tmp+basename + ".vtk";
filelist = filelist + " " + vtk_name;
ofstream vtk_os (vtk_name.c_str());
odiststream vtk (vtk_os);
if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
field_put_vtk (vtk, uh);
vtk.close();
//
// create python script
//
//
// 3) create python data file
//
string py_name = tmp+basename + "-cut.py";
string vtk_cut_name = tmp+basename + "-cut.vtk";
filelist = filelist + " " + py_name + " " + vtk_cut_name;
ofstream py (py_name.c_str());
if (verbose) clog << "! file \"" << py_name << "\" created.\n";
py << setprecision(numeric_limits<T>::digits10)
<< "from paraview.simple import *" << endl
<< "reader = LegacyVTKReader(FileNames=['" << vtk_name << "'])" << endl
<< "slice = Slice(SliceType=\"Plane\")" << endl
<< "slice.SliceOffsetValues = [0.0]" << endl
<< "slice.SliceType.Origin = " << render_option::python(origin) << endl
<< "slice.SliceType.Normal = " << render_option::python(normal) << endl
<< "writer = paraview.simple.CreateWriter(\"" << vtk_cut_name << "\", slice)" << endl
<< "writer.FileType = 'Ascii'" << endl
<< "writer.UpdatePipeline()" << endl
;
py.close();
//
// 3) run pyton
//
int status = 0;
string command;
command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " pvbatch " + py_name;
if (verbose) clog << "! " << command << endl;
status = system (command.c_str());
//
// 4) load vtk cutted mesh
//
ifstream vtk_polydata (vtk_cut_name.c_str());
check_macro (vtk_polydata, "field: vtk polydata file \"" << vtk_cut_name << "\" not found.");
if (verbose) clog << "! load `" << vtk_cut_name << "'." << endl;
idiststream ips_vtk_polydata (vtk_polydata);
geo_basic<T,sequential> gamma_cut;
geo_get_vtk (ips_vtk_polydata, gamma_cut);
gamma_cut.set_name (basename + "-cut");
check_macro (gamma_cut.n_node() > 0, "empty mesh & plane intersection: HINT check normal and origin");
//
// clean
//
if (clean) {
command = "/bin/rm -f " + filelist;
if (verbose) clog << "! " << command << endl;
status = system (command.c_str());
}
//
// project geometry from 2d/3d on plane in 1d/2d
//
// 1D) x := dot(OM, tangent)
//
// 2D) x1 := dot(OM, t1)
// x2 := dot(OM, t2)
//
bool use_projection = true;
if (use_projection) {
switch (uh.get_geo().dimension()) {
case 2: {
// project in 1D space
gamma_cut.set_dimension (1);
point_basic<T> tangent = point_basic<T>(normal[1], -normal[0]);
disarray<node_type,sequential> node = gamma_cut.get_nodes();
for (size_type i = 0, n = node.size(); i < n; i++) {
node[i] = point_basic<T> (dot(node[i], tangent), 0);
}
gamma_cut.set_nodes (node);
break;
}
case 3: {
gamma_cut.set_dimension (2);
point_basic<T> t1 = point_basic<T>(1, 0, 0);
if (norm(vect(t1, normal)) == Float(0)) t1 = point_basic<T>(0, 1, 0);
t1 = t1 - dot(t1,normal)*normal;
t1 = t1 / norm(t1);
point_basic<T> t2 = vect(normal,t1);
disarray<node_type,sequential> node = gamma_cut.get_nodes();
for (size_type i = 0, n = node.size(); i < n; i++) {
node[i] = point_basic<T> (dot(node[i] - origin, t1), dot(node[i] - origin, t2), 0);
}
gamma_cut.set_nodes (node);
break;
}
default: error_macro ("unexpected dimension="<< uh.get_geo().dimension());
}
}
// -------------------------------
// get field values
// polydata header:
// POINT_DATA <n_node>
// SCALARS scalars float
// LOOKUP_TABLE default
// -------------------------------
string data_type;
while (true) {
vtk_polydata >> data_type;
if ((data_type == "POINT_DATA") ||
(data_type == "CELL_DATA" )) {
break;
}
}
scatch (vtk_polydata, "default");
string approx = (data_type == "POINT_DATA") ? "P1" : "P0";
space_basic<T,sequential> V_cut (gamma_cut, approx);
field_basic<T,sequential> u_cut (V_cut);
u_cut.set_u().get_values (ips_vtk_polydata);
return u_cut;
}
// --------------------------------------------------------------------
// isovalue surface : via vtk-paraview script
// --------------------------------------------------------------------
template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
template <class T>
geo_basic<T,sequential>
paraview_extract_isosurface (const field_basic<T,sequential>& uh)
{
//
// 1) prerequises
//
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 typename geo_basic<float_type,sequential>::node_type node_type;
typedef point_basic<size_type> ilat;
ostream& os = std::cout;
bool verbose = iorheo::getverbose(os);
bool clean = iorheo::getclean(os);
bool execute = iorheo::getexecute(os);
T isovalue = iorheo::getisovalue(os);
string basename = iorheo::getbasename(os);
string tmp = get_tmpdir() + "/";
if (!clean) tmp = "";
//
// 2) output data
//
string filelist;
string vtk_name = tmp+basename + ".vtk";
filelist = filelist + " " + vtk_name;
ofstream vtk_os (vtk_name.c_str());
odiststream vtk (vtk_os);
if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
field_put_vtk (vtk, uh);
vtk.close();
//
// create python script
//
//
// 3) create python data file
//
string py_name = tmp+basename + "-iso.py";
string vtk_iso_name = tmp+basename + "-iso.vtk";
filelist = filelist + " " + py_name + " " + vtk_iso_name;
ofstream py (py_name.c_str());
if (verbose) clog << "! file \"" << py_name << "\" created.\n";
py << setprecision(numeric_limits<T>::digits10)
<< "from paraview.simple import *" << endl
<< "reader = LegacyVTKReader(FileNames=['" << vtk_name << "'])" << endl
<< "iso_surface = Contour(PointMergeMethod=\"Uniform Binning\")" << endl
<< "iso_surface.ContourBy = ['POINTS', 'scalar']" << endl
<< "iso_surface.Isosurfaces = [" << isovalue << "]" << endl
<< "writer = paraview.simple.CreateWriter(\"" << vtk_iso_name << "\", iso_surface)" << endl
<< "writer.FileType = 'Ascii'" << endl
<< "writer.UpdatePipeline()" << endl
;
py.close();
//
// 3) run pyton
//
int status = 0;
string command;
command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " pvbatch " + py_name;
if (verbose) clog << "! " << command << endl;
status = system (command.c_str());
//
// 4) load vtk isosurface mesh
//
ifstream vtk_polydata (vtk_iso_name.c_str());
check_macro (vtk_polydata, "field: vtk polydata file \"" << vtk_iso_name << "\" not found.");
if (verbose) clog << "! load `" << vtk_iso_name << "'." << endl;
idiststream ips_vtk_polydata (vtk_polydata);
geo_basic<T,sequential> gamma_iso;
geo_get_vtk (ips_vtk_polydata, gamma_iso);
gamma_iso.set_name (basename + "-iso");
gamma_iso.set_dimension (uh.get_geo().dimension());
check_macro (gamma_iso.n_node() > 0, "empty mesh & plane intersection: HINT check normal and origin");
//
// clean
//
if (clean) {
command = "/bin/rm -f " + filelist;
if (verbose) clog << "! " << command << endl;
status = system (command.c_str());
}
return gamma_iso;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template odiststream& \
visu_vtk_paraview<Float> ( \
odiststream&, \
const field_basic<Float,sequential>&); \
template field_basic<Float,sequential> \
paraview_plane_cut ( \
const field_basic<Float,sequential>&, \
const point_basic<Float>&, \
const point_basic<Float>&); \
template geo_basic<Float,sequential> \
paraview_extract_isosurface ( \
const field_basic<Float,sequential>&);
_RHEOLEF_instanciation(Float)
#undef _RHEOLEF_instanciation
}// namespace rheolef
rheolef::geo_basic
generic mesh with rerefence counting
Definition: domain_indirect.h:64
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::render_option::f_max
Float f_max
Definition: render_option.h:47
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::render_option::origin
point origin
Definition: render_option.h:49
rheolef::render_option::grid
bool grid
Definition: render_option.h:44
rheolef::normal
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
Definition: field_expr_terminal.h:439
rheolef::render_option::black_and_white
bool black_and_white
Definition: render_option.h:45
rheolef::get_tmpdir
std::string get_tmpdir()
get_tmpdir: see the rheostream page for the full documentation
Definition: rheostream.cc:50
rheolef::geo_get_vtk
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
Definition: geo_seq_get_vtk.cc:39
rheolef::field_basic::float_type
typename float_traits< T >::type float_type
Definition: field.h:245
rheolef::render_option::color
bool color
Definition: render_option.h:45
rheolef::render_option::fill
bool fill
Definition: render_option.h:44
rheolef::python
std::string python(const point_basic< T > &x, size_t d=3)
rheolef::visu_vtk_paraview
odiststream & visu_vtk_paraview(odiststream &, const field_basic< T, sequential > &)
Definition: field_seq_visu_vtk_paraview.cc:47
rheolef::render_option::elevation
bool elevation
Definition: render_option.h:44
rheolef::render_option::f_min
Float f_min
Definition: render_option.h:47
rheolef::odiststream::os
std::ostream & os()
Definition: diststream.h:236
rheolef::render_option::view_map
bool view_map
Definition: render_option.h:45
rheolef::render_option::xmin
point xmin
Definition: render_option.h:49
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::render_option::gray
bool gray
Definition: render_option.h:45
rheolef::paraview_plane_cut
field_basic< T, sequential > paraview_plane_cut(const field_basic< T, sequential > &uh, const point_basic< T > &origin, const point_basic< T > &normal)
Definition: field_seq_visu_vtk_paraview.cc:180
rheolef::geo_basic< T, sequential >::set_nodes
void set_nodes(const disarray< node_type, sequential > &x)
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::geo_basic< T, sequential >::set_dimension
void set_dimension(size_type dim)
rheolef::geo_basic< T, sequential >
Definition: geo.h:1108
rheolef::field_basic::set_u
vec< T, M > & set_u()
Definition: field.h:312
vtk
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format vtk
Definition: iorheo-members.h:113
rheolef::render_option::n_isovalue
size_t n_isovalue
Definition: render_option.h:46
rheolef::field_basic::min
T min() const
Definition: field.h:683
rheolef::render_option::stereo
bool stereo
Definition: render_option.h:44
rheolef::interpolate
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
Definition: interpolate.cc:233
rheolef::geo_basic< T, sequential >::set_name
void set_name(std::string name)
mkgeo_ball.clean
clean
Definition: mkgeo_ball.sh:335
rheolef::render_option::cut
bool cut
Definition: render_option.h:44
rheolef::render_option::mark
std::string mark
Definition: render_option.h:50
rheolef::field_basic::max
T max() const
Definition: field.h:699
rheolef::scatch
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:52
rheolef::field_basic
Definition: field_expr_utilities.h:38
rheolef::render_option::scale
Float scale
Definition: render_option.h:47
rheolef::geo_basic< T, sequential >::get_nodes
const disarray< node_type, sequential > & get_nodes() const
Definition: geo.h:1177
rheolef::render_option::volume
bool volume
Definition: render_option.h:44
rheolef::render_option::high_order
bool high_order
Definition: render_option.h:45
rheolef::render_option::xmax
point xmax
Definition: render_option.h:49
rheolef::render_option::normal
point normal
Definition: render_option.h:49
mkgeo_couette.basename
basename
Definition: mkgeo_couette.sh:73
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::render_option::label
std::string label
Definition: render_option.h:50
rheolef::render_option
Definition: render_option.h:38
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
rheolef::paraview_extract_isosurface
geo_basic< T, sequential > paraview_extract_isosurface(const field_basic< T, sequential > &uh)
Definition: field_seq_visu_vtk_paraview.cc:329
Float
see the Float page for the full documentation
point
see the point page for the full documentation
rheolef::disarray< node_type, sequential >
rheolef::render_option::isovalue
Float isovalue
Definition: render_option.h:47
mkgeo_ball.verbose
verbose
Definition: mkgeo_ball.sh:133
rheolef::render_option::resolution
point_basic< size_t > resolution
Definition: render_option.h:48
rheolef::geo_basic< T, sequential >::n_node
size_type n_node() const
Definition: geo.h:1172
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::render_option::style
std::string style
Definition: render_option.h:50
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T)
Definition: field_seq_visu_vtk_paraview.cc:413
_RHEOLEF_PKGDATADIR
#define _RHEOLEF_PKGDATADIR
Definition: config.h:240
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::render_option::iso
bool iso
Definition: render_option.h:44
rheolef::render_option::have_opacity_bug
bool have_opacity_bug
Definition: render_option.h:45
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
mkgeo_ball.tmp
tmp
Definition: mkgeo_ball.sh:380
rheolef::space_basic< T, sequential >
Definition: space.h:361
mkgeo_ball.dim
dim
Definition: mkgeo_ball.sh:307
rheolef::render_option::showlabel
bool showlabel
Definition: render_option.h:44
rheolef::render_option::n_isovalue_negative
size_t n_isovalue_negative
Definition: render_option.h:46
rheolef::render_option::valued
std::string valued
Definition: render_option.h:50
rheolef::field_put_vtk
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
Definition: field_seq_put_vtk.cc:154
rheolef::std
Definition: vec_expr_v2.h:391
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
mkgeo_ball.command
command
Definition: mkgeo_ball.sh:136
rheolef::details::dot
rheolef::details::is_vec dot
rheolef::render_option::view_2d
bool view_2d
Definition: render_option.h:44
rheolef::vect
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:265
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::render_option::format
std::string format
Definition: render_option.h:50