NavierStokes_TRBDF2_DG.h

Go to the documentation of this file.

 
 * 
 * @endcode
 * 
 * We declare class that describes the boundary conditions and initial one for velocity:
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   class Velocity: public Function<dim> {
 *   public:
 *     Velocity(const double initial_time = 0.0);
 * 
 *     virtual double value(const Point<dim>&  p,
 *                          const unsigned int component = 0) const override;
 * 
 *     virtual void vector_value(const Point<dim>& p,
 *                               Vector<double>&   values) const override;
 *   };
 * 
 * 
 *   template<int dim>
 *   Velocity<dim>::Velocity(const double initial_time): Function<dim>(dim, initial_time) {}
 * 
 * 
 *   template<int dim>
 *   double Velocity<dim>::value(const Point<dim>& p, const unsigned int component) const {
 *     AssertIndexRange(component, 3);
 *     if(component == 0) {
 *       const double Um = 1.5;
 *       const double H  = 4.1;
 * 
 *       return 4.0*Um*p(1)*(H - p(1))/(H*H);
 *     }
 *     else
 *       return 0.0;
 *   }
 * 
 * 
 *   template<int dim>
 *   void Velocity<dim>::vector_value(const Point<dim>& p, Vector<double>& values) const {
 *     Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
 * 
 *     for(unsigned int i = 0; i < dim; ++i)
 *       values[i] = value(p, i);
 *   }
 * 
 * 
 * @endcode
 * 
 * We do the same for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   class Pressure: public Function<dim> {
 *   public:
 *     Pressure(const double initial_time = 0.0);
 * 
 *     virtual double value(const Point<dim>&  p,
 *                          const unsigned int component = 0) const override;
 *   };
 * 
 * 
 *   template<int dim>
 *   Pressure<dim>::Pressure(const double initial_time): Function<dim>(1, initial_time) {}
 * 
 * 
 *   template<int dim>
 *   double Pressure<dim>::value(const Point<dim>&  p, const unsigned int component) const {
 *     (void)component;
 *     AssertIndexRange(component, 1);
 * 
 *     return 22.0 - p(0);
 *   }
 * 
 * } // namespace EquationData
 * @endcode
 
 
<a name="ann-navier_stokes_TRBDF2_DG.cc"></a>
<h1>Annotated version of navier_stokes_TRBDF2_DG.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 * /* Author: Giuseppe Orlando, 2022. */
 * 
 * @endcode
 * 
 * We start by including all the necessary deal.II header files and some C++
 * related ones.
 * 
 
 * 
 * 
 * @code
 * #include <deal.II/base/quadrature_lib.h>
 * #include <deal.II/base/multithread_info.h>
 * #include <deal.II/base/thread_management.h>
 * #include <deal.II/base/work_stream.h>
 * #include <deal.II/base/parallel.h>
 * #include <deal.II/base/utilities.h>
 * #include <deal.II/base/conditional_ostream.h>
 * 
 * #include <deal.II/lac/vector.h>
 * #include <deal.II/lac/solver_cg.h>
 * #include <deal.II/lac/precondition.h>
 * #include <deal.II/lac/solver_gmres.h>
 * #include <deal.II/lac/affine_constraints.h>
 * 
 * #include <deal.II/grid/tria.h>
 * #include <deal.II/grid/grid_generator.h>
 * #include <deal.II/grid/grid_tools.h>
 * #include <deal.II/grid/grid_refinement.h>
 * #include <deal.II/grid/tria_accessor.h>
 * #include <deal.II/grid/tria_iterator.h>
 * #include <deal.II/distributed/grid_refinement.h>
 * 
 * #include <deal.II/dofs/dof_handler.h>
 * #include <deal.II/dofs/dof_accessor.h>
 * #include <deal.II/dofs/dof_tools.h>
 * 
 * #include <deal.II/fe/fe_q.h>
 * #include <deal.II/fe/fe_dgq.h>
 * #include <deal.II/fe/fe_values.h>
 * #include <deal.II/fe/fe_tools.h>
 * #include <deal.II/fe/fe_system.h>
 * 
 * #include <deal.II/numerics/matrix_tools.h>
 * #include <deal.II/numerics/vector_tools.h>
 * #include <deal.II/numerics/data_out.h>
 * 
 * #include <fstream>
 * #include <cmath>
 * #include <iostream>
 * 
 * #include <deal.II/matrix_free/matrix_free.h>
 * #include <deal.II/matrix_free/operators.h>
 * #include <deal.II/matrix_free/fe_evaluation.h>
 * #include <deal.II/fe/component_mask.h>
 * 
 * #include <deal.II/base/timer.h>
 * #include <deal.II/distributed/solution_transfer.h>
 * #include <deal.II/numerics/error_estimator.h>
 * 
 * #include <deal.II/multigrid/multigrid.h>
 * #include <deal.II/multigrid/mg_transfer_matrix_free.h>
 * #include <deal.II/multigrid/mg_tools.h>
 * #include <deal.II/multigrid/mg_coarse.h>
 * #include <deal.II/multigrid/mg_smoother.h>
 * #include <deal.II/multigrid/mg_matrix.h>
 * 
 * #include <deal.II/meshworker/mesh_loop.h>
 * 
 * #include "runtime_parameters.h"
 * #include "equation_data.h"
 * 
 * @endcode
 * 
 * We include the code in a suitable namespace:
 * 
 
 * 
 * 
 * @code
 * namespace NS_TRBDF2 {
 *   using namespace dealii;
 * 
 * @endcode
 * 
 * The following class is an auxiliary one for post-processing of the vorticity
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   class PostprocessorVorticity: public DataPostprocessor<dim> {
 *   public:
 *     virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector<dim>& inputs,
 *                                        std::vector<Vector<double>>&                computed_quantities) const override;
 * 
 *     virtual std::vector<std::string> get_names() const override;
 * 
 *     virtual std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *     get_data_component_interpretation() const override;
 * 
 *     virtual UpdateFlags get_needed_update_flags() const override;
 *   };
 * 
 * @endcode
 * 
 * This function evaluates the vorticty in both 2D and 3D cases
 *   
 
 * 
 * 
 * @code
 *   template <int dim>
 *   void PostprocessorVorticity<dim>::evaluate_vector_field(const DataPostprocessorInputs::Vector<dim>& inputs,
 *                                                           std::vector<Vector<double>>&                computed_quantities) const {
 *     const unsigned int n_quadrature_points = inputs.solution_values.size();
 * 
 *     /*--- Check the correctness of all data structres ---*/
 *     Assert(inputs.solution_gradients.size() == n_quadrature_points, ExcInternalError());
 *     Assert(computed_quantities.size() == n_quadrature_points, ExcInternalError());
 * 
 *     Assert(inputs.solution_values[0].size() == dim, ExcInternalError());
 * 
 *     if(dim == 2) {
 *       Assert(computed_quantities[0].size() == 1, ExcInternalError());
 *     }
 *     else {
 *       Assert(computed_quantities[0].size() == dim, ExcInternalError());
 *     }
 * 
 *     /*--- Compute the vorticty ---*/
 *     if(dim == 2) {
 *       for(unsigned int q = 0; q < n_quadrature_points; ++q)
 *         computed_quantities[q](0) = inputs.solution_gradients[q][1][0] - inputs.solution_gradients[q][0][1];
 *     }
 *     else {
 *       for(unsigned int q = 0; q < n_quadrature_points; ++q) {
 *         computed_quantities[q](0) = inputs.solution_gradients[q][2][1] - inputs.solution_gradients[q][1][2];
 *         computed_quantities[q](1) = inputs.solution_gradients[q][0][2] - inputs.solution_gradients[q][2][0];
 *         computed_quantities[q](2) = inputs.solution_gradients[q][1][0] - inputs.solution_gradients[q][0][1];
 *       }
 *     }
 *   }
 * 
 * @endcode
 * 
 * This auxiliary function is required by the base class DataProcessor and simply
 * sets the name for the output file
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   std::vector<std::string> PostprocessorVorticity<dim>::get_names() const {
 *     std::vector<std::string> names;
 *     names.emplace_back("vorticity");
 *     if(dim == 3) {
 *       names.emplace_back("vorticity");
 *       names.emplace_back("vorticity");
 *     }
 * 
 *     return names;
 *   }
 * 
 * @endcode
 * 
 * This auxiliary function is required by the base class DataProcessor and simply
 * specifies if the vorticity is a scalar (2D) or a vector (3D)
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *   PostprocessorVorticity<dim>::get_data_component_interpretation() const {
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation> interpretation;
 *     if(dim == 2)
 *       interpretation.push_back(DataComponentInterpretation::component_is_scalar);
 *     else {
 *       interpretation.push_back(DataComponentInterpretation::component_is_part_of_vector);
 *       interpretation.push_back(DataComponentInterpretation::component_is_part_of_vector);
 *       interpretation.push_back(DataComponentInterpretation::component_is_part_of_vector);
 *     }
 * 
 *     return interpretation;
 *   }
 * 
 * @endcode
 * 
 * This auxiliary function is required by the base class DataProcessor and simply
 * sets which variables have to updated (only the gradients)
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   UpdateFlags PostprocessorVorticity<dim>::get_needed_update_flags() const {
 *     return update_gradients;
 *   }
 * 
 * 
 * @endcode
 * 
 * The following structs are auxiliary objects for mesh refinement. ScratchData simply sets
 * the FEValues object
 *   
 
 * 
 * 
 * @code
 *   template <int dim>
 *   struct ScratchData {
 *     ScratchData(const FiniteElement<dim>& fe,
 *                 const unsigned int        quadrature_degree,
 *                 const UpdateFlags         update_flags): fe_values(fe, QGauss<dim>(quadrature_degree), update_flags) {}
 * 
 *     ScratchData(const ScratchData<dim>& scratch_data): fe_values(scratch_data.fe_values.get_fe(),
 *                                                                  scratch_data.fe_values.get_quadrature(),
 *                                                                  scratch_data.fe_values.get_update_flags()) {}
 *     FEValues<dim> fe_values;
 *   };
 * 
 * 
 * @endcode
 * 
 * CopyData simply sets the cell index
 *   
 
 * 
 * 
 * @code
 *   struct CopyData {
 *     CopyData() : cell_index(numbers::invalid_unsigned_int), value(0.0) {}
 * 
 *     CopyData(const CopyData &) = default;
 * 
 *     unsigned int cell_index;
 *     double       value;
 *   };
 * 
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{ <code>NavierStokesProjectionOperator::NavierStokesProjectionOperator</code> }
 * 
 
 * 
 * The following class sets effecively the weak formulation of the problems for the different stages
 * and for both velocity and pressure.
 * The template parameters are the dimnesion of the problem, the polynomial degree for the pressure,
 * the polynomial degree for the velocity, the number of quadrature points for integrals for the pressure step,
 * the number of quadrature points for integrals for the velocity step, the type of vector for storage and the type
 * of floating point data (in general double or float for preconditioners structures if desired).
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   class NavierStokesProjectionOperator: public MatrixFreeOperators::Base<dim, Vec> {
 *   public:
 *     NavierStokesProjectionOperator();
 * 
 *     NavierStokesProjectionOperator(RunTimeParameters::Data_Storage& data);
 * 
 *     void set_dt(const double time_step);
 * 
 *     void set_TR_BDF2_stage(const unsigned int stage);
 * 
 *     void set_NS_stage(const unsigned int stage);
 * 
 *     void set_u_extr(const Vec& src);
 * 
 *     void vmult_rhs_velocity(Vec& dst, const std::vector<Vec>& src) const;
 * 
 *     void vmult_rhs_pressure(Vec& dst, const std::vector<Vec>& src) const;
 * 
 *     void vmult_grad_p_projection(Vec& dst, const Vec& src) const;
 * 
 *     virtual void compute_diagonal() override;
 * 
 *   protected:
 *     double       Re;
 *     double       dt;
 * 
 *     /*--- Parameters of time-marching scheme ---*/
 *     double       gamma;
 *     double       a31;
 *     double       a32;
 *     double       a33;
 * 
 *     unsigned int TR_BDF2_stage; /*--- Flag to denote at which stage of the TR-BDF2 are ---*/
 *     unsigned int NS_stage;      /*--- Flag to denote at which stage of NS solution inside each TR-BDF2 stage we are
 *                                       (solution of the velocity or of the pressure)---*/
 * 
 *     virtual void apply_add(Vec& dst, const Vec& src) const override;
 * 
 *   private:
 *     /*--- Auxiliary variable for the TR stage
 *           (just to avoid to report a lot of 0.5 and for my personal choice to be coherent with the article) ---*/
 *     const double a21 = 0.5;
 *     const double a22 = 0.5;
 * 
 *     /*--- Penalty method parameters, theta = 1 means SIP, while C_p and C_u are the penalization coefficients ---*/
 *     const double theta_v = 1.0;
 *     const double theta_p = 1.0;
 *     const double C_p = 1.0*(fe_degree_p + 1)*(fe_degree_p + 1);
 *     const double C_u = 1.0*(fe_degree_v + 1)*(fe_degree_v + 1);
 * 
 *     Vec                          u_extr; /*--- Auxiliary variable to update the extrapolated velocity ---*/
 * 
 *     EquationData::Velocity<dim>  vel_boundary_inflow; /*--- Auxiliary variable to impose velocity boundary conditions ---*/
 * 
 *     /*--- The following functions basically assemble the linear and bilinear forms. Their syntax is due to
 *           the base class MatrixFreeOperators::Base ---*/
 *     void assemble_rhs_cell_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                          Vec&                                         dst,
 *                                          const std::vector<Vec>&                      src,
 *                                          const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_rhs_face_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                          Vec&                                         dst,
 *                                          const std::vector<Vec>&                      src,
 *                                          const std::pair<unsigned int, unsigned int>& face_range) const;
 *     void assemble_rhs_boundary_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                              Vec&                                         dst,
 *                                              const std::vector<Vec>&                      src,
 *                                              const std::pair<unsigned int, unsigned int>& face_range) const;
 * 
 *     void assemble_rhs_cell_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                          Vec&                                         dst,
 *                                          const std::vector<Vec>&                      src,
 *                                          const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_rhs_face_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                          Vec&                                         dst,
 *                                          const std::vector<Vec>&                      src,
 *                                          const std::pair<unsigned int, unsigned int>& face_range) const;
 *     void assemble_rhs_boundary_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                              Vec&                                         dst,
 *                                              const std::vector<Vec>&                      src,
 *                                              const std::pair<unsigned int, unsigned int>& face_range) const;
 * 
 *     void assemble_cell_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                      Vec&                                         dst,
 *                                      const Vec&                                   src,
 *                                      const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_face_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                      Vec&                                         dst,
 *                                      const Vec&                                   src,
 *                                      const std::pair<unsigned int, unsigned int>& face_range) const;
 *     void assemble_boundary_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                          Vec&                                         dst,
 *                                          const Vec&                                   src,
 *                                          const std::pair<unsigned int, unsigned int>& face_range) const;
 * 
 *     void assemble_cell_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                      Vec&                                         dst,
 *                                      const Vec&                                   src,
 *                                      const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_face_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                      Vec&                                         dst,
 *                                      const Vec&                                   src,
 *                                      const std::pair<unsigned int, unsigned int>& face_range) const;
 *     void assemble_boundary_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                          Vec&                                         dst,
 *                                          const Vec&                                   src,
 *                                          const std::pair<unsigned int, unsigned int>& face_range) const;
 * 
 *     void assemble_cell_term_projection_grad_p(const MatrixFree<dim, Number>&               data,
 *                                               Vec&                                         dst,
 *                                               const Vec&                                   src,
 *                                               const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_rhs_cell_term_projection_grad_p(const MatrixFree<dim, Number>&               data,
 *                                                   Vec&                                         dst,
 *                                                   const Vec&                                   src,
 *                                                   const std::pair<unsigned int, unsigned int>& cell_range) const;
 * 
 *     void assemble_diagonal_cell_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                               Vec&                                         dst,
 *                                               const unsigned int&                          src,
 *                                               const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_diagonal_face_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                               Vec&                                         dst,
 *                                               const unsigned int&                          src,
 *                                               const std::pair<unsigned int, unsigned int>& face_range) const;
 *     void assemble_diagonal_boundary_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                                   Vec&                                         dst,
 *                                                   const unsigned int&                          src,
 *                                                   const std::pair<unsigned int, unsigned int>& face_range) const;
 * 
 *     void assemble_diagonal_cell_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                               Vec&                                         dst,
 *                                               const unsigned int&                          src,
 *                                               const std::pair<unsigned int, unsigned int>& cell_range) const;
 *     void assemble_diagonal_face_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                               Vec&                                         dst,
 *                                               const unsigned int&                          src,
 *                                               const std::pair<unsigned int, unsigned int>& face_range) const;
 *     void assemble_diagonal_boundary_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                                   Vec&                                         dst,
 *                                                   const unsigned int&                          src,
 *                                                   const std::pair<unsigned int, unsigned int>& face_range) const;
 *   };
 * 
 * 
 * @endcode
 * 
 * We start with the default constructor. It is important for MultiGrid, so it is fundamental
 * to properly set the parameters of the time scheme.
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   NavierStokesProjectionOperator():
 *     MatrixFreeOperators::Base<dim, Vec>(), Re(), dt(), gamma(2.0 - std::sqrt(2.0)), a31((1.0 - gamma)/(2.0*(2.0 - gamma))),
 *                                            a32(a31), a33(1.0/(2.0 - gamma)), TR_BDF2_stage(1), NS_stage(1), u_extr() {}
 * 
 * 
 * @endcode
 * 
 * We focus now on the constructor with runtime parameters storage
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   NavierStokesProjectionOperator(RunTimeParameters::Data_Storage& data):
 *     MatrixFreeOperators::Base<dim, Vec>(), Re(data.Reynolds), dt(data.dt),
 *                                            gamma(2.0 - std::sqrt(2.0)), a31((1.0 - gamma)/(2.0*(2.0 - gamma))),
 *                                            a32(a31), a33(1.0/(2.0 - gamma)), TR_BDF2_stage(1), NS_stage(1), u_extr(),
 *                                            vel_boundary_inflow(data.initial_time) {}
 * 
 * 
 * @endcode
 * 
 * Setter of time-step (called by Multigrid and in case a smaller time-step towards the end is needed)
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   set_dt(const double time_step) {
 *     dt = time_step;
 *   }
 * 
 * 
 * @endcode
 * 
 * Setter of TR-BDF2 stage (this can be known only during the effective execution
 * and so it has to be demanded to the class that really solves the problem)
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   set_TR_BDF2_stage(const unsigned int stage) {
 *     AssertIndexRange(stage, 3);
 *     Assert(stage > 0, ExcInternalError());
 * 
 *     TR_BDF2_stage = stage;
 *   }
 * 
 * 
 * @endcode
 * 
 * Setter of NS stage (this can be known only during the effective execution
 * and so it has to be demanded to the class that really solves the problem)
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   set_NS_stage(const unsigned int stage) {
 *     AssertIndexRange(stage, 4);
 *     Assert(stage > 0, ExcInternalError());
 * 
 *     NS_stage = stage;
 *   }
 * 
 * 
 * @endcode
 * 
 * Setter of extrapolated velocity for different stages
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   set_u_extr(const Vec& src) {
 *     u_extr = src;
 *     u_extr.update_ghost_values();
 *   }
 * 
 * 
 * @endcode
 * 
 * We are in a DG-MatrixFree framework, so it is convenient to compute separately cell contribution,
 * internal faces contributions and boundary faces contributions. We start by
 * assembling the rhs cell term for the velocity.
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_cell_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                   Vec&                                         dst,
 *                                   const std::vector<Vec>&                      src,
 *                                   const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       /*--- We first start by declaring the suitable instances to read the old velocity, the
 *       extrapolated velocity and the old pressure. 'phi' will be used only to submit the result.
 *       The second argument specifies which dof handler has to be used (in this implementation 0 stands for
 *       velocity and 1 for pressure). ---*/
 *       FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0),
 *                                                                    phi_old(data, 0),
 *                                                                    phi_old_extr(data, 0);
 *       FEEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_old_press(data, 1);
 * 
 *       /*--- We loop over the cells in the range ---*/
 *       for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *         /*--- Now we need to assign the current cell to each FEEvaluation object and then to specify which src vector
 *         it has to read (the proper order is clearly delegated to the user, which has to pay attention in the function
 *         call to be coherent). ---*/
 *         phi_old.reinit(cell);
 *         phi_old.gather_evaluate(src[0], true, true); /*--- The 'gather_evaluate' function reads data from the vector.
 *                                                            The second and third parameter specifies if you want to read
 *                                                            values and/or derivative related quantities ---*/
 *         phi_old_extr.reinit(cell);
 *         phi_old_extr.gather_evaluate(src[1], true, false);
 *         phi_old_press.reinit(cell);
 *         phi_old_press.gather_evaluate(src[2], true, false);
 *         phi.reinit(cell);
 * 
 *         /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& u_n                = phi_old.get_value(q);
 *           const auto& grad_u_n           = phi_old.get_gradient(q);
 *           const auto& u_n_gamma_ov_2     = phi_old_extr.get_value(q);
 *           const auto& tensor_product_u_n = outer_product(u_n, u_n_gamma_ov_2);
 *           const auto& p_n                = phi_old_press.get_value(q);
 *           auto p_n_times_identity        = tensor_product_u_n;
 *           p_n_times_identity = 0;
 *           for(unsigned int d = 0; d < dim; ++d)
 *             p_n_times_identity[d][d] = p_n;
 * 
 *           phi.submit_value(1.0/(gamma*dt)*u_n, q); /*--- 'submit_value' contains quantites that we want to test against the
 *                                                           test function ---*/
 *           phi.submit_gradient(-a21/Re*grad_u_n + a21*tensor_product_u_n + p_n_times_identity, q);
 *           /*--- 'submit_gradient' contains quantites that we want to test against the gradient of test function ---*/
 *         }
 *         phi.integrate_scatter(true, true, dst); /*--- 'integrate_scatter' is the responsible of distributing into dst.
 *                                                 The first two boolean parameters specify if we are testing against
 *                                                 the test function and/or its gradient ---*/
 *       }
 *     }
 *     else {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0),
 *                                                                    phi_old(data, 0),
 *                                                                    phi_int(data, 0);
 *       FEEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_old_press(data, 1);
 * 
 *       /*--- We loop over the cells in the range ---*/
 *       for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *         phi_old.reinit(cell);
 *         phi_old.gather_evaluate(src[0], true, true);
 *         phi_int.reinit(cell);
 *         phi_int.gather_evaluate(src[1], true, true);
 *         phi_old_press.reinit(cell);
 *         phi_old_press.gather_evaluate(src[2], true, false);
 *         phi.reinit(cell);
 * 
 *         /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& u_n                      = phi_old.get_value(q);
 *           const auto& grad_u_n                 = phi_old.get_gradient(q);
 *           const auto& u_n_gamma                = phi_int.get_value(q);
 *           const auto& grad_u_n_gamma           = phi_int.get_gradient(q);
 *           const auto& tensor_product_u_n       = outer_product(u_n, u_n);
 *           const auto& tensor_product_u_n_gamma = outer_product(u_n_gamma, u_n_gamma);
 *           const auto& p_n                      = phi_old_press.get_value(q);
 *           auto p_n_times_identity              = tensor_product_u_n;
 *           p_n_times_identity = 0;
 *           for(unsigned int d = 0; d < dim; ++d)
 *             p_n_times_identity[d][d] = p_n;
 * 
 *           phi.submit_value(1.0/((1.0 - gamma)*dt)*u_n_gamma, q);
 *           phi.submit_gradient(a32*tensor_product_u_n_gamma + a31*tensor_product_u_n -
 *                               a32/Re*grad_u_n_gamma - a31/Re*grad_u_n + p_n_times_identity, q);
 *         }
 *         phi.integrate_scatter(true, true, dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The followinf function assembles rhs face term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_face_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                   Vec&                                         dst,
 *                                   const std::vector<Vec>&                      src,
 *                                   const std::pair<unsigned int, unsigned int>& face_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. In this case
 *       we are at the face between two elements and this is the reason of 'FEFaceEvaluation'. It contains an extra
 *       input argument, the second one, that specifies if it is from 'interior' or not---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi_p(data, true, 0),
 *                                                                        phi_m(data, false, 0),
 *                                                                        phi_old_p(data, true, 0),
 *                                                                        phi_old_m(data, false, 0),
 *                                                                        phi_old_extr_p(data, true, 0),
 *                                                                        phi_old_extr_m(data, false, 0);
 *       FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_old_press_p(data, true, 1),
 *                                                                        phi_old_press_m(data, false, 1);
 * 
 *       /*--- We loop over the faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_old_p.reinit(face);
 *         phi_old_p.gather_evaluate(src[0], true, true);
 *         phi_old_m.reinit(face);
 *         phi_old_m.gather_evaluate(src[0], true, true);
 *         phi_old_extr_p.reinit(face);
 *         phi_old_extr_p.gather_evaluate(src[1], true, false);
 *         phi_old_extr_m.reinit(face);
 *         phi_old_extr_m.gather_evaluate(src[1], true, false);
 *         phi_old_press_p.reinit(face);
 *         phi_old_press_p.gather_evaluate(src[2], true, false);
 *         phi_old_press_m.reinit(face);
 *         phi_old_press_m.gather_evaluate(src[2], true, false);
 *         phi_p.reinit(face);
 *         phi_m.reinit(face);
 * 
 *         /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *         for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *           const auto& n_plus                 = phi_p.get_normal_vector(q); /*--- The normal vector is the same
 *                                                                                  for both phi_p and phi_m. If the face is interior,
 *                                                                                  it correspond to the outer normal ---*/
 * 
 *           const auto& avg_grad_u_old         = 0.5*(phi_old_p.get_gradient(q) + phi_old_m.get_gradient(q));
 *           const auto& avg_tensor_product_u_n = 0.5*(outer_product(phi_old_p.get_value(q), phi_old_extr_p.get_value(q)) +
 *                                                     outer_product(phi_old_m.get_value(q), phi_old_extr_m.get_value(q)));
 *           const auto& avg_p_old              = 0.5*(phi_old_press_p.get_value(q) + phi_old_press_m.get_value(q));
 * 
 *           phi_p.submit_value((a21/Re*avg_grad_u_old - a21*avg_tensor_product_u_n)*n_plus - avg_p_old*n_plus, q);
 *           phi_m.submit_value(-(a21/Re*avg_grad_u_old - a21*avg_tensor_product_u_n)*n_plus + avg_p_old*n_plus, q);
 *         }
 *         phi_p.integrate_scatter(true, false, dst);
 *         phi_m.integrate_scatter(true, false, dst);
 *       }
 *     }
 *     else {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi_p(data, true, 0),
 *                                                                        phi_m(data, false, 0),
 *                                                                        phi_old_p(data, true, 0),
 *                                                                        phi_old_m(data, false, 0),
 *                                                                        phi_int_p(data, true, 0),
 *                                                                        phi_int_m(data, false, 0);
 *       FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_old_press_p(data, true, 1),
 *                                                                        phi_old_press_m(data, false, 1);
 * 
 *       /*--- We loop over the faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++ face) {
 *         phi_old_p.reinit(face);
 *         phi_old_p.gather_evaluate(src[0], true, true);
 *         phi_old_m.reinit(face);
 *         phi_old_m.gather_evaluate(src[0], true, true);
 *         phi_int_p.reinit(face);
 *         phi_int_p.gather_evaluate(src[1], true, true);
 *         phi_int_m.reinit(face);
 *         phi_int_m.gather_evaluate(src[1], true, true);
 *         phi_old_press_p.reinit(face);
 *         phi_old_press_p.gather_evaluate(src[2], true, false);
 *         phi_old_press_m.reinit(face);
 *         phi_old_press_m.gather_evaluate(src[2], true, false);
 *         phi_p.reinit(face);
 *         phi_m.reinit(face);
 * 
 *         /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *         for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *           const auto& n_plus                       = phi_p.get_normal_vector(q);
 * 
 *           const auto& avg_grad_u_old               = 0.5*(phi_old_p.get_gradient(q) + phi_old_m.get_gradient(q));
 *           const auto& avg_grad_u_int               = 0.5*(phi_int_p.get_gradient(q) + phi_int_m.get_gradient(q));
 *           const auto& avg_tensor_product_u_n       = 0.5*(outer_product(phi_old_p.get_value(q), phi_old_p.get_value(q)) +
 *                                                           outer_product(phi_old_m.get_value(q), phi_old_m.get_value(q)));
 *           const auto& avg_tensor_product_u_n_gamma = 0.5*(outer_product(phi_int_p.get_value(q), phi_int_p.get_value(q)) +
 *                                                           outer_product(phi_int_m.get_value(q), phi_int_m.get_value(q)));
 *           const auto& avg_p_old                    = 0.5*(phi_old_press_p.get_value(q) + phi_old_press_m.get_value(q));
 * 
 *           phi_p.submit_value((a31/Re*avg_grad_u_old + a32/Re*avg_grad_u_int -
 *                               a31*avg_tensor_product_u_n - a32*avg_tensor_product_u_n_gamma)*n_plus - avg_p_old*n_plus, q);
 *           phi_m.submit_value(-(a31/Re*avg_grad_u_old + a32/Re*avg_grad_u_int -
 *                                a31*avg_tensor_product_u_n - a32*avg_tensor_product_u_n_gamma)*n_plus + avg_p_old*n_plus, q);
 *         }
 *         phi_p.integrate_scatter(true, false, dst);
 *         phi_m.integrate_scatter(true, false, dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The followinf function assembles rhs boundary term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_boundary_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                       Vec&                                         dst,
 *                                       const std::vector<Vec>&                      src,
 *                                       const std::pair<unsigned int, unsigned int>& face_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. Clearly on the boundary
 *       the second argument has to be true. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, true, 0),
 *                                                                        phi_old(data, true, 0),
 *                                                                        phi_old_extr(data, true, 0);
 *       FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_old_press(data, true, 1);
 * 
 *       /*--- We loop over the faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_old.reinit(face);
 *         phi_old.gather_evaluate(src[0], true, true);
 *         phi_old_extr.reinit(face);
 *         phi_old_extr.gather_evaluate(src[1], true, false);
 *         phi_old_press.reinit(face);
 *         phi_old_press.gather_evaluate(src[2], true, false);
 *         phi.reinit(face);
 * 
 *         const auto boundary_id = data.get_boundary_id(face); /*--- Get the id in order to impose the proper boundary condition ---*/
 *         const auto coef_jump   = (boundary_id == 1) ? 0.0 : C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
 *         const double aux_coeff = (boundary_id == 1) ? 0.0 : 1.0;
 * 
 *         /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& n_plus             = phi.get_normal_vector(q);
 * 
 *           const auto& grad_u_old         = phi_old.get_gradient(q);
 *           const auto& tensor_product_u_n = outer_product(phi_old.get_value(q), phi_old_extr.get_value(q));
 *           const auto& p_old              = phi_old_press.get_value(q);
 *           const auto& point_vectorized   = phi.quadrature_point(q);
 *           auto u_int_m                   = Tensor<1, dim, VectorizedArray<Number>>();
 *           if(boundary_id == 0) {
 *             for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
 *               Point<dim> point; /*--- The point returned by the 'quadrature_point' function is not an instance of Point
 *                                       and so it is not ready to be directly used. We need to pay attention to the
 *                                       vectorization ---*/
 *               for(unsigned int d = 0; d < dim; ++d)
 *                 point[d] = point_vectorized[d][v];
 *               for(unsigned int d = 0; d < dim; ++d)
 *                 u_int_m[d][v] = vel_boundary_inflow.value(point, d);
 *             }
 *           }
 *           const auto tensor_product_u_int_m = outer_product(u_int_m, phi_old_extr.get_value(q));
 *           const auto lambda                 = (boundary_id == 1) ? 0.0 : std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
 * 
 *           phi.submit_value((a21/Re*grad_u_old - a21*tensor_product_u_n)*n_plus - p_old*n_plus +
 *                            a22/Re*2.0*coef_jump*u_int_m -
 *                            aux_coeff*a22*tensor_product_u_int_m*n_plus + a22*lambda*u_int_m, q);
 *           phi.submit_normal_derivative(-aux_coeff*theta_v*a22/Re*u_int_m, q); /*--- This is equivalent to multiply to the gradient
 *                                                                                     with outer product and use 'submit_gradient' ---*/
 *         }
 *         phi.integrate_scatter(true, true, dst);
 *       }
 *     }
 *     else {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, true, 0),
 *                                                                        phi_old(data, true, 0),
 *                                                                        phi_int(data, true, 0),
 *                                                                        phi_int_extr(data, true, 0);
 *       FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_old_press(data, true, 1);
 * 
 *       /*--- We loop over the faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++ face) {
 *         phi_old.reinit(face);
 *         phi_old.gather_evaluate(src[0], true, true);
 *         phi_int.reinit(face);
 *         phi_int.gather_evaluate(src[1], true, true);
 *         phi_old_press.reinit(face);
 *         phi_old_press.gather_evaluate(src[2], true, false);
 *         phi_int_extr.reinit(face);
 *         phi_int_extr.gather_evaluate(src[3], true, false);
 *         phi.reinit(face);
 * 
 *         const auto boundary_id = data.get_boundary_id(face);
 *         const auto coef_jump   = (boundary_id == 1) ? 0.0 : C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
 *         const double aux_coeff = (boundary_id == 1) ? 0.0 : 1.0;
 * 
 *         /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& n_plus                   = phi.get_normal_vector(q);
 * 
 *           const auto& grad_u_old               = phi_old.get_gradient(q);
 *           const auto& grad_u_int               = phi_int.get_gradient(q);
 *           const auto& tensor_product_u_n       = outer_product(phi_old.get_value(q), phi_old.get_value(q));
 *           const auto& tensor_product_u_n_gamma = outer_product(phi_int.get_value(q), phi_int.get_value(q));
 *           const auto& p_old                    = phi_old_press.get_value(q);
 *           const auto& point_vectorized         = phi.quadrature_point(q);
 *           auto u_m                             = Tensor<1, dim, VectorizedArray<Number>>();
 *           if(boundary_id == 0) {
 *             for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
 *               Point<dim> point;
 *               for(unsigned int d = 0; d < dim; ++d)
 *                 point[d] = point_vectorized[d][v];
 *               for(unsigned int d = 0; d < dim; ++d)
 *                 u_m[d][v] = vel_boundary_inflow.value(point, d);
 *             }
 *           }
 *           const auto tensor_product_u_m = outer_product(u_m, phi_int_extr.get_value(q));
 *           const auto lambda             = (boundary_id == 1) ? 0.0 : std::abs(scalar_product(phi_int_extr.get_value(q), n_plus));
 * 
 *           phi.submit_value((a31/Re*grad_u_old + a32/Re*grad_u_int -
 *                            a31*tensor_product_u_n - a32*tensor_product_u_n_gamma)*n_plus - p_old*n_plus +
 *                            a33/Re*2.0*coef_jump*u_m -
 *                            aux_coeff*a33*tensor_product_u_m*n_plus + a33*lambda*u_m, q);
 *           phi.submit_normal_derivative(-aux_coeff*theta_v*a33/Re*u_m, q);
 *         }
 *         phi.integrate_scatter(true, true, dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Put together all the previous steps for velocity. This is done automatically by the loop function of 'MatrixFree' class
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   vmult_rhs_velocity(Vec& dst, const std::vector<Vec>& src) const {
 *     for(unsigned int d = 0; d < src.size(); ++d)
 *       src[d].update_ghost_values();
 * 
 *     this->data->loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_velocity,
 *                      &NavierStokesProjectionOperator::assemble_rhs_face_term_velocity,
 *                      &NavierStokesProjectionOperator::assemble_rhs_boundary_term_velocity,
 *                      this, dst, src, true,
 *                      MatrixFree<dim, Number>::DataAccessOnFaces::unspecified,
 *                      MatrixFree<dim, Number>::DataAccessOnFaces::unspecified);
 *   }
 * 
 * 
 * @endcode
 * 
 * Now we focus on computing the rhs for the projection step for the pressure with the same ratio.
 * The following function assembles rhs cell term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_cell_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                   Vec&                                         dst,
 *                                   const std::vector<Vec>&                      src,
 *                                   const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     /*--- We first start by declaring the suitable instances to read already available quantities.
 *           The third parameter specifies that we want to use the second quadrature formula stored. ---*/
 *     FEEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number>   phi(data, 1, 1),
 *                                                                  phi_old(data, 1, 1);
 *     FEEvaluation<dim, fe_degree_v, n_q_points_1d_p, dim, Number> phi_proj(data, 0, 1);
 * 
 *     const double coeff   = (TR_BDF2_stage == 1) ? 1.0e6*gamma*dt*gamma*dt : 1.0e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
 * 
 *     const double coeff_2 = (TR_BDF2_stage == 1) ? gamma*dt : (1.0 - gamma)*dt;
 * 
 *     /*--- We loop over cells in the range ---*/
 *     for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *       phi_proj.reinit(cell);
 *       phi_proj.gather_evaluate(src[0], true, false);
 *       phi_old.reinit(cell);
 *       phi_old.gather_evaluate(src[1], true, false);
 *       phi.reinit(cell);
 * 
 *       /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *       for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *         const auto& u_star_star = phi_proj.get_value(q);
 *         const auto& p_old       = phi_old.get_value(q);
 * 
 *         phi.submit_value(1.0/coeff*p_old, q);
 *         phi.submit_gradient(1.0/coeff_2*u_star_star, q);
 *       }
 *       phi.integrate_scatter(true, true, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles rhs face term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_face_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                   Vec&                                         dst,
 *                                   const std::vector<Vec>&                      src,
 *                                   const std::pair<unsigned int, unsigned int>& face_range) const {
 *     /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *     FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number>   phi_p(data, true, 1, 1),
 *                                                                      phi_m(data, false, 1, 1);
 *     FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_p, dim, Number> phi_proj_p(data, true, 0, 1),
 *                                                                      phi_proj_m(data, false, 0, 1);
 * 
 *     const double coeff = (TR_BDF2_stage == 1) ? 1.0/(gamma*dt) : 1.0/((1.0 - gamma)*dt);
 * 
 *     /*--- We loop over faces in the range ---*/
 *     for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *       phi_proj_p.reinit(face);
 *       phi_proj_p.gather_evaluate(src[0], true, false);
 *       phi_proj_m.reinit(face);
 *       phi_proj_m.gather_evaluate(src[0], true, false);
 *       phi_p.reinit(face);
 *       phi_m.reinit(face);
 * 
 *       /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *       for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *         const auto& n_plus           = phi_p.get_normal_vector(q);
 *         const auto& avg_u_star_star  = 0.5*(phi_proj_p.get_value(q) + phi_proj_m.get_value(q));
 * 
 *         phi_p.submit_value(-coeff*scalar_product(avg_u_star_star, n_plus), q);
 *         phi_m.submit_value(coeff*scalar_product(avg_u_star_star, n_plus), q);
 *       }
 *       phi_p.integrate_scatter(true, false, dst);
 *       phi_m.integrate_scatter(true, false, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles rhs boundary term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_boundary_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                       Vec&                                         dst,
 *                                       const std::vector<Vec>&                      src,
 *                                       const std::pair<unsigned int, unsigned int>& face_range) const {
 *     /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *     FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number>   phi(data, true, 1, 1);
 *     FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_p, dim, Number> phi_proj(data, true, 0, 1);
 * 
 *     const double coeff = (TR_BDF2_stage == 1) ? 1.0/(gamma*dt) : 1.0/((1.0 - gamma)*dt);
 * 
 *     /*--- We loop over faces in the range ---*/
 *     for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *       phi_proj.reinit(face);
 *       phi_proj.gather_evaluate(src[0], true, false);
 *       phi.reinit(face);
 * 
 *       /*--- Now we loop over all the quadrature points to compute the integrals ---*/
 *       for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *         const auto& n_plus = phi.get_normal_vector(q);
 * 
 *         phi.submit_value(-coeff*scalar_product(phi_proj.get_value(q), n_plus), q);
 *       }
 *       phi.integrate_scatter(true, false, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Put together all the previous steps for pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   vmult_rhs_pressure(Vec& dst, const std::vector<Vec>& src) const {
 *     for(unsigned int d = 0; d < src.size(); ++d)
 *       src[d].update_ghost_values();
 * 
 *     this->data->loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_pressure,
 *                      &NavierStokesProjectionOperator::assemble_rhs_face_term_pressure,
 *                      &NavierStokesProjectionOperator::assemble_rhs_boundary_term_pressure,
 *                      this, dst, src, true,
 *                      MatrixFree<dim, Number>::DataAccessOnFaces::unspecified,
 *                      MatrixFree<dim, Number>::DataAccessOnFaces::unspecified);
 *   }
 * 
 * 
 * @endcode
 * 
 * Now we need to build the 'matrices', i.e. the bilinear forms. We start by
 * assembling the cell term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_cell_term_velocity(const MatrixFree<dim, Number>&               data,
 *                               Vec&                                         dst,
 *                               const Vec&                                   src,
 *                               const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. Moreover 'phi' in
 *       this case serves for a bilinear form and so it will not used only to submit but also to read the src ---*/
 *       FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0),
 *                                                                    phi_old_extr(data, 0);
 * 
 *       /*--- We loop over all cells in the range ---*/
 *       for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *         phi.reinit(cell);
 *         phi.gather_evaluate(src, true, true);
 *         phi_old_extr.reinit(cell);
 *         phi_old_extr.gather_evaluate(u_extr, true, false);
 * 
 *         /*--- Now we loop over all quadrature points ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& u_int                = phi.get_value(q);
 *           const auto& grad_u_int           = phi.get_gradient(q);
 *           const auto& u_n_gamma_ov_2       = phi_old_extr.get_value(q);
 *           const auto& tensor_product_u_int = outer_product(u_int, u_n_gamma_ov_2);
 * 
 *           phi.submit_value(1.0/(gamma*dt)*u_int, q);
 *           phi.submit_gradient(-a22*tensor_product_u_int + a22/Re*grad_u_int, q);
 *         }
 *         phi.integrate_scatter(true, true, dst);
 *       }
 *     }
 *     else {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0),
 *                                                                    phi_int_extr(data, 0);
 * 
 *       /*--- We loop over all cells in the range ---*/
 *       for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *         phi.reinit(cell);
 *         phi.gather_evaluate(src, true, true);
 *         phi_int_extr.reinit(cell);
 *         phi_int_extr.gather_evaluate(u_extr, true, false);
 * 
 *         /*--- Now we loop over all quadrature points ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& u_curr                   = phi.get_value(q);
 *           const auto& grad_u_curr              = phi.get_gradient(q);
 *           const auto& u_n1_int                 = phi_int_extr.get_value(q);
 *           const auto& tensor_product_u_curr    = outer_product(u_curr, u_n1_int);
 * 
 *           phi.submit_value(1.0/((1.0 - gamma)*dt)*u_curr, q);
 *           phi.submit_gradient(-a33*tensor_product_u_curr + a33/Re*grad_u_curr, q);
 *         }
 *         phi.integrate_scatter(true, true, dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles face term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_face_term_velocity(const MatrixFree<dim, Number>&               data,
 *                               Vec&                                         dst,
 *                               const Vec&                                   src,
 *                               const std::pair<unsigned int, unsigned int>& face_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi_p(data, true, 0),
 *                                                                        phi_m(data, false, 0),
 *                                                                        phi_old_extr_p(data, true, 0),
 *                                                                        phi_old_extr_m(data, false, 0);
 * 
 *       /*--- We loop over all faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_p.reinit(face);
 *         phi_p.gather_evaluate(src, true, true);
 *         phi_m.reinit(face);
 *         phi_m.gather_evaluate(src, true, true);
 *         phi_old_extr_p.reinit(face);
 *         phi_old_extr_p.gather_evaluate(u_extr, true, false);
 *         phi_old_extr_m.reinit(face);
 *         phi_old_extr_m.gather_evaluate(u_extr, true, false);
 * 
 *         const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
 *                                         std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
 * 
 *         /*--- Now we loop over all quadrature points ---*/
 *         for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *           const auto& n_plus                   = phi_p.get_normal_vector(q);
 * 
 *           const auto& avg_grad_u_int           = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
 *           const auto& jump_u_int               = phi_p.get_value(q) - phi_m.get_value(q);
 *           const auto& avg_tensor_product_u_int = 0.5*(outer_product(phi_p.get_value(q), phi_old_extr_p.get_value(q)) +
 *                                                       outer_product(phi_m.get_value(q), phi_old_extr_m.get_value(q)));
 *           const auto  lambda                   = std::max(std::abs(scalar_product(phi_old_extr_p.get_value(q), n_plus)),
 *                                                           std::abs(scalar_product(phi_old_extr_m.get_value(q), n_plus)));
 * 
 *           phi_p.submit_value(a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) +
 *                              a22*avg_tensor_product_u_int*n_plus + 0.5*a22*lambda*jump_u_int, q);
 *           phi_m.submit_value(-a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) -
 *                               a22*avg_tensor_product_u_int*n_plus - 0.5*a22*lambda*jump_u_int, q);
 *           phi_p.submit_normal_derivative(-theta_v*a22/Re*0.5*jump_u_int, q);
 *           phi_m.submit_normal_derivative(-theta_v*a22/Re*0.5*jump_u_int, q);
 *         }
 *         phi_p.integrate_scatter(true, true, dst);
 *         phi_m.integrate_scatter(true, true, dst);
 *       }
 *     }
 *     else {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi_p(data, true, 0),
 *                                                                        phi_m(data, false, 0),
 *                                                                        phi_extr_p(data, true, 0),
 *                                                                        phi_extr_m(data, false, 0);
 * 
 *       /*--- We loop over all faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_p.reinit(face);
 *         phi_p.gather_evaluate(src, true, true);
 *         phi_m.reinit(face);
 *         phi_m.gather_evaluate(src, true, true);
 *         phi_extr_p.reinit(face);
 *         phi_extr_p.gather_evaluate(u_extr, true, false);
 *         phi_extr_m.reinit(face);
 *         phi_extr_m.gather_evaluate(u_extr, true, false);
 * 
 *         const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
 *                                         std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
 * 
 *         /*--- Now we loop over all quadrature points ---*/
 *         for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *           const auto& n_plus               = phi_p.get_normal_vector(q);
 * 
 *           const auto& avg_grad_u           = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
 *           const auto& jump_u               = phi_p.get_value(q) - phi_m.get_value(q);
 *           const auto& avg_tensor_product_u = 0.5*(outer_product(phi_p.get_value(q), phi_extr_p.get_value(q)) +
 *                                                   outer_product(phi_m.get_value(q), phi_extr_m.get_value(q)));
 *           const auto  lambda               = std::max(std::abs(scalar_product(phi_extr_p.get_value(q), n_plus)),
 *                                                       std::abs(scalar_product(phi_extr_m.get_value(q), n_plus)));
 * 
 *           phi_p.submit_value(a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) +
 *                              a33*avg_tensor_product_u*n_plus + 0.5*a33*lambda*jump_u, q);
 *           phi_m.submit_value(-a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) -
 *                               a33*avg_tensor_product_u*n_plus - 0.5*a33*lambda*jump_u, q);
 *           phi_p.submit_normal_derivative(-theta_v*a33/Re*0.5*jump_u, q);
 *           phi_m.submit_normal_derivative(-theta_v*a33/Re*0.5*jump_u, q);
 *         }
 *         phi_p.integrate_scatter(true, true, dst);
 *         phi_m.integrate_scatter(true, true, dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles boundary term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_boundary_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                   Vec&                                         dst,
 *                                   const Vec&                                   src,
 *                                   const std::pair<unsigned int, unsigned int>& face_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, true, 0),
 *                                                                        phi_old_extr(data, true, 0);
 * 
 *       /*--- We loop over all faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi.reinit(face);
 *         phi.gather_evaluate(src, true, true);
 *         phi_old_extr.reinit(face);
 *         phi_old_extr.gather_evaluate(u_extr, true, false);
 * 
 *         const auto boundary_id = data.get_boundary_id(face);
 *         const auto coef_jump   = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
 * 
 *         /*--- The application of the mirror principle is not so trivial because we have a Dirichlet condition
 *               on a single component for the outflow; so we distinguish the two cases ---*/
 *         if(boundary_id != 1) {
 *           const double coef_trasp = 0.0;
 * 
 *           /*--- Now we loop over all quadrature points ---*/
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& n_plus               = phi.get_normal_vector(q);
 *             const auto& grad_u_int           = phi.get_gradient(q);
 *             const auto& u_int                = phi.get_value(q);
 *             const auto& tensor_product_u_int = outer_product(phi.get_value(q), phi_old_extr.get_value(q));
 *             const auto& lambda               = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
 * 
 *             phi.submit_value(a22/Re*(-grad_u_int*n_plus + 2.0*coef_jump*u_int) +
 *                              a22*coef_trasp*tensor_product_u_int*n_plus + a22*lambda*u_int, q);
 *             phi.submit_normal_derivative(-theta_v*a22/Re*u_int, q);
 *           }
 *           phi.integrate_scatter(true, true, dst);
 *         }
 *         else {
 *           /*--- Now we loop over all quadrature points ---*/
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& n_plus               = phi.get_normal_vector(q);
 *             const auto& grad_u_int           = phi.get_gradient(q);
 *             const auto& u_int                = phi.get_value(q);
 *             const auto& lambda               = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
 * 
 *             const auto& point_vectorized     = phi.quadrature_point(q);
 *             auto u_int_m                     = u_int;
 *             auto grad_u_int_m                = grad_u_int;
 *             for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
 *               Point<dim> point;
 *               for(unsigned int d = 0; d < dim; ++d)
 *                 point[d] = point_vectorized[d][v];
 * 
 *               u_int_m[1][v] = -u_int_m[1][v];
 * 
 *               grad_u_int_m[0][0][v] = -grad_u_int_m[0][0][v];
 *               grad_u_int_m[0][1][v] = -grad_u_int_m[0][1][v];
 *             }
 * 
 *             phi.submit_value(a22/Re*(-(0.5*(grad_u_int + grad_u_int_m))*n_plus + coef_jump*(u_int - u_int_m)) +
 *                              a22*outer_product(0.5*(u_int + u_int_m), phi_old_extr.get_value(q))*n_plus +
 *                              a22*0.5*lambda*(u_int - u_int_m), q);
 *             phi.submit_normal_derivative(-theta_v*a22/Re*(u_int - u_int_m), q);
 *           }
 *           phi.integrate_scatter(true, true, dst);
 *         }
 *       }
 *     }
 *     else {
 *       /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, true, 0),
 *                                                                        phi_extr(data, true, 0);
 * 
 *       /*--- We loop over all faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi.reinit(face);
 *         phi.gather_evaluate(src, true, true);
 *         phi_extr.reinit(face);
 *         phi_extr.gather_evaluate(u_extr, true, false);
 * 
 *         const auto boundary_id = data.get_boundary_id(face);
 *         const auto coef_jump   = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
 * 
 *         if(boundary_id != 1) {
 *           const double coef_trasp = 0.0;
 * 
 *           /*--- Now we loop over all quadrature points ---*/
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& n_plus           = phi.get_normal_vector(q);
 *             const auto& grad_u           = phi.get_gradient(q);
 *             const auto& u                = phi.get_value(q);
 *             const auto& tensor_product_u = outer_product(phi.get_value(q), phi_extr.get_value(q));
 *             const auto& lambda           = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
 * 
 *             phi.submit_value(a33/Re*(-grad_u*n_plus + 2.0*coef_jump*u) +
 *                              a33*coef_trasp*tensor_product_u*n_plus + a33*lambda*u, q);
 *             phi.submit_normal_derivative(-theta_v*a33/Re*u, q);
 *           }
 *           phi.integrate_scatter(true, true, dst);
 *         }
 *         else {
 *           /*--- Now we loop over all quadrature points ---*/
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& n_plus           = phi.get_normal_vector(q);
 *             const auto& grad_u           = phi.get_gradient(q);
 *             const auto& u                = phi.get_value(q);
 *             const auto& lambda           = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
 * 
 *             const auto& point_vectorized = phi.quadrature_point(q);
 *             auto u_m                     = u;
 *             auto grad_u_m                = grad_u;
 *             for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
 *               Point<dim> point;
 *               for(unsigned int d = 0; d < dim; ++d)
 *                 point[d] = point_vectorized[d][v];
 * 
 *               u_m[1][v] = -u_m[1][v];
 * 
 *               grad_u_m[0][0][v] = -grad_u_m[0][0][v];
 *               grad_u_m[0][1][v] = -grad_u_m[0][1][v];
 *             }
 * 
 *             phi.submit_value(a33/Re*(-(0.5*(grad_u + grad_u_m))*n_plus + coef_jump*(u - u_m)) +
 *                              a33*outer_product(0.5*(u + u_m), phi_extr.get_value(q))*n_plus + a33*0.5*lambda*(u - u_m), q);
 *             phi.submit_normal_derivative(-theta_v*a33/Re*(u - u_m), q);
 *           }
 *           phi.integrate_scatter(true, true, dst);
 *         }
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Next, we focus on 'matrices' to compute the pressure. We first assemble cell term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_cell_term_pressure(const MatrixFree<dim, Number>&               data,
 *                               Vec&                                         dst,
 *                               const Vec&                                   src,
 *                               const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *     FEEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number> phi(data, 1, 1);
 * 
 *     const double coeff = (TR_BDF2_stage == 1) ? 1.0e6*gamma*dt*gamma*dt : 1.0e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
 * 
 *     /*--- Loop over all cells in the range ---*/
 *     for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *       phi.reinit(cell);
 *       phi.gather_evaluate(src, true, true);
 * 
 *       /*--- Now we loop over all quadrature points ---*/
 *       for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *         phi.submit_gradient(phi.get_gradient(q), q);
 *         phi.submit_value(1.0/coeff*phi.get_value(q), q);
 *       }
 * 
 *       phi.integrate_scatter(true, true, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles face term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_face_term_pressure(const MatrixFree<dim, Number>&               data,
 *                               Vec&                                         dst,
 *                               const Vec&                                   src,
 *                               const std::pair<unsigned int, unsigned int>& face_range) const {
 *     /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *     FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number> phi_p(data, true, 1, 1),
 *                                                                    phi_m(data, false, 1, 1);
 * 
 *     /*--- Loop over all faces in the range ---*/
 *     for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *       phi_p.reinit(face);
 *       phi_p.gather_evaluate(src, true, true);
 *       phi_m.reinit(face);
 *       phi_m.gather_evaluate(src, true, true);
 * 
 *       const auto coef_jump = C_p*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
 *                                       std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
 * 
 *       /*--- Loop over quadrature points ---*/
 *       for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *         const auto& n_plus        = phi_p.get_normal_vector(q);
 * 
 *         const auto& avg_grad_pres = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
 *         const auto& jump_pres     = phi_p.get_value(q) - phi_m.get_value(q);
 * 
 *         phi_p.submit_value(-scalar_product(avg_grad_pres, n_plus) + coef_jump*jump_pres, q);
 *         phi_m.submit_value(scalar_product(avg_grad_pres, n_plus) - coef_jump*jump_pres, q);
 *         phi_p.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
 *         phi_m.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
 *       }
 *       phi_p.integrate_scatter(true, true, dst);
 *       phi_m.integrate_scatter(true, true, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles boundary term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_boundary_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                   Vec&                                         dst,
 *                                   const Vec&                                   src,
 *                                   const std::pair<unsigned int, unsigned int>& face_range) const {
 *     FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number> phi(data, true, 1, 1);
 * 
 *     for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *       phi.reinit(face);
 *       phi.gather_evaluate(src, true, true);
 * 
 *       const auto coef_jump = C_p*std::abs((phi.get_normal_vector(0)*phi.inverse_jacobian(0))[dim - 1]);
 * 
 *       const auto boundary_id = data.get_boundary_id(face);
 * 
 *       if(boundary_id == 1) {
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           const auto& n_plus    = phi.get_normal_vector(q);
 * 
 *           const auto& grad_pres = phi.get_gradient(q);
 *           const auto& pres      = phi.get_value(q);
 * 
 *           phi.submit_value(-scalar_product(grad_pres, n_plus) + coef_jump*pres , q);
 *           phi.submit_normal_derivative(-theta_p*pres, q);
 *         }
 *         phi.integrate_scatter(true, true, dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Before coding the 'apply_add' function, which is the one that will perform the loop, we focus on
 * the linear system that arises to project the gradient of the pressure into the velocity space.
 * The following function assembles rhs cell term for the projection of gradient of pressure. Since no
 * integration by parts is performed, only a cell term contribution is present.
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_rhs_cell_term_projection_grad_p(const MatrixFree<dim, Number>&               data,
 *                                            Vec&                                         dst,
 *                                            const Vec&                                   src,
 *                                            const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
 *     FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0);
 *     FEEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number>   phi_pres(data, 1);
 * 
 *     /*--- Loop over all cells in the range ---*/
 *     for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *       phi_pres.reinit(cell);
 *       phi_pres.gather_evaluate(src, false, true);
 *       phi.reinit(cell);
 * 
 *       /*--- Loop over quadrature points ---*/
 *       for(unsigned int q = 0; q < phi.n_q_points; ++q)
 *         phi.submit_value(phi_pres.get_gradient(q), q);
 * 
 *       phi.integrate_scatter(true, false, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Put together all the previous steps for porjection of pressure gradient. Here we loop only over cells
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   vmult_grad_p_projection(Vec& dst, const Vec& src) const {
 *     this->data->cell_loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_projection_grad_p,
 *                           this, dst, src, true);
 *   }
 * 
 * 
 * @endcode
 * 
 * Assemble now cell term for the projection of gradient of pressure. This is nothing but a mass matrix
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_cell_term_projection_grad_p(const MatrixFree<dim, Number>&               data,
 *                                        Vec&                                         dst,
 *                                        const Vec&                                   src,
 *                                        const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0);
 * 
 *     /*--- Loop over all cells in the range ---*/
 *     for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *       phi.reinit(cell);
 *       phi.gather_evaluate(src, true, false);
 * 
 *       /*--- Loop over quadrature points ---*/
 *       for(unsigned int q = 0; q < phi.n_q_points; ++q)
 *         phi.submit_value(phi.get_value(q), q);
 * 
 *       phi.integrate_scatter(true, false, dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Put together all previous steps. This is the overriden function that effectively performs the
 * matrix-vector multiplication.
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   apply_add(Vec& dst, const Vec& src) const {
 *     if(NS_stage == 1) {
 *       this->data->loop(&NavierStokesProjectionOperator::assemble_cell_term_velocity,
 *                        &NavierStokesProjectionOperator::assemble_face_term_velocity,
 *                        &NavierStokesProjectionOperator::assemble_boundary_term_velocity,
 *                        this, dst, src, false,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified);
 *     }
 *     else if(NS_stage == 2) {
 *       this->data->loop(&NavierStokesProjectionOperator::assemble_cell_term_pressure,
 *                        &NavierStokesProjectionOperator::assemble_face_term_pressure,
 *                        &NavierStokesProjectionOperator::assemble_boundary_term_pressure,
 *                        this, dst, src, false,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified);
 *     }
 *     else if(NS_stage == 3) {
 *       this->data->cell_loop(&NavierStokesProjectionOperator::assemble_cell_term_projection_grad_p,
 *                             this, dst, src, false); /*--- Since we have only a cell term contribution, we use cell_loop ---*/
 *     }
 *     else
 *       Assert(false, ExcNotImplemented());
 *   }
 * 
 * 
 * @endcode
 * 
 * Finally, we focus on computing the diagonal for preconditioners and we start by assembling
 * the diagonal cell term for the velocity. Since we do not have access to the entries of the matrix,
 * in order to compute the element i, we test the matrix against a vector which is equal to 1 in position i and 0 elsewhere.
 * This is why 'src' will result as unused.
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_diagonal_cell_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                        Vec&                                         dst,
 *                                        const unsigned int&                          ,
 *                                        const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0),
 *                                                                    phi_old_extr(data, 0);
 * 
 *       AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal(phi.dofs_per_component);
 *       /*--- Build a vector of ones to be tested (here we will see the velocity as a whole vector, since
 *                                                  dof_handler_velocity is vectorial and so the dof values are vectors). ---*/
 *       Tensor<1, dim, VectorizedArray<Number>> tmp;
 *       for(unsigned int d = 0; d < dim; ++d)
 *         tmp[d] = make_vectorized_array<Number>(1.0);
 * 
 *       /*--- Loop over cells in the range ---*/
 *       for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *         phi_old_extr.reinit(cell);
 *         phi_old_extr.gather_evaluate(u_extr, true, false);
 *         phi.reinit(cell);
 * 
 *         /*--- Loop over dofs ---*/
 *         for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *           for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *             phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j); /*--- Set all dofs to zero ---*/
 *           phi.submit_dof_value(tmp, i); /*--- Set dof i equal to one ---*/
 *           phi.evaluate(true, true);
 * 
 *           /*--- Loop over quadrature points ---*/
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& u_int                = phi.get_value(q);
 *             const auto& grad_u_int           = phi.get_gradient(q);
 *             const auto& u_n_gamma_ov_2       = phi_old_extr.get_value(q);
 *             const auto& tensor_product_u_int = outer_product(u_int, u_n_gamma_ov_2);
 * 
 *             phi.submit_value(1.0/(gamma*dt)*u_int, q);
 *             phi.submit_gradient(-a22*tensor_product_u_int + a22/Re*grad_u_int, q);
 *           }
 *           phi.integrate(true, true);
 *           diagonal[i] = phi.get_dof_value(i);
 *         }
 *         for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *           phi.submit_dof_value(diagonal[i], i);
 *         phi.distribute_local_to_global(dst);
 *       }
 *     }
 *     else {
 *       FEEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, 0),
 *                                                                    phi_int_extr(data, 0);
 * 
 *       AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal(phi.dofs_per_component);
 *       Tensor<1, dim, VectorizedArray<Number>> tmp;
 *       for(unsigned int d = 0; d < dim; ++d)
 *         tmp[d] = make_vectorized_array<Number>(1.0);
 * 
 *       /*--- Loop over cells in the range ---*/
 *       for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *         phi_int_extr.reinit(cell);
 *         phi_int_extr.gather_evaluate(u_extr, true, false);
 *         phi.reinit(cell);
 * 
 *         /*--- Loop over dofs ---*/
 *         for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *           for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *             phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *           phi.submit_dof_value(tmp, i);
 *           phi.evaluate(true, true);
 * 
 *           /*--- Loop over quadrature points ---*/
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& u_curr                   = phi.get_value(q);
 *             const auto& grad_u_curr              = phi.get_gradient(q);
 *             const auto& u_n1_int                 = phi_int_extr.get_value(q);
 *             const auto& tensor_product_u_curr    = outer_product(u_curr, u_n1_int);
 * 
 *             phi.submit_value(1.0/((1.0 - gamma)*dt)*u_curr, q);
 *             phi.submit_gradient(-a33*tensor_product_u_curr + a33/Re*grad_u_curr, q);
 *           }
 *           phi.integrate(true, true);
 *           diagonal[i] = phi.get_dof_value(i);
 *         }
 *         for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *           phi.submit_dof_value(diagonal[i], i);
 *         phi.distribute_local_to_global(dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles diagonal face term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_diagonal_face_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                        Vec&                                         dst,
 *                                        const unsigned int&                          ,
 *                                        const std::pair<unsigned int, unsigned int>& face_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi_p(data, true, 0),
 *                                                                        phi_m(data, false, 0),
 *                                                                        phi_old_extr_p(data, true, 0),
 *                                                                        phi_old_extr_m(data, false, 0);
 * 
 *       AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component); /*--- We just assert for safety that dimension match,
 *                                                                                 in the sense that we have selected the proper
 *                                                                                 space ---*/
 *       AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal_p(phi_p.dofs_per_component),
 *                                                              diagonal_m(phi_m.dofs_per_component);
 *       Tensor<1, dim, VectorizedArray<Number>> tmp;
 *       for(unsigned int d = 0; d < dim; ++d)
 *         tmp[d] = make_vectorized_array<Number>(1.0); /*--- We build the usal vector of ones that we will use as dof value ---*/
 * 
 *       /*--- Now we loop over faces ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_old_extr_p.reinit(face);
 *         phi_old_extr_p.gather_evaluate(u_extr, true, false);
 *         phi_old_extr_m.reinit(face);
 *         phi_old_extr_m.gather_evaluate(u_extr, true, false);
 *         phi_p.reinit(face);
 *         phi_m.reinit(face);
 * 
 *         const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
 *                                         std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
 * 
 *         /*--- Loop over dofs. We will set all equal to zero apart from the current one ---*/
 *         for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
 *           for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
 *             phi_p.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *             phi_m.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *           }
 *           phi_p.submit_dof_value(tmp, i);
 *           phi_p.evaluate(true, true);
 *           phi_m.submit_dof_value(tmp, i);
 *           phi_m.evaluate(true, true);
 * 
 *           /*--- Loop over quadrature points to compute the integral ---*/
 *           for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *             const auto& n_plus                   = phi_p.get_normal_vector(q);
 *             const auto& avg_grad_u_int           = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
 *             const auto& jump_u_int               = phi_p.get_value(q) - phi_m.get_value(q);
 *             const auto& avg_tensor_product_u_int = 0.5*(outer_product(phi_p.get_value(q), phi_old_extr_p.get_value(q)) +
 *                                                         outer_product(phi_m.get_value(q), phi_old_extr_m.get_value(q)));
 *             const auto  lambda                   = std::max(std::abs(scalar_product(phi_old_extr_p.get_value(q), n_plus)),
 *                                                             std::abs(scalar_product(phi_old_extr_m.get_value(q), n_plus)));
 * 
 *             phi_p.submit_value(a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) +
 *                                a22*avg_tensor_product_u_int*n_plus + 0.5*a22*lambda*jump_u_int , q);
 *             phi_m.submit_value(-a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) -
 *                                a22*avg_tensor_product_u_int*n_plus - 0.5*a22*lambda*jump_u_int, q);
 *             phi_p.submit_normal_derivative(-theta_v*0.5*a22/Re*jump_u_int, q);
 *             phi_m.submit_normal_derivative(-theta_v*0.5*a22/Re*jump_u_int, q);
 *           }
 *           phi_p.integrate(true, true);
 *           diagonal_p[i] = phi_p.get_dof_value(i);
 *           phi_m.integrate(true, true);
 *           diagonal_m[i] = phi_m.get_dof_value(i);
 *         }
 *         for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
 *           phi_p.submit_dof_value(diagonal_p[i], i);
 *           phi_m.submit_dof_value(diagonal_m[i], i);
 *         }
 *         phi_p.distribute_local_to_global(dst);
 *         phi_m.distribute_local_to_global(dst);
 *       }
 *     }
 *     else {
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi_p(data, true, 0),
 *                                                                        phi_m(data, false, 0),
 *                                                                        phi_extr_p(data, true, 0),
 *                                                                        phi_extr_m(data, false, 0);
 * 
 *       AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component);
 *       AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal_p(phi_p.dofs_per_component),
 *                                                              diagonal_m(phi_m.dofs_per_component);
 *       Tensor<1, dim, VectorizedArray<Number>> tmp;
 *       for(unsigned int d = 0; d < dim; ++d)
 *         tmp[d] = make_vectorized_array<Number>(1.0);
 * 
 *       /*--- Now we loop over faces ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_extr_p.reinit(face);
 *         phi_extr_p.gather_evaluate(u_extr, true, false);
 *         phi_extr_m.reinit(face);
 *         phi_extr_m.gather_evaluate(u_extr, true, false);
 *         phi_p.reinit(face);
 *         phi_m.reinit(face);
 * 
 *         const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
 *                                         std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
 * 
 *         /*--- Loop over dofs. We will set all equal to zero apart from the current one ---*/
 *         for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
 *           for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
 *             phi_p.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *             phi_m.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *           }
 *           phi_p.submit_dof_value(tmp, i);
 *           phi_p.evaluate(true, true);
 *           phi_m.submit_dof_value(tmp, i);
 *           phi_m.evaluate(true, true);
 * 
 *           /*--- Loop over quadrature points to compute the integral ---*/
 *           for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *             const auto& n_plus               = phi_p.get_normal_vector(q);
 *             const auto& avg_grad_u           = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
 *             const auto& jump_u               = phi_p.get_value(q) - phi_m.get_value(q);
 *             const auto& avg_tensor_product_u = 0.5*(outer_product(phi_p.get_value(q), phi_extr_p.get_value(q)) +
 *                                                     outer_product(phi_m.get_value(q), phi_extr_m.get_value(q)));
 *             const auto  lambda               = std::max(std::abs(scalar_product(phi_extr_p.get_value(q), n_plus)),
 *                                                         std::abs(scalar_product(phi_extr_m.get_value(q), n_plus)));
 * 
 *             phi_p.submit_value(a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) +
 *                                a33*avg_tensor_product_u*n_plus + 0.5*a33*lambda*jump_u, q);
 *             phi_m.submit_value(-a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) -
 *                                a33*avg_tensor_product_u*n_plus - 0.5*a33*lambda*jump_u, q);
 *             phi_p.submit_normal_derivative(-theta_v*0.5*a33/Re*jump_u, q);
 *             phi_m.submit_normal_derivative(-theta_v*0.5*a33/Re*jump_u, q);
 *           }
 *           phi_p.integrate(true, true);
 *           diagonal_p[i] = phi_p.get_dof_value(i);
 *           phi_m.integrate(true, true);
 *           diagonal_m[i] = phi_m.get_dof_value(i);
 *         }
 *         for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
 *           phi_p.submit_dof_value(diagonal_p[i], i);
 *           phi_m.submit_dof_value(diagonal_m[i], i);
 *         }
 *         phi_p.distribute_local_to_global(dst);
 *         phi_m.distribute_local_to_global(dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles boundary term for the velocity
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_diagonal_boundary_term_velocity(const MatrixFree<dim, Number>&               data,
 *                                            Vec&                                         dst,
 *                                            const unsigned int&                          ,
 *                                            const std::pair<unsigned int, unsigned int>& face_range) const {
 *     if(TR_BDF2_stage == 1) {
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, true, 0),
 *                                                                        phi_old_extr(data, true, 0);
 * 
 *       AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal(phi.dofs_per_component);
 *       Tensor<1, dim, VectorizedArray<Number>> tmp;
 *       for(unsigned int d = 0; d < dim; ++d)
 *         tmp[d] = make_vectorized_array<Number>(1.0);
 * 
 *       /*--- Loop over all faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_old_extr.reinit(face);
 *         phi_old_extr.gather_evaluate(u_extr, true, false);
 *         phi.reinit(face);
 * 
 *         const auto boundary_id = data.get_boundary_id(face);
 *         const auto coef_jump   = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
 * 
 *         if(boundary_id != 1) {
 *           const double coef_trasp = 0.0;
 * 
 *           /*--- Loop over all dofs ---*/
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *             for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *               phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *             phi.submit_dof_value(tmp, i);
 *             phi.evaluate(true, true);
 * 
 *             /*--- Loop over quadrature points to compute the integral ---*/
 *             for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *               const auto& n_plus               = phi.get_normal_vector(q);
 *               const auto& grad_u_int           = phi.get_gradient(q);
 *               const auto& u_int                = phi.get_value(q);
 *               const auto& tensor_product_u_int = outer_product(phi.get_value(q), phi_old_extr.get_value(q));
 *               const auto& lambda               = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
 * 
 *               phi.submit_value(a22/Re*(-grad_u_int*n_plus + 2.0*coef_jump*u_int) +
 *                                a22*coef_trasp*tensor_product_u_int*n_plus + a22*lambda*u_int, q);
 *               phi.submit_normal_derivative(-theta_v*a22/Re*u_int, q);
 *             }
 *             phi.integrate(true, true);
 *             diagonal[i] = phi.get_dof_value(i);
 *           }
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *             phi.submit_dof_value(diagonal[i], i);
 *           phi.distribute_local_to_global(dst);
 *         }
 *         else {
 *           /*--- Loop over all dofs ---*/
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *             for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *               phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *             phi.submit_dof_value(tmp, i);
 *             phi.evaluate(true, true);
 * 
 *             /*--- Loop over quadrature points to compute the integral ---*/
 *             for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *               const auto& n_plus               = phi.get_normal_vector(q);
 *               const auto& grad_u_int           = phi.get_gradient(q);
 *               const auto& u_int                = phi.get_value(q);
 *               const auto& lambda               = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
 * 
 *               const auto& point_vectorized     = phi.quadrature_point(q);
 *               auto u_int_m                     = u_int;
 *               auto grad_u_int_m                = grad_u_int;
 *               for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
 *                 Point<dim> point;
 *                 for(unsigned int d = 0; d < dim; ++d)
 *                   point[d] = point_vectorized[d][v];
 * 
 *                 u_int_m[1][v] = -u_int_m[1][v];
 * 
 *                 grad_u_int_m[0][0][v] = -grad_u_int_m[0][0][v];
 *                 grad_u_int_m[0][1][v] = -grad_u_int_m[0][1][v];
 *               }
 * 
 *               phi.submit_value(a22/Re*(-(0.5*(grad_u_int + grad_u_int_m))*n_plus + coef_jump*(u_int - u_int_m)) +
 *                                a22*outer_product(0.5*(u_int + u_int_m), phi_old_extr.get_value(q))*n_plus +
 *                                a22*0.5*lambda*(u_int - u_int_m), q);
 *               phi.submit_normal_derivative(-theta_v*a22/Re*(u_int - u_int_m), q);
 *             }
 *             phi.integrate(true, true);
 *             diagonal[i] = phi.get_dof_value(i);
 *           }
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *             phi.submit_dof_value(diagonal[i], i);
 *           phi.distribute_local_to_global(dst);
 *         }
 *       }
 *     }
 *     else {
 *       FEFaceEvaluation<dim, fe_degree_v, n_q_points_1d_v, dim, Number> phi(data, true, 0),
 *                                                                        phi_extr(data, true, 0);
 * 
 *       AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal(phi.dofs_per_component);
 *       Tensor<1, dim, VectorizedArray<Number>> tmp;
 *       for(unsigned int d = 0; d < dim; ++d)
 *         tmp[d] = make_vectorized_array<Number>(1.0);
 * 
 *       /*--- Loop over all faces in the range ---*/
 *       for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *         phi_extr.reinit(face);
 *         phi_extr.gather_evaluate(u_extr, true, false);
 *         phi.reinit(face);
 * 
 *         const auto boundary_id = data.get_boundary_id(face);
 *         const auto coef_jump   = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
 * 
 *         if(boundary_id != 1) {
 *           const double coef_trasp = 0.0;
 * 
 *           /*--- Loop over all dofs ---*/
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *             for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *               phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *             phi.submit_dof_value(tmp, i);
 *             phi.evaluate(true, true);
 * 
 *             /*--- Loop over quadrature points to compute the integral ---*/
 *             for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *               const auto& n_plus           = phi.get_normal_vector(q);
 *               const auto& grad_u           = phi.get_gradient(q);
 *               const auto& u                = phi.get_value(q);
 *               const auto& tensor_product_u = outer_product(phi.get_value(q), phi_extr.get_value(q));
 *               const auto& lambda           = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
 * 
 *               phi.submit_value(a33/Re*(-grad_u*n_plus + 2.0*coef_jump*u) +
 *                                a33*coef_trasp*tensor_product_u*n_plus + a33*lambda*u, q);
 *               phi.submit_normal_derivative(-theta_v*a33/Re*u, q);
 *             }
 *             phi.integrate(true, true);
 *             diagonal[i] = phi.get_dof_value(i);
 *           }
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *             phi.submit_dof_value(diagonal[i], i);
 *           phi.distribute_local_to_global(dst);
 *         }
 *         else {
 *           /*--- Loop over all dofs ---*/
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *             for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *               phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
 *             phi.submit_dof_value(tmp, i);
 *             phi.evaluate(true, true);
 * 
 *             /*--- Loop over quadrature points to compute the integral ---*/
 *             for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *               const auto& n_plus           = phi.get_normal_vector(q);
 *               const auto& grad_u           = phi.get_gradient(q);
 *               const auto& u                = phi.get_value(q);
 *               const auto& lambda           = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
 * 
 *               const auto& point_vectorized = phi.quadrature_point(q);
 *               auto u_m                     = u;
 *               auto grad_u_m                = grad_u;
 *               for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
 *                 Point<dim> point;
 *                 for(unsigned int d = 0; d < dim; ++d)
 *                   point[d] = point_vectorized[d][v];
 * 
 *                 u_m[1][v] = -u_m[1][v];
 * 
 *                 grad_u_m[0][0][v] = -grad_u_m[0][0][v];
 *                 grad_u_m[0][1][v] = -grad_u_m[0][1][v];
 *               }
 * 
 *               phi.submit_value(a33/Re*(-(0.5*(grad_u + grad_u_m))*n_plus + coef_jump*(u - u_m)) +
 *                                a33*outer_product(0.5*(u + u_m), phi_extr.get_value(q))*n_plus +
 *                                a33*0.5*lambda*(u - u_m), q);
 *               phi.submit_normal_derivative(-theta_v*a33/Re*(u - u_m), q);
 *             }
 *             phi.integrate(true, true);
 *             diagonal[i] = phi.get_dof_value(i);
 *           }
 *           for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *             phi.submit_dof_value(diagonal[i], i);
 *           phi.distribute_local_to_global(dst);
 *         }
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Now we consider the pressure related bilinear forms. We first assemble diagonal cell term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_diagonal_cell_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                        Vec&                                         dst,
 *                                        const unsigned int&                          ,
 *                                        const std::pair<unsigned int, unsigned int>& cell_range) const {
 *     FEEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number> phi(data, 1, 1);
 * 
 *     AlignedVector<VectorizedArray<Number>> diagonal(phi.dofs_per_component); /*--- Here we are using dofs_per_component but
 *                                                                                    it coincides with dofs_per_cell since it is
 *                                                                                    scalar finite element space ---*/
 * 
 *     const double coeff = (TR_BDF2_stage == 1) ? 1e6*gamma*dt*gamma*dt : 1e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
 * 
 *     /*--- Loop over all cells in the range ---*/
 *     for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
 *       phi.reinit(cell);
 * 
 *       /*--- Loop over all dofs ---*/
 *       for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *         for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *           phi.submit_dof_value(VectorizedArray<Number>(), j); /*--- We set all dofs to zero ---*/
 *         phi.submit_dof_value(make_vectorized_array<Number>(1.0), i); /*--- Now we set the current one to 1; since it is scalar,
 *                                                                            we can directly use 'make_vectorized_array' without
 *                                                                            relying on 'Tensor' ---*/
 *         phi.evaluate(true, true);
 * 
 *         /*--- Loop over quadrature points ---*/
 *         for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *           phi.submit_value(1.0/coeff*phi.get_value(q), q);
 *           phi.submit_gradient(phi.get_gradient(q), q);
 *         }
 *         phi.integrate(true, true);
 *         diagonal[i] = phi.get_dof_value(i);
 *       }
 *       for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *         phi.submit_dof_value(diagonal[i], i);
 * 
 *       phi.distribute_local_to_global(dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function assembles diagonal face term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_diagonal_face_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                        Vec&                                         dst,
 *                                        const unsigned int&                          ,
 *                                        const std::pair<unsigned int, unsigned int>& face_range) const {
 *     FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number> phi_p(data, true, 1, 1),
 *                                                                    phi_m(data, false, 1, 1);
 * 
 *     AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component);
 *     AlignedVector<VectorizedArray<Number>> diagonal_p(phi_p.dofs_per_component),
 *                                            diagonal_m(phi_m.dofs_per_component); /*--- Again, we just assert for safety that dimension
 *                                                                                        match, in the sense that we have selected
 *                                                                                        the proper space ---*/
 * 
 *     /*--- Loop over all faces ---*/
 *     for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *       phi_p.reinit(face);
 *       phi_m.reinit(face);
 * 
 *       const auto coef_jump = C_p*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
 *                                       std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
 * 
 *       /*--- Loop over all dofs ---*/
 *       for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
 *         for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
 *           phi_p.submit_dof_value(VectorizedArray<Number>(), j);
 *           phi_m.submit_dof_value(VectorizedArray<Number>(), j);
 *         }
 *         phi_p.submit_dof_value(make_vectorized_array<Number>(1.0), i);
 *         phi_m.submit_dof_value(make_vectorized_array<Number>(1.0), i);
 *         phi_p.evaluate(true, true);
 *         phi_m.evaluate(true, true);
 * 
 *         /*--- Loop over all quadrature points to compute the integral ---*/
 *         for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
 *           const auto& n_plus        = phi_p.get_normal_vector(q);
 * 
 *           const auto& avg_grad_pres = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
 *           const auto& jump_pres     = phi_p.get_value(q) - phi_m.get_value(q);
 * 
 *           phi_p.submit_value(-scalar_product(avg_grad_pres, n_plus) + coef_jump*jump_pres, q);
 *           phi_m.submit_value(scalar_product(avg_grad_pres, n_plus) - coef_jump*jump_pres, q);
 *           phi_p.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
 *           phi_m.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
 *         }
 *         phi_p.integrate(true, true);
 *         diagonal_p[i] = phi_p.get_dof_value(i);
 *         phi_m.integrate(true, true);
 *         diagonal_m[i] = phi_m.get_dof_value(i);
 *       }
 *       for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
 *         phi_p.submit_dof_value(diagonal_p[i], i);
 *         phi_m.submit_dof_value(diagonal_m[i], i);
 *       }
 *       phi_p.distribute_local_to_global(dst);
 *       phi_m.distribute_local_to_global(dst);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Eventually, we assemble diagonal boundary term for the pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   assemble_diagonal_boundary_term_pressure(const MatrixFree<dim, Number>&               data,
 *                                            Vec&                                         dst,
 *                                            const unsigned int&                          ,
 *                                            const std::pair<unsigned int, unsigned int>& face_range) const {
 *     FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_p, 1, Number> phi(data, true, 1, 1);
 * 
 *     AlignedVector<VectorizedArray<Number>> diagonal(phi.dofs_per_component);
 * 
 *     for(unsigned int face = face_range.first; face < face_range.second; ++face) {
 *       phi.reinit(face);
 * 
 *       const auto coef_jump = C_p*std::abs((phi.get_normal_vector(0)*phi.inverse_jacobian(0))[dim - 1]);
 * 
 *       const auto boundary_id = data.get_boundary_id(face);
 * 
 *       if(boundary_id == 1) {
 *         for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
 *           for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
 *             phi.submit_dof_value(VectorizedArray<Number>(), j);
 *           phi.submit_dof_value(make_vectorized_array<Number>(1.0), i);
 *           phi.evaluate(true, true);
 * 
 *           for(unsigned int q = 0; q < phi.n_q_points; ++q) {
 *             const auto& n_plus    = phi.get_normal_vector(q);
 * 
 *             const auto& grad_pres = phi.get_gradient(q);
 *             const auto& pres      = phi.get_value(q);
 * 
 *             phi.submit_value(-scalar_product(grad_pres, n_plus) + 2.0*coef_jump*pres , q);
 *             phi.submit_normal_derivative(-theta_p*pres, q);
 *           }
 *           phi.integrate(true, true);
 *           diagonal[i] = phi.get_dof_value(i);
 *         }
 *         for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
 *           phi.submit_dof_value(diagonal[i], i);
 *         phi.distribute_local_to_global(dst);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * Put together all previous steps. We create a dummy auxliary vector that serves for the src input argument in
 * the previous functions that as we have seen before is unused. Then everything is done by the 'loop' function
 * and it is saved in the field 'inverse_diagonal_entries' already present in the base class. Anyway since there is
 * only one field, we need to resize properly depending on whether we are considering the velocity or the pressure.
 *   
 
 * 
 * 
 * @code
 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec, typename Number>
 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec, Number>::
 *   compute_diagonal() {
 *     Assert(NS_stage == 1 || NS_stage == 2, ExcInternalError());
 *     if(NS_stage == 1) {
 *       this->inverse_diagonal_entries.reset(new DiagonalMatrix<Vec>());
 *       auto& inverse_diagonal = this->inverse_diagonal_entries->get_vector();
 *       this->data->initialize_dof_vector(inverse_diagonal, 0);
 *       const unsigned int dummy = 0;
 * 
 *       this->data->loop(&NavierStokesProjectionOperator::assemble_diagonal_cell_term_velocity,
 *                        &NavierStokesProjectionOperator::assemble_diagonal_face_term_velocity,
 *                        &NavierStokesProjectionOperator::assemble_diagonal_boundary_term_velocity,
 *                        this, inverse_diagonal, dummy, false,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified);
 * 
 *       for(unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i) {
 *         Assert(inverse_diagonal.local_element(i) != 0.0,
 *                ExcMessage("No diagonal entry in a definite operator should be zero"));
 *         inverse_diagonal.local_element(i) = 1.0/inverse_diagonal.local_element(i);
 *       }
 *     }
 *     else if(NS_stage == 2) {
 *       this->inverse_diagonal_entries.reset(new DiagonalMatrix<Vec>());
 *       auto& inverse_diagonal = this->inverse_diagonal_entries->get_vector();
 *       this->data->initialize_dof_vector(inverse_diagonal, 1);
 *       const unsigned int dummy = 0;
 * 
 *       this->data->loop(&NavierStokesProjectionOperator::assemble_diagonal_cell_term_pressure,
 *                        &NavierStokesProjectionOperator::assemble_diagonal_face_term_pressure,
 *                        &NavierStokesProjectionOperator::assemble_diagonal_boundary_term_pressure,
 *                        this, inverse_diagonal, dummy, false,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified,
 *                        MatrixFree<dim, Number>::DataAccessOnFaces::unspecified);
 * 
 *       for(unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i) {
 *         Assert(inverse_diagonal.local_element(i) != 0.0,
 *                ExcMessage("No diagonal entry in a definite operator should be zero"));
 *         inverse_diagonal.local_element(i) = 1.0/inverse_diagonal.local_element(i);
 *       }
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{The <code>NavierStokesProjection</code> class}
 * 
 
 * 
 * Now we are ready for the main class of the program. It implements the calls to the various steps
 * of the projection method for Navier-Stokes equations.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   class NavierStokesProjection {
 *   public:
 *     NavierStokesProjection(RunTimeParameters::Data_Storage& data);
 * 
 *     void run(const bool verbose = false, const unsigned int output_interval = 10);
 * 
 *   protected:
 *     const double t_0;
 *     const double T;
 *     const double gamma;         //--- TR-BDF2 parameter
 *     unsigned int TR_BDF2_stage; //--- Flag to check at which current stage of TR-BDF2 are
 *     const double Re;
 *     double       dt;
 * 
 *     EquationData::Velocity<dim> vel_init;
 *     EquationData::Pressure<dim> pres_init; /*--- Instance of 'Velocity' and 'Pressure' classes to initialize. ---*/
 * 
 *     parallel::distributed::Triangulation<dim> triangulation;
 * 
 *     /*--- Finite Element spaces ---*/
 *     FESystem<dim> fe_velocity;
 *     FESystem<dim> fe_pressure;
 * 
 *     /*--- Handler for dofs ---*/
 *     DoFHandler<dim> dof_handler_velocity;
 *     DoFHandler<dim> dof_handler_pressure;
 * 
 *     /*--- Quadrature formulas for velocity and pressure, respectively ---*/
 *     QGauss<dim> quadrature_pressure;
 *     QGauss<dim> quadrature_velocity;
 * 
 *     /*--- Now we define all the vectors for the solution. We start from the pressure
 *           with p^n, p^(n+gamma) and a vector for rhs ---*/
 *     LinearAlgebra::distributed::Vector<double> pres_n;
 *     LinearAlgebra::distributed::Vector<double> pres_int;
 *     LinearAlgebra::distributed::Vector<double> rhs_p;
 * 
 *     /*--- Next, we move to the velocity, with u^n, u^(n-1), u^(n+gamma/2),
 *           u^(n+gamma) and other two auxiliary vectors as well as the rhs ---*/
 *     LinearAlgebra::distributed::Vector<double> u_n;
 *     LinearAlgebra::distributed::Vector<double> u_n_minus_1;
 *     LinearAlgebra::distributed::Vector<double> u_extr;
 *     LinearAlgebra::distributed::Vector<double> u_n_gamma;
 *     LinearAlgebra::distributed::Vector<double> u_star;
 *     LinearAlgebra::distributed::Vector<double> u_tmp;
 *     LinearAlgebra::distributed::Vector<double> rhs_u;
 *     LinearAlgebra::distributed::Vector<double> grad_pres_int;
 * 
 *     Vector<double> Linfty_error_per_cell_vel;
 * 
 *     DeclException2(ExcInvalidTimeStep,
 *                    double,
 *                    double,
 *                    << " The time step " << arg1 << " is out of range."
 *                    << std::endl
 *                    << " The permitted range is (0," << arg2 << "]");
 * 
 *     void create_triangulation(const unsigned int n_refines);
 * 
 *     void setup_dofs();
 * 
 *     void initialize();
 * 
 *     void interpolate_velocity();
 * 
 *     void diffusion_step();
 * 
 *     void projection_step();
 * 
 *     void project_grad(const unsigned int flag);
 * 
 *     double get_maximal_velocity();
 * 
 *     double get_maximal_difference();
 * 
 *     void output_results(const unsigned int step);
 * 
 *     void refine_mesh();
 * 
 *     void interpolate_max_res(const unsigned int level);
 * 
 *     void save_max_res();
 * 
 *   private:
 *     void compute_lift_and_drag();
 * 
 *     /*--- Technical member to handle the various steps ---*/
 *     std::shared_ptr<MatrixFree<dim, double>> matrix_free_storage;
 * 
 *     /*--- Now we need an instance of the class implemented before with the weak form ---*/
 *     NavierStokesProjectionOperator<dim, EquationData::degree_p, EquationData::degree_p + 1,
 *                                    EquationData::degree_p + 1, EquationData::degree_p + 2,
 *                                    LinearAlgebra::distributed::Vector<double>, double> navier_stokes_matrix;
 * 
 *     /*--- This is an instance for geometric multigrid preconditioner ---*/
 *     MGLevelObject<NavierStokesProjectionOperator<dim, EquationData::degree_p, EquationData::degree_p + 1,
 *                                                  EquationData::degree_p + 1, EquationData::degree_p + 2,
 *                                                  LinearAlgebra::distributed::Vector<float>, float>> mg_matrices;
 * 
 *     /*--- Here we define two 'AffineConstraints' instance, one for each finite element space.
 *           This is just a technical issue, due to MatrixFree requirements. In general
 *           this class is used to impose boundary conditions (or any kind of constraints), but in this case, since
 *           we are using a weak imposition of bcs, everything is already in the weak forms and so these instances
 *           will be default constructed ---*/
 *     AffineConstraints<double> constraints_velocity,
 *                               constraints_pressure;
 * 
 *     /*--- Now a bunch of variables handled by 'ParamHandler' introduced at the beginning of the code ---*/
 *     unsigned int max_its;
 *     double       eps;
 * 
 *     unsigned int max_loc_refinements;
 *     unsigned int min_loc_refinements;
 *     unsigned int refinement_iterations;
 * 
 *     std::string saving_dir;
 * 
 *     /*--- Finally, some output related streams ---*/
 *     ConditionalOStream pcout;
 * 
 *     std::ofstream      time_out;
 *     ConditionalOStream ptime_out;
 *     TimerOutput        time_table;
 * 
 *     std::ofstream output_n_dofs_velocity;
 *     std::ofstream output_n_dofs_pressure;
 * 
 *     std::ofstream output_lift;
 *     std::ofstream output_drag;
 *   };
 * 
 * 
 * @endcode
 * 
 * In the constructor, we just read all the data from the
 * <code>Data_Storage</code> object that is passed as an argument, verify that
 * the data we read are reasonable and, finally, create the triangulation and
 * load the initial data.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   NavierStokesProjection<dim>::NavierStokesProjection(RunTimeParameters::Data_Storage& data):
 *     t_0(data.initial_time),
 *     T(data.final_time),
 *     gamma(2.0 - std::sqrt(2.0)),  //--- Save also in the NavierStokes class the TR-BDF2 parameter value
 *     TR_BDF2_stage(1),             //--- Initialize the flag for the TR_BDF2 stage
 *     Re(data.Reynolds),
 *     dt(data.dt),
 *     vel_init(data.initial_time),
 *     pres_init(data.initial_time),
 *     triangulation(MPI_COMM_WORLD, parallel::distributed::Triangulation<dim>::limit_level_difference_at_vertices,
 *                   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
 *     fe_velocity(FE_DGQ<dim>(EquationData::degree_p + 1), dim),
 *     fe_pressure(FE_DGQ<dim>(EquationData::degree_p), 1),
 *     dof_handler_velocity(triangulation),
 *     dof_handler_pressure(triangulation),
 *     quadrature_pressure(EquationData::degree_p + 1),
 *     quadrature_velocity(EquationData::degree_p + 2),
 *     navier_stokes_matrix(data),
 *     max_its(data.max_iterations),
 *     eps(data.eps),
 *     max_loc_refinements(data.max_loc_refinements),
 *     min_loc_refinements(data.min_loc_refinements),
 *     refinement_iterations(data.refinement_iterations),
 *     saving_dir(data.dir),
 *     pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
 *     time_out("./" + data.dir + "/time_analysis_" +
 *              Utilities::int_to_string(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)) + "proc.dat"),
 *     ptime_out(time_out, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
 *     time_table(ptime_out, TimerOutput::summary, TimerOutput::cpu_and_wall_times),
 *     output_n_dofs_velocity("./" + data.dir + "/n_dofs_velocity.dat", std::ofstream::out),
 *     output_n_dofs_pressure("./" + data.dir + "/n_dofs_pressure.dat", std::ofstream::out),
 *     output_lift("./" + data.dir + "/lift.dat", std::ofstream::out),
 *     output_drag("./" + data.dir + "/drag.dat", std::ofstream::out) {
 *       if(EquationData::degree_p < 1) {
 *         pcout
 *         << " WARNING: The chosen pair of finite element spaces is not stable."
 *         << std::endl
 *         << " The obtained results will be nonsense" << std::endl;
 *       }
 * 
 *       AssertThrow(!((dt <= 0.0) || (dt > 0.5*T)), ExcInvalidTimeStep(dt, 0.5*T));
 * 
 *       matrix_free_storage = std::make_shared<MatrixFree<dim, double>>();
 * 
 *       create_triangulation(data.n_refines);
 *       setup_dofs();
 *       initialize();
 *   }
 * 
 * 
 * @endcode
 * 
 * The method that creates the triangulation and refines it the needed number
 * of times.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::create_triangulation(const unsigned int n_refines) {
 *     TimerOutput::Scope t(time_table, "Create triangulation");
 * 
 *     GridGenerator::plate_with_a_hole(triangulation, 0.5, 1.0, 1.0, 1.1, 1.0, 19.0, Point<2>(2.0, 2.0), 0, 1, 1.0, 2, true);
 *     /*--- We strongly advice to check the documentation to verify the meaning of all input parameters. ---*/
 * 
 *     pcout << "Number of refines = " << n_refines << std::endl;
 *     triangulation.refine_global(n_refines);
 *   }
 * 
 * 
 * @endcode
 * 
 * After creating the triangulation, it creates the mesh dependent
 * data, i.e. it distributes degrees of freedom, and
 * initializes the vectors that we will use.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::setup_dofs() {
 *     pcout << "Number of active cells: " << triangulation.n_global_active_cells() << std::endl;
 *     pcout << "Number of levels: "       << triangulation.n_global_levels()       << std::endl;
 * 
 *     /*--- Distribute dofs and prepare for multigrid ---*/
 *     dof_handler_velocity.distribute_dofs(fe_velocity);
 *     dof_handler_pressure.distribute_dofs(fe_pressure);
 * 
 *     pcout << "dim (X_h) = " << dof_handler_velocity.n_dofs()
 *           << std::endl
 *           << "dim (M_h) = " << dof_handler_pressure.n_dofs()
 *           << std::endl
 *           << "Re        = " << Re << std::endl
 *           << std::endl;
 * 
 *     if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
 *       output_n_dofs_velocity << dof_handler_velocity.n_dofs() << std::endl;
 *       output_n_dofs_pressure << dof_handler_pressure.n_dofs() << std::endl;
 *     }
 * 
 *     typename MatrixFree<dim, double>::AdditionalData additional_data;
 *     additional_data.mapping_update_flags                = (update_gradients | update_JxW_values |
 *                                                            update_quadrature_points | update_values);
 *     additional_data.mapping_update_flags_inner_faces    = (update_gradients | update_JxW_values | update_quadrature_points |
 *                                                            update_normal_vectors | update_values);
 *     additional_data.mapping_update_flags_boundary_faces = (update_gradients | update_JxW_values | update_quadrature_points |
 *                                                            update_normal_vectors | update_values);
 *     additional_data.tasks_parallel_scheme               = MatrixFree<dim, double>::AdditionalData::none;
 * 
 *     std::vector<const DoFHandler<dim>*> dof_handlers; /*--- Vector of dof_handlers to feed the 'MatrixFree'. Here the order
 *                                                             counts and enters into the game as parameter of FEEvaluation and
 *                                                             FEFaceEvaluation in the previous class ---*/
 *     dof_handlers.push_back(&dof_handler_velocity);
 *     dof_handlers.push_back(&dof_handler_pressure);
 * 
 *     constraints_velocity.clear();
 *     constraints_velocity.close();
 *     constraints_pressure.clear();
 *     constraints_pressure.close();
 *     std::vector<const AffineConstraints<double>*> constraints;
 *     constraints.push_back(&constraints_velocity);
 *     constraints.push_back(&constraints_pressure);
 * 
 *     std::vector<QGauss<1>> quadratures; /*--- We cannot directly use 'quadrature_velocity' and 'quadrature_pressure',
 *                                               because the 'MatrixFree' structure wants a quadrature formula for 1D
 *                                               (this is way the template parameter of the previous class was called 'n_q_points_1d_p'
 *                                                and 'n_q_points_1d_v' and the reason of '1' as QGauss template parameter). ---*/
 *     quadratures.push_back(QGauss<1>(EquationData::degree_p + 2));
 *     quadratures.push_back(QGauss<1>(EquationData::degree_p + 1));
 * 
 *     /*--- Initialize the matrix-free structure and size properly the vectors. Here again the
 *           second input argument of the 'initialize_dof_vector' method depends on the order of 'dof_handlers' ---*/
 *     matrix_free_storage->reinit(MappingQ1<dim>(),dof_handlers, constraints, quadratures, additional_data);
 *     matrix_free_storage->initialize_dof_vector(u_star, 0);
 *     matrix_free_storage->initialize_dof_vector(rhs_u, 0);
 *     matrix_free_storage->initialize_dof_vector(u_n, 0);
 *     matrix_free_storage->initialize_dof_vector(u_extr, 0);
 *     matrix_free_storage->initialize_dof_vector(u_n_minus_1, 0);
 *     matrix_free_storage->initialize_dof_vector(u_n_gamma, 0);
 *     matrix_free_storage->initialize_dof_vector(u_tmp, 0);
 *     matrix_free_storage->initialize_dof_vector(grad_pres_int, 0);
 * 
 *     matrix_free_storage->initialize_dof_vector(pres_int, 1);
 *     matrix_free_storage->initialize_dof_vector(pres_n, 1);
 *     matrix_free_storage->initialize_dof_vector(rhs_p, 1);
 * 
 *     /*--- Initialize the multigrid structure. We dedicate ad hoc 'dof_handlers_mg' and 'constraints_mg' because
 *           we use float as type. Moreover we can initialize already with the index of the finite element of the pressure;
 *           anyway we need by requirement to declare also structures for the velocity for coherence (basically because
 *           the index of finite element space has to be the same, so the pressure has to be the second).---*/
 *     mg_matrices.clear_elements();
 *     dof_handler_velocity.distribute_mg_dofs();
 *     dof_handler_pressure.distribute_mg_dofs();
 * 
 *     const unsigned int nlevels = triangulation.n_global_levels();
 *     mg_matrices.resize(0, nlevels - 1);
 *     for(unsigned int level = 0; level < nlevels; ++level) {
 *       typename MatrixFree<dim, float>::AdditionalData additional_data_mg;
 *       additional_data_mg.tasks_parallel_scheme               = MatrixFree<dim, float>::AdditionalData::none;
 *       additional_data_mg.mapping_update_flags                = (update_gradients | update_JxW_values);
 *       additional_data_mg.mapping_update_flags_inner_faces    = (update_gradients | update_JxW_values);
 *       additional_data_mg.mapping_update_flags_boundary_faces = (update_gradients | update_JxW_values);
 *       additional_data_mg.mg_level = level;
 * 
 *       std::vector<const DoFHandler<dim>*> dof_handlers_mg;
 *       dof_handlers_mg.push_back(&dof_handler_velocity);
 *       dof_handlers_mg.push_back(&dof_handler_pressure);
 *       std::vector<const AffineConstraints<float>*> constraints_mg;
 *       AffineConstraints<float> constraints_velocity_mg;
 *       constraints_velocity_mg.clear();
 *       constraints_velocity_mg.close();
 *       constraints_mg.push_back(&constraints_velocity_mg);
 *       AffineConstraints<float> constraints_pressure_mg;
 *       constraints_pressure_mg.clear();
 *       constraints_pressure_mg.close();
 *       constraints_mg.push_back(&constraints_pressure_mg);
 * 
 *       std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(new MatrixFree<dim, float>());
 *       mg_mf_storage_level->reinit(MappingQ1<dim>(),dof_handlers_mg, constraints_mg, quadratures, additional_data_mg);
 *       const std::vector<unsigned int> tmp = {1};
 *       mg_matrices[level].initialize(mg_mf_storage_level, tmp, tmp);
 *       mg_matrices[level].set_dt(dt);
 *       mg_matrices[level].set_NS_stage(2);
 *     }
 * 
 *     Linfty_error_per_cell_vel.reinit(triangulation.n_active_cells());
 *   }
 * 
 * 
 * @endcode
 * 
 * This method loads the initial data. It simply uses the class <code>Pressure</code> instance for the pressure
 * and the class <code>Velocity</code> instance for the velocity.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::initialize() {
 *     TimerOutput::Scope t(time_table, "Initialize pressure and velocity");
 * 
 *     VectorTools::interpolate(dof_handler_pressure, pres_init, pres_n);
 * 
 *     VectorTools::interpolate(dof_handler_velocity, vel_init, u_n_minus_1);
 *     VectorTools::interpolate(dof_handler_velocity, vel_init, u_n);
 *   }
 * 
 * 
 * @endcode
 * 
 * This function computes the extrapolated velocity to be used in the momentum predictor
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::interpolate_velocity() {
 *     TimerOutput::Scope t(time_table, "Interpolate velocity");
 * 
 * @endcode
 * 
 * --- TR-BDF2 first step
 * 
 * @code
 *     if(TR_BDF2_stage == 1) {
 *       u_extr.equ(1.0 + gamma/(2.0*(1.0 - gamma)), u_n);
 *       u_tmp.equ(gamma/(2.0*(1.0 - gamma)), u_n_minus_1);
 *       u_extr -= u_tmp;
 *     }
 * @endcode
 * 
 * --- TR-BDF2 second step
 * 
 * @code
 *     else {
 *       u_extr.equ(1.0 + (1.0 - gamma)/gamma, u_n_gamma);
 *       u_tmp.equ((1.0 - gamma)/gamma, u_n);
 *       u_extr -= u_tmp;
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * We are finally ready to solve the diffusion step.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::diffusion_step() {
 *     TimerOutput::Scope t(time_table, "Diffusion step");
 * 
 *     /*--- We first speicify that we want to deal with velocity dof_handler (index 0, since it is the first one
 *           in the 'dof_handlers' vector) ---*/
 *     const std::vector<unsigned int> tmp = {0};
 *     navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
 * 
 *     /*--- Next, we specify at we are at stage 1, namely the diffusion step ---*/
 *     navier_stokes_matrix.set_NS_stage(1);
 * 
 *     /*--- Now, we compute the right-hand side and we set the convective velocity. The necessity of 'set_u_extr' is
 *           that this quantity is required in the bilinear forms and we can't use a vector of src like on the right-hand side,
 *           so it has to be available ---*/
 *     if(TR_BDF2_stage == 1) {
 *       navier_stokes_matrix.vmult_rhs_velocity(rhs_u, {u_n, u_extr, pres_n});
 *       navier_stokes_matrix.set_u_extr(u_extr);
 *       u_star = u_extr;
 *     }
 *     else {
 *       navier_stokes_matrix.vmult_rhs_velocity(rhs_u, {u_n, u_n_gamma, pres_int, u_extr});
 *       navier_stokes_matrix.set_u_extr(u_extr);
 *       u_star = u_extr;
 *     }
 * 
 *     /*--- Build the linear solver; in this case we specifiy the maximum number of iterations and residual ---*/
 *     SolverControl solver_control(max_its, eps*rhs_u.l2_norm());
 *     SolverGMRES<LinearAlgebra::distributed::Vector<double>> gmres(solver_control);
 * 
 *     /*--- Build a Jacobi preconditioner and solve ---*/
 *     PreconditionJacobi<NavierStokesProjectionOperator<dim,
 *                                                       EquationData::degree_p,
 *                                                       EquationData::degree_p + 1,
 *                                                       EquationData::degree_p + 1,
 *                                                       EquationData::degree_p + 2,
 *                                                       LinearAlgebra::distributed::Vector<double>,
 *                                                       double>> preconditioner;
 *     navier_stokes_matrix.compute_diagonal();
 *     preconditioner.initialize(navier_stokes_matrix);
 * 
 *     gmres.solve(navier_stokes_matrix, u_star, rhs_u, preconditioner);
 *   }
 * 
 * 
 * @endcode
 * 
 * Next, we solve the projection step.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::projection_step() {
 *     TimerOutput::Scope t(time_table, "Projection step pressure");
 * 
 *     /*--- We start in the same way of 'diffusion_step': we first reinitialize with the index of FE space,
 *           we specify that this is the second stage and we compute the right-hand side ---*/
 *     const std::vector<unsigned int> tmp = {1};
 *     navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
 * 
 *     navier_stokes_matrix.set_NS_stage(2);
 * 
 *     if(TR_BDF2_stage == 1)
 *       navier_stokes_matrix.vmult_rhs_pressure(rhs_p, {u_star, pres_n});
 *     else
 *       navier_stokes_matrix.vmult_rhs_pressure(rhs_p, {u_star, pres_int});
 * 
 *     /*--- Build the linear solver (Conjugate Gradient in this case) ---*/
 *     SolverControl solver_control(max_its, eps*rhs_p.l2_norm());
 *     SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
 * 
 *     /*--- Build the preconditioner (as in @ref step_37 "step-37") ---*/
 *     MGTransferMatrixFree<dim, float> mg_transfer;
 *     mg_transfer.build(dof_handler_pressure);
 * 
 *     using SmootherType = PreconditionChebyshev<NavierStokesProjectionOperator<dim,
 *                                                                               EquationData::degree_p,
 *                                                                               EquationData::degree_p + 1,
 *                                                                               EquationData::degree_p + 1,
 *                                                                               EquationData::degree_p + 2,
 *                                                                               LinearAlgebra::distributed::Vector<float>,
 *                                                                               float>,
 *                                                LinearAlgebra::distributed::Vector<float>>;
 *     mg::SmootherRelaxation<SmootherType, LinearAlgebra::distributed::Vector<float>> mg_smoother;
 *     MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
 *     smoother_data.resize(0, triangulation.n_global_levels() - 1);
 *     for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level) {
 *       if(level > 0) {
 *         smoother_data[level].smoothing_range     = 15.0;
 *         smoother_data[level].degree              = 3;
 *         smoother_data[level].eig_cg_n_iterations = 10;
 *       }
 *       else {
 *         smoother_data[0].smoothing_range     = 2e-2;
 *         smoother_data[0].degree              = numbers::invalid_unsigned_int;
 *         smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
 *       }
 *       mg_matrices[level].compute_diagonal();
 *       smoother_data[level].preconditioner = mg_matrices[level].get_matrix_diagonal_inverse();
 *     }
 *     mg_smoother.initialize(mg_matrices, smoother_data);
 * 
 *     PreconditionIdentity                                identity;
 *     SolverCG<LinearAlgebra::distributed::Vector<float>> cg_mg(solver_control);
 *     MGCoarseGridIterativeSolver<LinearAlgebra::distributed::Vector<float>,
 *                                 SolverCG<LinearAlgebra::distributed::Vector<float>>,
 *                                 NavierStokesProjectionOperator<dim,
 *                                                                EquationData::degree_p,
 *                                                                EquationData::degree_p + 1,
 *                                                                EquationData::degree_p + 1,
 *                                                                EquationData::degree_p + 2,
 *                                                                LinearAlgebra::distributed::Vector<float>,
 *                                                                float>,
 *                                 PreconditionIdentity> mg_coarse(cg_mg, mg_matrices[0], identity);
 * 
 *     mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(mg_matrices);
 * 
 *     Multigrid<LinearAlgebra::distributed::Vector<float>> mg(mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
 * 
 *     PreconditionMG<dim,
 *                    LinearAlgebra::distributed::Vector<float>,
 *                    MGTransferMatrixFree<dim, float>> preconditioner(dof_handler_pressure, mg, mg_transfer);
 * 
 *     /*--- Solve the linear system ---*/
 *     if(TR_BDF2_stage == 1) {
 *       pres_int = pres_n;
 *       cg.solve(navier_stokes_matrix, pres_int, rhs_p, preconditioner);
 *     }
 *     else {
 *       pres_n = pres_int;
 *       cg.solve(navier_stokes_matrix, pres_n, rhs_p, preconditioner);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * This implements the projection step for the gradient of pressure
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::project_grad(const unsigned int flag) {
 *     TimerOutput::Scope t(time_table, "Gradient of pressure projection");
 * 
 *     /*--- The input parameter flag is used just to specify where we want to save the result ---*/
 *     AssertIndexRange(flag, 3);
 *     Assert(flag > 0, ExcInternalError());
 * 
 *     /*--- We need to select the dof handler related to the velocity since the result lives there ---*/
 *     const std::vector<unsigned int> tmp = {0};
 *     navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
 * 
 *     if(flag == 1)
 *       navier_stokes_matrix.vmult_grad_p_projection(rhs_u, pres_n);
 *     else if(flag == 2)
 *       navier_stokes_matrix.vmult_grad_p_projection(rhs_u, pres_int);
 * 
 *     /*--- We conventionally decide that the this corresponds to third stage ---*/
 *     navier_stokes_matrix.set_NS_stage(3);
 * 
 *     /*--- Solve the system ---*/
 *     SolverControl solver_control(max_its, 1e-12*rhs_u.l2_norm());
 *     SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
 *     cg.solve(navier_stokes_matrix, u_tmp, rhs_u, PreconditionIdentity());
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function is used in determining the maximal velocity
 * in order to compute the Courant number.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   double NavierStokesProjection<dim>::get_maximal_velocity() {
 *     VectorTools::integrate_difference(dof_handler_velocity, u_n, ZeroFunction<dim>(dim),
 *                                       Linfty_error_per_cell_vel, quadrature_velocity, VectorTools::Linfty_norm);
 *     const double res = VectorTools::compute_global_error(triangulation, Linfty_error_per_cell_vel, VectorTools::Linfty_norm);
 * 
 *     return res;
 *   }
 * 
 * 
 * @endcode
 * 
 * The following function is used in determining the maximal nodal difference
 * in order to see if we have reched steady-state. We simply use integrate_difference testing
 * u_n - u_n_minus_1 against the zero function.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   double NavierStokesProjection<dim>::get_maximal_difference() {
 *     u_tmp = u_n;
 *     u_tmp -= u_n_minus_1;
 * 
 *     VectorTools::integrate_difference(dof_handler_velocity, u_tmp, ZeroFunction<dim>(dim),
 *                                       Linfty_error_per_cell_vel, quadrature_velocity, VectorTools::Linfty_norm);
 *     const double res = VectorTools::compute_global_error(triangulation, Linfty_error_per_cell_vel, VectorTools::Linfty_norm);
 *     pcout << "Maximum nodal difference = " << res <<std::endl;
 * 
 *     return res;
 *   }
 * 
 * 
 * @endcode
 * 
 * This method plots the current solution. The main difficulty is that we want
 * to create a single output file that contains the data for all velocity
 * components and the pressure. On the other hand, velocities and the pressure
 * live on separate DoFHandler objects, so we need to pay attention when we use
 * 'add_data_vector' to select the proper space.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::output_results(const unsigned int step) {
 *     TimerOutput::Scope t(time_table, "Output results");
 * 
 *     DataOut<dim> data_out;
 * 
 *     std::vector<std::string> velocity_names(dim, "v");
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *     component_interpretation_velocity(dim, DataComponentInterpretation::component_is_part_of_vector);
 *     u_n.update_ghost_values();
 *     data_out.add_data_vector(dof_handler_velocity, u_n, velocity_names, component_interpretation_velocity);
 *     pres_n.update_ghost_values();
 *     data_out.add_data_vector(dof_handler_pressure, pres_n, "p", {DataComponentInterpretation::component_is_scalar});
 * 
 *     std::vector<std::string> velocity_names_old(dim, "v_old");
 *     u_n_minus_1.update_ghost_values();
 *     data_out.add_data_vector(dof_handler_velocity, u_n_minus_1, velocity_names_old, component_interpretation_velocity);
 * 
 *     /*--- Here we rely on the postprocessor we have built ---*/
 *     PostprocessorVorticity<dim> postprocessor;
 *     data_out.add_data_vector(dof_handler_velocity, u_n, postprocessor);
 * 
 *     data_out.build_patches(MappingQ1<dim>(), 1, DataOut<dim>::curved_inner_cells);
 * 
 *     const std::string output = "./" + saving_dir + "/solution-" + Utilities::int_to_string(step, 5) + ".vtu";
 *     data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{<code>NavierStokesProjection::compute_lift_and_drag</code>}
 * 
 
 * 
 * This routine computes the lift and the drag forces in a non-dimensional framework
 * (so basically for the classical coefficients, it is necessary to multiply by a factor 2).
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::compute_lift_and_drag() {
 *     QGauss<dim - 1> face_quadrature_formula(EquationData::degree_p + 2);
 *     const int n_q_points = face_quadrature_formula.size();
 * 
 *     std::vector<double>                      pressure_values(n_q_points);
 *     std::vector<std::vector<Tensor<1, dim>>> velocity_gradients(n_q_points, std::vector<Tensor<1, dim>>(dim));
 * 
 *     Tensor<1, dim> normal_vector;
 *     Tensor<2, dim> fluid_stress;
 *     Tensor<2, dim> fluid_pressure;
 *     Tensor<1, dim> forces;
 * 
 *     /*--- We need to compute the integral over the cylinder boundary, so we need to use 'FEFaceValues' instances.
 *           For the velocity we need the gradients, for the pressure the values. ---*/
 *     FEFaceValues<dim> fe_face_values_velocity(fe_velocity, face_quadrature_formula,
 *                                               update_quadrature_points | update_gradients |
 *                                               update_JxW_values | update_normal_vectors);
 *     FEFaceValues<dim> fe_face_values_pressure(fe_pressure, face_quadrature_formula, update_values);
 * 
 *     double local_drag = 0.0;
 *     double local_lift = 0.0;
 * 
 *     /*--- We need to perform a unique loop because the whole stress tensor takes into account contributions of
 *           velocity and pressure obviously. However, the two dof_handlers are different, so we neede to create an ad-hoc
 *           iterator for the pressure that we update manually. It is guaranteed that the cells are visited in the same order
 *           (see the documentation) ---*/
 *     auto tmp_cell = dof_handler_pressure.begin_active();
 *     for(const auto& cell : dof_handler_velocity.active_cell_iterators()) {
 *       if(cell->is_locally_owned()) {
 *         for(unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
 *           if(cell->face(face)->at_boundary() && cell->face(face)->boundary_id() == 4) {
 *             fe_face_values_velocity.reinit(cell, face);
 *             fe_face_values_pressure.reinit(tmp_cell, face);
 * 
 *             fe_face_values_velocity.get_function_gradients(u_n, velocity_gradients); /*--- velocity gradients ---*/
 *             fe_face_values_pressure.get_function_values(pres_n, pressure_values); /*--- pressure values ---*/
 * 
 *             for(int q = 0; q < n_q_points; q++) {
 *               normal_vector = -fe_face_values_velocity.normal_vector(q);
 * 
 *               for(unsigned int d = 0; d < dim; ++ d) {
 *                 fluid_pressure[d][d] = pressure_values[q];
 *                 for(unsigned int k = 0; k < dim; ++k)
 *                   fluid_stress[d][k] = 1.0/Re*velocity_gradients[q][d][k];
 *               }
 *               fluid_stress = fluid_stress - fluid_pressure;
 * 
 *               forces = fluid_stress*normal_vector*fe_face_values_velocity.JxW(q);
 * 
 *               local_drag += forces[0];
 *               local_lift += forces[1];
 *             }
 *           }
 *         }
 *       }
 *       ++tmp_cell;
 *     }
 * 
 *     /*--- At the end, each processor has computed the contribution to the boundary cells it owns and, therefore,
 *           we need to sum up all the contributions. ---*/
 *     double lift = Utilities::MPI::sum(local_lift, MPI_COMM_WORLD);
 *     double drag = Utilities::MPI::sum(local_drag, MPI_COMM_WORLD);
 *     if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
 *       output_lift << lift << std::endl;
 *       output_drag << drag << std::endl;
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{ <code>NavierStokesProjection::refine_mesh</code>}
 * 
 
 * 
 * After finding a good initial guess on the coarse mesh, we hope to
 * decrease the error through refining the mesh. We also need to transfer the current solution to the
 * next mesh using the SolutionTransfer class.
 *   
 
 * 
 * 
 * @code
 *   template <int dim>
 *   void NavierStokesProjection<dim>::refine_mesh() {
 *     TimerOutput::Scope t(time_table, "Refine mesh");
 * 
 *     /*--- We first create a proper vector for computing estimator ---*/
 *     IndexSet locally_relevant_dofs;
 *     DoFTools::extract_locally_relevant_dofs(dof_handler_velocity, locally_relevant_dofs);
 *     LinearAlgebra::distributed::Vector<double> tmp_velocity;
 *     tmp_velocity.reinit(dof_handler_velocity.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
 *     tmp_velocity = u_n;
 *     tmp_velocity.update_ghost_values();
 * 
 *     using Iterator = typename DoFHandler<dim>::active_cell_iterator;
 *     Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 * 
 *     /*--- This is basically the indicator per cell computation (see @ref step_50 "step-50"). Since it is not so complciated
 *           we implement it through a lambda expression ---*/
 *     auto cell_worker = [&](const Iterator&   cell,
 *                            ScratchData<dim>& scratch_data,
 *                            CopyData&         copy_data) {
 *       FEValues<dim>& fe_values = scratch_data.fe_values; /*--- Here we finally use the 'FEValues' inside ScratchData ---*/
 *       fe_values.reinit(cell);
 * 
 *       /*--- Compute the gradients for all quadrature points ---*/
 *       std::vector<std::vector<Tensor<1, dim>>> gradients(fe_values.n_quadrature_points, std::vector<Tensor<1, dim>>(dim));
 *       fe_values.get_function_gradients(tmp_velocity, gradients);
 *       copy_data.cell_index = cell->active_cell_index();
 *       double vorticity_norm_square = 0.0;
 *       /*--- Loop over quadrature points and evaluate the integral multiplying the vorticty
 *             by the weights and the determinant of the Jacobian (which are included in 'JxW') ---*/
 *       for(unsigned k = 0; k < fe_values.n_quadrature_points; ++k) {
 *         const double vorticity = gradients[k][1][0] - gradients[k][0][1];
 *         vorticity_norm_square += vorticity*vorticity*fe_values.JxW(k);
 *       }
 *       copy_data.value = cell->diameter()*cell->diameter()*vorticity_norm_square;
 *     };
 * 
 *     const UpdateFlags cell_flags = update_gradients | update_quadrature_points | update_JxW_values;
 * 
 *     auto copier = [&](const CopyData &copy_data) {
 *       if(copy_data.cell_index != numbers::invalid_unsigned_int)
 *         estimated_error_per_cell[copy_data.cell_index] += copy_data.value;
 *     };
 * 
 *     /*--- Now everything is 'automagically' handled by 'mesh_loop' ---*/
 *     ScratchData scratch_data(fe_velocity, EquationData::degree_p + 2, cell_flags);
 *     CopyData copy_data;
 *     MeshWorker::mesh_loop(dof_handler_velocity.begin_active(),
 *                           dof_handler_velocity.end(),
 *                           cell_worker,
 *                           copier,
 *                           scratch_data,
 *                           copy_data,
 *                           MeshWorker::assemble_own_cells);
 * 
 *     /*--- Refine grid. In case the refinement level is above a certain value (or the coarsening level is below)
 *           we clear the flags. ---*/
 *     parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(triangulation, estimated_error_per_cell, 0.01, 0.3);
 *     for(const auto& cell: triangulation.active_cell_iterators()) {
 *       if(cell->refine_flag_set() && static_cast<unsigned int>(cell->level()) == max_loc_refinements)
 *         cell->clear_refine_flag();
 *       if(cell->coarsen_flag_set() && static_cast<unsigned int>(cell->level()) == min_loc_refinements)
 *         cell->clear_coarsen_flag();
 *     }
 *     triangulation.prepare_coarsening_and_refinement();
 * 
 *     /*--- Now we prepare the object for transfering, basically saving the old quantities using SolutionTransfer.
 *           Since the 'prepare_for_coarsening_and_refinement' method can be called only once, but we have two vectors
 *           for dof_handler_velocity, we need to put them in an auxiliary vector. ---*/
 *     std::vector<const LinearAlgebra::distributed::Vector<double>*> velocities;
 *     velocities.push_back(&u_n);
 *     velocities.push_back(&u_n_minus_1);
 *     parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
 *     solution_transfer_velocity(dof_handler_velocity);
 *     solution_transfer_velocity.prepare_for_coarsening_and_refinement(velocities);
 *     parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
 *     solution_transfer_pressure(dof_handler_pressure);
 *     solution_transfer_pressure.prepare_for_coarsening_and_refinement(pres_n);
 * 
 *     triangulation.execute_coarsening_and_refinement(); /*--- Effectively perform the remeshing ---*/
 * 
 *     /*--- First DoFHandler objects are set up within the new grid ----*/
 *     setup_dofs();
 * 
 *     /*--- Interpolate current solutions to new mesh. This is done using auxliary vectors just for safety,
 *           but the new u_n or pres_n could be used. Again, the only point is that the function 'interpolate'
 *           can be called once and so the vectors related to 'dof_handler_velocity' have to collected in an auxiliary vector. ---*/
 *     LinearAlgebra::distributed::Vector<double> transfer_velocity,
 *                                                transfer_velocity_minus_1,
 *                                                transfer_pressure;
 *     transfer_velocity.reinit(u_n);
 *     transfer_velocity.zero_out_ghost_values();
 *     transfer_velocity_minus_1.reinit(u_n_minus_1);
 *     transfer_velocity_minus_1.zero_out_ghost_values();
 *     transfer_pressure.reinit(pres_n);
 *     transfer_pressure.zero_out_ghost_values();
 * 
 *     std::vector<LinearAlgebra::distributed::Vector<double>*> transfer_velocities;
 *     transfer_velocities.push_back(&transfer_velocity);
 *     transfer_velocities.push_back(&transfer_velocity_minus_1);
 *     solution_transfer_velocity.interpolate(transfer_velocities);
 *     transfer_velocity.update_ghost_values();
 *     transfer_velocity_minus_1.update_ghost_values();
 *     solution_transfer_pressure.interpolate(transfer_pressure);
 *     transfer_pressure.update_ghost_values();
 * 
 *     u_n         = transfer_velocity;
 *     u_n_minus_1 = transfer_velocity_minus_1;
 *     pres_n      = transfer_pressure;
 *   }
 * 
 * 
 * @endcode
 * 
 * Interpolate the locally refined solution to a mesh with maximal resolution
 * and transfer velocity and pressure.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::interpolate_max_res(const unsigned int level) {
 *     parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
 *     solution_transfer_velocity(dof_handler_velocity);
 *     std::vector<const LinearAlgebra::distributed::Vector<double>*> velocities;
 *     velocities.push_back(&u_n);
 *     velocities.push_back(&u_n_minus_1);
 *     solution_transfer_velocity.prepare_for_coarsening_and_refinement(velocities);
 * 
 *     parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
 *     solution_transfer_pressure(dof_handler_pressure);
 *     solution_transfer_pressure.prepare_for_coarsening_and_refinement(pres_n);
 * 
 *     for(const auto& cell: triangulation.active_cell_iterators_on_level(level)) {
 *       if(cell->is_locally_owned())
 *         cell->set_refine_flag();
 *     }
 *     triangulation.execute_coarsening_and_refinement();
 * 
 *     setup_dofs();
 * 
 *     LinearAlgebra::distributed::Vector<double> transfer_velocity, transfer_velocity_minus_1,
 *                                                transfer_pressure;
 * 
 *     transfer_velocity.reinit(u_n);
 *     transfer_velocity.zero_out_ghost_values();
 *     transfer_velocity_minus_1.reinit(u_n_minus_1);
 *     transfer_velocity_minus_1.zero_out_ghost_values();
 * 
 *     transfer_pressure.reinit(pres_n);
 *     transfer_pressure.zero_out_ghost_values();
 * 
 *     std::vector<LinearAlgebra::distributed::Vector<double>*> transfer_velocities;
 * 
 *     transfer_velocities.push_back(&transfer_velocity);
 *     transfer_velocities.push_back(&transfer_velocity_minus_1);
 *     solution_transfer_velocity.interpolate(transfer_velocities);
 *     transfer_velocity.update_ghost_values();
 *     transfer_velocity_minus_1.update_ghost_values();
 * 
 *     solution_transfer_pressure.interpolate(transfer_pressure);
 *     transfer_pressure.update_ghost_values();
 * 
 *     u_n         = transfer_velocity;
 *     u_n_minus_1 = transfer_velocity_minus_1;
 *     pres_n      = transfer_pressure;
 *   }
 * 
 * 
 * @endcode
 * 
 * Save maximum resolution to a mesh adapted.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::save_max_res() {
 *     parallel::distributed::Triangulation<dim> triangulation_tmp(MPI_COMM_WORLD);
 *     GridGenerator::plate_with_a_hole(triangulation_tmp, 0.5, 1.0, 1.0, 1.1, 1.0, 19.0, Point<2>(2.0, 2.0), 0, 1, 1.0, 2, true);
 *     triangulation_tmp.refine_global(triangulation.n_global_levels() - 1);
 * 
 *     DoFHandler<dim> dof_handler_velocity_tmp(triangulation_tmp);
 *     DoFHandler<dim> dof_handler_pressure_tmp(triangulation_tmp);
 *     dof_handler_velocity_tmp.distribute_dofs(fe_velocity);
 *     dof_handler_pressure_tmp.distribute_dofs(fe_pressure);
 * 
 *     LinearAlgebra::distributed::Vector<double> u_n_tmp,
 *                                                pres_n_tmp;
 *     u_n_tmp.reinit(dof_handler_velocity_tmp.n_dofs());
 *     pres_n_tmp.reinit(dof_handler_pressure_tmp.n_dofs());
 * 
 *     DataOut<dim> data_out;
 *     std::vector<std::string> velocity_names(dim, "v");
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *     component_interpretation_velocity(dim, DataComponentInterpretation::component_is_part_of_vector);
 *     VectorTools::interpolate_to_different_mesh(dof_handler_velocity, u_n, dof_handler_velocity_tmp, u_n_tmp);
 *     u_n_tmp.update_ghost_values();
 *     data_out.add_data_vector(dof_handler_velocity_tmp, u_n_tmp, velocity_names, component_interpretation_velocity);
 *     VectorTools::interpolate_to_different_mesh(dof_handler_pressure, pres_n, dof_handler_pressure_tmp, pres_n_tmp);
 *     pres_n_tmp.update_ghost_values();
 *     data_out.add_data_vector(dof_handler_pressure_tmp, pres_n_tmp, "p", {DataComponentInterpretation::component_is_scalar});
 *     PostprocessorVorticity<dim> postprocessor;
 *     data_out.add_data_vector(dof_handler_velocity_tmp, u_n_tmp, postprocessor);
 * 
 *     data_out.build_patches(MappingQ1<dim>(), 1, DataOut<dim>::curved_inner_cells);
 *     const std::string output = "./" + saving_dir + "/solution_max_res_end.vtu";
 *     data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{ <code>NavierStokesProjection::run</code> }
 * 
 
 * 
 * This is the time marching function, which starting at <code>t_0</code>
 * advances in time using the projection method with time step <code>dt</code>
 * until <code>T</code>.
 *   
 
 * 
 * Its second parameter, <code>verbose</code> indicates whether the function
 * should output information what it is doing at any given moment:
 * we use the ConditionalOStream class to do that for us.
 *   
 
 * 
 * 
 * @code
 *   template<int dim>
 *   void NavierStokesProjection<dim>::run(const bool verbose, const unsigned int output_interval) {
 *     ConditionalOStream verbose_cout(std::cout, verbose && Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0);
 * 
 *     output_results(1);
 *     double time = t_0 + dt;
 *     unsigned int n = 1;
 *     while(std::abs(T - time) > 1e-10) {
 *       time += dt;
 *       n++;
 *       pcout << "Step = " << n << " Time = " << time << std::endl;
 * 
 *       /*--- First stage of TR-BDF2 and we start by setting the proper flag ---*/
 *       TR_BDF2_stage = 1;
 *       navier_stokes_matrix.set_TR_BDF2_stage(TR_BDF2_stage);
 *       for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
 *         mg_matrices[level].set_TR_BDF2_stage(TR_BDF2_stage);
 * 
 *       verbose_cout << "  Interpolating the velocity stage 1" << std::endl;
 *       interpolate_velocity();
 * 
 *       verbose_cout << "  Diffusion Step stage 1 " << std::endl;
 *       diffusion_step();
 * 
 *       verbose_cout << "  Projection Step stage 1" << std::endl;
 *       project_grad(1);
 *       u_tmp.equ(gamma*dt, u_tmp);
 *       u_star += u_tmp; /*--- In the rhs of the projection step we need u_star + gamma*dt*grad(pres_n) and we save it into u_star ---*/
 *       projection_step();
 * 
 *       verbose_cout << "  Updating the Velocity stage 1" << std::endl;
 *       u_n_gamma.equ(1.0, u_star);
 *       project_grad(2);
 *       grad_pres_int.equ(1.0, u_tmp); /*--- We save grad(pres_int), because we will need it soon ---*/
 *       u_tmp.equ(-gamma*dt, u_tmp);
 *       u_n_gamma += u_tmp; /*--- u_n_gamma = u_star - gamma*dt*grad(pres_int) ---*/
 *       u_n_minus_1 = u_n;
 * 
 *       /*--- Second stage of TR-BDF2 ---*/
 *       TR_BDF2_stage = 2;
 *       for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
 *         mg_matrices[level].set_TR_BDF2_stage(TR_BDF2_stage);
 *       navier_stokes_matrix.set_TR_BDF2_stage(TR_BDF2_stage);
 * 
 *       verbose_cout << "  Interpolating the velocity stage 2" << std::endl;
 *       interpolate_velocity();
 * 
 *       verbose_cout << "  Diffusion Step stage 2 " << std::endl;
 *       diffusion_step();
 * 
 *       verbose_cout << "  Projection Step stage 2" << std::endl;
 *       u_tmp.equ((1.0 - gamma)*dt, grad_pres_int);
 *       u_star += u_tmp;  /*--- In the rhs of the projection step we need u_star + (1 - gamma)*dt*grad(pres_int) ---*/
 *       projection_step();
 * 
 *       verbose_cout << "  Updating the Velocity stage 2" << std::endl;
 *       u_n.equ(1.0, u_star);
 *       project_grad(1);
 *       u_tmp.equ((gamma - 1.0)*dt, u_tmp);
 *       u_n += u_tmp;  /*--- u_n = u_star - (1 - gamma)*dt*grad(pres_n) ---*/
 * 
 *       const double max_vel = get_maximal_velocity();
 *       pcout<< "Maximal velocity = " << max_vel << std::endl;
 *       /*--- The Courant number is computed taking into account the polynomial degree for the velocity ---*/
 *       pcout << "CFL = " << dt*max_vel*(EquationData::degree_p + 1)*
 *                            std::sqrt(dim)/GridTools::minimal_cell_diameter(triangulation) << std::endl;
 *       compute_lift_and_drag();
 *       if(n % output_interval == 0) {
 *         verbose_cout << "Plotting Solution final" << std::endl;
 *         output_results(n);
 *       }
 *       /*--- In case dt is not a multiple of T, we reduce dt in order to end up at T ---*/
 *       if(T - time < dt && T - time > 1e-10) {
 *         dt = T - time;
 *         navier_stokes_matrix.set_dt(dt);
 *         for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
 *           mg_matrices[level].set_dt(dt);
 *       }
 *       /*--- Perform the refinement if desired ---*/
 *       if(refinement_iterations > 0 && n % refinement_iterations == 0) {
 *         verbose_cout << "Refining mesh" << std::endl;
 *         refine_mesh();
 *       }
 *     }
 *     if(n % output_interval != 0) {
 *       verbose_cout << "Plotting Solution final" << std::endl;
 *       output_results(n);
 *     }
 *     if(refinement_iterations > 0) {
 *       for(unsigned int lev = 0; lev < triangulation.n_global_levels() - 1; ++ lev)
 *         interpolate_max_res(lev);
 *       save_max_res();
 *     }
 *   }
 * 
 * } // namespace NS_TRBDF2
 * 
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{ The main function }
 * 
 
 * 
 * The main function looks very much like in all the other tutorial programs. We first initialize MPI,
 * we initialize the class 'NavierStokesProjection' with the dimension as template parameter and then
 * let the method 'run' do the job.
 * 
 
 * 
 * 
 * @code
 * int main(int argc, char *argv[]) {
 *   try {
 *     using namespace NS_TRBDF2;
 * 
 *     RunTimeParameters::Data_Storage data;
 *     data.read_data("parameter-file.prm");
 * 
 *     Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, -1);
 * 
 *     const auto& curr_rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
 *     deallog.depth_console(data.verbose && curr_rank == 0 ? 2 : 0);
 * 
 *     NavierStokesProjection<2> test(data);
 *     test.run(data.verbose, data.output_interval);
 * 
 *     if(curr_rank == 0)
 *       std::cout << "----------------------------------------------------"
 *                 << std::endl
 *                 << "Apparently everything went fine!" << std::endl
 *                 << "Don't forget to brush your teeth :-)" << std::endl
 *                 << std::endl;
 * 
 *     return 0;
 *   }
 *   catch(std::exception &exc) {
 *     std::cerr << std::endl
 *               << std::endl
 *               << "----------------------------------------------------"
 *               << std::endl;
 *     std::cerr << "Exception on processing: " << std::endl
 *               << exc.what() << std::endl
 *               << "Aborting!" << std::endl
 *               << "----------------------------------------------------"
 *               << std::endl;
 *     return 1;
 *   }
 *   catch(...) {
 *     std::cerr << std::endl
 *               << std::endl
 *               << "----------------------------------------------------"
 *               << std::endl;
 *     std::cerr << "Unknown exception!" << std::endl
 *               << "Aborting!" << std::endl
 *               << "----------------------------------------------------"
 *               << std::endl;
 *     return 1;
 *   }
 * 
 * }
 * @endcode
 
 
<a name="ann-runtime_parameters.h"></a>
<h1>Annotated version of runtime_parameters.h</h1>
 * 
 * 
 * 
 * We start by including all the necessary deal.II header files
 * 
 
 * 
 * 
 * @code
 * #include <deal.II/base/parameter_handler.h>
 * 
 * @endcode
 * 
 * 
 * <a name=""></a> 
 * @sect{Run time parameters}
 * 
 
 * 
 * Since our method has several parameters that can be fine-tuned we put them
 * into an external file, so that they can be determined at run-time.
 * 
 
 * 
 * 
 * @code
 * namespace RunTimeParameters {
 *   using namespace dealii;
 * 
 *   class Data_Storage {
 *   public:
 *     Data_Storage();
 * 
 *     void read_data(const std::string& filename);
 * 
 *     double initial_time;
 *     double final_time;
 * 
 *     double Reynolds;
 *     double dt;
 * 
 *     unsigned int n_refines;            /*--- Number of refinements ---*/
 *     unsigned int max_loc_refinements; /*--- Number of maximum local refinements allowed ---*/
 *     unsigned int min_loc_refinements; /*--- Number of minimum local refinements allowed
 *                                             once reached that level ---*/
 * 
 *     /*--- Parameters related to the linear solver ---*/
 *     unsigned int max_iterations;
 *     double       eps;
 * 
 *     bool         verbose;
 *     unsigned int output_interval;
 * 
 *     std::string dir; /*--- Auxiliary string variable for output storage ---*/
 * 
 *     unsigned int refinement_iterations; /*--- Auxiliary variable about how many steps perform remeshing ---*/
 * 
 *   protected:
 *     ParameterHandler prm;
 *   };
 * 
 * @endcode
 * 
 * In the constructor of this class we declare all the parameters in suitable (but arbitrary) subsections.
 *   
 
 * 
 * 
 * @code
 *   Data_Storage::Data_Storage(): initial_time(0.0),
 *                                 final_time(1.0),
 *                                 Reynolds(1.0),
 *                                 dt(5e-4),
 *                                 n_refines(0),
 *                                 max_loc_refinements(0),
 *                                 min_loc_refinements(0),
 *                                 max_iterations(1000),
 *                                 eps(1e-12),
 *                                 verbose(true),
 *                                 output_interval(15),
 *                                 refinement_iterations(0) {
 *     prm.enter_subsection("Physical data");
 *     {
 *       prm.declare_entry("initial_time",
 *                         "0.0",
 *                         Patterns::Double(0.0),
 *                         " The initial time of the simulation. ");
 *       prm.declare_entry("final_time",
 *                         "1.0",
 *                         Patterns::Double(0.0),
 *                         " The final time of the simulation. ");
 *       prm.declare_entry("Reynolds",
 *                         "1.0",
 *                         Patterns::Double(0.0),
 *                         " The Reynolds number. ");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.enter_subsection("Time step data");
 *     {
 *       prm.declare_entry("dt",
 *                         "5e-4",
 *                         Patterns::Double(0.0),
 *                         " The time step size. ");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.enter_subsection("Space discretization");
 *     {
 *       prm.declare_entry("n_of_refines",
 *                         "100",
 *                         Patterns::Integer(0, 1500),
 *                         " The number of cells we want on each direction of the mesh. ");
 *       prm.declare_entry("max_loc_refinements",
 *                         "4",
 *                          Patterns::Integer(0, 10),
 *                          " The number of maximum local refinements. ");
 *       prm.declare_entry("min_loc_refinements",
 *                         "2",
 *                          Patterns::Integer(0, 10),
 *                          " The number of minimum local refinements. ");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.enter_subsection("Data solve");
 *     {
 *       prm.declare_entry("max_iterations",
 *                         "1000",
 *                         Patterns::Integer(1, 30000),
 *                         " The maximal number of iterations linear solvers must make. ");
 *       prm.declare_entry("eps",
 *                         "1e-12",
 *                         Patterns::Double(0.0),
 *                         " The stopping criterion. ");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.declare_entry("refinement_iterations",
 *                       "0",
 *                       Patterns::Integer(0),
 *                       " This number indicates how often we need to "
 *                       "refine the mesh");
 * 
 *     prm.declare_entry("saving directory", "SimTest");
 * 
 *     prm.declare_entry("verbose",
 *                       "true",
 *                       Patterns::Bool(),
 *                       " This indicates whether the output of the solution "
 *                       "process should be verbose. ");
 * 
 *     prm.declare_entry("output_interval",
 *                       "1",
 *                       Patterns::Integer(1),
 *                       " This indicates between how many time steps we print "
 *                       "the solution. ");
 *   }
 * 
 * @endcode
 * 
 * We need now a routine to read all declared parameters in the constructor
 *   
 
 * 
 * 
 * @code
 *   void Data_Storage::read_data(const std::string& filename) {
 *     std::ifstream file(filename);
 *     AssertThrow(file, ExcFileNotOpen(filename));
 * 
 *     prm.parse_input(file);
 * 
 *     prm.enter_subsection("Physical data");
 *     {
 *       initial_time = prm.get_double("initial_time");
 *       final_time   = prm.get_double("final_time");
 *       Reynolds     = prm.get_double("Reynolds");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.enter_subsection("Time step data");
 *     {
 *       dt = prm.get_double("dt");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.enter_subsection("Space discretization");
 *     {
 *       n_refines           = prm.get_integer("n_of_refines");
 *       max_loc_refinements = prm.get_integer("max_loc_refinements");
 *       min_loc_refinements = prm.get_integer("min_loc_refinements");
 *     }
 *     prm.leave_subsection();
 * 
 *     prm.enter_subsection("Data solve");
 *     {
 *       max_iterations = prm.get_integer("max_iterations");
 *       eps            = prm.get_double("eps");
 *     }
 *     prm.leave_subsection();
 * 
 *     dir = prm.get("saving directory");
 * 
 *     refinement_iterations = prm.get_integer("refinement_iterations");
 * 
 *     verbose = prm.get_bool("verbose");
 * 
 *     output_interval = prm.get_integer("output_interval");
 *   }
 * 
 * } // namespace RunTimeParameters
 * @endcode
 
 
*/