Rheolef  7.1
an efficient C++ finite element environment
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 
25 namespace rheolef {
332 } // namespace rheolef
333 
334 # include "rheolef.h"
335 # include "rheolef/iofem.h"
336 using namespace rheolef;
337 using namespace std;
338 void 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:
376 namespace rheolef {
378  const field_basic<T,sequential>& uh);
379  template <class T> field_basic<T,sequential> paraview_plane_cut (
380  const field_basic<T,sequential>& uh,
381  const point_basic<T>& origin,
382  const point_basic<T>& normal);
383 } // namespace rheolef
384 
385 typedef enum {
395 
396 typedef enum {
401 
402 typedef enum {
407 } show_type;
408 
409 int 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;
425  vector_style_type vector_style = deformation_style;
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()) {
691  const space_basic<Float,sequential>& Uh = uh.get_space();
692  size_t k = Uh.degree();
693  if (k == 0) k++;
694  std::string approx = (use_proj_approx == "") ? "P" + itos(k) : 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()) {
727  case space_constant::vector: {
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:491
void usage()
Definition: field.cc:338
int main(int argc, char **argv)
Definition: field.cc:409
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 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
Definition: environment.h:115
const space_type & get_space() const
Definition: field.h:300
valued_type valued_tag() const
Definition: field.h:304
std::string get_approx() const
Definition: field.h:303
const geo_type & get_geo() const
Definition: field.h:301
const std::string & valued() const
Definition: field.h:305
T min() const
Definition: field.h:683
vec< T, M > & set_u()
Definition: field.h:312
T max() const
Definition: field.h:699
generic mesh with rerefence counting
Definition: geo.h:1089
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
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
#define error_macro(message)
Definition: dis_macros.h:49
#define warning_macro(message)
Definition: dis_macros.h:53
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
Definition: rheostream.cc:222
geo_basic< T, sequential > paraview_extract_isosurface(const field_basic< T, sequential > &uh)
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
field_basic< T, sequential > paraview_plane_cut(const field_basic< T, sequential > &uh, const point_basic< T > &origin, const point_basic< T > &normal)
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
Definition: interpolate.cc:233
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
std::enable_if< details::is_field_expr_affine_homogeneous< Expr >::value,field_basic< typename Expr::scalar_type,typename Expr::memory_type >>::type round(const Expr &expr, const T2 &prec)
Definition: round.h:37
Float to_float(const string &s)
to_float: see the rheostream page for the full documentation
Definition: rheostream.cc:494
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
rheolef - reference manual
Definition: sphere.icc:25
Definition: leveque.h:25
Float u(const point &x)
Float epsilon