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

« back to all changes in this revision

Viewing changes to nfem/pbin/branch.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
* 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:
 
1
 
 
2
///
 
3
/// This file is part of Rheolef.
 
4
///
 
5
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
6
///
 
7
/// Rheolef is free software; you can redistribute it and/or modify
 
8
/// it under the terms of the GNU General Public License as published by
 
9
/// the Free Software Foundation; either version 2 of the License, or
 
10
/// (at your option) any later version.
 
11
///
 
12
/// Rheolef is distributed in the hope that it will be useful,
 
13
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
14
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
15
/// GNU General Public License for more details.
 
16
///
 
17
/// You should have received a copy of the GNU General Public License
 
18
/// along with Rheolef; if not, write to the Free Software
 
19
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
20
///
 
21
/// =========================================================================
 
22
//
 
23
// TODO: filter options
 
24
// ------------
 
25
//  -extract N
 
26
// ------------
 
27
//  -mask  u
 
28
//  -unmask  u
 
29
//  -mask-all
 
30
//  -unmask-all
 
31
//
 
32
//  -select u <==> -mask-all -unmask u
 
33
//
 
34
//  -show-names
 
35
// ------------
 
36
//  -norm {-l2|-linf|-h1}
 
37
//      -> normes in direct via gnuplot via branch (geo with dim=0)
 
38
//
 
39
//  Q. how to combine and synchronize images u(t) and |u(t)| ?
 
40
//     -> field n(t,x), mesh with dim=0, one point p0, value n(t,p0) = |u(t)| 
 
41
// ------------
 
42
// -energy-norm
 
43
//  sqrt(|u^{n+1}-u^n|/(t^{n+1}-t^n)) : mesure a(u,u), l'energie (norme l2 pour u)
 
44
// ------------
 
45
// graphic:
 
46
//   -height      : with -topography, interprets scalar as height 
 
47
//   -yield value : mask all scalars lesser than
 
48
// menu:
 
49
//   stop/start
 
50
//   start number
 
51
//   transparency : of the height
 
52
//   edit color map for height : can put uniform color
 
53
//   stop/start
 
54
 
 
55
/*Prog:branch
 
56
NAME: @code{branch} -- handle a family of fields
 
57
 
 
58
@cindex animation
 
59
@cindex continuation methods
 
60
@cindex time-dependent problems
 
61
@pindex branch
 
62
@cindex RHEOPATH environment variable
 
63
@fiindex @file{.branch} family of fields 
 
64
 
 
65
SYNOPSIS:
 
66
  @example
 
67
        branch [@var{options}] @var{filename}
 
68
  @end example
 
69
 
 
70
EXAMPLE:
 
71
  Generates vtk file colection for visualization with paraview:
 
72
  @example
 
73
        branch output.branch -paraview
 
74
  @end example
 
75
DESCRIPTION:      
 
76
  Read and output a branch of finite element fields from file,
 
77
  in field text file format.
 
78
INPUT FILE SPECIFICATION:
 
79
    @table @code
 
80
 
 
81
    @itemx -I@var{dir}
 
82
        add @var{dir} to the RHEOPATH search path.
 
83
        See also @ref{geo class} for RHEOPATH mechanism. 
 
84
 
 
85
    @itemx @var{filename}
 
86
        specifies the name of the file containing
 
87
        the input field.
 
88
 
 
89
    @itemx -
 
90
        read field on standard input instead on a file.
 
91
 
 
92
    @itemx -ndigit @var{int}
 
93
        Number of digits used to print floating point values
 
94
        when using the @samp{-geo} option.
 
95
        Default depends upon the machine precision associated to the
 
96
        @code{Float} type.
 
97
 
 
98
    @end table
 
99
 
 
100
OUTPUT AND RENDER SPECIFICATION:
 
101
    @table @code
 
102
 
 
103
    @itemx -extract @var{int}
 
104
        Extract the i-th record in the file. The output
 
105
        is a field or multi-field file format.
 
106
 
 
107
    @itemx -branch
 
108
        Output on stdout in @file{.branch} format.
 
109
        This is the default.
 
110
 
 
111
@toindex gnuplot
 
112
@toindex plotmtv
 
113
@toindex paraview
 
114
@toindex vtk
 
115
 
 
116
    @itemx -paraview
 
117
        Generate a collection of @code{vtk} files
 
118
        for using @code{paraview}.
 
119
 
 
120
    @itemx -vtk
 
121
        Generate a single @code{vtk} file
 
122
        with numbered fields.
 
123
 
 
124
    @itemx -gnuplot
 
125
        Run 1d animation using @code{gnuplot}.
 
126
 
 
127
    @itemx -plotmtv
 
128
        This driver is unsupported for animations.
 
129
 
 
130
    @end table
 
131
 
 
132
OTHER OPTIONS:
 
133
    @table @code
 
134
 
 
135
    @itemx -umin @var{float}
 
136
    @itemx -umax @var{float}
 
137
        set the solution range for the @code{gnuplot} driver.
 
138
        By default this range is computed from the first field
 
139
        of the branch, and this could be problematic when this
 
140
        field is initialy zero.
 
141
 
 
142
@cindex topography
 
143
    @itemx -topography @var{filename}[.field[.gz]]
 
144
        performs a tridimensionnal elevation view based
 
145
        on the topographic data.
 
146
 
 
147
@cindex projection
 
148
@apindex P0
 
149
@apindex P1
 
150
    @itemx -proj
 
151
        performs a @code{P1} projection on the fly.
 
152
        This option is useful when rendering @code{P0} data
 
153
        while @code{vtk} render requieres @code{P1} description.
 
154
 
 
155
@cindex elevation
 
156
    @itemx -elevation
 
157
        For  two  dimensional  field, represent values as elevation in the third dimension.
 
158
        This is the default.
 
159
 
 
160
    @itemx -noelevation
 
161
        Prevent from the elevation representation. 
 
162
 
 
163
    @itemx -scale @var{float}
 
164
        applies a multiplicative factor to the field.
 
165
        This is useful e.g. in conjonction with the @code{elevation} option.
 
166
        The default value is 1.
 
167
 
 
168
    @itemx -verbose
 
169
        print messages related to graphic files created and
 
170
       command system calls (this is the default).
 
171
 
 
172
    @itemx -noverbose
 
173
        does not print previous messages.
 
174
 
 
175
    @itemx -clean
 
176
        clear temporary graphic files (this is the default).
 
177
   
 
178
    @itemx -noclean
 
179
        does not clear temporary graphic files.
 
180
 
 
181
    @itemx -execute
 
182
        execute graphic command (this is the default).
 
183
   
 
184
    @itemx -noexecute
 
185
        does not execute graphic command. Generates only
 
186
        graphic files. This is usefull in conjuction with the
 
187
        @code{-noclean} command.
 
188
 
 
189
    @end table
 
190
 
 
191
BRANCH FILE FORMAT:
 
192
 
 
193
@pindex field
 
194
@pindex mfield
 
195
@fiindex @file{.field} field 
 
196
 
 
197
   The @file{.branch} file format bases on the @file{.field} one:
 
198
   @example
 
199
       EXAMPLE          GENERAL FORM
 
200
 
 
201
        #!branch        #!branch
 
202
        branch          branch
 
203
          1 1 11        <version> <nfield=1> <nvalue=N>
 
204
          time u        <key> <field name>
 
205
 
 
206
          #time 3.14    #<key> <key value 1>
 
207
          #u            #<field name>
 
208
          field         <field 1>
 
209
          .....         ....
 
210
 
 
211
          .....         ....
 
212
          #time 6.28    #<key> <key value N>
 
213
          #u            #<field name>
 
214
          field         <field N>
 
215
          .....         ....
 
216
   @end example
 
217
   The key name is here @code{time}, but could be any string (without spaces).
 
218
   The previous example contains one @code{field} at each time step.
 
219
   Labels appears all along the file to facilitate direct jumps and field and step skips.
 
220
 
 
221
   The format supports several fields, such as (t,u(t),p(t)), where u could
 
222
   be a multi-component (e.g. a vector) field:
 
223
   @example
 
224
        #!branch
 
225
        branch
 
226
          1 2 11
 
227
          time u p
 
228
 
 
229
          #time 3.14
 
230
          #u
 
231
          mfield
 
232
          1 2
 
233
          #u0
 
234
          field
 
235
          ...
 
236
          #u1
 
237
          field
 
238
          ...
 
239
          #p
 
240
 
 
241
          #time 6.28
 
242
          ...
 
243
   @end example
 
244
 
 
245
END:*/
 
246
// -------------------------------------------------------------
 
247
// program
 
248
// -------------------------------------------------------------
 
249
#include <rheolef.h>
 
250
#include <rheolef/iofem.h>
 
251
using namespace rheolef;
 
252
using namespace std;
 
253
 
 
254
void usage()
 
255
{
 
256
      cerr << "branch: usage: branch"
 
257
           << " [-ndigit int]"
 
258
           << " {-Igeodir}*"
 
259
           << " [-extract int]"
 
260
           << " [-[no]elevation]"
 
261
           << " [-scale float]"
 
262
           << " [-umin float]"
 
263
           << " [-umax float]"
 
264
           << " [-topography filename]"
 
265
           << " [-proj]"
 
266
           << " [-paraview|-vtk|-gnuplot|-plotmtv]"
 
267
           << " [-[no]verbose|-[no]clean|-[no]execute]"
 
268
           << " -|filename"
 
269
           << endl;
 
270
      exit (1);
 
271
}
 
272
typedef field::size_type size_type;
 
273
 
 
274
void
 
275
extract (idiststream& in, odiststream& out, bool do_proj, size_type extract_id,const Float& scale_value)
 
276
{
 
277
    branch_basic<Float,sequential> event;
 
278
    in  >> event.header();
 
279
    out << event.header();
 
280
    size_t i = 0;
 
281
    do {
 
282
        in >> catchmark (event.parameter_name());
 
283
        i++;
 
284
    } while (i <= extract_id);
 
285
    if (!in) return;
 
286
 
 
287
    Float t;
 
288
    in >> t;
 
289
    if (!in) return;
 
290
    event.set_parameter (t);
 
291
    for (size_t i = 0; in && i < event.size(); i++) {
 
292
       in >> catchmark (event[i].first) >> event[i].second;
 
293
    }
 
294
    field_basic<Float,sequential> u = event[0].second;
 
295
    if (do_proj) {
 
296
            fatal_macro ("-proj: not yet -- sorry");
 
297
#ifdef TODO
 
298
            space_basic<Float,sequential> V = u.get_space();
 
299
            space_basic<Float,sequential> V1 (u.get_geo(), "P1");
 
300
            field_basic<Float,sequential> one (V1, 1);
 
301
            form_basic<Float,sequential> proj (V, V1, "mass");
 
302
            form_diag_basic<Float,sequential> m (V1, "mass");
 
303
            field_basic<Float,sequential> u1(V1);
 
304
            u1.u = (1./m.uu)*(proj.uu*u.u);
 
305
            u = u1;
 
306
#endif // TODO
 
307
    }
 
308
    if (scale_value != Float(1)) {
 
309
      u *= scale_value;
 
310
    }
 
311
    out << event(t,u);
 
312
    out << event.finalize();
 
313
}
 
314
void
 
315
put (
 
316
        idiststream&                  in,
 
317
        odiststream&                  out,
 
318
        bool                          do_proj,
 
319
        size_type                     extract_id,
 
320
        const Float&                  scale_value,
 
321
        const std::pair<Float,Float>& u_range)
 
322
{
 
323
    if (extract_id != numeric_limits<size_type>::max()) {
 
324
      extract (in, out, do_proj, extract_id, scale_value);
 
325
      return;
 
326
    }
 
327
    Float t = 0;
 
328
    field_basic<Float,sequential> u;
 
329
    branch_basic<Float,sequential> event;
 
330
    event.set_range (u_range);
 
331
    size_type n = 0;
 
332
    in  >> event.header();
 
333
    out << event.header();
 
334
    while (in >> event) {
 
335
      for (size_t i = 0; i < event.size(); i++) {
 
336
        u = event[i].second;
 
337
        if (do_proj) {
 
338
            fatal_macro ("-proj: not yet -- sorry");
 
339
#ifdef TODO
 
340
            space_basic<Float,sequential> V = u.get_space();
 
341
            size_t n_comp = u.n_component();
 
342
            space_basic<Float,sequential> V_new_i (u.get_geo(), "P1");
 
343
            space_basic<Float,sequential> V_new   = pow(V_new_i, n_comp);
 
344
            field_basic<Float,sequential> u_new(V_new);
 
345
            field_basic<Float,sequential> one (V_new_i, 1);
 
346
            for (size_t i = 0; i < n_comp; i++) {
 
347
              space_basic<Float,sequential> Vi = V[i];
 
348
              form_basic<Float,sequential> proj (Vi, V_new_i, "mass");
 
349
              form_diag_basic<Float,sequential> m (V_new_i, "mass");
 
350
              field_basic<Float,sequential> u_new_i(V_new_i);
 
351
              u_new_i.u = (1./m.uu)*(proj.uu*u[i].u);
 
352
              u_new[i] = u_new_i;
 
353
            }
 
354
            u = u_new;
 
355
#endif // TODO
 
356
        }
 
357
        if (scale_value != Float(1)) {
 
358
          u *= scale_value;
 
359
        }
 
360
        event[i].second = u;
 
361
      }
 
362
      out << event;
 
363
      n++;
 
364
    }
 
365
    out << event.finalize();
 
366
}
 
367
int main(int argc, char**argv)
 
368
{
 
369
    environment distributed(argc, argv);
 
370
    if (argc <= 1) usage();
 
371
    clog << verbose;
 
372
    dout.os() << elevation;
 
373
    bool on_stdin = false;
 
374
    bool do_proj  = false;
 
375
    int digits10 = numeric_limits<Float>::digits10;
 
376
    typedef enum { text_render, paraview_render, vtk_render, plotmtv_render, gnuplot_render } render_type;
 
377
    render_type render = text_render;
 
378
    size_type extract_id = numeric_limits<size_type>::max();
 
379
    Float scale_value = 1;
 
380
    string file_name;
 
381
    std::pair<Float,Float> u_range;
 
382
    u_range.first  = std::numeric_limits<Float>::min();
 
383
    u_range.second = std::numeric_limits<Float>::max();
 
384
 
 
385
    for (int i = 1; i < argc; i++) {
 
386
 
 
387
        if      (strcmp (argv[i], "-ndigit") == 0)    { digits10 = atoi(argv[++i]); }
 
388
        else if (strcmp (argv[i], "-extract") == 0)   { extract_id = atoi(argv[++i]); }
 
389
        else if (strcmp (argv[i], "-branch") == 0)    { render = text_render; }
 
390
        else if (strcmp (argv[i], "-paraview") == 0)  { render = paraview_render; dout.os() << paraview; }
 
391
        else if (strcmp (argv[i], "-vtk") == 0)       { render = vtk_render; dout.os() << vtk; }
 
392
        else if (strcmp (argv[i], "-plotmtv") == 0)   { render = plotmtv_render; dout.os() << plotmtv; }
 
393
        else if (strcmp (argv[i], "-gnuplot") == 0)   { render = gnuplot_render; dout.os() << gnuplot; }
 
394
        else if (strcmp (argv[i], "-proj") == 0)      { do_proj = true; }
 
395
        else if (strcmp (argv[i], "-elevation") == 0) { dout.os() << elevation; }
 
396
        else if (strcmp (argv[i], "-noelevation") == 0) { dout.os() << noelevation; }
 
397
        else if (strcmp (argv[i], "-umin") == 0)   {
 
398
            if (i+1 == argc || !is_float(argv[i+1])) usage();
 
399
            u_range.first = to_float (argv[++i]);
 
400
        } else if (strcmp (argv[i], "-umax") == 0)   {
 
401
            if (i+1 == argc || !is_float(argv[i+1])) usage();
 
402
            u_range.second = to_float (argv[++i]);
 
403
        } else if (strcmp (argv[i], "-scale") == 0)   {
 
404
            if (i+1 == argc || !is_float(argv[i+1])) usage();
 
405
            scale_value = to_float (argv[++i]);
 
406
        } else if (strcmp (argv[i], "-topography") == 0)   {
 
407
 
 
408
            if (i+1 == argc) usage();
 
409
            idiststream zin (argv[++i]);
 
410
            field_basic<Float,sequential> z;
 
411
            zin >> z;
 
412
            dout.os() << settopography(z);
 
413
        }
 
414
        else if (strcmp (argv[i], "-noclean") == 0)   clog << noclean;
 
415
        else if (strcmp (argv[i], "-clean") == 0)     clog << clean;
 
416
        else if (strcmp (argv[i], "-noexecute") == 0) clog << noexecute;
 
417
        else if (strcmp (argv[i], "-execute") == 0)   clog << execute;
 
418
        else if (strcmp (argv[i], "-verbose") == 0)   clog << verbose;
 
419
        else if (strcmp (argv[i], "-noverbose") == 0) clog << noverbose;
 
420
 
 
421
        else if (argv [i][0] == '-' && argv [i][1] == 'I') {
 
422
 
 
423
            append_dir_to_rheo_path (argv[i]+2);
 
424
        }
 
425
        else if (strcmp (argv [i], "-") == 0) {
 
426
            
 
427
            on_stdin = true;
 
428
            dout.os() << setbasename("output") << reader_on_stdin;
 
429
            file_name = "output";
 
430
        }
 
431
        else {
 
432
 
 
433
            // input on file
 
434
            string dir_name = get_dirname(argv[i]);
 
435
            prepend_dir_to_rheo_path (dir_name);
 
436
            file_name = get_basename(delete_suffix (delete_suffix (argv[i], "gz"), "branch"));
 
437
        }
 
438
    }
 
439
    if (!on_stdin && file_name == "") {
 
440
        cerr << "branch: no input specified" << endl;
 
441
        usage();
 
442
    }
 
443
    dout.os() << setbasename(file_name)
 
444
              << setprecision(digits10);
 
445
 
 
446
    if (on_stdin) {
 
447
        put(din,dout, do_proj, extract_id, scale_value, u_range);
 
448
    } else {
 
449
        idiststream in (file_name, "branch");
 
450
        check_macro(in.good(), "\"" << file_name << "[.branch[.gz]]\" not found.");
 
451
        put(in, dout, do_proj, extract_id, scale_value, u_range);
 
452
    }
 
453
}