Updated deal.II in Spack (markdown)
[dealii.wiki.git] / 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)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.