Function-Plotting-Tool.md

   1 # A program to convert an analytically known function into a visualization
   2
   3 <!-- No auto-Table of Contents support! -->
   4
   5 ## Introduction
   6
   7 This is a little program that reads the formula of a function from an input file (along with other information) and generates a graphical representation of this function in one of the formats supported by the deal.II DataOutBase class for plotting in one of the widely used visualization programs.
   8 To use this program, you have to cut and paste the `plot.cc` and `Makefile` to a directory of your choice.
   9
  10 Change the Makefile to reflect your installation of deal.II. Then say
  11 ```
  12 make
  13 ```
  14 Run the executable (default is `plot`). This will create a file `parameters.par` that you can edit to change the function, change the output format, etc.
  15
  16
  17 If you call the executable with a command line option, e.g.
  18 ```
  19 test@somewhere.com> ./plot myparameters.par
  20
  21 0.0689890:DEAL:SetUpLogging::Log file output.log
  22 0.0709890:DEAL::Finite Element: FE_Q<2>(1)
  23 0.0709890:DEAL::Generated a 1 valued FunctionParser.
  24 0.0969850:DEAL:CreateMesh::Active Cells: 256
  25 0.0969850:DEAL:CreateMesh::Distributing DOFS
  26 0.263959:DEAL:cg::Starting value 0.0818334
  27 0.269958:DEAL:cg::Convergence step 18 value 3.86807e-17
  28 0.269958:DEAL:OutputSolution::Writing system solution: output.gmv
  29 ```
  30 the program will read `myparameters.par` and write on whatever you specified in this file.
  31
  32
  33
  34 ## The source file
  35
  36 ```
  37 //----------------------------  plot.cc  ---------------------------
  38 //
  39 //    Copyright (C) 2005 by Luca Heltai
  40 //
  41 //    This file is subject to QPL and may not be  distributed
  42 //    without copyright and license information. Please refer
  43 //    to the file deal.II/doc/license.html for the  text  and
  44 //    further information on this license.
  45 //
  46 //----------------------------  plot.cc  ---------------------------
  47
  48 // Base libraries
  49 #include <base/function_parser.h>
  50 #include <base/parameter_handler.h>
  51 #include <base/logstream.h>
  52 #include <base/quadrature_lib.h>
  53
  54 // Grid libraries
  55 #include <grid/tria.h>
  56 #include <grid/grid_in.h>
  57 #include <grid/grid_refinement.h>
  58 #include <grid/grid_generator.h>
  59
  60 // Dofs libraries
  61 #include <dofs/dof_handler.h>
  62 #include <dofs/dof_constraints.h>
  63
  64 // Finite Element libraries.
  65 #include <fe/fe_system.h>
  66 #include <fe/mapping_q.h>
  67 #include <fe/fe_tools.h>
  68
  69 // Linear Algebra Class libraries.
  70 #include <lac/vector.h>
  71
  72 // Numerics libraries
  73 #include <numerics/vectors.h>
  74 #include <numerics/data_out.h>
  75
  76 // C++ libraries
  77 #include <fstream>
  78
  79
  80 #ifdef HAVE_STD_STRINGSTREAM
  81 #  include <sstream>
  82 #  define MY_STR str().c_str()
  83 #  define MY_OSTRSTREAM std::ostringstream
  84 #else
  85 #  include <strstream>
  86 #  define MY_STR str()
  87 #  define MY_OSTRSTREAM  std::ostrstream
  88 #endif
  89
  90 class Parameters : public ParameterHandler
  91 {
  92  public:
  93   Parameters(int, char**);
  94   ~Parameters();
  95
  96   /** Reads the parameters. This method reads the input file
  97       "parameters.par" and set all the variable of this class to
  98       be used later in the problem. */
  99   void set_parameters();
 100
 101   /** Log file name. Used by the deallog class. */
 102   std::string   log_file;
 103   /** Verbosity of console logging. 0 disables it. */
 104   unsigned int  log_console_level;
 105   /** Verbosity of file logging. 0 disables it. */
 106   unsigned int  log_file_level;
 107   /** Write time information. Decides wether to write time information
 108       in the log file. */
 109   bool          log_record_time;
 110   /** Use differential time when loggin.*/
 111   bool          log_differential_time;
 112
 113   /** Read the domain mesh from a file. If set to false, the domain is
 114       the hypercube [bool read_mesh_from_file;
 115
 116   /** The global initial refinement.*/
 117   unsigned int  grid_refinement;
 118
 119   /** The Finite Element to use.*/
 120   std::string     finite_element;
 121
 122   /** The file to read. If there is a command line argument, this file
 123       is read, if it is there, else it is filled with default
 124       values.*/
 125   std::string     file_parameters;
 126
 127   /** File of Solution.*/
 128   std::string     file_solution;
 129
 130   /** Format of the solution. Supported formats:
 131       dx|ucd|gmv|gpl|eps|pov|tec|tecbin|vtk*/
 132   std::string     file_solution_format;
 133
 134   /** Input file for the domain mesh. Supported format: gmsh.*/
 135   std::string     file_domain_mesh;
 136
 137   /** Variables to use in our function. This is the string that
 138       contains the name of the variables used in the function to
 139       plot. */
 140   std::string   variables;
 141   /** Expression of the function to plot. */
 142   std::string   expression;
 143
 144   /** Number of Quadrature points used to project. It has to be
 145       enough.*/
 146   unsigned int n_q_points;
 147 };
 148
 149
 150
 151 /**
 152    Output Object. This class is used for all output related things.
 153 */
 154 template <int dim,
 155           typename VECTOR=Vector<double> >
 156 class OutputObject
 157 {
 158   public:
 159
 160   OutputObject (Parameters * some_par) {
 161     par = some_par;
 162   };
 163
 164   /** Output system Solution. The output format is defined in
 165       according to values set in "parameters.par". \see Parameters
 166       for available options. */
 167   void output_solution (const DoFHandler<dim> &dof_handler,
 168                         const VECTOR &solution,
 169                         unsigned int n_components);
 170
 171   /** Setup Logging.*/
 172   void set_log_file ();
 173
 174   private:
 175   /**
 176       Holds the problem parameters.
 177   */
 178   Parameters * par;
 179
 180   /** Output File. */
 181   std::ofstream log_file;
 182 };
 183
 184
 185 /**
 186  - This class implements a drawer. We make use of a few classes from
 187  - the deal.II library, to show how to output a given function
 188  - interpolated on a given finite element space.
 189  - @author Luca Heltai, 2005
 190  */
 191 template <int dim, typename VECTOR=Vector<double> >
 192 class FunctionDraw {
 193   public:
 194   FunctionDraw(int, char**);
 195   ~FunctionDraw();
 196
 197   void run();
 198
 199   private:
 200   void initialize();
 201   void create_mesh();
 202   void output_function();
 203
 204   Parameters par;
 205
 206   Triangulation<dim> triangulation;
 207
 208   DoFHandler<dim> dof_handler;
 209
 210   OutputObject<dim> out;
 211
 212   VECTOR solution;
 213
 214   SmartPointer<FiniteElement<dim> > fe;
 215
 216   SmartPointer<FunctionParser<dim> > function;
 217
 218 };
 219
 220 /** The initialization of the parameters follows the general structure of the
 221     \class ParameterHandler, with structured input.
 222 */
 223 Parameters::Parameters(int argc, char ** argv)
 224 {
 225   enter_subsection ("Logging");
 226
 227   declare_entry("Log file", "output.log", Patterns::Anything());
 228   declare_entry("Log console level", "5", Patterns::Integer());
 229   declare_entry("Log file level", "3", Patterns::Integer());
 230   declare_entry("Write time information", "true", Patterns::Bool());
 231   declare_entry("Use differential time", "false", Patterns::Bool());
 232
 233   leave_subsection();
 234
 235   enter_subsection ("I/O Options");
 236
 237   declare_entry("Write parameter file", "true", Patterns::Bool());
 238   declare_entry("Parameter file", "parameters.par", Patterns::Anything());
 239
 240   declare_entry("Solution file", "output",  Patterns::Anything());
 241   declare_entry("Solution format", "gmv",
 242                 Patterns::Selection("dx|ucd|gmv|gpl|eps|pov|tec|tecbin|vtk"));
 243   declare_entry("Read domain mesh from file", "false", Patterns::Bool());
 244   declare_entry("Domain coarse mesh gmsh file", "square.msh",
 245                 Patterns::Anything());
 246
 247   leave_subsection();
 248
 249   enter_subsection("Numerical Parameters");
 250   declare_entry("Refinement", "4", Patterns::Integer());
 251   declare_entry ("Finite element", "FE_Q(1)", Patterns::Anything());
 252   declare_entry ("Number of quadrature points", "5", Patterns::Integer());
 253   leave_subsection();
 254
 255   enter_subsection ("Problem Data");
 256   declare_entry ("Variables", "x,y", Patterns::Anything());
 257   declare_entry ("Function to plot", "0", Patterns::Anything());
 258   leave_subsection ();
 259
 260   if(argc > 1)
 261     file_parameters = argv[1](0,1]^dim.*/);
 262   else
 263     file_parameters = "parameters.par";
 264 }
 265
 266 Parameters::~Parameters()
 267 {}
 268
 269 void Parameters::set_parameters ()
 270 {
 271   read_input(file_parameters);
 272
 273   enter_subsection ("Logging");
 274
 275   log_file              = get("Log file");
 276   log_console_level     = get_integer("Log console level");
 277   log_file_level        = get_integer("Log file level");
 278   log_record_time       = get_bool("Write time information");
 279   log_differential_time = get_bool("Use differential time");
 280
 281   leave_subsection();
 282
 283   enter_subsection ("I/O Options");
 284
 285
 286   file_parameters       = get ("Parameter file");
 287
 288   file_solution         = get ("Solution file");
 289   file_solution_format  = get ("Solution format");
 290
 291   read_mesh_from_file   = get_bool("Read domain mesh from file");
 292   file_domain_mesh      = get ("Domain coarse mesh gmsh file");
 293
 294   leave_subsection();
 295
 296   enter_subsection ("Numerical Parameters");
 297   grid_refinement       = get_integer ("Refinement");
 298   finite_element        = get("Finite element");
 299   n_q_points            = get_integer("Number of quadrature points");
 300
 301   leave_subsection();
 302
 303
 304   enter_subsection ("Problem Data");
 305   variables = get("Variables");
 306   expression = get("Function to plot");
 307   leave_subsection ();
 308
 309   std::ofstream out (file_parameters.c_str());
 310   print_parameters(out, Text);
 311
 312 }
 313
 314 template <int dim, typename VECTOR>
 315 void OutputObject<dim,VECTOR>::set_log_file () {
 316    // Initialize Logging to output file.
 317   log_file.open(par->log_file.c_str());
 318   deallog.attach(log_file);
 319   deallog.depth_console(par->log_console_level);
 320   deallog.depth_file(par->log_file_level);
 321   deallog.log_execution_time (par->log_record_time);
 322   deallog.log_time_differences (par->log_differential_time);
 323   deallog.push("SetUpLogging");
 324   deallog << "Log file "
 325           << par->log_file << std::endl;
 326   deallog.pop();
 327 }
 328
 329 template <int dim, typename VECTOR>
 330 void OutputObject<dim,VECTOR>::output_solution(const DoFHandler<dim> &dof_handler,
 331                                                const VECTOR &solution,
 332                                                unsigned int n_components)
 333 {
 334   deallog.push("OutputSolution");
 335
 336   char   filename[sprintf(filename,"%s.%s", par->file_solution.c_str(),
 337           par->file_solution_format.c_str());
 338
 339   deallog << "Writing system solution: " << filename
 340           << std::endl;
 341
 342
 343   std::ofstream output (filename );
 344
 345   DataOut<dim> data_out;
 346
 347   std::vector<std::string> solution_names;
 348
 349   for (unsigned int i=1; i<= n_components; ++i) {
 350     MY_OSTRSTREAM tmp;
 351     tmp << "v_" << i;
 352     solution_names.push_back( tmp.MY_STR );
 353   }
 354
 355   data_out.attach_dof_handler (dof_handler);
 356   data_out.add_data_vector (solution, solution_names);
 357
 358   data_out.build_patches (1);
 359   if (par->file_solution_format == "dx") data_out.write_dx (output);
 360   else if (par->file_solution_format == "ucd" ) data_out.write_ucd (output);
 361   else if (par->file_solution_format == "gpl" ) data_out.write_gnuplot (output);
 362   else if (par->file_solution_format == "eps" ) data_out.write_eps (output);
 363   else if (par->file_solution_format == "pov" ) data_out.write_povray (output);
 364   else if (par->file_solution_format == "tec" ) data_out.write_tecplot (output);
 365   else if (par->file_solution_format == "tecbin" ) data_out.write_tecplot_binary (output);
 366   else if (par->file_solution_format == "vtk" ) data_out.write_vtk (output);
 367   else if (par->file_solution_format == "gmv" ) data_out.write_gmv (output);
 368   else {
 369     Assert(false, ExcInternalError());
 370   }
 371   deallog.pop();
 372 }
 373
 374 template <int dim, typename VECTOR>
 375 FunctionDraw<dim,VECTOR>::FunctionDraw(int argc, char ** argv) :
 376   par(argc, argv),
 377   dof_handler(triangulation),
 378   out(&par)
 379 {}
 380
 381 template <int dim, typename VECTOR>
 382 FunctionDraw<dim,VECTOR>::~FunctionDraw()
 383 {
 384   dof_handler.clear();
 385
 386   FunctionParser<dim> * function_real = function;
 387   function = 0;
 388   delete function_real;
 389
 390   FiniteElement<dim> * fe_real = fe;
 391   fe = 0;
 392   delete fe_real;
 393 }
 394
 395
 396 template <int dim, typename VECTOR>
 397 void FunctionDraw<dim,VECTOR>::initialize()
 398 {
 399   par.set_parameters();
 400   out.set_log_file();
 401
 402   // Creates a new pointer to a FESystem
 403   FiniteElement<dim> *fe_real =
 404     FETools::get_fe_from_name<dim>(par.finite_element);
 405   // Associates a smart pointer to the newly created pointer.
 406   fe = fe_real;
 407
 408   deallog << "Finite Element: " << fe->get_name() << std::endl;
 409
 410   // generate a new function parser
 411   FunctionParser<dim> * function_real =
 412     new FunctionParser<dim>(fe->n_components());
 413   // assign it to our function
 414   function = function_real;
 415
 416   deallog << "Generated a " << fe->n_components()
 417           << " valued FunctionParser." << std::endl;
 418
 419
 420   // Constants we want to pass to the function parser
 421   std::map<std::string, double> constants;
 422   constants["pi"](100];) = M_PI;
 423
 424   function->initialize(par.variables,
 425                        par.expression,
 426                        constants);
 427
 428 }
 429
 430 template <int dim, typename VECTOR>
 431 void FunctionDraw<dim,VECTOR>::run()
 432 {
 433   initialize();
 434
 435   create_mesh();
 436
 437   output_function();
 438 }
 439
 440 template <int dim, typename VECTOR>
 441 void FunctionDraw<dim,VECTOR>::create_mesh()
 442 {
 443   deallog.push("CreateMesh");
 444
 445   if(par.read_mesh_from_file) {
 446     GridIn<dim> grid_in;
 447     grid_in.attach_triangulation (triangulation);
 448     std::ifstream input_file(par.file_domain_mesh.c_str());
 449     grid_in.read_msh(input_file);
 450   } else {
 451     GridGenerator::hyper_cube (triangulation, 0, 1);
 452   }
 453
 454   triangulation.refine_global (par.grid_refinement);
 455
 456   deallog << "Active Cells: "
 457           << triangulation.n_active_cells()
 458           << std::endl;
 459
 460   deallog << "Distributing DOFS" << std::endl;
 461
 462   dof_handler.distribute_dofs (*fe);
 463   solution.reinit(dof_handler.n_dofs());
 464
 465   deallog.pop();
 466 }
 467
 468 template <int dim, typename VECTOR>
 469 void FunctionDraw<dim,VECTOR>::output_function()
 470 {
 471   QGauss< dim > quadrature(par.n_q_points);
 472   ConstraintMatrix constraints;
 473   constraints.close();
 474   MappingQ<dim> mapping(1);
 475
 476   VectorTools::project(mapping,
 477                        dof_handler,
 478                        constraints,
 479                        quadrature,
 480                        *function,
 481                        solution);
 482
 483   out.output_solution(dof_handler,
 484                       solution,
 485                       fe->n_components() );
 486
 487 }
 488
 489 int main (int argc, char ** argv)
 490 {
 491   try
 492     {
 493       FunctionDraw<2> draw(argc, argv);
 494       draw.run();
 495     }
 496   catch (std::exception &exc)
 497     {
 498       std::cerr << std::endl << std::endl
 499                 << "----------------------------------------------------"
 500                 << std::endl;
 501       std::cerr << "Exception on processing: " << std::endl
 502                 << exc.what() << std::endl
 503                 << "Aborting!" << std::endl
 504                 << "----------------------------------------------------"
 505                 << std::endl;
 506
 507       return 1;
 508     }
 509   catch (...)
 510     {
 511       std::cerr << std::endl << std::endl
 512                 << "----------------------------------------------------"
 513                 << std::endl;
 514       std::cerr << "Unknown exception!" << std::endl
 515                 << "Aborting!" << std::endl
 516                 << "----------------------------------------------------"
 517                 << std::endl;
 518       return 1;
 519     };
 520
 521   return 0;
 522 }
 523 ```
 524
 525
 526 ## The Makefile
 527
 528 ```
 529
 530 target = plot
 531
 532 debug-mode = on
 533
 534
 535 D = /usr/local/deal.II
 536
 537
 538 clean-up-files = *gmv *gnuplot *gpl *eps **pov
 539
 540
 541
 542
 543 include $D/common/Make.global_options
 544
 545
 546 libs.g   = $(lib-deal2-2d.g) \
 547            $(lib-lac.g)      \
 548            $(lib-base.g)
 549 libs.o   = $(lib-deal2-2d.o) \
 550            $(lib-lac.o)      \
 551            $(lib-base.o)
 552
 553
 554 ifeq ($(debug-mode),on)
 555   libraries = $(target).g.$(OBJEXT) $(libs.g)
 556 else
 557   libraries = $(target).$(OBJEXT) $(libs.o)
 558 endif
 559
 560
 561 $(target) : $(libraries)
 562         @echo ============================ Linking $@
 563         @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
 564
 565
 566 run: $(target)
 567         @echo ============================ Running $<
 568         @./$(target)$(EXEEXT)
 569
 570
 571 clean:
 572         -rm -f **.$(OBJEXT) **~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
 573
 574
 575 ./%.g.$(OBJEXT) :
 576         @echo ==============debug========= $(<F)
 577         @$(CXX) $(CXXFLAGS.g) -c $< -o $@
 578 ./%.$(OBJEXT) :
 579         @echo ==============optimized===== $(<F)
 580         @$(CXX) $(CXXFLAGS.o) -c $< -o $@
 581
 582
 583 .PHONY: run clean
 584
 585
 586 Makefile.dep: $(target).cc Makefile \
 587               $(shell echo $D/**/include/**/**.h)
 588         @echo ============================ Remaking $@
 589         @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
 590                 > Makefile.dep
 591         @if test -s $@ ; then : else rm $@ ; fi
 592
 593
 594 include Makefile.dep
 595
 596 ```
 597
 598
 599 ## An example parameter file
 600
 601 ```
 602 # Listing of Parameters
 603 # ---------------------
 604 subsection I/O Options
 605   set Domain coarse mesh gmsh file = square.msh
 606   set Parameter file               = parameters.par
 607   set Read domain mesh from file   = false
 608   set Solution file                = output
 609   set Solution format              = gmv
 610   set Write parameter file         = true
 611 end
 612
 613
 614 subsection Logging
 615   set Log console level      = 5
 616   set Log file               = output.log
 617   set Log file level         = 3
 618   set Use differential time  = false
 619   set Write time information = true
 620 end
 621
 622
 623 subsection Numerical Parameters
 624   set Finite element              = FE_Q(1)
 625   set Number of quadrature points = 5
 626   set Refinement                  = 4
 627 end
 628
 629
 630 subsection Problem Data
 631   set Function to plot = exp(x*y) # default: 0
 632   set Variables        = x,y
 633 end
 634 </pre>
 635
 636 ## Another example parameter file
 637
 638 <pre>
 639 # Listing of Parameters
 640 # ---------------------
 641 subsection I/O Options
 642   set Domain coarse mesh gmsh file = square.msh
 643   set Parameter file               = parameters.par
 644   set Read domain mesh from file   = false
 645   set Solution file                = output
 646   set Solution format              = gmv
 647   set Write parameter file         = true
 648 end
 649
 650
 651 subsection Logging
 652   set Log console level      = 5
 653   set Log file               = output.log
 654   set Log file level         = 3
 655   set Use differential time  = false
 656   set Write time information = true
 657 end
 658
 659
 660 subsection Numerical Parameters
 661   set Finite element              = FESystem[set Number of quadrature points = 5
 662   set Refinement                  = 4
 663 end
 664
 665
 666 subsection Problem Data
 667   set Function to plot = sin(pi*y);sin(pi*x) # default: 0
 668   set Variables        = x,y
 669 end
 670 ```
 671
 672 --[[User:Luca|Luca]](FE_Q(1)^2]) 08:39, 5 Apr 2005 (CEST)