Loading...
Searching...
No Matches
geo_seq_visu_gnuplot.cc
Go to the documentation of this file.
//
// gnuplot geo visualisation
//
// author: Pierre.Saramito@imag.fr
//
// date: 16 sept 2011
//
#include "rheolef/geo.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/reference_element.h"
#include "rheolef/piola_util.h"
namespace rheolef {
// ----------------------------------------------------------------------------
// element puts
// ----------------------------------------------------------------------------
template<class T>
static
void
put_vertex (std::ostream& gdat, const geo_element& P, const geo_basic<T,sequential>& omega)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
size_type dim = max (omega.dimension(), size_type(2));
omega.node (P[0]).put (gdat, dim); gdat << endl;
}
// simple edge, could be curved (order > 1 => n_node = order+1 > 2)
template<class T>
static
void
put_edge (std::ostream& gdat, const geo_element& E, const geo_basic<T,sequential>& omega,
const basis_on_pointset<T>& pointset, size_t subdivide)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
size_type dim = max (omega.dimension(), size_type(2));
std::vector<size_type> dis_inod;
omega.dis_inod (E, dis_inod);
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x;
piola_transformation (omega, pointset, E.variant(), dis_inod, x);
for (size_type i = 0; i <= subdivide; i++) {
size_type loc_inod = reference_element_e::ilat2loc_inod (subdivide, ilat(i));
x[loc_inod].put (gdat, dim); gdat << endl;
}
gdat << endl << endl;
}
// 2d triangle or quadrangle
template<class T>
static
void
put_2d_face (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
const basis_on_pointset<T>& pointset, size_t subdivide)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
size_type dim = max (omega.dimension(), size_type(2));
std::vector<size_type> dis_inod;
omega.dis_inod (F, dis_inod);
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x;
piola_transformation (omega, pointset, F.variant(), dis_inod, x);
if (F.variant() == reference_element::q) {
// 2d quadrangle: draw lines of the lattice
// ICI
for (size_type j = 0; j < subdivide; j++) {
for (size_type i = 0; i < subdivide; i++) {
size_type loc_inod00 = reference_element_q::ilat2loc_inod (subdivide, ilat(i, j));
size_type loc_inod10 = reference_element_q::ilat2loc_inod (subdivide, ilat(i+1, j));
size_type loc_inod01 = reference_element_q::ilat2loc_inod (subdivide, ilat(i, j+1));
size_type loc_inod11 = reference_element_q::ilat2loc_inod (subdivide, ilat(i+1, j+1));
const point_basic<T>& x00 = x[loc_inod00];
const point_basic<T>& x10 = x[loc_inod10];
const point_basic<T>& x01 = x[loc_inod01];
const point_basic<T>& x11 = x[loc_inod11];
if (i+1 == subdivide) {
x11.put (gdat, dim); gdat << endl;
}
x10.put (gdat, dim); gdat << endl;
x00.put (gdat, dim); gdat << endl;
x01.put (gdat, dim); gdat << endl;
if (j+1 == subdivide) {
x11.put (gdat, dim); gdat << endl;
}
gdat << endl;
}
}
gdat << endl;
} else {
// 2d triangle: draw lines of the lattice
for (size_type i = 0; i < subdivide; i++) {
for (size_type j = 0; j < subdivide - i; j++) {
size_type loc_inod00 = reference_element_t::ilat2loc_inod (subdivide, ilat(i, j));
size_type loc_inod10 = reference_element_t::ilat2loc_inod (subdivide, ilat(i+1, j));
size_type loc_inod01 = reference_element_t::ilat2loc_inod (subdivide, ilat(i, j+1));
const point_basic<T>& x00 = x[loc_inod00];
const point_basic<T>& x10 = x[loc_inod10];
const point_basic<T>& x01 = x[loc_inod01];
x10.put (gdat, dim); gdat << endl;
x00.put (gdat, dim); gdat << endl;
x01.put (gdat, dim); gdat << endl;
if (i+j+1 == subdivide) {
x10.put (gdat, dim); gdat << endl;
}
gdat << endl;
}
}
gdat << endl;
}
}
// 3d triangle or quadrangle: splot lattice data, for hidden faces removal
template<class T>
static
void
put_3d_face (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
size_type dim = max (omega.dimension(), size_type(2));
size_type order = omega.order();
std::vector<size_type> inod;
omega.dis_inod (F, inod);
gdat << "# F"<<F.dis_ie()<<endl;
for (size_type j = 0; j <= order; j++) {
for (size_type i = 0; i <= order; i++) {
size_type loc_inod;
if (F.variant() == reference_element::q) { // TODO: virtual in geo_elment_v2::loc_ilat2loc_inod
loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
} else {
loc_inod = reference_element_t::ilat2loc_inod (order, ilat(min(i,order-j), j));
}
omega.node (inod[loc_inod]).put (gdat, dim); gdat << endl;
}
gdat << endl;
}
gdat << endl;
}
template<class T>
static
void
put_face (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
const basis_on_pointset<T>& pointset, size_t subdivide)
{
switch (omega.dimension()) {
case 2: put_2d_face (gdat, F, omega, pointset, subdivide); break;
case 3: put_3d_face (gdat, F, omega); break;
default: break;
}
}
template<class T>
static
void
put_volume (std::ostream& gdat, const geo_element& K, const geo_basic<T,sequential>& omega)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
size_type dim = max (omega.dimension(), size_type(2));
size_type order = omega.order();
std::vector<size_type> inod;
omega.dis_inod (K, inod);
gdat << "# K"<<K.dis_ie()<<endl;
if (K.variant() == reference_element::T) {
// 3d tetra: draw lattice
for (size_type i = 0; i < order; i++) {
for (size_type j = 0; j < order - i; j++) {
for (size_type k = 0; k < order - i - j; k++) {
size_type loc_inod000 = reference_element_T::ilat2loc_inod (order, ilat(i, j, k));
size_type loc_inod100 = reference_element_T::ilat2loc_inod (order, ilat(i+1, j, k));
size_type loc_inod010 = reference_element_T::ilat2loc_inod (order, ilat(i, j+1, k));
size_type loc_inod001 = reference_element_T::ilat2loc_inod (order, ilat(i, j, k+1));
gdat << omega.node (inod[loc_inod100]) << endl
<< omega.node (inod[loc_inod000]) << endl
<< omega.node (inod[loc_inod010]) << endl << endl
<< omega.node (inod[loc_inod000]) << endl
<< omega.node (inod[loc_inod001]) << endl << endl;
if (i+j+k+1 == order) {
gdat << omega.node (inod[loc_inod100]) << endl
<< omega.node (inod[loc_inod010]) << endl
<< omega.node (inod[loc_inod001]) << endl
<< omega.node (inod[loc_inod100]) << endl << endl;
}
}
}
}
} else if (K.variant() == reference_element::H) {
// 3d hexa: draw lattice
for (size_type k = 0; k < order; k++) {
for (size_type j = 0; j < order; j++) {
for (size_type i = 0; i < order; i++) {
size_type loc_inod000 = reference_element_H::ilat2loc_inod (order, ilat(i, j, k));
size_type loc_inod100 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j, k));
size_type loc_inod110 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j+1, k));
size_type loc_inod010 = reference_element_H::ilat2loc_inod (order, ilat(i, j+1, k));
size_type loc_inod001 = reference_element_H::ilat2loc_inod (order, ilat(i, j, k+1));
size_type loc_inod101 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j, k+1));
size_type loc_inod111 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j+1, k+1));
size_type loc_inod011 = reference_element_H::ilat2loc_inod (order, ilat(i, j+1, k+1));
gdat << omega.node (inod[loc_inod100]) << endl
<< omega.node (inod[loc_inod000]) << endl
<< omega.node (inod[loc_inod010]) << endl << endl
<< omega.node (inod[loc_inod000]) << endl
<< omega.node (inod[loc_inod001]) << endl << endl;
if (i+1 == order) {
gdat << omega.node (inod[loc_inod101]) << endl
<< omega.node (inod[loc_inod100]) << endl
<< omega.node (inod[loc_inod110]) << endl << endl;
}
if (j+1 == order) {
gdat << omega.node (inod[loc_inod011]) << endl
<< omega.node (inod[loc_inod010]) << endl
<< omega.node (inod[loc_inod110]) << endl << endl;
}
if (k+1 == order) {
gdat << omega.node (inod[loc_inod101]) << endl
<< omega.node (inod[loc_inod001]) << endl
<< omega.node (inod[loc_inod011]) << endl << endl;
}
if (i+1 == order && j+1 == order) {
gdat << omega.node (inod[loc_inod110]) << endl
<< omega.node (inod[loc_inod111]) << endl << endl;
}
if (i+1 == order && k+1 == order) {
gdat << omega.node (inod[loc_inod101]) << endl
<< omega.node (inod[loc_inod111]) << endl << endl;
}
if (j+1 == order && k+1 == order) {
gdat << omega.node (inod[loc_inod011]) << endl
<< omega.node (inod[loc_inod111]) << endl << endl;
}
}
}
}
} else { // prism:
error_macro ("gnuplot volume '" << K.name() << "': not yet!");
}
}
template<class T>
static
void
put (std::ostream& gdat, const geo_element& K, const geo_basic<T,sequential>& omega,
const basis_on_pointset<T>& pointset, size_t subdivide)
{
switch (K.dimension()) {
case 0: put_vertex (gdat, K, omega); break;
case 1: put_edge (gdat, K, omega, pointset, subdivide); break;
case 2: put_face (gdat, K, omega, pointset, subdivide); break;
case 3: put_volume (gdat, K, omega); break;
default: break;
}
}
// ----------------------------------------------------------------------------
// mesh puts
// ----------------------------------------------------------------------------
template <class T>
odiststream&
visu_gnuplot (odiststream& ops, const geo_basic<T,sequential>& omega)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
ostream& os = ops.os();
bool verbose = iorheo::getverbose(os);
bool clean = iorheo::getclean(os);
bool execute = iorheo::getexecute(os);
string basename = iorheo::getbasename(os);
bool full = iorheo::getfull(os); // all edges & !faces in 3d
bool lattice = iorheo::getlattice(os); // face & volume lattices
bool reader_on_stdin = iorheo::getreader_on_stdin(os);
bool color = iorheo::getcolor(os);
string format = iorheo::getimage_format(os);
size_type subdivide = iorheo::getsubdivide(os);
if (basename.length() == 0) basename = "output";
string filelist;
size_type dim = max (omega.dimension(), size_type(2));
size_type map_dim = omega.map_dimension();
size_type order = omega.order();
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();
if (order == 1) lattice = false;
bool show_volumes = (lattice && full && map_dim == 3);
bool show_faces = ((!full || lattice) && map_dim == 3) || (lattice && map_dim < 3);
bool show_edges = ((full && map_dim == 3) || (map_dim < 3)) || (order > 1);
bool show_vertices = (map_dim <= 1);
bool show_domains = true;
subdivide = std::max(omega.order(), subdivide);
string outfile_fmt = "";
string tmp = get_tmpdir() + "/";
if (!clean) tmp = "";
string filename = tmp+basename + ".plot";
filelist = filelist + " " + filename;
ofstream plot (filename.c_str());
ofstream gdat;
if (verbose) clog << "! file \"" << filename << "\" created." << endl;
basis_basic<T> subdivide_pointset ("P"+std::to_string(subdivide));
basis_on_pointset<T> pointset (subdivide_pointset, omega.get_piola_basis());
plot << "#!gnuplot" << endl
<< setprecision(numeric_limits<T>::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;
}
if (format == "") {
plot << "set title \"" << basename << ": " << ne << " elements, " << nv << " vertices\"" << endl;
}
check_macro (omega.dimension() > 0, "unsupported 0d geo gnuplot output");
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;
plot << "set xrange [" << xmin[0] << ":" << xmax[0] << "]" << endl;
if (omega.dimension() == 1) {
plot << "set yrange [-1:1]" << endl;
}
if (omega.dimension() >= 2) {
plot << "set yrange [" << xmin[1] << ":" << xmax[1] << "]" << endl;
}
if (omega.dimension() == 2) {
plot << "set size ratio -1 # equal scales" << endl
<< "#set key left Right at graph 1,1" << endl;
}
if (omega.dimension() == 3) {
plot << "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;
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 (!full) {
plot << "set hidden3d nooffset" << endl;
}
}
if (format != "") {
plot << "set nokey" << endl
<< "set noborder" << endl
<< "set notics" << endl;
}
if (omega.dimension() <= 2) {
plot << "plot \\" << endl;
} else {
plot << "splot \\" << endl;
}
bool first_plot_line = true;
//
// plot internal volume, as lattice, for high order elements
//
if (show_volumes) {
filename = tmp+basename + "-vol.gdat";
if (!first_plot_line) plot << ",\\" << endl;
first_plot_line = false;
plot << " \"" << filename << "\" u 1:2:3:(0.0) title \"volumes\" with l lc 3 lw 1";
filelist = filelist + " " + filename;
gdat.open (filename.c_str());
if (verbose) clog << "! file \"" << filename << "\" created." << endl;
gdat << setprecision(numeric_limits<T>::digits10);
for (size_type ivol = 0; ivol < nvol; ivol++) {
const geo_element& K = omega.get_geo_element(3,ivol);
put_volume (gdat, K, omega);
}
gdat.close();
}
//
// plot faces
//
if (show_faces) {
filename = tmp+basename + "-fac.gdat";
if (!first_plot_line) plot << ",\\" << endl;
first_plot_line = false;
plot << " \"" << filename << "\" title \"faces\" with l lc 3 lw 1";
filelist = filelist + " " + filename;
gdat.open (filename.c_str());
if (verbose) clog << "! file \"" << filename << "\" created." << endl;
gdat << setprecision(numeric_limits<T>::digits10);
for (size_type ifac = 0; ifac < nfac; ifac++) {
const geo_element& F = omega.get_geo_element(2,ifac);
put_face (gdat, F, omega, pointset, subdivide);
}
gdat.close();
}
//
// plot edges
//
if (show_edges) {
filename = tmp+basename + "-edg.gdat";
if (!first_plot_line) plot << ",\\" << endl;
first_plot_line = false;
plot << " \"" << filename << "\" title \"edges\" with l lc 1 lw 1.5";
filelist = filelist + " " + filename;
gdat.open (filename.c_str());
if (verbose) clog << "! file \"" << filename << "\" created." << endl;
gdat << setprecision(numeric_limits<T>::digits10);
for (size_type iedg = 0; iedg < nedg; iedg++) {
const geo_element& E = omega.get_geo_element(1,iedg);
put_edge (gdat, E, omega, pointset, subdivide);
}
gdat.close();
}
//
// plot vertices
//
if (show_vertices) {
if (!first_plot_line) plot << ",\\" << endl;
first_plot_line = false;
plot << " \"" << tmp+basename << "-ver.gdat\" title \"vertices\" with p lc 0";
filename = tmp+basename + "-ver.gdat";
filelist = filelist + " " + filename;
gdat.open (filename.c_str());
if (verbose) clog << "! file \"" << filename << "\" created." << endl;
gdat << setprecision(numeric_limits<T>::digits10);
for (size_type iv = 0; iv < nv; iv++) {
const geo_element& P = omega.get_geo_element(0,iv);
put_vertex (gdat, P, omega);
}
gdat.close();
}
//
// plot domains, by decreasing dimension order
// 3d: faces opaques si (!full) ou edges sinon : couleur = code du domaine
// 2d: edges
// 1d: vertex
if (show_domains) {
size_type line_color = 3;
for (size_type dim = omega.map_dimension(); dim+1 >= 1; dim--) {
for (size_type idom = 0; idom < omega.n_domain_indirect(); idom++) {
const domain_indirect_basic<sequential>& dom = omega.get_domain_indirect (idom);
if (dom.map_dimension() != dim) continue;
if (dom.size() == 0) continue;
if (dim == 3) continue; // TODO: only d=0,1,2 domains yet
filename = tmp+basename + "-dom-" + dom.name() + ".gdat";
if (!first_plot_line) plot << ",\\" << endl;
first_plot_line = false;
plot << " \"" << filename + "\" title \"" << dom.name() << "\"";
if (dom.map_dimension() == 0) {
plot << " with p";
} else if (dom.map_dimension() == 1 && omega.dimension() == 1) {
plot << " with lp";
} else {
plot << " with l";
}
plot << " lc " << line_color << " lw 2";
line_color++;
filelist = filelist + " " + filename;
gdat.open (filename.c_str());
if (verbose) clog << "! file \"" << filename << "\" created." << endl;
gdat << setprecision(numeric_limits<T>::digits10);
for (size_type ioige = 0; ioige < dom.size(); ioige++) {
size_type ige = dom.oige(ioige).index();
const geo_element& S = omega.get_geo_element (dim, ige);
put (gdat, S, omega, pointset, subdivide);
}
gdat.close();
}
}
}
//
// end of plot
//
plot << endl;
if (format == "" && !reader_on_stdin) {
plot << "pause -1 \"<return>\"" << endl;
}
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 ops;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template odiststream& visu_gnuplot (odiststream& ops, const geo_basic<Float,sequential>& omega);
} // rheolef namespace
rheolef::basis_basic
Definition basis.h:566
rheolef::basis_on_pointset
Definition basis_on_pointset.h:192
rheolef::domain_indirect_basic
the finite element boundary domain
Definition domain_indirect.h:335
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
rheolef::odiststream::os
std::ostream & os()
Definition diststream.h:247
rheolef::point_basic
Definition point.h:86
rheolef::point_basic::put
std::ostream & put(std::ostream &s, int d=3) const
Definition point.h:388
rheolef::pointset
Definition pointset.h:79
rheolef::reference_element_H::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:790
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:582
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:350
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:481
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:413
rheolef::reference_element::H
static const variant_type H
Definition reference_element.h:81
rheolef::reference_element::q
static const variant_type q
Definition reference_element.h:78
rheolef::reference_element::T
static const variant_type T
Definition reference_element.h:79
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.dim
int dim
Definition mkgeo_ball.sh:307
mkgeo_ball.order
int order
Definition mkgeo_ball.sh:343
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition space_numbering.cc:177
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::piola_transformation
void piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x)
Definition piola_util.cc:43
rheolef::get_tmpdir
std::string get_tmpdir()
get_tmpdir: see the rheostream page for the full documentation
Definition rheostream.cc:54
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
rheolef::file_exists
bool file_exists(const std::string &filename)
file_exists: see the rheostream page for the full documentation
Definition scatch.icc:34
rheolef::visu_gnuplot
odiststream & visu_gnuplot(odiststream &, const field_basic< T, sequential > &)
Definition field_seq_visu_gnuplot.cc:461
std
STL namespace.