Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field.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 field unix command
22// author: Pierre.Saramito@imag.fr
23// last change: 11 oct 2011. initial version: 7 july 1997
24
25namespace rheolef {
332} // namespace rheolef
333
334# include "rheolef.h"
335# include "rheolef/iofem.h"
336using namespace rheolef;
337using namespace std;
338void usage() {
339 cerr << "field: usage:" << endl
340 << "field "
341 << "{-|file[.field[.gz]]}"
342 << "[-Idir|-I dir] "
343 << "[-min|-max|-get-geo] "
344 << "[-name string] "
345 << "[-[catch]mark string] "
346 << "[-comp string] "
347 << "[-domain string] "
348 << "[-round [float]] "
349 << "[-proj [string]] "
350 << "[-lumped-proj] "
351 << "[-paraview|-gnuplot] "
352 << "[-field|-text|-gmsh|-gmsh-pos|-bamg-bb] "
353 << "[-[no]clean] [-[no]execute] [-[no]verbose] "
354 << "[-color|-gray|-black-and-white|-bw] "
355 << "[-[no]elevation] "
356 << "[-[no]volume] "
357 << "[-[no]fill] "
358 << "[-[no]showlabel] "
359 << "[-label string] "
360 << "[-[no]stereo] "
361 << "[-[no]cut] "
362 << "[-normal x [y [z]]] "
363 << "[-origin x [y [z]]] "
364 << "[-scale float] "
365 << "[-iso[value] float|-noiso[value]] "
366 << "[-n-iso int] "
367 << "[-n-iso-negative int] "
368 << "[-velocity|-deformation] "
369 << "[-image-format string] "
370 << "[-resolution int [int]] "
371 << "[-subdivide int] "
372 << endl;
373 exit (1);
374}
375// extern:
376namespace rheolef {
378 const field_basic<T,sequential>& uh);
381 const point_basic<T>& origin,
382 const point_basic<T>& normal);
383} // namespace rheolef
384
395
401
408
409int main(int argc, char**argv) {
410 //environment_option_type eopt;
411 //eopt.thread_level = environment_option_type::no_thread;
412 //environment rheolef (argc, argv, eopt);
413 environment rheolef (argc, argv);
414 check_macro (communicator().size() == 1, "field: command may be used as mono-process only");
415 // ----------------------------
416 // default options
417 // ----------------------------
418 string filename = "";
419 string name = "output";
420 dout.os() << setbasename(name);
421 string format = "";
422 dout.os() << verbose; bool bverbose = true;
423 std::string mark = "";
424 render_type render = no_render;
427 bool def_fill_opt = false;
428 bool do_iso3d = false;
429 std::string i_comp_name = "";
430 bool do_proj = false;
431 string use_proj_approx = "";
432 bool do_lumped_mass = false;
433 bool do_cut = false;
434 bool do_round = false;
435 std::string reduce_to_domain = "";
436 Float round_prec = sqrt(std::numeric_limits<Float>::epsilon());
437 show_type show = show_none;
438 dout.os() << showlabel;
439 std::string label;
440 // this normal is not so bad for the dirichlet.cc demo on the cube:
441 cout << setnormal(point(-0.015940197423022637, -0.9771157601293953, -0.21211011624358989));
442 cout << setorigin(point(std::numeric_limits<Float>::max()));
443 // ----------------------------
444 // scan the command line
445 // ----------------------------
446 for (int i = 1; i < argc; i++) {
447
448 // general options:
449 if (strcmp (argv[i], "-clean") == 0) dout.os() << clean;
450 else if (strcmp (argv[i], "-noclean") == 0) dout.os() << noclean;
451 else if (strcmp (argv[i], "-execute") == 0) dout.os() << execute;
452 else if (strcmp (argv[i], "-noexecute") == 0) dout.os() << noexecute;
453 else if (strcmp (argv[i], "-verbose") == 0) { bverbose = true; dout.os() << verbose; }
454 else if (strcmp (argv[i], "-noverbose") == 0) { bverbose = false; dout.os() << noverbose; }
455 else if (strcmp (argv[i], "-I") == 0) {
456 if (i+1 == argc) { cerr << "geo -I: option argument missing" << endl; usage(); }
457 append_dir_to_rheo_path (argv[++i]);
458 }
459 else if (argv [i][0] == '-' && argv [i][1] == 'I') { append_dir_to_rheo_path (argv[i]+2); }
460 // output file option:
461 else if (strcmp (argv[i], "-field") == 0) { dout.os() << rheo; render = file_render; }
462 else if (strcmp (argv[i], "-text") == 0) { dout.os() << rheo; render = file_render; }
463 else if (strcmp (argv[i], "-gmsh") == 0) { dout.os() << gmsh; render = file_render; }
464 else if (strcmp (argv[i], "-gmsh-pos") == 0) { dout.os() << gmsh_pos; render = file_render; }
465 else if (strcmp (argv[i], "-bamg-bb") == 0) { dout.os() << bamg; render = file_render; }
466 else if (strcmp (argv[i], "-name") == 0) {
467 if (i+1 == argc) { std::cerr << "field -name: option argument missing" << std::endl; usage(); }
468 name = argv[++i];
469 }
470 else if (strcmp (argv[i], "-label") == 0) {
471 if (i+1 == argc) { std::cerr << "field -label: option argument missing" << std::endl; usage(); }
472 label = argv[++i];
473 dout.os() << setlabel(label);
474 }
475 // render spec:
476 else if (strcmp (argv[i], "-gnuplot") == 0) { dout.os() << gnuplot; render = gnuplot_render; }
477 else if (strcmp (argv[i], "-paraview") == 0) { dout.os() << paraview; render = paraview_render; }
478
479 // filter
480 else if (strcmp (argv[i], "-proj") == 0) {
481 do_proj = true;
482 if (i+1 < argc && argv[i+1][0] != '-') {
483 use_proj_approx = argv[++i];
484 }
485 }
486 else if (strcmp (argv[i], "-lumped-proj") == 0) {
487 do_proj = true;
488 do_lumped_mass = true;
489 use_proj_approx = "P1";
490 }
491 // inquire option
492 else if (strcmp (argv[i], "-min") == 0) { show = show_min; }
493 else if (strcmp (argv[i], "-max") == 0) { show = show_max; }
494 else if (strcmp (argv[i], "-get-geo") == 0) { show = show_geo; }
495
496 // render options:
497 else if (strcmp (argv[i], "-velocity") == 0) { dout.os() << velocity; vector_style = velocity_style; }
498 else if (strcmp (argv[i], "-deformation") == 0) { dout.os() << deformation; vector_style = deformation_style; }
499 else if (strcmp (argv[i], "-fill") == 0) { dout.os() << fill; def_fill_opt = true; }
500 else if (strcmp (argv[i], "-nofill") == 0) { dout.os() << nofill; def_fill_opt = false; }
501 else if (strcmp (argv[i], "-elevation") == 0) { dout.os() << elevation; }
502 else if (strcmp (argv[i], "-noelevation") == 0) { dout.os() << noelevation; }
503 else if (strcmp (argv[i], "-color") == 0) { dout.os() << color; }
504 else if (strcmp (argv[i], "-gray") == 0) { dout.os() << gray; }
505 else if (strcmp (argv[i], "-black-and-white") == 0) { dout.os() << black_and_white; }
506 else if (strcmp (argv[i], "-bw") == 0) { dout.os() << black_and_white; }
507 else if (strcmp (argv[i], "-showlabel") == 0) { dout.os() << showlabel; }
508 else if (strcmp (argv[i], "-noshowlabel") == 0) { dout.os() << noshowlabel; }
509 else if (strcmp (argv[i], "-stereo") == 0) { dout.os() << stereo;
510 if (render != paraview_render) {
511 dout.os() << paraview;
512 render = paraview_render;
513 }
514 }
515 else if (strcmp (argv[i], "-nostereo") == 0) { dout.os() << nostereo; }
516 else if (strcmp (argv[i], "-volume") == 0) { dout.os() << paraview << volume;
517 render = paraview_render; }
518 else if (strcmp (argv[i], "-novolume") == 0) { dout.os() << novolume; }
519 else if (strcmp (argv[i], "-cut") == 0) { do_cut = true; }
520 else if (strcmp (argv[i], "-nocut") == 0) { do_cut = false; }
521 else if (strcmp (argv[i], "-noisovalue") == 0) { dout.os() << noiso; do_iso3d = false; }
522 else if (strcmp (argv[i], "-isovalue") == 0 || strcmp (argv[i], "-iso") == 0) {
523
524 do_iso3d = true;
525 dout.os() << iso;
526 if (i+1 < argc && is_float(argv[i+1])) {
527 Float iso_value = to_float (argv[++i]);
528 dout.os() << setisovalue(iso_value);
529 }
530 } else if (strcmp (argv[i], "-n-iso") == 0) {
531
532 if (i+1 == argc || !isdigit(argv[i+1][0])) usage();
533 size_t idx = atoi (argv[++i]);
534 dout.os() << setn_isovalue(idx);
535
536 } else if (strcmp (argv[i], "-n-iso-negative") == 0) {
537
538 if (i+1 == argc || !isdigit(argv[i+1][0])) usage();
539 size_t idx = atoi (argv[++i]);
540 dout.os() << setn_isovalue_negative(idx);
541
542 } else if (strcmp (argv[i], "-scale") == 0) {
543
544 if (i+1 == argc || !is_float(argv[i+1])) usage();
545 Float scale = to_float (argv[++i]);
546 cout << setvectorscale (scale);
547
548 } else if (strcmp (argv[i], "-image-format") == 0) {
549 if (i == argc-1) {
550 cerr << "field -image-format: option argument missing" << endl;
551 usage();
552 }
553 string format = argv[++i];
554 if (format == "jpeg") format = "jpg";
555 if (format == "postscript") format = "ps";
556 dout.os() << setimage_format(format);
557 }
558 else if (strcmp (argv[i], "-resolution") == 0) {
559 if (i == argc-1 || !isdigit(argv[i+1][0])) { std::cerr << "geo -resolution: option argument missing" << std::endl; usage(); }
560 size_t nx = atoi(argv[++i]);
561 size_t ny = (i < argc-1 && isdigit(argv[i+1][0])) ? atoi(argv[++i]) : nx;
562 dout.os() << setresolution(point_basic<size_t>(nx,ny));
563 }
564 else if (strcmp (argv[i], "-mark") == 0 || strcmp (argv[i], "-catch") == 0 || strcmp (argv[i], "-catchmark") == 0) {
565 if (i == argc-1) {
566 cerr << "field -mark: option argument missing" << endl;
567 usage();
568 }
569 mark = argv[++i];
570 }
571 else if (strcmp (argv[i], "-subdivide") == 0) {
572 if (i == argc-1) { cerr << "field -subdivide: option argument missing" << endl; usage(); }
573 size_t nsub = atoi(argv[++i]);
574 dout.os() << setsubdivide (nsub);
575 }
576 else if (strcmp (argv[i], "-comp") == 0) {
577
578 if (i+1 == argc || !isdigit(argv[i+1][0])) usage();
579 i_comp_name = argv[++i];
580 }
581 else if (strcmp (argv[i], "-domain") == 0) {
582
583 if (i+1 == argc) usage();
584 reduce_to_domain = argv[++i];
585 }
586 else if (strcmp (argv[i], "-round") == 0) {
587
588 do_round = true;
589 if (i+1 < argc && is_float(argv[i+1])) {
590 round_prec = to_float (argv[++i]);
591 }
592 }
593 else if ((strcmp (argv[i], "-origin") == 0) || (strcmp (argv[i], "-normal") == 0)) {
594
595 point x;
596 unsigned int io = i;
597 if (i+1 == argc || !is_float(argv[i+1])) {
598 warning_macro ("invalid argument to `" << argv[i] << "'");
599 usage();
600 }
601 x[0] = to_float (argv[++i]);
602 if (i+1 < argc && is_float(argv[i+1])) {
603 x[1] = to_float (argv[++i]);
604 if (i+1 < argc && is_float(argv[i+1])) {
605 x[2] = to_float (argv[++i]);
606 }
607 }
608 if (strcmp (argv[io], "-origin") == 0) {
609 cout << setorigin(x);
610 } else {
611 cout << setnormal(x);
612 }
613 }
614 // input options:
615 else if (strcmp (argv[i], "-") == 0) {
616
617 // field on stdin
618 filename = "-";
619 }
620 else if (argv[i][0] != '-') {
621 // input field on file
622 filename = argv[i];
623 }
624 else {
625 cerr << "field: unknown option `" << argv[i] << "'" << endl;
626 usage();
627 }
628 }
629 // ----------------------------
630 // field treatment
631 // ----------------------------
632 if (filename == "") {
633 cerr << "field: no input file specified" << endl;
634 usage();
635 } else if (filename == "-") {
636 // field on stdin
637 std::string thename;
638 if (name != "output") thename = name;
639 std::cin >> setbasename(thename);
640 if (mark != "") din >> catchmark(mark);
641 dout.os() << setbasename(name) << reader_on_stdin;
642 din >> uh;
643 } else {
644 idiststream ids;
645 ids.open (filename, "field");
646 check_macro(ids.good(), "\"" << filename << "[.field[.gz]]\" not found.");
647 if (mark != "") ids >> catchmark(mark);
648 ids >> uh;
649 std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "field");
650 name = get_basename (root_name);
651 dout.os() << setbasename(name);
652 }
653 if (uh.valued_tag() == space_constant::vector) {
654 if (vector_style == velocity_style && render != paraview_render) {
655 // gnuplot does not support yet velocity style: requieres paraview
656 dout.os() << paraview; render = paraview_render;
657 }
658 }
659 if (render == no_render) {
660 // try to choose the best render from dimension
661 if (uh.get_geo().map_dimension() >= 2 ||
662 uh.get_geo().dimension() == 3) {
663 dout.os() << paraview;
664 } else {
665 dout.os() << gnuplot;
666 }
667 }
668 if (uh.get_geo().map_dimension() == 3) {
669 // default 3D is iso+cut and nofill; default 2D is nocut and fill...
670 if (!def_fill_opt) dout.os() << nofill;
671 }
672 if (reduce_to_domain != "") {
673 geo_basic<Float,sequential> dom = uh.get_geo()[reduce_to_domain];
675 if (uh.get_space().get_basis().is_discontinuous() && dom.map_dimension()+1 == uh.get_geo().map_dimension()) {
676 // brdy trace of a discontinuous field
678 uh.get_geo().neighbour_guard();
679 vh = interpolate (Wh, uh);
680 } else {
681 vh = uh[reduce_to_domain];
682 }
683 uh = vh;
684 if (render == file_render) {
685 // put field on dout: save also the cutted geo:
686 odiststream geo_out (uh.get_geo().name(), "geo");
687 geo_out << uh.get_geo();
688 }
689 }
690 if (do_proj && uh.get_space().get_basis().have_compact_support_inside_element()) {
692 size_t k = Uh.degree();
693 if (k == 0) k++;
694 std::string approx = (do_lumped_mass || use_proj_approx == "") ? "P1" : use_proj_approx;
695 space_basic<Float,sequential> Vh (uh.get_geo(), approx, uh.valued());
698 integrate_option fopt;
699 fopt.lump = do_lumped_mass;
701 switch (Uh.valued_tag()) {
703 m = integrate(v*vt, fopt);
704 p = integrate(u*vt);
705 break;
707 m = integrate(dot(v,vt), fopt);
708 p = integrate(dot(u,vt));
709 break;
712 m = integrate(ddot(v,vt), fopt);
713 p = integrate(ddot(u,vt));
714 break;
715 default:
716 error_macro ("proj: unexpected valued field: " << Uh.valued());
717 }
720 vh.set_u() = sm.solve((p*uh).u());
721 uh = vh;
722 }
723 if (i_comp_name != "") {
724 // compute uh component
726 switch (uh.valued_tag()) {
728 size_t i_comp = atoi (i_comp_name.c_str());
729 vh = uh[i_comp];
730 break;
731 }
734 // string as "00" "21" "22" etc
735 check_macro (i_comp_name.size() == 2, "invalid component `"<<i_comp_name<<"'");
736 size_t i = i_comp_name[0] - '0';
737 size_t j = i_comp_name[1] - '0';
738 vh = uh(i,j);
739 break;
740 }
741 default: error_macro ("component: unexpected "<<uh.valued()<<"-valued field");
742 }
743 uh = vh;
744 }
745 if (do_cut) {
746 point origin = iofem::getorigin(cout);
747 point normal = iofem::getnormal(cout);
748 uh = paraview_plane_cut (uh, origin, normal);
749 if (render == file_render) {
750 // put field on dout: save also the cutted geo:
751 odiststream geo_out (uh.get_geo().name(), "geo");
752 geo_out << uh.get_geo();
753 }
754 }
755 if (do_round) {
756 uh = round (uh, round_prec);
757 }
758 if (show == show_min) {
759 cout << std::setprecision(std::numeric_limits<Float>::digits10)
760 << uh.min() << endl;
761 return 0;
762 }
763 if (show == show_max) {
764 cout << std::setprecision(std::numeric_limits<Float>::digits10)
765 << uh.max() << endl;
766 return 0;
767 }
768 if (show == show_geo) {
769 cout << uh.get_geo().name() << endl;
770 return 0;
771 }
772 if (do_iso3d && (render == file_render || render == gnuplot_render)) {
773 // explicitely build the isosurface: for file or gnuplot
775 dout << uh_iso_surface;
776 return 0;
777 }
778 // generates label:
779 if (label == "") {
780 string root_filename = iorheo::getbasename(dout.os());
781 label = (root_filename == "output") ? "" : root_filename;
782 if (mark != "") label = mark;
783 if (i_comp_name != "") {
784 label = (label == "") ? "value" : label;
785 label = label + "[" + i_comp_name + "]";
786 }
787 if (reduce_to_domain != "") {
788 label = (label == "") ? "value" : label;
789 label = label + "[" + reduce_to_domain + "]";
790 }
791 dout.os() << setlabel (label);
792 }
793 // final output:
794 dout << uh;
795}
render_type
Definition branch.cc:432
void usage()
Definition field.cc:338
show_type
Definition field.cc:402
@ show_geo
Definition field.cc:406
@ show_max
Definition field.cc:405
@ show_min
Definition field.cc:404
@ show_none
Definition field.cc:403
vector_style_type
Definition field.cc:396
@ velocity_style
Definition field.cc:398
@ deformation_style
Definition field.cc:397
@ undef_vector_style
Definition field.cc:399
render_type
Definition field.cc:385
@ file_render
Definition field.cc:393
@ gnuplot_render
Definition field.cc:387
@ no_render
Definition field.cc:386
@ paraview_render
Definition field.cc:388
@ vtk_render
Definition field.cc:390
@ x3d_render
Definition field.cc:392
@ plotmtv_render
Definition field.cc:389
@ atom_render
Definition field.cc:391
see the Float page for the full documentation
see the communicator page for the full documentation
see the point page for the full documentation
see the catchmark page for the full documentation
Definition catchmark.h:67
see the environment page for the full documentation
valued_type valued_tag() const
Definition field.h:273
std::string get_approx() const
Definition field.h:272
vec< T, M > & set_u()
Definition field.h:284
T min() const
Definition field.h:786
T max() const
Definition field.h:802
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
generic mesh with rerefence counting
Definition geo.h:1089
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
void open(std::string filename, std::string suffix="", const communicator &comm=communicator())
This routine opens a physical input file.
Definition diststream.cc:85
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 bamg
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 gmsh_pos
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 color format format format format format format gmsh
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 delete_suffix(const string &name, const string &suffix)
delete_suffix: see the rheostream page for the full documentation
std::enable_if< details::is_field_expr_affine_homogeneous< Expr >::value, field_basic< typenameExpr::scalar_type, typenameExpr::memory_type > >::type round(const Expr &expr, const T2 &prec)
Definition round.h:37
geo_basic< T, sequential > paraview_extract_isosurface(const field_basic< T, sequential > &uh)
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
field_basic< T, sequential > paraview_plane_cut(const field_basic< T, sequential > &uh, const point_basic< T > &origin, const point_basic< T > &normal)
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
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)
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
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