← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/geo_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/geo_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/geo.h"
 
29
#include "rheolef/rheostream.h"
 
30
#include "rheolef/iorheo.h"
 
31
#include "rheolef/reference_element.h"
 
32
#include "rheolef/piola.h"
 
33
 
 
34
namespace rheolef {
 
35
 
 
36
// ----------------------------------------------------------------------------
 
37
// element puts
 
38
// ----------------------------------------------------------------------------
 
39
template<class T>
 
40
static
 
41
void
 
42
put_vertex (std::ostream& gdat, const geo_element& P, const geo_basic<T,sequential>& omega)
 
43
{
 
44
  using namespace std;
 
45
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
46
  size_type dim = max (omega.dimension(), size_type(2));
 
47
  omega.node (P[0]).put (gdat, dim); gdat << endl;
 
48
}
 
49
// simple edge, could be curved (order > 1 => n_node = order+1 > 2)
 
50
template<class T>
 
51
static
 
52
void
 
53
put_edge (std::ostream& gdat, const geo_element& E, const geo_basic<T,sequential>& omega,
 
54
        const basis_on_nodal_basis& pointset)
 
55
{
 
56
  using namespace std;
 
57
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
58
  typedef point_basic<size_type>                      ilat;
 
59
  size_type dim = max (omega.dimension(), size_type(2));
 
60
  std::vector<size_type> dis_inod;
 
61
  size_type order = pointset.pointset_order();
 
62
  omega.dis_inod (E, dis_inod);
 
63
  for (size_type i = 0; i <= order; i++) {
 
64
    size_type loc_inod = reference_element_e::ilat2loc_inod (order, ilat(i));
 
65
    point_basic<T> xi = piola_transformation (omega, pointset, E.variant(), dis_inod, loc_inod);
 
66
    xi.put (gdat, dim); gdat << endl;
 
67
  }
 
68
  gdat << endl << endl;
 
69
}
 
70
// 2d triangle or quadrangle
 
71
template<class T>
 
72
static
 
73
void
 
74
put_2d_face (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
 
75
        const basis_on_nodal_basis& pointset)
 
76
{
 
77
  using namespace std;
 
78
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
79
  typedef point_basic<size_type>                      ilat;
 
80
  size_type dim = max (omega.dimension(), size_type(2));
 
81
  size_type order = pointset.pointset_order();
 
82
  std::vector<size_type> dis_inod;
 
83
  omega.dis_inod (F, dis_inod);
 
84
  if (F.variant() == reference_element::q) {
 
85
    // 2d quadrangle: draw lines of the lattice
 
86
    // ICI
 
87
    for (size_type j = 0; j < order; j++) {
 
88
      for (size_type i = 0; i < order; i++) {
 
89
        size_type loc_inod00 = reference_element_q::ilat2loc_inod (order, ilat(i,   j));
 
90
        size_type loc_inod10 = reference_element_q::ilat2loc_inod (order, ilat(i+1, j));
 
91
        size_type loc_inod01 = reference_element_q::ilat2loc_inod (order, ilat(i,   j+1));
 
92
        size_type loc_inod11 = reference_element_q::ilat2loc_inod (order, ilat(i+1, j+1));
 
93
        point_basic<T> x00 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod00);
 
94
        point_basic<T> x10 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod10);
 
95
        point_basic<T> x01 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod01);
 
96
        point_basic<T> x11 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod11);
 
97
        if (i+1 == order) {
 
98
          x11.put (gdat, dim); gdat << endl;
 
99
        }
 
100
        x10.put (gdat, dim); gdat << endl;
 
101
        x00.put (gdat, dim); gdat << endl;
 
102
        x01.put (gdat, dim); gdat << endl;
 
103
        if (j+1 == order) {
 
104
          x11.put (gdat, dim); gdat << endl;
 
105
        }
 
106
        gdat << endl;
 
107
      }
 
108
    }
 
109
    gdat << endl;
 
110
  } else {
 
111
    // 2d triangle: draw lines of the lattice
 
112
    for (size_type i = 0; i < order; i++) {
 
113
      for (size_type j = 0; j < order - i; j++) {
 
114
        size_type loc_inod00 = reference_element_t::ilat2loc_inod (order, ilat(i,   j));
 
115
        size_type loc_inod10 = reference_element_t::ilat2loc_inod (order, ilat(i+1, j));
 
116
        size_type loc_inod01 = reference_element_t::ilat2loc_inod (order, ilat(i,   j+1));
 
117
        point_basic<T> x00 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod00);
 
118
        point_basic<T> x10 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod10);
 
119
        point_basic<T> x01 = piola_transformation (omega, pointset, F.variant(), dis_inod, loc_inod01);
 
120
        x10.put (gdat, dim); gdat << endl;
 
121
        x00.put (gdat, dim); gdat << endl;
 
122
        x01.put (gdat, dim); gdat << endl;
 
123
        if (i+j+1 == order) {
 
124
          x10.put (gdat, dim); gdat << endl;
 
125
        }
 
126
        gdat << endl;
 
127
      }
 
128
    }
 
129
    gdat << endl;
 
130
  }
 
131
}
 
132
// 3d triangle or quadrangle: splot lattice data, for hidden faces removal
 
133
template<class T>
 
134
static
 
135
void
 
136
put_3d_face (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega)
 
137
{
 
138
  using namespace std;
 
139
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
140
  typedef point_basic<size_type>                      ilat;
 
141
  size_type dim = max (omega.dimension(), size_type(2));
 
142
  size_type order = omega.order();
 
143
  std::vector<size_type> inod;
 
144
  omega.dis_inod (F, inod);
 
145
  gdat << "# F"<<F.dis_ie()<<endl;
 
146
  for (size_type j = 0; j <= order; j++) {
 
147
    for (size_type i = 0; i <= order; i++) {
 
148
      size_type loc_inod;
 
149
      if (F.variant() == reference_element::q) { // TODO: virtual in geo_elment_v2::loc_ilat2loc_inod
 
150
        loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
 
151
      } else {
 
152
        loc_inod = reference_element_t::ilat2loc_inod (order, ilat(min(i,order-j), j));
 
153
      }
 
154
      omega.node (inod[loc_inod]).put (gdat, dim); gdat << endl;
 
155
    }
 
156
    gdat << endl;
 
157
  }
 
158
  gdat << endl;
 
159
}
 
160
template<class T>
 
161
static
 
162
void
 
163
put_face (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
 
164
        const basis_on_nodal_basis& pointset)
 
165
{
 
166
  switch (omega.dimension()) {
 
167
    case 2: put_2d_face (gdat, F, omega, pointset); break;
 
168
    case 3: put_3d_face (gdat, F, omega); break;
 
169
    default: break;
 
170
  }
 
171
}
 
172
template<class T>
 
173
static
 
174
void
 
175
put_volume (std::ostream& gdat, const geo_element& K, const geo_basic<T,sequential>& omega)
 
176
{
 
177
  using namespace std;
 
178
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
179
  typedef point_basic<size_type>                      ilat;
 
180
  size_type dim = max (omega.dimension(), size_type(2));
 
181
  size_type order = omega.order();
 
182
  std::vector<size_type> inod;
 
183
  omega.dis_inod (K, inod);
 
184
  gdat << "# K"<<K.dis_ie()<<endl;
 
185
  if (K.variant() == reference_element::T) {
 
186
        // 3d tetra: draw lattice
 
187
        for (size_type i = 0; i < order; i++) {
 
188
          for (size_type j = 0; j < order - i; j++) {
 
189
            for (size_type k = 0; k < order - i - j; k++) {
 
190
              size_type loc_inod000 = reference_element_T::ilat2loc_inod (order, ilat(i,   j,   k));
 
191
              size_type loc_inod100 = reference_element_T::ilat2loc_inod (order, ilat(i+1, j,   k));
 
192
              size_type loc_inod010 = reference_element_T::ilat2loc_inod (order, ilat(i,   j+1, k));
 
193
              size_type loc_inod001 = reference_element_T::ilat2loc_inod (order, ilat(i,   j,   k+1));
 
194
              gdat << omega.node (inod[loc_inod100]) << endl
 
195
                   << omega.node (inod[loc_inod000]) << endl
 
196
                   << omega.node (inod[loc_inod010]) << endl << endl
 
197
                   << omega.node (inod[loc_inod000]) << endl
 
198
                   << omega.node (inod[loc_inod001]) << endl << endl;
 
199
              if (i+j+k+1 == order) {
 
200
                gdat << omega.node (inod[loc_inod100]) << endl
 
201
                     << omega.node (inod[loc_inod010]) << endl
 
202
                     << omega.node (inod[loc_inod001]) << endl
 
203
                     << omega.node (inod[loc_inod100]) << endl << endl;
 
204
              }
 
205
            }
 
206
          }
 
207
        }
 
208
  } else if (K.variant() == reference_element::H) {
 
209
        // 3d hexa: draw lattice
 
210
        for (size_type k = 0; k < order; k++) {
 
211
          for (size_type j = 0; j < order; j++) {
 
212
            for (size_type i = 0; i < order; i++) {
 
213
              size_type loc_inod000 = reference_element_H::ilat2loc_inod (order, ilat(i,   j,   k));
 
214
              size_type loc_inod100 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j,   k));
 
215
              size_type loc_inod110 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j+1, k));
 
216
              size_type loc_inod010 = reference_element_H::ilat2loc_inod (order, ilat(i,   j+1, k));
 
217
              size_type loc_inod001 = reference_element_H::ilat2loc_inod (order, ilat(i,   j,   k+1));
 
218
              size_type loc_inod101 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j,   k+1));
 
219
              size_type loc_inod111 = reference_element_H::ilat2loc_inod (order, ilat(i+1, j+1, k+1));
 
220
              size_type loc_inod011 = reference_element_H::ilat2loc_inod (order, ilat(i,   j+1, k+1));
 
221
              gdat << omega.node (inod[loc_inod100]) << endl
 
222
                   << omega.node (inod[loc_inod000]) << endl
 
223
                   << omega.node (inod[loc_inod010]) << endl << endl
 
224
                   << omega.node (inod[loc_inod000]) << endl
 
225
                   << omega.node (inod[loc_inod001]) << endl << endl;
 
226
              if (i+1 == order) {
 
227
                gdat << omega.node (inod[loc_inod101]) << endl
 
228
                     << omega.node (inod[loc_inod100]) << endl
 
229
                     << omega.node (inod[loc_inod110]) << endl << endl;
 
230
              }
 
231
              if (j+1 == order) {
 
232
                gdat << omega.node (inod[loc_inod011]) << endl
 
233
                     << omega.node (inod[loc_inod010]) << endl
 
234
                     << omega.node (inod[loc_inod110]) << endl << endl;
 
235
              }
 
236
              if (k+1 == order) {
 
237
                gdat << omega.node (inod[loc_inod101]) << endl
 
238
                     << omega.node (inod[loc_inod001]) << endl
 
239
                     << omega.node (inod[loc_inod011]) << endl << endl;
 
240
              }
 
241
              if (i+1 == order && j+1 == order) {
 
242
                gdat << omega.node (inod[loc_inod110]) << endl
 
243
                     << omega.node (inod[loc_inod111]) << endl << endl;
 
244
              }
 
245
              if (i+1 == order && k+1 == order) {
 
246
                gdat << omega.node (inod[loc_inod101]) << endl
 
247
                     << omega.node (inod[loc_inod111]) << endl << endl;
 
248
              }
 
249
              if (j+1 == order && k+1 == order) {
 
250
                gdat << omega.node (inod[loc_inod011]) << endl
 
251
                     << omega.node (inod[loc_inod111]) << endl << endl;
 
252
              }
 
253
            }
 
254
          }
 
255
        }
 
256
  } else { // prism:
 
257
        error_macro ("gnuplot volume '" << K.name() << "': not yet!");
 
258
  }
 
259
}
 
260
template<class T>
 
261
static
 
262
void
 
263
put (std::ostream& gdat, const geo_element& K, const geo_basic<T,sequential>& omega,
 
264
        const basis_on_nodal_basis& pointset)
 
265
{
 
266
  switch (K.dimension()) {
 
267
    case 0:     put_vertex (gdat, K, omega); break;
 
268
    case 1:     put_edge   (gdat, K, omega, pointset); break;
 
269
    case 2:     put_face   (gdat, K, omega, pointset); break;
 
270
    case 3:     put_volume (gdat, K, omega); break;
 
271
    default: break;
 
272
  }
 
273
}
 
274
// ----------------------------------------------------------------------------
 
275
// mesh puts
 
276
// ----------------------------------------------------------------------------
 
277
template <class T>
 
278
odiststream&
 
279
visu_gnuplot (odiststream& ops, const geo_basic<T,sequential>& omega)
 
280
{
 
281
  using namespace std;
 
282
  typedef typename geo_basic<T,sequential>::size_type size_type;
 
283
  typedef point_basic<size_type>                      ilat;
 
284
  ostream& os = ops.os();
 
285
  bool   verbose  = iorheo::getverbose(os);
 
286
  bool   clean    = iorheo::getclean(os);
 
287
  bool   execute  = iorheo::getexecute(os);
 
288
  string basename = iorheo::getbasename(os);
 
289
  bool   full     = iorheo::getfull(os);    // all edges & !faces in 3d
 
290
  bool   lattice  = iorheo::getlattice(os); // face & volume lattices
 
291
  bool   reader_on_stdin  = iorheo::getreader_on_stdin(os);
 
292
  string format   = iorheo::getimage_format(os);
 
293
  size_type subdivide = iorheo::getsubdivide(os);
 
294
  if (basename.length() == 0) basename = "output";
 
295
  string filelist;
 
296
 
 
297
  size_type dim     = max (omega.dimension(), size_type(2));
 
298
  size_type map_dim = omega.map_dimension();
 
299
  size_type order   = omega.order();
 
300
  size_type nv      = omega.sizes().ownership_by_dimension[0].size();
 
301
  size_type nedg    = omega.sizes().ownership_by_dimension[1].size();
 
302
  size_type nfac    = omega.sizes().ownership_by_dimension[2].size();
 
303
  size_type nvol    = omega.sizes().ownership_by_dimension[3].size();
 
304
  size_type ne      = omega.sizes().ownership_by_dimension[map_dim].size();
 
305
 
 
306
  if (order == 1) lattice = false;
 
307
  bool show_volumes  = (lattice && full && map_dim == 3);
 
308
  bool show_faces    = ((!full || lattice) && map_dim == 3) || (lattice && map_dim < 3);
 
309
  bool show_edges    = ((full && map_dim == 3) || (map_dim < 3)) || (order > 1);
 
310
  bool show_vertices = (map_dim <= 1);
 
311
  bool show_domains  = true;
 
312
  subdivide = std::max(omega.order(), subdivide);
 
313
 
 
314
  string outfile_fmt = "";
 
315
  string tmp = get_tmpdir() + "/";
 
316
  if (!clean) tmp = "";
 
317
  string filename = tmp+basename + ".plot";
 
318
  filelist = filelist + " " + filename;
 
319
  ofstream plot (filename.c_str());
 
320
  ofstream gdat;
 
321
  if (verbose) clog << "! file \"" << filename << "\" created." << endl;
 
322
 
 
323
  basis subdivide_pointset ("P"+itos(subdivide));
 
324
  basis_on_nodal_basis pointset (subdivide_pointset, omega.get_piola_basis());
 
325
 
 
326
  plot << "#!gnuplot" << endl
 
327
       << setprecision(numeric_limits<T>::digits10);
 
328
  if (format != "") {
 
329
    outfile_fmt = basename + "." + format;
 
330
    string terminal = format;
 
331
    if (terminal == "ps")  {
 
332
      terminal = "postscript eps";
 
333
      if (color) terminal += " color";
 
334
    }
 
335
    if (terminal == "jpg") terminal = "jpeg";
 
336
    if (terminal == "jpeg" || terminal == "png" || terminal == "gif") {
 
337
      terminal += " crop";
 
338
    }
 
339
    plot << "set terminal " << terminal    << endl
 
340
         << "set output \"" << outfile_fmt << "\"" << endl;
 
341
  }
 
342
  if (format == "") {
 
343
    plot << "set title \"" << basename << ": " << ne << " elements, " << nv << " vertices\"" << endl;
 
344
  }
 
345
  check_macro (omega.dimension() > 0, "unsupported 0d geo gnuplot output");
 
346
  point_basic<T> dx = 0.1*(omega.xmax() - omega.xmin());
 
347
  T dx_max = max(dx[0],max(dx[1],dx[2]));
 
348
  if (dx_max == 0) dx_max = 0.1;
 
349
  dx[0] = max(dx[0],dx_max);
 
350
  if (omega.dimension() >= 2) dx[1] = max(dx[1],dx_max);
 
351
  if (omega.dimension() == 3) dx[2] = max(dx[2],dx_max);
 
352
  point_basic<T> xmin = omega.xmin() - dx;
 
353
  point_basic<T> xmax = omega.xmax() + dx;
 
354
  plot << "set xrange [" << xmin[0] << ":" << xmax[0] << "]" << endl;
 
355
  if (omega.dimension() == 1) {
 
356
      plot << "set yrange [-1:1]" << endl;
 
357
  } 
 
358
  if (omega.dimension() >= 2) {
 
359
    plot << "set yrange [" << xmin[1] << ":" << xmax[1] << "]" << endl;
 
360
  }
 
361
  if (omega.dimension() == 2) {
 
362
    plot << "set size ratio -1  # equal scales" << endl
 
363
         << "set key left Right at graph 1,1" << endl;
 
364
  }
 
365
  if (omega.dimension() == 3) {
 
366
    plot << "set zrange [" << xmin[2] << ":" << xmax[2] << "]" << endl
 
367
         << "set xyplane at " << xmin[2] << endl
 
368
         << "set view equal xyz # equal scales" << endl
 
369
         << "set view 70,120" << endl
 
370
         << "set xlabel \"x\"" << endl
 
371
         << "set ylabel \"y\"" << endl
 
372
         << "set zlabel \"z\"" << endl;
 
373
    if (!full) {
 
374
      plot << "set hidden3d nooffset" << endl;
 
375
    }
 
376
  }
 
377
  if (format != "") {
 
378
    plot << "set nokey" << endl
 
379
         << "set noborder" << endl
 
380
         << "set notics" << endl;
 
381
  }
 
382
  if (omega.dimension() <= 2) {
 
383
      plot << "plot \\" << endl;
 
384
  } else {
 
385
      plot << "splot \\" << endl;
 
386
  }
 
387
  bool first_plot_line = true;
 
388
  //
 
389
  // plot internal volume, as lattice, for high order elements
 
390
  //
 
391
  if (show_volumes) {
 
392
    filename = tmp+basename + "-vol.gdat";
 
393
    if (!first_plot_line) plot << ",\\" << endl; first_plot_line = false;
 
394
    plot << "  \"" << filename << "\" u 1:2:3:(0.0) title \"volumes\" with l lc 3 lw 1";
 
395
 
 
396
    filelist = filelist + " " + filename;
 
397
    gdat.open (filename.c_str());
 
398
    if (verbose) clog << "! file \"" << filename << "\" created." << endl;
 
399
    gdat << setprecision(numeric_limits<T>::digits10);
 
400
    for (size_type ivol = 0; ivol < nvol; ivol++) {
 
401
      const geo_element& K = omega.get_geo_element(3,ivol);
 
402
      put_volume (gdat, K, omega);
 
403
    }
 
404
    gdat.close();
 
405
  }
 
406
  //
 
407
  // plot faces
 
408
  //
 
409
  if (show_faces) {
 
410
    filename = tmp+basename + "-fac.gdat";
 
411
    if (!first_plot_line) plot << ",\\" << endl; first_plot_line = false;
 
412
    plot << "  \"" << filename << "\" title \"faces\" with l lc 3 lw 1";
 
413
 
 
414
    filelist = filelist + " " + filename;
 
415
    gdat.open (filename.c_str());
 
416
    if (verbose) clog << "! file \"" << filename << "\" created." << endl;
 
417
    gdat << setprecision(numeric_limits<T>::digits10);
 
418
    for (size_type ifac = 0; ifac < nfac; ifac++) {
 
419
      const geo_element& F = omega.get_geo_element(2,ifac);
 
420
      put_face (gdat, F, omega, pointset);
 
421
    }
 
422
    gdat.close();
 
423
  }
 
424
  //
 
425
  // plot edges
 
426
  //
 
427
  if (show_edges) {
 
428
    filename = tmp+basename + "-edg.gdat";
 
429
    if (!first_plot_line) plot << ",\\" << endl; first_plot_line = false;
 
430
    plot << "  \"" << filename << "\" title \"edges\" with l lc 1 lw 1.5";
 
431
 
 
432
    filelist = filelist + " " + filename;
 
433
    gdat.open (filename.c_str());
 
434
    if (verbose) clog << "! file \"" << filename << "\" created." << endl;
 
435
    gdat << setprecision(numeric_limits<T>::digits10);
 
436
    for (size_type iedg = 0; iedg < nedg; iedg++) {
 
437
      const geo_element& E = omega.get_geo_element(1,iedg);
 
438
      put_edge (gdat, E, omega, pointset);
 
439
    }
 
440
    gdat.close();
 
441
  }
 
442
  //
 
443
  // plot vertices
 
444
  //
 
445
  if (show_vertices) {
 
446
    if (!first_plot_line) plot << ",\\" << endl; first_plot_line = false;
 
447
    plot << "  \"" << tmp+basename << "-ver.gdat\" title \"vertices\" with p lc 0";
 
448
 
 
449
    filename = tmp+basename + "-ver.gdat";
 
450
    filelist = filelist + " " + filename;
 
451
    gdat.open (filename.c_str());
 
452
    if (verbose) clog << "! file \"" << filename << "\" created." << endl;
 
453
    gdat << setprecision(numeric_limits<T>::digits10);
 
454
    for (size_type iv = 0; iv < nv; iv++) {
 
455
      const geo_element& P = omega.get_geo_element(0,iv);
 
456
      put_vertex (gdat, P, omega);
 
457
    }
 
458
    gdat.close();
 
459
  }
 
460
  //
 
461
  // plot domains, by decreasing dimension order
 
462
  //  3d: faces opaques si (!full) ou edges sinon : couleur = code du domaine
 
463
  //  2d: edges
 
464
  //  1d: vertex
 
465
  if (show_domains) {
 
466
    size_type line_color = 3;
 
467
    for (size_type dim = omega.map_dimension(); dim+1 >= 1; dim--) {
 
468
      for (size_type idom = 0; idom < omega.n_domain_indirect(); idom++) {
 
469
        const domain_indirect_basic<sequential>& dom = omega.get_domain_indirect (idom);
 
470
        if (dom.map_dimension() != dim) continue;
 
471
        if (dom.size() == 0) continue;
 
472
        if (dim == 3) continue; // TODO: only d=0,1,2 domains yet
 
473
  
 
474
        filename = tmp+basename + "-dom-" + dom.name() + ".gdat";
 
475
        if (!first_plot_line) plot << ",\\" << endl; first_plot_line = false;
 
476
        plot << "  \"" << filename + "\" title \"" << dom.name() << "\"";
 
477
        if (dom.map_dimension() == 0) {
 
478
          plot << " with p";
 
479
        } else if (dom.map_dimension() == 1 && omega.dimension() == 1) {
 
480
          plot << " with lp";
 
481
        } else {
 
482
          plot << " with l";
 
483
        }
 
484
        plot << " lc " << line_color << " lw 2";
 
485
        line_color++;
 
486
  
 
487
        filelist = filelist + " " + filename;
 
488
        gdat.open (filename.c_str());
 
489
        if (verbose) clog << "! file \"" << filename << "\" created." << endl;
 
490
        gdat << setprecision(numeric_limits<T>::digits10);
 
491
        for (size_type ioige = 0; ioige < dom.size(); ioige++) {
 
492
        size_type ige = dom.oige(ioige).index();
 
493
          const geo_element& S = omega.get_geo_element (dim, ige);
 
494
          put (gdat, S, omega, pointset);
 
495
        }
 
496
        gdat.close();
 
497
      }
 
498
    }
 
499
  }
 
500
  //
 
501
  // end of plot
 
502
  //
 
503
  plot << endl;
 
504
  if (format == "" && !reader_on_stdin) {
 
505
    plot << "pause -1 \"<return>\"" << endl;
 
506
  }
 
507
  plot.close();
 
508
  //
 
509
  // run gnuplot
 
510
  //
 
511
  int status = 0;
 
512
  string command;
 
513
  if (execute) {
 
514
      command = "gnuplot ";
 
515
      if (reader_on_stdin) command += "-persist ";
 
516
      command += tmp + basename + ".plot";
 
517
      if (verbose) clog << "! " << command << endl;
 
518
      cin.sync();
 
519
      status = system (command.c_str());
 
520
      if (format != "") {
 
521
        check_macro (file_exists (outfile_fmt), "! file \"" << outfile_fmt << "\" creation failed");
 
522
        if (verbose) clog << "! file \"" << outfile_fmt << "\" created" << endl;
 
523
      }
 
524
  }
 
525
  //
 
526
  // clear gnuplot data
 
527
  //
 
528
  if (clean) {
 
529
      command = "/bin/rm -f " + filelist;
 
530
      if (verbose) clog << "! " << command << endl;
 
531
      status = system (command.c_str());
 
532
  }
 
533
  return ops;
 
534
}
 
535
// ----------------------------------------------------------------------------
 
536
// instanciation in library
 
537
// ----------------------------------------------------------------------------
 
538
template odiststream& visu_gnuplot (odiststream& ops, const geo_basic<Float,sequential>& omega);
 
539
 
 
540
} // rheolef namespace