28#include "rheolef/geo.h"
29#include "rheolef/rheostream.h"
30#include "rheolef/iorheo.h"
31#include "rheolef/reference_element.h"
32#include "rheolef/piola_util.h"
42put_vertex (std::ostream& gdat,
const geo_element& P,
const geo_basic<T,sequential>& omega)
45 typedef typename geo_basic<T,sequential>::size_type
size_type;
47 omega.node (P[0]).put (gdat, dim); gdat << endl;
53put_edge (std::ostream& gdat,
const geo_element& E,
const geo_basic<T,sequential>& omega,
54 const basis_on_pointset<T>& pointset,
size_t subdivide)
57 typedef typename geo_basic<T,sequential>::size_type
size_type;
58 typedef point_basic<size_type> ilat;
61 omega.dis_inod (E, dis_inod);
62 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x;
64 for (
size_type i = 0; i <= subdivide; i++) {
66 x[loc_inod].put (gdat, dim); gdat << endl;
74put_2d_face (std::ostream& gdat,
const geo_element& F,
const geo_basic<T,sequential>& omega,
75 const basis_on_pointset<T>& pointset,
size_t subdivide)
78 typedef typename geo_basic<T,sequential>::size_type
size_type;
79 typedef point_basic<size_type> ilat;
82 omega.dis_inod (F, dis_inod);
83 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x;
88 for (
size_type j = 0; j < subdivide; j++) {
89 for (
size_type i = 0; i < subdivide; i++) {
94 const point_basic<T>& x00 = x[loc_inod00];
95 const point_basic<T>& x10 = x[loc_inod10];
96 const point_basic<T>& x01 = x[loc_inod01];
97 const point_basic<T>& x11 = x[loc_inod11];
98 if (i+1 == subdivide) {
99 x11.put (gdat, dim); gdat << endl;
101 x10.put (gdat, dim); gdat << endl;
102 x00.put (gdat, dim); gdat << endl;
103 x01.put (gdat, dim); gdat << endl;
104 if (j+1 == subdivide) {
105 x11.put (gdat, dim); gdat << endl;
113 for (
size_type i = 0; i < subdivide; i++) {
114 for (
size_type j = 0; j < subdivide - i; j++) {
118 const point_basic<T>& x00 = x[loc_inod00];
119 const point_basic<T>& x10 = x[loc_inod10];
120 const point_basic<T>& x01 = x[loc_inod01];
121 x10.put (gdat, dim); gdat << endl;
122 x00.put (gdat, dim); gdat << endl;
123 x01.put (gdat, dim); gdat << endl;
124 if (i+j+1 == subdivide) {
125 x10.put (gdat, dim); gdat << endl;
137put_3d_face (std::ostream& gdat,
const geo_element& F,
const geo_basic<T,sequential>& omega)
140 typedef typename geo_basic<T,sequential>::size_type
size_type;
141 typedef point_basic<size_type> ilat;
144 std::vector<size_type> inod;
145 omega.dis_inod (F, inod);
146 gdat <<
"# F"<<F.dis_ie()<<endl;
155 omega.node (inod[loc_inod]).put (gdat, dim); gdat << endl;
164put_face (std::ostream& gdat,
const geo_element& F,
const geo_basic<T,sequential>& omega,
165 const basis_on_pointset<T>& pointset,
size_t subdivide)
167 switch (omega.dimension()) {
168 case 2: put_2d_face (gdat, F, omega, pointset, subdivide);
break;
169 case 3: put_3d_face (gdat, F, omega);
break;
176put_volume (std::ostream& gdat,
const geo_element& K,
const geo_basic<T,sequential>& omega)
179 typedef typename geo_basic<T,sequential>::size_type
size_type;
180 typedef point_basic<size_type> ilat;
183 std::vector<size_type> inod;
184 omega.dis_inod (K, inod);
185 gdat <<
"# K"<<K.dis_ie()<<endl;
195 gdat << omega.node (inod[loc_inod100]) << endl
196 << omega.node (inod[loc_inod000]) << endl
197 << omega.node (inod[loc_inod010]) << endl << endl
198 << omega.node (inod[loc_inod000]) << endl
199 << omega.node (inod[loc_inod001]) << endl << endl;
200 if (i+j+k+1 == order) {
201 gdat << omega.node (inod[loc_inod100]) << endl
202 << omega.node (inod[loc_inod010]) << endl
203 << omega.node (inod[loc_inod001]) << endl
204 << omega.node (inod[loc_inod100]) << endl << endl;
222 gdat << omega.node (inod[loc_inod100]) << endl
223 << omega.node (inod[loc_inod000]) << endl
224 << omega.node (inod[loc_inod010]) << endl << endl
225 << omega.node (inod[loc_inod000]) << endl
226 << omega.node (inod[loc_inod001]) << endl << endl;
228 gdat << omega.node (inod[loc_inod101]) << endl
229 << omega.node (inod[loc_inod100]) << endl
230 << omega.node (inod[loc_inod110]) << endl << endl;
233 gdat << omega.node (inod[loc_inod011]) << endl
234 << omega.node (inod[loc_inod010]) << endl
235 << omega.node (inod[loc_inod110]) << endl << endl;
238 gdat << omega.node (inod[loc_inod101]) << endl
239 << omega.node (inod[loc_inod001]) << endl
240 << omega.node (inod[loc_inod011]) << endl << endl;
242 if (i+1 == order && j+1 == order) {
243 gdat << omega.node (inod[loc_inod110]) << endl
244 << omega.node (inod[loc_inod111]) << endl << endl;
246 if (i+1 == order && k+1 == order) {
247 gdat << omega.node (inod[loc_inod101]) << endl
248 << omega.node (inod[loc_inod111]) << endl << endl;
250 if (j+1 == order && k+1 == order) {
251 gdat << omega.node (inod[loc_inod011]) << endl
252 << omega.node (inod[loc_inod111]) << endl << endl;
258 error_macro (
"gnuplot volume '" << K.name() <<
"': not yet!");
264put (std::ostream& gdat,
const geo_element& K,
const geo_basic<T,sequential>& omega,
265 const basis_on_pointset<T>& pointset,
size_t subdivide)
267 switch (K.dimension()) {
268 case 0: put_vertex (gdat, K, omega);
break;
269 case 1: put_edge (gdat, K, omega, pointset, subdivide);
break;
270 case 2: put_face (gdat, K, omega, pointset, subdivide);
break;
271 case 3: put_volume (gdat, K, omega);
break;
285 ostream& os = ops.
os();
286 bool verbose = iorheo::getverbose(os);
287 bool clean = iorheo::getclean(os);
288 bool execute = iorheo::getexecute(os);
289 string basename = iorheo::getbasename(os);
290 bool full = iorheo::getfull(os);
291 bool lattice = iorheo::getlattice(os);
292 bool reader_on_stdin = iorheo::getreader_on_stdin(os);
293 bool color = iorheo::getcolor(os);
294 string format = iorheo::getimage_format(os);
295 size_type subdivide = iorheo::getsubdivide(os);
296 if (basename.length() == 0) basename =
"output";
300 size_type map_dim = omega.map_dimension();
302 size_type nv = omega.sizes().ownership_by_dimension[0].size();
303 size_type nedg = omega.sizes().ownership_by_dimension[1].size();
304 size_type nfac = omega.sizes().ownership_by_dimension[2].size();
305 size_type nvol = omega.sizes().ownership_by_dimension[3].size();
306 size_type ne = omega.sizes().ownership_by_dimension[map_dim].size();
308 if (order == 1) lattice =
false;
309 bool show_volumes = (lattice && full && map_dim == 3);
310 bool show_faces = ((!full || lattice) && map_dim == 3) || (lattice && map_dim < 3);
311 bool show_edges = ((full && map_dim == 3) || (map_dim < 3)) || (order > 1);
312 bool show_vertices = (map_dim <= 1);
313 bool show_domains =
true;
314 subdivide = std::max(omega.order(), subdivide);
316 string outfile_fmt =
"";
318 if (!clean) tmp =
"";
319 string filename = tmp+basename +
".plot";
320 filelist = filelist +
" " + filename;
321 ofstream plot (filename.c_str());
323 if (verbose) clog <<
"! file \"" << filename <<
"\" created." << endl;
325 basis_basic<T> subdivide_pointset (
"P"+std::to_string(subdivide));
328 plot <<
"#!gnuplot" << endl
329 << setprecision(numeric_limits<T>::digits10);
331 outfile_fmt = basename +
"." + format;
332 string terminal = format;
333 if (terminal ==
"ps") {
334 terminal =
"postscript eps";
335 if (color) terminal +=
" color";
337 if (terminal ==
"jpg") terminal =
"jpeg";
338 if (terminal ==
"jpeg" || terminal ==
"png" || terminal ==
"gif") {
341 plot <<
"set terminal " << terminal << endl
342 <<
"set output \"" << outfile_fmt <<
"\"" << endl;
345 plot <<
"set title \"" << basename <<
": " << ne <<
" elements, " << nv <<
" vertices\"" << endl;
347 check_macro (omega.dimension() > 0,
"unsupported 0d geo gnuplot output");
349 T dx_max = max(dx[0],max(dx[1],dx[2]));
350 if (dx_max == 0) dx_max = 0.1;
351 dx[0] = max(dx[0],dx_max);
352 if (omega.dimension() >= 2) dx[1] = max(dx[1],dx_max);
353 if (omega.dimension() == 3) dx[2] = max(dx[2],dx_max);
356 plot <<
"set xrange [" << xmin[0] <<
":" << xmax[0] <<
"]" << endl;
357 if (omega.dimension() == 1) {
358 plot <<
"set yrange [-1:1]" << endl;
360 if (omega.dimension() >= 2) {
361 plot <<
"set yrange [" << xmin[1] <<
":" << xmax[1] <<
"]" << endl;
363 if (omega.dimension() == 2) {
364 plot <<
"set size ratio -1 # equal scales" << endl
365 <<
"#set key left Right at graph 1,1" << endl;
367 if (omega.dimension() == 3) {
368 plot <<
"set zrange [" << xmin[2] <<
":" << xmax[2] <<
"]" << endl
369 <<
"set xyplane at " << xmin[2] << endl
370 <<
"set view equal xyz # equal scales" << endl
371 <<
"set view 70,120" << endl;
373 plot <<
"set noxlabel" << endl
374 <<
"set noylabel" << endl
375 <<
"set nozlabel" << endl;
377 plot <<
"set xlabel \"x\"" << endl
378 <<
"set ylabel \"y\"" << endl
379 <<
"set zlabel \"z\"" << endl;
382 plot <<
"set hidden3d nooffset" << endl;
386 plot <<
"set nokey" << endl
387 <<
"set noborder" << endl
388 <<
"set notics" << endl;
390 if (omega.dimension() <= 2) {
391 plot <<
"plot \\" << endl;
393 plot <<
"splot \\" << endl;
395 bool first_plot_line =
true;
400 filename = tmp+basename +
"-vol.gdat";
401 if (!first_plot_line) plot <<
",\\" << endl;
402 first_plot_line =
false;
403 plot <<
" \"" << filename <<
"\" u 1:2:3:(0.0) title \"volumes\" with l lc 3 lw 1";
405 filelist = filelist +
" " + filename;
406 gdat.open (filename.c_str());
407 if (verbose) clog <<
"! file \"" << filename <<
"\" created." << endl;
408 gdat << setprecision(numeric_limits<T>::digits10);
409 for (
size_type ivol = 0; ivol < nvol; ivol++) {
410 const geo_element& K = omega.get_geo_element(3,ivol);
411 put_volume (gdat, K, omega);
419 filename = tmp+basename +
"-fac.gdat";
420 if (!first_plot_line) plot <<
",\\" << endl;
421 first_plot_line =
false;
422 plot <<
" \"" << filename <<
"\" title \"faces\" with l lc 3 lw 1";
424 filelist = filelist +
" " + filename;
425 gdat.open (filename.c_str());
426 if (verbose) clog <<
"! file \"" << filename <<
"\" created." << endl;
427 gdat << setprecision(numeric_limits<T>::digits10);
428 for (
size_type ifac = 0; ifac < nfac; ifac++) {
429 const geo_element& F = omega.get_geo_element(2,ifac);
430 put_face (gdat, F, omega,
pointset, subdivide);
438 filename = tmp+basename +
"-edg.gdat";
439 if (!first_plot_line) plot <<
",\\" << endl;
440 first_plot_line =
false;
441 plot <<
" \"" << filename <<
"\" title \"edges\" with l lc 1 lw 1.5";
443 filelist = filelist +
" " + filename;
444 gdat.open (filename.c_str());
445 if (verbose) clog <<
"! file \"" << filename <<
"\" created." << endl;
446 gdat << setprecision(numeric_limits<T>::digits10);
447 for (
size_type iedg = 0; iedg < nedg; iedg++) {
448 const geo_element& E = omega.get_geo_element(1,iedg);
449 put_edge (gdat, E, omega,
pointset, subdivide);
457 if (!first_plot_line) plot <<
",\\" << endl;
458 first_plot_line =
false;
459 plot <<
" \"" << tmp+basename <<
"-ver.gdat\" title \"vertices\" with p lc 0";
461 filename = tmp+basename +
"-ver.gdat";
462 filelist = filelist +
" " + filename;
463 gdat.open (filename.c_str());
464 if (verbose) clog <<
"! file \"" << filename <<
"\" created." << endl;
465 gdat << setprecision(numeric_limits<T>::digits10);
467 const geo_element& P = omega.get_geo_element(0,iv);
468 put_vertex (gdat, P, omega);
479 for (
size_type dim = omega.map_dimension(); dim+1 >= 1; dim--) {
480 for (
size_type idom = 0; idom < omega.n_domain_indirect(); idom++) {
481 const domain_indirect_basic<sequential>& dom = omega.get_domain_indirect (idom);
482 if (dom.map_dimension() != dim)
continue;
483 if (dom.size() == 0)
continue;
484 if (dim == 3)
continue;
486 filename = tmp+basename +
"-dom-" + dom.name() +
".gdat";
487 if (!first_plot_line) plot <<
",\\" << endl;
488 first_plot_line =
false;
489 plot <<
" \"" << filename +
"\" title \"" << dom.name() <<
"\"";
490 if (dom.map_dimension() == 0) {
492 }
else if (dom.map_dimension() == 1 && omega.dimension() == 1) {
497 plot <<
" lc " << line_color <<
" lw 2";
500 filelist = filelist +
" " + filename;
501 gdat.open (filename.c_str());
502 if (verbose) clog <<
"! file \"" << filename <<
"\" created." << endl;
503 gdat << setprecision(numeric_limits<T>::digits10);
504 for (
size_type ioige = 0; ioige < dom.size(); ioige++) {
506 const geo_element& S = omega.get_geo_element (dim, ige);
517 if (format ==
"" && !reader_on_stdin) {
518 plot <<
"pause -1 \"<return>\"" << endl;
527 command =
"gnuplot ";
528 if (reader_on_stdin) command +=
"-persist ";
529 command += tmp + basename +
".plot";
530 if (verbose) clog <<
"! " << command << endl;
532 status = system (command.c_str());
535 if (verbose) clog <<
"! file \"" << outfile_fmt <<
"\" created" << endl;
542 command =
"/bin/rm -f " + filelist;
543 if (verbose) clog <<
"! " << command << endl;
544 status = system (command.c_str());
field::size_type size_type
generic mesh with rerefence counting
see the geo_element page for the full documentation
odiststream: see the diststream page for the full documentation
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static const variant_type H
static const variant_type q
static const variant_type T
#define error_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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)
This file is part of Rheolef.
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)
std::string get_tmpdir()
get_tmpdir: see the rheostream page for the full documentation
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
bool file_exists(const std::string &filename)
file_exists: see the rheostream page for the full documentation
odiststream & visu_gnuplot(odiststream &, const field_basic< T, sequential > &)