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