← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/field_seq_visu_gnuplot.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/plib/field_seq_visu_gnuplot.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
21
//
 
22
// gnuplot geo visualisation
 
23
//
 
24
// author: Pierre.Saramito@imag.fr
 
25
//
 
26
// date: 16 sept 2011
 
27
//
 
28
# include "rheolef/field.h"
 
29
# include "rheolef/field_indirect.h"
 
30
# include "rheolef/field_expr.h"
 
31
# include "rheolef/field_expr_ops.h"
 
32
# include "rheolef/field_evaluate.h"
 
33
# include "rheolef/piola.h"
 
34
# include "rheolef/rounder.h"
 
35
# include "rheolef/rheostream.h"
 
36
 
 
37
namespace rheolef {
 
38
 
 
39
template <class T>
 
40
odiststream& visu_gnuplot (odiststream& ops, const geo_basic<T,sequential>& omega);
 
41
 
 
42
// ----------------------------------------------------------------------------
 
43
// bounds: for curved geometries, need to commpute it on the fly
 
44
// ----------------------------------------------------------------------------
 
45
template <class T>
 
46
struct bound_type {
 
47
// data:
 
48
  point_basic<T> xmin, xmax;
 
49
  T              umin, umax;
 
50
// modifier:
 
51
  void update (const point_basic<T>& x, const T& u) {
 
52
    for (size_t i = 0; i < 3; i++) {
 
53
      xmin[i] = std::min(xmin[i], x[i]);
 
54
      xmax[i] = std::max(xmax[i], x[i]);
 
55
    }
 
56
    umin = std::min(umin, u);
 
57
    umax = std::max(umax, u);
 
58
  }
 
59
};
 
60
// ----------------------------------------------------------------------------
 
61
// puts for one element
 
62
// ----------------------------------------------------------------------------
 
63
template<class T>
 
64
static
 
65
void
 
66
put_edge (std::ostream& gdat, const geo_element& E, const geo_basic<T,sequential>& omega,
 
67
        const field_basic<T,sequential>& uh, 
 
68
        const basis_on_nodal_basis& geo_on_pointset,
 
69
        const basis_on_nodal_basis& field_on_pointset,
 
70
        bound_type<T>& bbox)
 
71
{
 
72
  using namespace std;
 
73
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
74
  typedef point_basic<size_type>                      ilat;
 
75
  size_type dim = omega.dimension();
 
76
  size_type my_order = geo_on_pointset.pointset_order();
 
77
  std::vector<size_type> dis_idof1;
 
78
  uh.get_space().dis_idof (E, dis_idof1);
 
79
  std::vector<size_type> dis_inod;
 
80
  omega.dis_inod (E, dis_inod);
 
81
  std::vector<point_basic<T> > hat_xi;
 
82
  uh.get_space().get_numbering().get_basis().hat_node (E.variant(), hat_xi);
 
83
  for (size_type i = 0; i <= my_order; i++) {
 
84
    size_type iloc = reference_element_e::ilat2loc_inod (my_order, ilat(i));
 
85
    point_basic<T> xi = piola_transformation (omega, geo_on_pointset,   E.variant(), dis_inod, iloc);
 
86
    T              ui = field_evaluate       (uh,    field_on_pointset, E.variant(), dis_idof1, iloc);
 
87
    xi.put (gdat, dim); gdat << " " << ui << endl;
 
88
    bbox.update (xi, ui);
 
89
  }
 
90
  gdat << endl;
 
91
}
 
92
template<class T>
 
93
static
 
94
void
 
95
put_triangle (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
 
96
        const field_basic<T,sequential>& uh,
 
97
        const basis_on_nodal_basis& geo_on_pointset,
 
98
        const basis_on_nodal_basis& field_on_pointset,
 
99
        bound_type<T>& bbox)
 
100
{
 
101
  using namespace std;
 
102
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
103
  typedef point_basic<size_type>                      ilat;
 
104
  size_type dim = omega.dimension();
 
105
  std::vector<size_type> dis_idof1;
 
106
  std::vector<size_type> dis_inod;
 
107
  omega.dis_inod (F, dis_inod);
 
108
  uh.get_space().dis_idof (F, dis_idof1);
 
109
  size_type nloc = geo_on_pointset.size(F);
 
110
  std::vector<point_basic<T> > xdof (nloc);
 
111
  std::vector<T>               udof (nloc);
 
112
  for (size_type iloc = 0; iloc < nloc; iloc++) {
 
113
    xdof[iloc] = piola_transformation (omega, geo_on_pointset,   F.variant(), dis_inod, iloc);
 
114
    udof[iloc] = field_evaluate       (uh,    field_on_pointset, F.variant(), dis_idof1, iloc);
 
115
    bbox.update (xdof[iloc], udof[iloc]);
 
116
  }
 
117
  size_type my_order = geo_on_pointset.pointset_order();
 
118
  for (size_type j = 0; j <= my_order; j++) {
 
119
    for (size_type i1 = 0; i1 <= my_order; i1++) {
 
120
      size_type i = std::min(i1, my_order-j);
 
121
      size_type iloc = reference_element_t::ilat2loc_inod (my_order, ilat(i, j));
 
122
      xdof [iloc].put (gdat, dim); gdat << " " << udof[iloc] << endl;
 
123
    }
 
124
    gdat << endl;
 
125
  }
 
126
  gdat << endl << endl;
 
127
}
 
128
template<class T>
 
129
static
 
130
void
 
131
put_quadrangle (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
 
132
        const field_basic<T,sequential>& uh, const basis_on_nodal_basis& b)
 
133
{
 
134
  using namespace std;
 
135
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
136
  typedef point_basic<size_type>                      ilat;
 
137
  size_type dim = omega.dimension();
 
138
  size_type degree = uh.get_space().get_numbering().get_basis().degree();
 
139
  std::vector<size_type> dis_idof1;
 
140
  uh.get_space().dis_idof (F, dis_idof1);
 
141
  std::vector<size_type> dis_inod;
 
142
  omega.dis_inod (F, dis_inod);
 
143
  std::vector<point_basic<T> > xdof (dis_idof1.size());
 
144
  for (size_type loc_idof = 0, loc_ndof = xdof.size(); loc_idof < loc_ndof; loc_idof++) {
 
145
    xdof[loc_idof] = piola_transformation (omega, b, F.variant(), dis_inod, loc_idof);
 
146
  }
 
147
  for (size_type j = 0; j < degree+1; j++) {
 
148
    for (size_type i = 0; i < degree+1; i++) {
 
149
      size_type loc_idof00 = reference_element_q::ilat2loc_inod (degree, ilat(i, j));
 
150
      xdof [loc_idof00].put (gdat, dim); gdat << " " << uh.dof(dis_idof1[loc_idof00]) << endl;
 
151
    }
 
152
    gdat << endl;
 
153
  }
 
154
  gdat << endl << endl;
 
155
}
 
156
template<class T>
 
157
static
 
158
void
 
159
put (std::ostream& gdat, const geo_element& K, const geo_basic<T,sequential>& omega,
 
160
        const field_basic<T,sequential>& uh,
 
161
        const basis_on_nodal_basis& geo_on_pointset,
 
162
        const basis_on_nodal_basis& field_on_pointset,
 
163
        bound_type<T>& bbox)
 
164
{
 
165
  switch (K.variant()) {
 
166
   case reference_element::e: put_edge       (gdat, K, omega, uh, geo_on_pointset, field_on_pointset, bbox); break;
 
167
   case reference_element::t: put_triangle   (gdat, K, omega, uh, geo_on_pointset, field_on_pointset, bbox); break;
 
168
   case reference_element::q: put_quadrangle (gdat, K, omega, uh, geo_on_pointset); break;
 
169
   default: error_macro ("unsupported element variant `"<<K.variant()<<"'");
 
170
  }
 
171
}
 
172
// ----------------------------------------------------------------------------
 
173
// scalar field puts
 
174
// ----------------------------------------------------------------------------
 
175
template <class T>
 
176
odiststream&
 
177
visu_gnuplot_scalar (odiststream& ods, const field_basic<T,sequential>& uh)
 
178
{
 
179
  using namespace std;
 
180
  typedef typename field_basic<T,sequential>::float_type float_type;
 
181
  typedef typename geo_basic<float_type,sequential>::size_type size_type;
 
182
  typedef point_basic<size_type>                      ilat;
 
183
  ostream& os = ods.os();
 
184
  bool verbose  = iorheo::getverbose(os);
 
185
  bool clean    = iorheo::getclean(os);
 
186
  bool execute  = iorheo::getexecute(os);
 
187
  bool fill     = iorheo::getfill(os);    // show grid or fill elements
 
188
  bool elevation = iorheo::getelevation(os);
 
189
  bool color   = iorheo::getcolor(os);
 
190
  bool gray    = iorheo::getgray(os);
 
191
  bool black_and_white = iorheo::getblack_and_white(os);
 
192
  bool reader_on_stdin = iorheo::getreader_on_stdin(os);
 
193
  string format   = iorheo::getimage_format(os);
 
194
  string basename = iorheo::getbasename(os);
 
195
  size_type subdivide = iorheo::getsubdivide(os);
 
196
  size_type n_isovalue = iorheo::getn_isovalue(os);
 
197
  size_type n_isovalue_negative = iorheo::getn_isovalue_negative(os);
 
198
  string outfile_fmt = "";
 
199
  string tmp = get_tmpdir() + "/";
 
200
  if (!clean) tmp = "";
 
201
 
 
202
  const geo_basic<float_type,sequential>& omega = uh.get_geo();
 
203
  size_type dim     = omega.dimension();
 
204
  size_type map_dim = omega.map_dimension();
 
205
  size_type nv      = omega.sizes().ownership_by_dimension[0].size();
 
206
  size_type nedg    = omega.sizes().ownership_by_dimension[1].size();
 
207
  size_type nfac    = omega.sizes().ownership_by_dimension[2].size();
 
208
  size_type nvol    = omega.sizes().ownership_by_dimension[3].size();
 
209
  size_type ne      = omega.sizes().ownership_by_dimension[map_dim].size();
 
210
 
 
211
  const numbering<float_type,sequential>& fem = uh.get_space().get_numbering();
 
212
  if (subdivide == 0) { // subdivide is unset: use default
 
213
    subdivide = std::max(omega.order(), subdivide);
 
214
    subdivide = std::max(fem.degree (), subdivide);
 
215
  }
 
216
  basis subdivide_pointset ("P"+itos(subdivide));
 
217
  basis_on_nodal_basis   geo_on_pointset (subdivide_pointset, omega.get_piola_basis());
 
218
  basis_on_nodal_basis field_on_pointset (subdivide_pointset, fem.get_basis());
 
219
 
 
220
  bound_type<T> bbox;
 
221
  bbox.xmin = omega.xmin();
 
222
  bbox.xmax = omega.xmax();
 
223
  bbox.umin = uh.min();
 
224
  bbox.umax = uh.max();
 
225
 
 
226
  std::vector<T> values;
 
227
  if (n_isovalue_negative != 0) {
 
228
    for (size_t i = 0; i <= n_isovalue_negative; i++) {
 
229
      values.push_back (bbox.umin*(n_isovalue_negative - i)/n_isovalue_negative);
 
230
    }
 
231
    for (size_t i = 1; i <= n_isovalue - n_isovalue_negative; i++) {
 
232
      values.push_back (bbox.umax*i/(n_isovalue - n_isovalue_negative));
 
233
    }
 
234
  }
 
235
  string filelist;
 
236
  //
 
237
  // output .gdat
 
238
  //
 
239
  string filename = tmp+basename + ".gdat";
 
240
  string gdatname = filename;
 
241
  filelist = filelist + " " + filename;
 
242
  ofstream gdat (filename.c_str());
 
243
  if (verbose) clog << "! file \"" << filename << "\" created.\n";
 
244
  gdat << setprecision(numeric_limits<float_type>::digits10);
 
245
  size_type used_dim = (fill ? map_dim : 1);
 
246
  for (size_type ie = 0, ne = omega.size(used_dim); ie < ne; ie++) {
 
247
    const geo_element& K = omega.get_geo_element(used_dim,ie);
 
248
    put (gdat, K, omega, uh, geo_on_pointset, field_on_pointset, bbox);
 
249
  }
 
250
  gdat.close();
 
251
  //
 
252
  // rounding bounds
 
253
  //
 
254
  Float eps = 1e-7;
 
255
  bbox.umin = floorer(eps) (bbox.umin);
 
256
  bbox.umax = ceiler(eps)  (bbox.umax);
 
257
  //
 
258
  // output .plot
 
259
  //
 
260
  filename = tmp+basename + ".plot";
 
261
  filelist = filelist + " " + filename;
 
262
  ofstream plot (filename.c_str());
 
263
  if (verbose) clog << "! file \"" << filename << "\" created.\n";
 
264
 
 
265
  plot << "#!gnuplot" << endl
 
266
       << setprecision(numeric_limits<float_type>::digits10);
 
267
  if (format != "") {
 
268
    outfile_fmt = basename + "." + format;
 
269
    string terminal = format;
 
270
    if (terminal == "ps")  {
 
271
      terminal = "postscript eps";
 
272
      if (color) terminal += " color";
 
273
    }
 
274
    if (terminal == "jpg") terminal = "jpeg";
 
275
    if (terminal == "jpeg" || terminal == "png" || terminal == "gif") {
 
276
      terminal += " crop";
 
277
    }
 
278
    plot << "set terminal " << terminal    << endl
 
279
         << "set output \"" << outfile_fmt << "\"" << endl;
 
280
  }
 
281
  if (!black_and_white) {
 
282
    plot << "set noborder" << endl;
 
283
  }
 
284
  if (dim == 2) {
 
285
    plot << "umin = " << bbox.umin << endl
 
286
         << "umax = " << bbox.umax << endl;
 
287
    if (bbox.xmin[0] >= bbox.xmax[0]) plot << "#";
 
288
    plot << "set xrange [" << bbox.xmin[0] << ":" << bbox.xmax[0] << "]" << endl;
 
289
    if (bbox.xmin[1] >= bbox.xmax[1]) plot << "#";
 
290
    plot << "set yrange [" << bbox.xmin[1] << ":" << bbox.xmax[1] << "]" << endl;
 
291
    if (bbox.umin >= bbox.umax) plot << "#";
 
292
    plot << "set zrange [umin:umax]" << endl;
 
293
    if (bbox.xmin[0] >= bbox.xmax[0] || bbox.xmin[1] >= bbox.xmax[1]) {
 
294
      plot << "set size square" << endl;
 
295
    } else {
 
296
      plot << "set size ratio -1 # equal scales" << endl
 
297
           << "set noxtics" << endl
 
298
           << "set noytics" << endl;
 
299
    }
 
300
    if (elevation) {
 
301
      plot << "set xyplane at umin-0.1*(umax-umin)" << endl;
 
302
    } else {
 
303
      plot << "set view map" << endl;
 
304
    }
 
305
    if (color || gray) {
 
306
      plot << "set pm3d interpolate 10,10 corners2color mean" << endl;
 
307
      if (gray) {
 
308
        plot << "set palette gray" << endl;
 
309
      } else {
 
310
        plot << "set palette rgbformulae 33,13,-4" << endl;
 
311
      }
 
312
    } else { // black-and-white
 
313
      if (values.size() != 0) {
 
314
        plot << "set cntrparam levels discrete ";
 
315
        for (size_t i = 0, n = values.size(); i < n; i++) {
 
316
          plot << values[i];
 
317
          if (i != n-1) plot << ", ";
 
318
        }
 
319
        plot << endl;
 
320
      } else {
 
321
        plot << "set cntrparam levels incremental umin,(umax-umin)/10.0,umax"<< endl;
 
322
      }
 
323
      plot << "unset surface" << endl
 
324
           << "set contour base" << endl
 
325
           << "set nokey" << endl
 
326
           << "set palette rgbformulae 0,0,0" << endl
 
327
           << "unset colorbox" << endl;
 
328
    }
 
329
  }
 
330
  if (dim == 1) {
 
331
    plot << "plot \"" << gdatname << "\" notitle with lines";
 
332
    if (black_and_white) {
 
333
      plot << " lc 0" << endl;
 
334
    }
 
335
    plot << endl;
 
336
  } else {
 
337
    if (!fill && dim == 2) {
 
338
      plot << "plot";
 
339
    } else {
 
340
      plot << "splot";
 
341
    } 
 
342
    plot << " \"" << gdatname << "\" notitle";
 
343
    if (map_dim == 2) {
 
344
      if (!black_and_white && fill) {
 
345
        plot << " with pm3d" << endl;
 
346
      } else {
 
347
        plot << " with lines palette" << endl;
 
348
      }
 
349
    } else { // a 2d line
 
350
      plot << " with lines palette lw 2" << endl;
 
351
    }
 
352
  }
 
353
  //
 
354
  // end of plot
 
355
  //
 
356
  if (format == "" && !reader_on_stdin) {
 
357
    plot << "pause -1 \"<return>\"\n";
 
358
  }
 
359
  plot.close();
 
360
  //
 
361
  // run gnuplot
 
362
  //
 
363
  int status = 0;
 
364
  string command;
 
365
  if (execute) {
 
366
      command = "gnuplot ";
 
367
      if (reader_on_stdin) command += "-persist ";
 
368
      command += tmp + basename + ".plot";
 
369
      if (verbose) clog << "! " << command << endl;
 
370
      cin.sync();
 
371
      status = system (command.c_str());
 
372
      if (format != "") {
 
373
        check_macro (file_exists (outfile_fmt), "! file \"" << outfile_fmt << "\" creation failed");
 
374
        if (verbose) clog << "! file \"" << outfile_fmt << "\" created" << endl;
 
375
      }
 
376
  }
 
377
  //
 
378
  // clear gnuplot data
 
379
  //
 
380
  if (clean) {
 
381
      command = "/bin/rm -f " + filelist;
 
382
      if (verbose) clog << "! " << command << endl;
 
383
      status = system (command.c_str());
 
384
  }
 
385
  return ods;
 
386
}
 
387
// ----------------------------------------------------------------------------
 
388
// vector field puts
 
389
// ----------------------------------------------------------------------------
 
390
template <class T>
 
391
odiststream&
 
392
visu_gnuplot_vector (odiststream& ods, const field_basic<T,sequential>& uh)
 
393
{
 
394
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
395
  const geo_basic<T,sequential>& omega = uh.get_geo();
 
396
  check_macro (omega.get_piola_basis().name() == uh.get_space().get_numbering().name(),
 
397
        "gnuplot vector: unsupported non-isoparametric approx " << uh.get_space().get_numbering().name());
 
398
  size_type n_elt = omega.size();
 
399
  size_type d     = omega.dimension();
 
400
  Float diam_omega = norm (omega.xmax() - omega.xmin()); // characteristic length in omega
 
401
  Float h_moy = diam_omega/pow(n_elt,1./d);
 
402
  field_basic<T,sequential> norm_uh = norm(uh);
 
403
  Float norm_max_uh = norm_uh.max_abs();      // get max vector length
 
404
  if (norm_max_uh + 1 == 1) norm_max_uh = 1;
 
405
  Float scale = h_moy/norm_max_uh;
 
406
  size_type n_comp = uh.size();
 
407
  std::vector<field_component_const<T,sequential> > uh_comp (n_comp);
 
408
  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
 
409
    uh_comp[i_comp] = uh[i_comp];
 
410
  }
 
411
  array<point_basic<T>, sequential> x = omega.get_nodes();
 
412
  for (size_type inod = 0, nnod = x.size(); inod < nnod; inod++) {
 
413
    point vi;
 
414
    for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
 
415
      vi[i_comp] = uh[i_comp].dof(inod); // TODO: non-isoparam => interpolate
 
416
    }
 
417
    x[inod] += vi;
 
418
  }
 
419
  geo_basic<T,sequential> deformed_omega = omega;
 
420
  deformed_omega.set_nodes(x);
 
421
  space_basic<T,sequential> deformed_Vh (deformed_omega, norm_uh.get_space().get_numbering().name());
 
422
  field_basic<T,sequential> deformed_norm_uh (deformed_Vh);
 
423
  std::copy (norm_uh.begin_dof(), norm_uh.end_dof(), deformed_norm_uh.begin_dof());
 
424
  visu_gnuplot (ods, deformed_norm_uh);
 
425
  return ods;
 
426
}
 
427
// ----------------------------------------------------------------------------
 
428
// switch
 
429
// ----------------------------------------------------------------------------
 
430
template <class T>
 
431
odiststream&
 
432
visu_gnuplot (odiststream& ods, const field_basic<T,sequential>& uh)
 
433
{
 
434
  switch (uh.get_space().valued_tag()) {
 
435
    case space_constant::scalar: visu_gnuplot_scalar (ods, uh); break;
 
436
    case space_constant::vector: visu_gnuplot_vector (ods, uh); break;
 
437
    default: error_macro ("do not known how to print " << uh.valued() << "-valued field");
 
438
  }
 
439
  return ods;
 
440
}
 
441
// ----------------------------------------------------------------------------
 
442
// instanciation in library
 
443
// ----------------------------------------------------------------------------
 
444
template odiststream& visu_gnuplot<Float> (odiststream&, const field_basic<Float,sequential>&);
 
445
 
 
446
} // rheolef namespace