Rheolef  7.1
an efficient C++ finite element environment
branch.cc
Go to the documentation of this file.
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 // the branch unix command
22 // author: Pierre.Saramito@imag.fr
23 //
24 
25 namespace rheolef {
345 } // namespace rheolef
346 
347 //
348 // TODO: filter options
349 // ------------
350 // -extract N
351 // ------------
352 // -mask u
353 // -unmask u
354 // -mask-all
355 // -unmask-all
356 //
357 // -select u <==> -mask-all -unmask u
358 //
359 // -show-names
360 // ------------
361 // -norm {-l2|-linf|-h1}
362 // -> normes in direct via gnuplot via branch (geo with dim=0)
363 //
364 // Q. how to combine and synchronize images u(t) and |u(t)| ?
365 // -> field n(t,x), mesh with dim=0, one point p0, value n(t,p0) = |u(t)|
366 // ------------
367 // -energy-norm
368 // sqrt(|u^{n+1}-u^n|/(t^{n+1}-t^n)) : mesure a(u,u), l'energie (norme l2 pour u)
369 // ------------
370 // graphic:
371 // -height : with -topography, interprets scalar as height
372 // -yield value : mask all scalars lesser than
373 // menu:
374 // stop/start
375 // start number
376 // transparency : of the height
377 // edit color map for height : can put uniform color
378 // stop/start
379 
380 // -------------------------------------------------------------
381 // program
382 // -------------------------------------------------------------
383 #include <rheolef.h>
384 #include <rheolef/iofem.h>
385 using namespace rheolef;
386 using namespace std;
387 
388 void usage()
389 {
390  cerr << "branch: usage: branch "
391  << "-|file[.branch[.gz]]"
392  << "[-toc] "
393  << "[-Idir|-I dir] "
394  << "[-name string] "
395  << "[-index int|-extract int] "
396  << "[-[catch]mark string] "
397  << "[-ndigit int] "
398  << "[-proj] "
399  << "[-lumped-proj] "
400  << "[-if {branch,vtk}] "
401  << "[-paraview|-gnuplot] "
402  << "[-vtk|-branch] "
403  << "[-[no]verbose|-[no]clean|-[no]execute] "
404  << "[-color|-gray|-black-and-white|-bw] "
405  << "[-[no]elevation] "
406  << "[-[no]volume] "
407  << "[-[no]fill] "
408  << "[-[no]showlabel] "
409  << "[-label string] "
410  << "[-[no]stereo] "
411  << "[-[no]cut] "
412  << "[-normal x [y [z]]] "
413  << "[-origin x [y [z]]] "
414  << "[-scale float] "
415  << "[-iso[value] float|-noiso[value]] "
416  << "[-n-iso int] "
417  << "[-n-iso-negative int] "
418  << "[-topography filename] "
419  << "[-umin float] "
420  << "[-umax float] "
421  << "[-subdivide int] "
422  << endl;
423  exit (1);
424 }
426 
428 proj (const field_basic<Float,sequential>& uh, bool do_lumped_mass) {
429  const space_basic<Float,sequential>& Uh = uh.get_space();
430  size_t k = Uh.degree();
431  if (k == 0) k++;
432  std::string approx = "P" + itos(k);
433  space_basic<Float,sequential> Vh (uh.get_geo(), approx, uh.valued());
436  integrate_option fopt;
437  fopt.lump = do_lumped_mass;
439  switch (Uh.valued_tag()) {
441  m = integrate(v*vt, fopt);
442  p = integrate(u*vt);
443  break;
445  m = integrate(dot(v,vt), fopt);
446  p = integrate(dot(u,vt));
447  break;
450  m = integrate(ddot(v,vt), fopt);
451  p = integrate(ddot(u,vt));
452  break;
453  default:
454  error_macro ("proj: unexpected valued field: " << Uh.valued());
455  }
458  vh.set_u() = sm.solve((p*uh).u());
459  return vh;
460 }
461 void
462 extract (idiststream& in, odiststream& out, bool do_proj, bool do_lumped_mass, size_type extract_id,const Float& scale_value)
463 {
465  in >> event.header();
466  out << event.header();
467  size_t i = 0;
468  do {
469  in >> catchmark (event.parameter_name());
470  i++;
471  } while (i <= extract_id);
472  if (!in) return;
473 
474  Float t;
475  in >> t;
476  if (!in) return;
477  event.set_parameter (t);
478  for (size_t i = 0; in && i < event.size(); i++) {
479  in >> catchmark (event[i].first) >> event[i].second;
480  }
481  field_basic<Float,sequential> uh = event[0].second;
482  if (do_proj && uh.get_space().get_basis().have_compact_support_inside_element()) {
483  uh = proj (uh, do_lumped_mass);
484  }
485  if (scale_value != Float(1)) {
486  uh *= scale_value;
487  }
488  out << event(t,uh);
489  out << event.finalize();
490 }
491 typedef enum {
497  toc_render
499 void
501  idiststream& in,
502  odiststream& out,
503  bool do_proj,
504  bool do_lumped_mass,
505  bool def_fill_opt,
506  size_type extract_id,
507  const Float& scale_value,
508  const std::pair<Float,Float>& u_range,
509  render_type render)
510 {
511  if (extract_id != numeric_limits<size_type>::max()) {
512  extract (in, out, do_proj, do_lumped_mass, extract_id, scale_value);
513  return;
514  }
515  Float t = 0;
517  if (u_range.first != std::numeric_limits<Float>::max() ||
518  u_range.second != -std::numeric_limits<Float>::max()) {
519  event.set_range (u_range);
520  }
521  in >> event.header();
522  if (render == toc_render) {
523  in.is() >> noverbose;
524  out << setprecision(numeric_limits<Float>::digits10)
525  << "# i " << event.parameter_name() << endl;
526  Float param;
527  for (size_t n = 0; in >> catchmark (event.parameter_name()) >> param; ++n) {
528  out << n << " " << param << endl;
529  }
530  return;
531  }
533  size_type n = 0;
534  while (in >> event) {
535  for (size_t i = 0; i < event.size(); i++) {
536  uh = event[i].second;
537  if (n == 0) {
538  // default 1D graphic render is gnuplot
539  if (uh.get_geo().dimension() == 1 && render == paraview_render) {
540  dout.os() << gnuplot;
541  }
542  // default 3D is iso+cut and nofill; default 2D is nocut and fill...
543  if (uh.get_geo().map_dimension() == 3) {
544 #ifdef TODO
545  if (!def_plane_cut_opt) dout.os() << cut;
546 #endif // TODO
547  if (!def_fill_opt) dout.os() << nofill;
548  } else {
549 #ifdef TODO
550  if (!def_plane_cut_opt) dout.os() << nocut;
551 #endif // TODO
552  if (!def_fill_opt) dout.os() << fill;
553  }
554  out << event.header();
555  }
556  if (do_proj && uh.get_space().get_basis().have_compact_support_inside_element()) {
557  uh = proj (uh, do_lumped_mass); // TODO: avoid recomputing/inverting the mass at all step !
558  }
559  if (scale_value != Float(1)) {
560  uh *= scale_value;
561  }
562  event[i].second = uh;
563  }
564  out << event;
565  n++;
566  }
567  out << event.finalize();
568 }
569 void set_input_format (idiststream& in, std::string input_format)
570 {
571  if (input_format == "vtk") in.is() >> vtk;
572  else if (input_format == "branch") in.is() >> rheo;
573  else {
574  std::cerr << "branch: invalid input format \""<<input_format<<"\"" << std::endl;
575  usage();
576  }
577 }
578 int main(int argc, char**argv)
579 {
580  environment distributed(argc, argv);
581  if (argc <= 1) usage();
582  clog << verbose;
583  dout.os() << noelevation;
584  bool on_stdin = false;
585  bool do_proj = false;
586  bool do_lumped_mass = false;
587  bool do_cut = false;
588  bool def_fill_opt = false;
589  int digits10 = numeric_limits<Float>::digits10;
590  render_type render = paraview_render; dout.os() << paraview;
591  size_type extract_id = numeric_limits<size_type>::max();
592  Float scale_value = 1;
593  string file_name, name, input_format = "branch";
594  std::pair<Float,Float> u_range;
595  u_range.first = std::numeric_limits<Float>::max();
596  u_range.second = -std::numeric_limits<Float>::max();
597  dout.os() << showlabel;
598  // this normal is not so bad for the dirichlet.cc demo on the cube:
599  cout << setnormal(point(-0.015940197423022637, -0.9771157601293953, -0.21211011624358989));
600  cout << setorigin(point(std::numeric_limits<Float>::max()));
601 
602  for (int i = 1; i < argc; i++) {
603 
604  if (strcmp (argv[i], "-ndigit") == 0) { digits10 = atoi(argv[++i]); }
605  else if (strcmp (argv[i], "-toc") == 0) { render = toc_render; }
606  else if (strcmp (argv[i], "-index") == 0 || strcmp (argv[i], "-extract") == 0)
607  { extract_id = atoi(argv[++i]); render = text_render; dout.os() << rheo; }
608  else if (strcmp (argv[i], "-branch") == 0) { render = text_render; dout.os() << rheo; }
609  else if (strcmp (argv[i], "-vtk") == 0) { render = vtk_render; dout.os() << vtk; }
610  else if (strcmp (argv[i], "-gnuplot") == 0) { render = gnuplot_render; dout.os() << gnuplot; }
611  else if (strcmp (argv[i], "-paraview") == 0) { render = paraview_render; dout.os() << paraview; }
612  else if (strcmp (argv[i], "-skipvtk") == 0) { dout.os() << skipvtk; }
613  else if (strcmp (argv[i], "-proj") == 0) { do_proj = true; }
614  else if (strcmp (argv[i], "-lumped-proj") == 0){ do_proj = do_lumped_mass = true; }
615  else if (strcmp (argv[i], "-elevation") == 0) { dout.os() << elevation; }
616  else if (strcmp (argv[i], "-noelevation") == 0) { dout.os() << noelevation; }
617  else if (strcmp (argv[i], "-color") == 0) { dout.os() << color; }
618  else if (strcmp (argv[i], "-gray") == 0) { dout.os() << gray; }
619  else if (strcmp (argv[i], "-black-and-white") == 0) { dout.os() << black_and_white; }
620  else if (strcmp (argv[i], "-bw") == 0) { dout.os() << black_and_white; }
621  else if (strcmp (argv[i], "-showlabel") == 0) { dout.os() << showlabel; }
622  else if (strcmp (argv[i], "-noshowlabel") == 0) { dout.os() << noshowlabel; }
623  else if (strcmp (argv[i], "-fill") == 0) { dout.os() << fill; def_fill_opt = true; }
624  else if (strcmp (argv[i], "-nofill") == 0) { dout.os() << nofill; def_fill_opt = true; }
625  else if (strcmp (argv[i], "-stereo") == 0) { dout.os() << stereo;
626  if (render != paraview_render) {
627  dout.os() << paraview;
628  render = paraview_render;
629  }
630  }
631  else if (strcmp (argv[i], "-nostereo") == 0) { dout.os() << nostereo; }
632  else if (strcmp (argv[i], "-volume") == 0) { dout.os() << paraview << volume;
633  render = paraview_render; }
634  else if (strcmp (argv[i], "-novolume") == 0) { dout.os() << novolume; }
635  else if (strcmp (argv[i], "-cut") == 0) { do_cut = true; }
636  else if (strcmp (argv[i], "-nocut") == 0) { do_cut = false; }
637  else if (strcmp (argv[i], "-umin") == 0) {
638  if (i+1 == argc || !is_float(argv[i+1])) usage();
639  u_range.first = to_float (argv[++i]);
640  } else if (strcmp (argv[i], "-umax") == 0) {
641  if (i+1 == argc || !is_float(argv[i+1])) usage();
642  u_range.second = to_float (argv[++i]);
643  } else if (strcmp (argv[i], "-scale") == 0) {
644  if (i+1 == argc || !is_float(argv[i+1])) usage();
645  scale_value = to_float (argv[++i]);
646  dout.os() << setvectorscale (scale_value);
647  } else if (strcmp (argv[i], "-noisovalue") == 0) {
648  dout.os() << noiso;
649  } else if (strcmp (argv[i], "-isovalue") == 0 || strcmp (argv[i], "-iso") == 0) {
650 
651  dout.os() << iso;
652  if (i+1 < argc && is_float(argv[i+1])) {
653  Float iso_value = to_float (argv[++i]);
654  dout.os() << setisovalue(iso_value);
655  }
656  } else if (strcmp (argv[i], "-n-iso") == 0) {
657 
658  if (i+1 == argc || !isdigit(argv[i+1][0])) usage();
659  size_t idx = atoi (argv[++i]);
660  dout.os() << setn_isovalue(idx);
661 
662  } else if (strcmp (argv[i], "-n-iso-negative") == 0) {
663 
664  if (i+1 == argc || !isdigit(argv[i+1][0])) usage();
665  size_t idx = atoi (argv[++i]);
666  dout.os() << setn_isovalue_negative(idx);
667 
668  } else if (strcmp (argv[i], "-subdivide") == 0) {
669  if (i == argc-1) { cerr << "branch -subdivide: option argument missing" << endl; usage(); }
670  size_t nsub = atoi(argv[++i]);
671  dout.os() << setsubdivide (nsub);
672  } else if (strcmp (argv[i], "-topography") == 0) {
673 
674  if (i+1 == argc) usage();
675  idiststream zin (argv[++i]);
677  zin >> z;
678  dout.os() << settopography(z);
679  }
680  else if (strcmp (argv[i], "-I") == 0) {
681  if (i+1 == argc) { cerr << "geo -I: option argument missing" << endl; usage(); }
682  append_dir_to_rheo_path (argv[++i]);
683  }
684  else if (argv [i][0] == '-' && argv [i][1] == 'I') { append_dir_to_rheo_path (argv[i]+2); }
685  else if (strcmp (argv[i], "-noclean") == 0) clog << noclean;
686  else if (strcmp (argv[i], "-clean") == 0) clog << clean;
687  else if (strcmp (argv[i], "-noexecute") == 0) clog << noexecute;
688  else if (strcmp (argv[i], "-execute") == 0) clog << execute;
689  else if (strcmp (argv[i], "-verbose") == 0) clog << verbose;
690  else if (strcmp (argv[i], "-noverbose") == 0) clog << noverbose;
691  else if ((strcmp(argv[i], "-origin") == 0) || (strcmp (argv[i], "-normal") == 0)) {
692 
693  point x;
694  unsigned int io = i;
695  if (i+1 == argc || !is_float(argv[i+1])) {
696  warning_macro ("invalid argument to `" << argv[i] << "'");
697  usage();
698  }
699  x[0] = to_float (argv[++i]);
700  if (i+1 < argc && is_float(argv[i+1])) {
701  x[1] = to_float (argv[++i]);
702  if (i+1 < argc && is_float(argv[i+1])) {
703  x[2] = to_float (argv[++i]);
704  }
705  }
706  if (strcmp (argv[io], "-origin") == 0) {
707  cout << setorigin(x);
708  } else {
709  cout << setnormal(x);
710 warning_macro("normal="<<x);
711  }
712  } else if (strcmp (argv[i], "-image-format") == 0) {
713  if (i == argc-1) {
714  cerr << "field -image-format: option argument missing" << endl;
715  usage();
716  }
717  string format = argv[++i];
718  dout.os() << setimage_format(format);
719  }
720  else if (strcmp (argv[i], "-resolution") == 0) {
721  if (i == argc-1 || !isdigit(argv[i+1][0])) { std::cerr << "geo -resolution: option argument missing" << std::endl; usage(); }
722  size_t nx = atoi(argv[++i]);
723  size_t ny = (i < argc-1 && isdigit(argv[i+1][0])) ? atoi(argv[++i]) : nx;
724  dout.os() << setresolution(point_basic<size_t>(nx,ny));
725  }
726  else if (argv [i][0] == '-' && argv [i][1] == 'I') {
727 
728  append_dir_to_rheo_path (argv[i]+2);
729  }
730  else if (strcmp (argv[i], "-name") == 0) {
731  if (i+1 == argc) { std::cerr << "field -name: option argument missing" << std::endl; usage(); }
732  name = argv[++i];
733  }
734  else if (strcmp (argv[i], "-label") == 0) {
735  if (i+1 == argc) { std::cerr << "field -label: option argument missing" << std::endl; usage(); }
736  string label = argv[++i];
737  dout.os() << setlabel(label);
738  }
739  else if (strcmp (argv[i], "-catchmark") == 0 || strcmp (argv[i], "-mark") == 0) {
740  if (i+1 == argc) { std::cerr << "field -mark: option argument missing" << std::endl; usage(); }
741  string mark = argv[++i];
742  dout.os() << setmark(mark);
743  }
744  else if (strcmp (argv[i], "-if") == 0 ||
745  strcmp (argv[i], "-input-format") == 0) {
746  if (i == argc-1) { std::cerr << "branch "<<argv[i]<<": option argument missing" << std::endl; usage(); }
747  input_format = argv[++i];
748  }
749  else if (strcmp (argv [i], "-") == 0) {
750 
751  on_stdin = true;
752  dout.os() << setbasename("output") << reader_on_stdin;
753  file_name = "output";
754  }
755  else if (argv [i][0] == '-') {
756  cerr << "branch: invalid option `" << argv[i] << "'" << endl;
757  usage();
758  }
759  else {
760 
761  // input on file
762  string dir_name = get_dirname(argv[i]);
763  prepend_dir_to_rheo_path (dir_name);
764  file_name = get_basename(delete_suffix (delete_suffix (argv[i], "gz"), "branch"));
765  }
766  }
767  if (!on_stdin && file_name == "") {
768  cerr << "branch: no input specified" << endl;
769  usage();
770  }
771  string basename = (name != "") ? name : file_name;
772  dout.os() << setbasename(basename)
773  << setprecision(digits10);
774 
775  if (on_stdin) {
776  set_input_format (din, input_format);
777  put(din,dout, do_proj, do_lumped_mass, def_fill_opt, extract_id, scale_value, u_range, render);
778  } else {
779  idiststream in (file_name, "branch");
780  check_macro(in.good(), "\"" << file_name << "[.branch[.gz]]\" not found.");
781  set_input_format (in, input_format);
782  put(in, dout, do_proj, do_lumped_mass, def_fill_opt, extract_id, scale_value, u_range, render);
783  }
784 }
void usage()
Definition: branch.cc:388
void extract(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, size_type extract_id, const Float &scale_value)
Definition: branch.cc:462
int main(int argc, char **argv)
Definition: branch.cc:578
void set_input_format(idiststream &in, std::string input_format)
Definition: branch.cc:569
field::size_type size_type
Definition: branch.cc:425
render_type
Definition: branch.cc:491
@ gnuplot_render
Definition: branch.cc:496
@ paraview_render
Definition: branch.cc:493
@ vtk_render
Definition: branch.cc:494
@ plotmtv_render
Definition: branch.cc:495
@ toc_render
Definition: branch.cc:497
@ text_render
Definition: branch.cc:492
see the Float page for the full documentation
see the point page for the full documentation
void set_range(const std::pair< T, T > &u_range)
Definition: branch.h:371
const std::string & parameter_name() const
Definition: branch.h:333
__branch_header< T, M > header()
Definition: branch.h:404
see the catchmark page for the full documentation
Definition: catchmark.h:67
see the environment page for the full documentation
Definition: environment.h:104
const space_type & get_space() const
Definition: field.h:300
const geo_type & get_geo() const
Definition: field.h:301
const std::string & valued() const
Definition: field.h:305
vec< T, M > & set_u()
Definition: field.h:312
see the integrate_option page for the full documentation
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
std::ostream & os()
Definition: diststream.h:236
vec< T, M > solve(const vec< T, M > &b) const
Definition: solver.h:345
the finite element space
Definition: space.h:352
point_basic< Float > point
Definition: point.h:164
size_t size_type
Definition: basis_get.cc:76
double Float
see the Float page for the full documentation
Definition: Float.h:143
idiststream din
see the diststream page for the full documentation
Definition: diststream.h:427
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format format format format paraview
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format format gnuplot
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format vtk
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color rheo
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin gray
This file is part of Rheolef.
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:155
string delete_suffix(const string &name, const string &suffix)
delete_suffix: see the rheostream page for the full documentation
Definition: rheostream.cc:222
void prepend_dir_to_rheo_path(const string &dir)
prepend_dir_to_rheo_path: see the rheostream page for the full documentation
Definition: rheostream.cc:339
field_basic< T, M > proj(const field_basic< T, M > &uh, const std::string &approx="P1")
Definition: adapt.cc:39
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
string get_basename(const string &name)
get_basename: see the rheostream page for the full documentation
Definition: rheostream.cc:254
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&! is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result dummy=Result())
see the integrate page for the full documentation
Definition: integrate.h:202
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
bool is_float(const string &s)
is_float: see the rheostream page for the full documentation
Definition: rheostream.cc:476
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
void append_dir_to_rheo_path(const string &dir)
append_dir_to_rheo_path: see the rheostream page for the full documentation
Definition: rheostream.cc:330
string get_dirname(const string &name)
get_dirname: see the rheostream page for the full documentation
Definition: rheostream.cc:263
Float to_float(const string &s)
to_float: see the rheostream page for the full documentation
Definition: rheostream.cc:494
rheolef - reference manual
Definition: sphere.icc:25
Definition: leveque.h:25
Float u(const point &x)