Updated Notes on visualizing high order output (markdown)
[dealii.wiki.git] / Notes-on-visualizing-high-order-output.md
1 ### Motivation
2
3 As of version 9.1, deal.II allows one to generate VTK and VTU output using high-order Lagrange [cells](https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/). As their name suggests, these cells are described by a set of Lagrange points with numerical quantities attached to them. This becomes useful when you are working with high-order elements and/or high-order mappings, since these objects can be represented more accurately. 
4
5 This page provides instructions on how these high-order meshes can be visualized using ParaView (note: this feature was implemented in the version 5.5, older versions will not be able to visualize these meshes).
6
7 ### Creating high-order output
8
9 First, we produced 2D and 3D VTU output of a spherical shell mesh. Internally, the mesh was generated using the `GridGenerator::hyper_shell` function, which attaches a `SphericalManifold` to all cells. Subsequently, the mesh was attached to a `DoFHandler` object with a `FE_Q` element of the order four. A trigonometric function was interpolated to the underlying finite-element space using the `MappingQGeneric` mapping of order four. Finally, we write the mesh along with the finite-element representation of the function to a VTU file using `DataOut::write_vtu` with high-order curved cells enabled.
10
11 ### Visualization
12
13 The produced VTU output was loaded into ParaView v5.7. Even though the mesh is high-order, ParaView uses (bi-/tri-)linear interpolation between the Lagrange points by default. To enable the high-order interpolation, you need to locate the property "Nonlinear Subdivision Level" in the object properties (the easiest way to do this is to search for this field by its name) and set the value to a number larger than one, according to your needs:
14
15 ![Nonlinear Subdivision Level](https://github.com/agrayver/dealii_wiki_imgs/blob/master/subdivision_option.png)
16
17 For our example, the 2D result looks like this:
18
19 ![2D shell](https://github.com/agrayver/dealii_wiki_imgs/blob/master/shell_2d.png)
20
21 Accordingly, here is the 3D mesh:
22
23 ![3D shell](https://github.com/agrayver/dealii_wiki_imgs/blob/master/shell_3d.png)
24
25 You can see now that cells have curved shapes and the quantity within each cell varies as a high-order (specifically, 4-th order) function.  
26
27 ### Extra details
28
29 For more information on deal.II implementation of this feature, see the documentation of the [DataOut::build_patches](https://www.dealii.org/current/doxygen/deal.II/classDataOut.html#a5eb51872b8736849bb7e8d2007fae086) method. Additionally, step-53 and step-65 provide some more examples and details about high-order mappings, their internal representation and output. 
30
31 ### Code
32
33 The images shown above were generated using the VTU files produced by the program below:
34
35 ```c++
36 #include <deal.II/grid/tria.h>
37 #include <deal.II/grid/grid_generator.h>
38 #include <deal.II/numerics/data_out.h>
39 #include <deal.II/dofs/dof_handler.h>
40 #include <deal.II/numerics/vector_tools.h>
41 #include <deal.II/fe/mapping_q_generic.h>
42
43 #include <iostream>
44 #include <fstream>
45 #include <sstream>
46 #include <cmath>
47
48 using namespace dealii;
49
50 template<int dim>
51 class TestFunction: public Function<dim>
52 {
53 public:
54   double value(const Point<dim> &p, const unsigned int component = 0) const
55   {
56     double v = 0;
57     for(int d = 0; d < dim; ++d)
58     {
59       v += cos(2 * numbers::PI * p[d]);
60     }
61
62     return v;
63   }
64 };
65
66 template<int dim>
67 void shell_grid(unsigned fe_order,
68                 unsigned mapping_order,
69                 unsigned output_mesh_order)
70 {
71   FE_Q<dim> fe(fe_order);
72   MappingQGeneric<dim> mapping(mapping_order);
73
74   Triangulation<dim> triangulation;
75   GridGenerator::hyper_shell(triangulation, Point<dim>(), 0.5, 1., 0, true);
76
77   for(unsigned n = 0; n < 2; ++n)
78   {
79     for(auto cell: triangulation.active_cell_iterators())
80     {
81       for(unsigned face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
82         if(cell->face(face)->at_boundary() &&
83            cell->face(face)->boundary_id() == 1)
84           cell->set_refine_flag();
85     }
86     triangulation.execute_coarsening_and_refinement();
87   }
88
89   DoFHandler<dim> dof_handler(triangulation);
90   dof_handler.distribute_dofs(fe);
91
92   TestFunction<dim> function;
93   Vector<double> vec(dof_handler.n_dofs());
94   VectorTools::interpolate(mapping, dof_handler, function, vec);
95
96   DataOutBase::VtkFlags flags;
97   flags.write_higher_order_cells = true;
98
99   std::stringstream ss;
100   ss << "shell_dim=" << dim
101      << "_p=" << fe_order
102      << "_mapping=" << mapping_order
103      << "_n=" << output_mesh_order
104      << ".vtu";
105
106   std::ofstream out(ss.str());
107   DataOut<dim> data_out;
108   data_out.set_flags(flags);
109   data_out.attach_dof_handler(dof_handler);
110   data_out.add_data_vector(vec, "vec");
111   data_out.build_patches(mapping, output_mesh_order, DataOut<dim>::curved_inner_cells);
112   data_out.write_vtu(out);
113 }
114
115 int main()
116 {
117   const unsigned fe_order = 4;
118   const unsigned mapping_order = 4;
119   const unsigned output_mesh_order = 4;
120
121   shell_grid<2>(fe_order, mapping_order, output_mesh_order);
122   shell_grid<3>(fe_order, mapping_order, output_mesh_order);
123 }
124 ```

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.