Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
NavierStokes_TRBDF2_DG.h
Go to the documentation of this file.
1
165 *  
166 * @endcode
167 *
168 * We start by including all the necessary deal.II header files.
169 *
170
171 *
172 *
173 * @code
174 *   #include <deal.II/base/point.h>
175 *   #include <deal.II/base/function.h>
176 *  
177 * @endcode
178 *
179 *
180 * <a name="equation_data.h-"></a>
181 * @sect{Equation data}
182 *
183
184 *
185 * In the next namespace, we declare and implement suitable functions that may be used for the initial and boundary conditions
186 *
187
188 *
189 *
190 * @code
191 *   namespace EquationData {
192 *   using namespace dealii;
193 *  
194 *   static const unsigned int degree_p = 1; /*--- Polynomial degree for the pressure. The one for the velocity
195 *   will be equal to degree_p + 1, but its value can be easily changed
196 *   in the template parameter list of the classes with weak form ---*/
197 *  
198 * @endcode
199 *
200 * We declare class that describes the boundary conditions and initial one for velocity:
201 *
202
203 *
204 *
205 * @code
206 *   template<int dim>
207 *   class Velocity: public Function<dim> {
208 *   public:
209 *   Velocity(const double initial_time = 0.0);
210 *  
211 *   virtual double value(const Point<dim>& p,
212 *   const unsigned int component = 0) const override;
213 *  
214 *   virtual void vector_value(const Point<dim>& p,
215 *   Vector<double>& values) const override;
216 *   };
217 *  
218 *  
219 *   template<int dim>
220 *   Velocity<dim>::Velocity(const double initial_time): Function<dim>(dim, initial_time) {}
221 *  
222 *  
223 *   template<int dim>
224 *   double Velocity<dim>::value(const Point<dim>& p, const unsigned int component) const {
225 *   AssertIndexRange(component, 3);
226 *   if(component == 0) {
227 *   const double Um = 1.5;
228 *   const double H = 4.1;
229 *  
230 *   return 4.0*Um*p(1)*(H - p(1))/(H*H);
231 *   }
232 *   else
233 *   return 0.0;
234 *   }
235 *  
236 *  
237 *   template<int dim>
238 *   void Velocity<dim>::vector_value(const Point<dim>& p, Vector<double>& values) const {
239 *   Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
240 *  
241 *   for(unsigned int i = 0; i < dim; ++i)
242 *   values[i] = value(p, i);
243 *   }
244 *  
245 *  
246 * @endcode
247 *
248 * We do the same for the pressure
249 *
250
251 *
252 *
253 * @code
254 *   template<int dim>
255 *   class Pressure: public Function<dim> {
256 *   public:
257 *   Pressure(const double initial_time = 0.0);
258 *  
259 *   virtual double value(const Point<dim>& p,
260 *   const unsigned int component = 0) const override;
261 *   };
262 *  
263 *  
264 *   template<int dim>
265 *   Pressure<dim>::Pressure(const double initial_time): Function<dim>(1, initial_time) {}
266 *  
267 *  
268 *   template<int dim>
269 *   double Pressure<dim>::value(const Point<dim>& p, const unsigned int component) const {
270 *   (void)component;
271 *   AssertIndexRange(component, 1);
272 *  
273 *   return 22.0 - p(0);
274 *   }
275 *  
276 *   } // namespace EquationData
277 * @endcode
278
279
280<a name="ann-navier_stokes_TRBDF2_DG.cc"></a>
281<h1>Annotated version of navier_stokes_TRBDF2_DG.cc</h1>
282 *
283 *
284 *
285 *
286 * @code
287 *   /* -----------------------------------------------------------------------------
288 *   *
289 *   * SPDX-License-Identifier: LGPL-2.1-or-later
290 *   * Copyright (C) 2022 by Giuseppe Orlando
291 *   *
292 *   * This file is part of the deal.II code gallery.
293 *   *
294 *   * -----------------------------------------------------------------------------
295 *   *
296 *   * Author: Giuseppe Orlando, 2022.
297 *   */
298 *  
299 * @endcode
300 *
301 * We start by including all the necessary deal.II header files and some C++
302 * related ones.
303 *
304
305 *
306 *
307 * @code
308 *   #include <deal.II/base/quadrature_lib.h>
309 *   #include <deal.II/base/multithread_info.h>
310 *   #include <deal.II/base/thread_management.h>
311 *   #include <deal.II/base/work_stream.h>
312 *   #include <deal.II/base/parallel.h>
313 *   #include <deal.II/base/utilities.h>
314 *   #include <deal.II/base/conditional_ostream.h>
315 *  
316 *   #include <deal.II/lac/vector.h>
317 *   #include <deal.II/lac/solver_cg.h>
318 *   #include <deal.II/lac/precondition.h>
319 *   #include <deal.II/lac/solver_gmres.h>
320 *   #include <deal.II/lac/affine_constraints.h>
321 *  
322 *   #include <deal.II/grid/tria.h>
323 *   #include <deal.II/grid/grid_generator.h>
324 *   #include <deal.II/grid/grid_tools.h>
325 *   #include <deal.II/grid/grid_refinement.h>
326 *   #include <deal.II/grid/tria_accessor.h>
327 *   #include <deal.II/grid/tria_iterator.h>
328 *   #include <deal.II/distributed/grid_refinement.h>
329 *  
330 *   #include <deal.II/dofs/dof_handler.h>
331 *   #include <deal.II/dofs/dof_accessor.h>
332 *   #include <deal.II/dofs/dof_tools.h>
333 *  
334 *   #include <deal.II/fe/fe_q.h>
335 *   #include <deal.II/fe/fe_dgq.h>
336 *   #include <deal.II/fe/fe_values.h>
337 *   #include <deal.II/fe/fe_tools.h>
338 *   #include <deal.II/fe/fe_system.h>
339 *  
340 *   #include <deal.II/numerics/matrix_tools.h>
341 *   #include <deal.II/numerics/vector_tools.h>
342 *   #include <deal.II/numerics/data_out.h>
343 *  
344 *   #include <fstream>
345 *   #include <cmath>
346 *   #include <iostream>
347 *  
348 *   #include <deal.II/matrix_free/matrix_free.h>
349 *   #include <deal.II/matrix_free/operators.h>
350 *   #include <deal.II/matrix_free/fe_evaluation.h>
351 *   #include <deal.II/fe/component_mask.h>
352 *  
353 *   #include <deal.II/base/timer.h>
354 *   #include <deal.II/distributed/solution_transfer.h>
355 *   #include <deal.II/numerics/error_estimator.h>
356 *  
357 *   #include <deal.II/multigrid/multigrid.h>
358 *   #include <deal.II/multigrid/mg_transfer_matrix_free.h>
359 *   #include <deal.II/multigrid/mg_tools.h>
360 *   #include <deal.II/multigrid/mg_coarse.h>
361 *   #include <deal.II/multigrid/mg_smoother.h>
362 *   #include <deal.II/multigrid/mg_matrix.h>
363 *  
364 *   #include <deal.II/meshworker/mesh_loop.h>
365 *  
366 *   #include "runtime_parameters.h"
367 *   #include "equation_data.h"
368 *  
369 *   namespace MatrixFreeTools {
370 *   using namespace dealii;
371 *  
372 *   template<int dim, typename Number, typename VectorizedArrayType>
374 *   LinearAlgebra::distributed::Vector<Number>& diagonal_global,
375 *   const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType>&,
377 *   const unsigned int&,
378 *   const std::pair<unsigned int, unsigned int>&)>& cell_operation,
379 *   const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType>&,
381 *   const unsigned int&,
382 *   const std::pair<unsigned int, unsigned int>&)>& face_operation,
383 *   const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType>&,
385 *   const unsigned int&,
386 *   const std::pair<unsigned int, unsigned int>&)>& boundary_operation,
387 *   const unsigned int dof_no = 0) {
388 * @endcode
389 *
390 * initialize vector
391 *
392 * @code
393 *   matrix_free.initialize_dof_vector(diagonal_global, dof_no);
394 *  
395 *   const unsigned int dummy = 0;
396 *  
397 *   matrix_free.loop(cell_operation, face_operation, boundary_operation,
398 *   diagonal_global, dummy, false,
401 *   }
402 *   }
403 *  
404 * @endcode
405 *
406 * We include the code in a suitable namespace:
407 *
408
409 *
410 *
411 * @code
412 *   namespace NS_TRBDF2 {
413 *   using namespace dealii;
414 *  
415 * @endcode
416 *
417 * The following class is an auxiliary one for post-processing of the vorticity
418 *
419
420 *
421 *
422 * @code
423 *   template<int dim>
424 *   class PostprocessorVorticity: public DataPostprocessor<dim> {
425 *   public:
426 *   virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector<dim>& inputs,
427 *   std::vector<Vector<double>>& computed_quantities) const override;
428 *  
429 *   virtual std::vector<std::string> get_names() const override;
430 *  
431 *   virtual std::vector<DataComponentInterpretation::DataComponentInterpretation>
432 *   get_data_component_interpretation() const override;
433 *  
434 *   virtual UpdateFlags get_needed_update_flags() const override;
435 *   };
436 *  
437 * @endcode
438 *
439 * This function evaluates the vorticty in both 2D and 3D cases
440 *
441
442 *
443 *
444 * @code
445 *   template <int dim>
446 *   void PostprocessorVorticity<dim>::evaluate_vector_field(const DataPostprocessorInputs::Vector<dim>& inputs,
447 *   std::vector<Vector<double>>& computed_quantities) const {
448 *   const unsigned int n_quadrature_points = inputs.solution_values.size();
449 *  
450 *   /*--- Check the correctness of all data structures ---*/
451 *   Assert(inputs.solution_gradients.size() == n_quadrature_points, ExcInternalError());
452 *   Assert(computed_quantities.size() == n_quadrature_points, ExcInternalError());
453 *  
454 *   Assert(inputs.solution_values[0].size() == dim, ExcInternalError());
455 *  
456 *   if(dim == 2) {
457 *   Assert(computed_quantities[0].size() == 1, ExcInternalError());
458 *   }
459 *   else {
460 *   Assert(computed_quantities[0].size() == dim, ExcInternalError());
461 *   }
462 *  
463 *   /*--- Compute the vorticty ---*/
464 *   if(dim == 2) {
465 *   for(unsigned int q = 0; q < n_quadrature_points; ++q)
466 *   computed_quantities[q](0) = inputs.solution_gradients[q][1][0] - inputs.solution_gradients[q][0][1];
467 *   }
468 *   else {
469 *   for(unsigned int q = 0; q < n_quadrature_points; ++q) {
470 *   computed_quantities[q](0) = inputs.solution_gradients[q][2][1] - inputs.solution_gradients[q][1][2];
471 *   computed_quantities[q](1) = inputs.solution_gradients[q][0][2] - inputs.solution_gradients[q][2][0];
472 *   computed_quantities[q](2) = inputs.solution_gradients[q][1][0] - inputs.solution_gradients[q][0][1];
473 *   }
474 *   }
475 *   }
476 *  
477 * @endcode
478 *
479 * This auxiliary function is required by the base class DataProcessor and simply
480 * sets the name for the output file
481 *
482
483 *
484 *
485 * @code
486 *   template<int dim>
487 *   std::vector<std::string> PostprocessorVorticity<dim>::get_names() const {
488 *   std::vector<std::string> names;
489 *   names.emplace_back("vorticity");
490 *   if(dim == 3) {
491 *   names.emplace_back("vorticity");
492 *   names.emplace_back("vorticity");
493 *   }
494 *  
495 *   return names;
496 *   }
497 *  
498 * @endcode
499 *
500 * This auxiliary function is required by the base class DataProcessor and simply
501 * specifies if the vorticity is a scalar (2D) or a vector (3D)
502 *
503
504 *
505 *
506 * @code
507 *   template<int dim>
508 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
509 *   PostprocessorVorticity<dim>::get_data_component_interpretation() const {
510 *   std::vector<DataComponentInterpretation::DataComponentInterpretation> interpretation;
511 *   if(dim == 2)
512 *   interpretation.push_back(DataComponentInterpretation::component_is_scalar);
513 *   else {
517 *   }
518 *  
519 *   return interpretation;
520 *   }
521 *  
522 * @endcode
523 *
524 * This auxiliary function is required by the base class DataProcessor and simply
525 * sets which variables have to updated (only the gradients)
526 *
527
528 *
529 *
530 * @code
531 *   template<int dim>
532 *   UpdateFlags PostprocessorVorticity<dim>::get_needed_update_flags() const {
533 *   return update_gradients;
534 *   }
535 *  
536 *  
537 * @endcode
538 *
539 * The following structs are auxiliary objects for mesh refinement. ScratchData simply sets
540 * the FEValues object
541 *
542
543 *
544 *
545 * @code
546 *   template <int dim>
547 *   struct ScratchData {
548 *   ScratchData(const FiniteElement<dim>& fe,
549 *   const unsigned int quadrature_degree,
550 *   const UpdateFlags update_flags): fe_values(fe, QGauss<dim>(quadrature_degree), update_flags) {}
551 *  
552 *   ScratchData(const ScratchData<dim>& scratch_data): fe_values(scratch_data.fe_values.get_fe(),
553 *   scratch_data.fe_values.get_quadrature(),
554 *   scratch_data.fe_values.get_update_flags()) {}
555 *   FEValues<dim> fe_values;
556 *   };
557 *  
558 *  
559 * @endcode
560 *
561 * CopyData simply sets the cell index
562 *
563
564 *
565 *
566 * @code
567 *   struct CopyData {
568 *   CopyData() : cell_index(numbers::invalid_unsigned_int), value(0.0) {}
569 *  
570 *   CopyData(const CopyData &) = default;
571 *  
572 *   unsigned int cell_index;
573 *   double value;
574 *   };
575 *  
576 *  
577 * @endcode
578 *
579 *
580 * <a name="navier_stokes_TRBDF2_DG.cc-"></a>
581 * @sect{ <code>NavierStokesProjectionOperator::NavierStokesProjectionOperator</code> }
582 *
583
584 *
585 * The following class sets effectively the weak formulation of the problems for the different stages
586 * and for both velocity and pressure.
587 * The template parameters are the dimnesion of the problem, the polynomial degree for the pressure,
588 * the polynomial degree for the velocity, the number of quadrature points for integrals for the pressure step,
589 * the number of quadrature points for integrals for the velocity step, the type of vector for storage and the type
590 * of floating point data (in general double or float for preconditioners structures if desired).
591 *
592
593 *
594 *
595 * @code
596 *   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>
597 *   class NavierStokesProjectionOperator: public MatrixFreeOperators::Base<dim, Vec> {
598 *   public:
599 *   using Number = typename Vec::value_type;
600 *  
601 *   NavierStokesProjectionOperator();
602 *  
603 *   NavierStokesProjectionOperator(RunTimeParameters::Data_Storage& data);
604 *  
605 *   void set_dt(const double time_step);
606 *  
607 *   void set_TR_BDF2_stage(const unsigned int stage);
608 *  
609 *   void set_NS_stage(const unsigned int stage);
610 *  
611 *   void set_u_extr(const Vec& src);
612 *  
613 *   void vmult_rhs_velocity(Vec& dst, const std::vector<Vec>& src) const;
614 *  
615 *   void vmult_rhs_pressure(Vec& dst, const std::vector<Vec>& src) const;
616 *  
617 *   void vmult_grad_p_projection(Vec& dst, const Vec& src) const;
618 *  
619 *   virtual void compute_diagonal() override;
620 *  
621 *   protected:
622 *   double Re;
623 *   double dt;
624 *  
625 *   /*--- Parameters of time-marching scheme ---*/
626 *   double gamma;
627 *   double a31;
628 *   double a32;
629 *   double a33;
630 *  
631 *   unsigned int TR_BDF2_stage; /*--- Flag to denote at which stage of the TR-BDF2 are ---*/
632 *   unsigned int NS_stage; /*--- Flag to denote at which stage of NS solution inside each TR-BDF2 stage we are
633 *   (solution of the velocity or of the pressure)---*/
634 *  
635 *   virtual void apply_add(Vec& dst, const Vec& src) const override;
636 *  
637 *   private:
638 *   /*--- Auxiliary variable for the TR stage
639 *   (just to avoid to report a lot of 0.5 and for my personal choice to be coherent with the article) ---*/
640 *   const double a21 = 0.5;
641 *   const double a22 = 0.5;
642 *  
643 *   /*--- Penalty method parameters, theta = 1 means SIP, while C_p and C_u are the penalization coefficients ---*/
644 *   const double theta_v = 1.0;
645 *   const double theta_p = 1.0;
646 *   const double C_p = 1.0*(fe_degree_p + 1)*(fe_degree_p + 1);
647 *   const double C_u = 1.0*(fe_degree_v + 1)*(fe_degree_v + 1);
648 *  
649 *   Vec u_extr; /*--- Auxiliary variable to update the extrapolated velocity ---*/
650 *  
651 *   EquationData::Velocity<dim> vel_boundary_inflow; /*--- Auxiliary variable to impose velocity boundary conditions ---*/
652 *  
653 *   /*--- The following functions basically assemble the linear and bilinear forms. Their syntax is due to
654 *   the base class MatrixFreeOperators::Base ---*/
655 *   void assemble_rhs_cell_term_velocity(const MatrixFree<dim, Number>& data,
656 *   Vec& dst,
657 *   const std::vector<Vec>& src,
658 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
659 *   void assemble_rhs_face_term_velocity(const MatrixFree<dim, Number>& data,
660 *   Vec& dst,
661 *   const std::vector<Vec>& src,
662 *   const std::pair<unsigned int, unsigned int>& face_range) const;
663 *   void assemble_rhs_boundary_term_velocity(const MatrixFree<dim, Number>& data,
664 *   Vec& dst,
665 *   const std::vector<Vec>& src,
666 *   const std::pair<unsigned int, unsigned int>& face_range) const;
667 *  
668 *   void assemble_rhs_cell_term_pressure(const MatrixFree<dim, Number>& data,
669 *   Vec& dst,
670 *   const std::vector<Vec>& src,
671 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
672 *   void assemble_rhs_face_term_pressure(const MatrixFree<dim, Number>& data,
673 *   Vec& dst,
674 *   const std::vector<Vec>& src,
675 *   const std::pair<unsigned int, unsigned int>& face_range) const;
676 *   void assemble_rhs_boundary_term_pressure(const MatrixFree<dim, Number>& data,
677 *   Vec& dst,
678 *   const std::vector<Vec>& src,
679 *   const std::pair<unsigned int, unsigned int>& face_range) const;
680 *  
681 *   void assemble_cell_term_velocity(const MatrixFree<dim, Number>& data,
682 *   Vec& dst,
683 *   const Vec& src,
684 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
685 *   void assemble_face_term_velocity(const MatrixFree<dim, Number>& data,
686 *   Vec& dst,
687 *   const Vec& src,
688 *   const std::pair<unsigned int, unsigned int>& face_range) const;
689 *   void assemble_boundary_term_velocity(const MatrixFree<dim, Number>& data,
690 *   Vec& dst,
691 *   const Vec& src,
692 *   const std::pair<unsigned int, unsigned int>& face_range) const;
693 *  
694 *   void assemble_cell_term_pressure(const MatrixFree<dim, Number>& data,
695 *   Vec& dst,
696 *   const Vec& src,
697 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
698 *   void assemble_face_term_pressure(const MatrixFree<dim, Number>& data,
699 *   Vec& dst,
700 *   const Vec& src,
701 *   const std::pair<unsigned int, unsigned int>& face_range) const;
702 *   void assemble_boundary_term_pressure(const MatrixFree<dim, Number>& data,
703 *   Vec& dst,
704 *   const Vec& src,
705 *   const std::pair<unsigned int, unsigned int>& face_range) const;
706 *  
707 *   void assemble_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
708 *   Vec& dst,
709 *   const Vec& src,
710 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
711 *   void assemble_rhs_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
712 *   Vec& dst,
713 *   const Vec& src,
714 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
715 *  
716 *   void assemble_diagonal_cell_term_velocity(const MatrixFree<dim, Number>& data,
717 *   Vec& dst,
718 *   const unsigned int& src,
719 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
720 *   void assemble_diagonal_face_term_velocity(const MatrixFree<dim, Number>& data,
721 *   Vec& dst,
722 *   const unsigned int& src,
723 *   const std::pair<unsigned int, unsigned int>& face_range) const;
724 *   void assemble_diagonal_boundary_term_velocity(const MatrixFree<dim, Number>& data,
725 *   Vec& dst,
726 *   const unsigned int& src,
727 *   const std::pair<unsigned int, unsigned int>& face_range) const;
728 *  
729 *   void assemble_diagonal_cell_term_pressure(const MatrixFree<dim, Number>& data,
730 *   Vec& dst,
731 *   const unsigned int& src,
732 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
733 *   void assemble_diagonal_face_term_pressure(const MatrixFree<dim, Number>& data,
734 *   Vec& dst,
735 *   const unsigned int& src,
736 *   const std::pair<unsigned int, unsigned int>& face_range) const;
737 *   void assemble_diagonal_boundary_term_pressure(const MatrixFree<dim, Number>& data,
738 *   Vec& dst,
739 *   const unsigned int& src,
740 *   const std::pair<unsigned int, unsigned int>& face_range) const;
741 *   };
742 *  
743 *  
744 * @endcode
745 *
746 * We start with the default constructor. It is important for MultiGrid, so it is fundamental
747 * to properly set the parameters of the time scheme.
748 *
749
750 *
751 *
752 * @code
753 *   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>
754 *   NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
755 *   NavierStokesProjectionOperator():
756 *   MatrixFreeOperators::Base<dim, Vec>(), Re(), dt(), gamma(2.0 - std::sqrt(2.0)), a31((1.0 - gamma)/(2.0*(2.0 - gamma))),
757 *   a32(a31), a33(1.0/(2.0 - gamma)), TR_BDF2_stage(1), NS_stage(1), u_extr() {}
758 *  
759 *  
760 * @endcode
761 *
762 * We focus now on the constructor with runtime parameters storage
763 *
764
765 *
766 *
767 * @code
768 *   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>
769 *   NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
770 *   NavierStokesProjectionOperator(RunTimeParameters::Data_Storage& data):
771 *   MatrixFreeOperators::Base<dim, Vec>(), Re(data.Reynolds), dt(data.dt),
772 *   gamma(2.0 - std::sqrt(2.0)), a31((1.0 - gamma)/(2.0*(2.0 - gamma))),
773 *   a32(a31), a33(1.0/(2.0 - gamma)), TR_BDF2_stage(1), NS_stage(1), u_extr(),
774 *   vel_boundary_inflow(data.initial_time) {}
775 *  
776 *  
777 * @endcode
778 *
779 * Setter of time-step (called by Multigrid and in case a smaller time-step towards the end is needed)
780 *
781
782 *
783 *
784 * @code
785 *   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>
786 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
787 *   set_dt(const double time_step) {
788 *   dt = time_step;
789 *   }
790 *  
791 *  
792 * @endcode
793 *
794 * Setter of TR-BDF2 stage (this can be known only during the effective execution
795 * and so it has to be demanded to the class that really solves the problem)
796 *
797
798 *
799 *
800 * @code
801 *   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>
802 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
803 *   set_TR_BDF2_stage(const unsigned int stage) {
804 *   AssertIndexRange(stage, 3);
805 *   Assert(stage > 0, ExcInternalError());
806 *  
807 *   TR_BDF2_stage = stage;
808 *   }
809 *  
810 *  
811 * @endcode
812 *
813 * Setter of NS stage (this can be known only during the effective execution
814 * and so it has to be demanded to the class that really solves the problem)
815 *
816
817 *
818 *
819 * @code
820 *   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>
821 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
822 *   set_NS_stage(const unsigned int stage) {
823 *   AssertIndexRange(stage, 4);
824 *   Assert(stage > 0, ExcInternalError());
825 *  
826 *   NS_stage = stage;
827 *   }
828 *  
829 *  
830 * @endcode
831 *
832 * Setter of extrapolated velocity for different stages
833 *
834
835 *
836 *
837 * @code
838 *   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>
839 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
840 *   set_u_extr(const Vec& src) {
841 *   u_extr = src;
842 *   u_extr.update_ghost_values();
843 *   }
844 *  
845 *  
846 * @endcode
847 *
848 * We are in a DG-MatrixFree framework, so it is convenient to compute separately cell contribution,
849 * internal faces contributions and boundary faces contributions. We start by
850 * assembling the rhs cell term for the velocity.
851 *
852
853 *
854 *
855 * @code
856 *   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>
857 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
858 *   assemble_rhs_cell_term_velocity(const MatrixFree<dim, Number>& data,
859 *   Vec& dst,
860 *   const std::vector<Vec>& src,
861 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
862 *   if(TR_BDF2_stage == 1) {
863 *   /*--- We first start by declaring the suitable instances to read the old velocity, the
864 *   extrapolated velocity and the old pressure. 'phi' will be used only to submit the result.
865 *   The second argument specifies which dof handler has to be used (in this implementation 0 stands for
866 *   velocity and 1 for pressure). ---*/
868 *   phi_old(data, 0),
869 *   phi_old_extr(data, 0);
871 *  
872 *   /*--- We loop over the cells in the range ---*/
873 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
874 *   /*--- Now we need to assign the current cell to each FEEvaluation object and then to specify which src vector
875 *   it has to read (the proper order is clearly delegated to the user, which has to pay attention in the function
876 *   call to be coherent). ---*/
877 *   phi_old.reinit(cell);
878 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
879 *   /*--- The 'gather_evaluate' function reads data from the vector.
880 *   The second and third parameter specifies if you want to read
881 *   values and/or derivative related quantities ---*/
882 *   phi_old_extr.reinit(cell);
883 *   phi_old_extr.gather_evaluate(src[1], EvaluationFlags::values);
884 *   phi_old_press.reinit(cell);
885 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
886 *   phi.reinit(cell);
887 *  
888 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
889 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
890 *   const auto& u_n = phi_old.get_value(q);
891 *   const auto& grad_u_n = phi_old.get_gradient(q);
892 *   const auto& u_n_gamma_ov_2 = phi_old_extr.get_value(q);
893 *   const auto& tensor_product_u_n = outer_product(u_n, u_n_gamma_ov_2);
894 *   const auto& p_n = phi_old_press.get_value(q);
895 *   auto p_n_times_identity = tensor_product_u_n;
896 *   p_n_times_identity = 0;
897 *   for(unsigned int d = 0; d < dim; ++d)
898 *   p_n_times_identity[d][d] = p_n;
899 *  
900 *   phi.submit_value(1.0/(gamma*dt)*u_n, q); /*--- 'submit_value' contains quantites that we want to test against the
901 *   test function ---*/
902 *   phi.submit_gradient(-a21/Re*grad_u_n + a21*tensor_product_u_n + p_n_times_identity, q);
903 *   /*--- 'submit_gradient' contains quantites that we want to test against the gradient of test function ---*/
904 *   }
905 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
906 *   /*--- 'integrate_scatter' is the responsible of distributing into dst.
907 *   The flag parameter specifies if we are testing against the test function and/or its gradient ---*/
908 *   }
909 *   }
910 *   else {
911 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
913 *   phi_old(data, 0),
914 *   phi_int(data, 0);
916 *  
917 *   /*--- We loop over the cells in the range ---*/
918 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
919 *   phi_old.reinit(cell);
920 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
921 *   phi_int.reinit(cell);
922 *   phi_int.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
923 *   phi_old_press.reinit(cell);
924 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
925 *   phi.reinit(cell);
926 *  
927 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
928 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
929 *   const auto& u_n = phi_old.get_value(q);
930 *   const auto& grad_u_n = phi_old.get_gradient(q);
931 *   const auto& u_n_gamma = phi_int.get_value(q);
932 *   const auto& grad_u_n_gamma = phi_int.get_gradient(q);
933 *   const auto& tensor_product_u_n = outer_product(u_n, u_n);
934 *   const auto& tensor_product_u_n_gamma = outer_product(u_n_gamma, u_n_gamma);
935 *   const auto& p_n = phi_old_press.get_value(q);
936 *   auto p_n_times_identity = tensor_product_u_n;
937 *   p_n_times_identity = 0;
938 *   for(unsigned int d = 0; d < dim; ++d)
939 *   p_n_times_identity[d][d] = p_n;
940 *  
941 *   phi.submit_value(1.0/((1.0 - gamma)*dt)*u_n_gamma, q);
942 *   phi.submit_gradient(a32*tensor_product_u_n_gamma + a31*tensor_product_u_n -
943 *   a32/Re*grad_u_n_gamma - a31/Re*grad_u_n + p_n_times_identity, q);
944 *   }
945 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
946 *   }
947 *   }
948 *   }
949 *  
950 *  
951 * @endcode
952 *
953 * The followinf function assembles rhs face term for the velocity
954 *
955
956 *
957 *
958 * @code
959 *   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>
960 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
961 *   assemble_rhs_face_term_velocity(const MatrixFree<dim, Number>& data,
962 *   Vec& dst,
963 *   const std::vector<Vec>& src,
964 *   const std::pair<unsigned int, unsigned int>& face_range) const {
965 *   if(TR_BDF2_stage == 1) {
966 *   /*--- We first start by declaring the suitable instances to read already available quantities. In this case
967 *   we are at the face between two elements and this is the reason of 'FEFaceEvaluation'. It contains an extra
968 *   input argument, the second one, that specifies if it is from 'interior' or not---*/
970 *   phi_m(data, false, 0),
971 *   phi_old_p(data, true, 0),
972 *   phi_old_m(data, false, 0),
973 *   phi_old_extr_p(data, true, 0),
974 *   phi_old_extr_m(data, false, 0);
975 *   FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number> phi_old_press_p(data, true, 1),
976 *   phi_old_press_m(data, false, 1);
977 *  
978 *   /*--- We loop over the faces in the range ---*/
979 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
980 *   phi_old_p.reinit(face);
981 *   phi_old_p.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
982 *   phi_old_m.reinit(face);
983 *   phi_old_m.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
984 *   phi_old_extr_p.reinit(face);
985 *   phi_old_extr_p.gather_evaluate(src[1], EvaluationFlags::values);
986 *   phi_old_extr_m.reinit(face);
987 *   phi_old_extr_m.gather_evaluate(src[1], EvaluationFlags::values);
988 *   phi_old_press_p.reinit(face);
989 *   phi_old_press_p.gather_evaluate(src[2], EvaluationFlags::values);
990 *   phi_old_press_m.reinit(face);
991 *   phi_old_press_m.gather_evaluate(src[2], EvaluationFlags::values);
992 *   phi_p.reinit(face);
993 *   phi_m.reinit(face);
994 *  
995 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
996 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
997 *   const auto& n_plus = phi_p.get_normal_vector(q); /*--- The normal vector is the same
998 *   for both phi_p and phi_m. If the face is interior,
999 *   it correspond to the outer normal ---*/
1000 *  
1001 *   const auto& avg_grad_u_old = 0.5*(phi_old_p.get_gradient(q) + phi_old_m.get_gradient(q));
1002 *   const auto& avg_tensor_product_u_n = 0.5*(outer_product(phi_old_p.get_value(q), phi_old_extr_p.get_value(q)) +
1003 *   outer_product(phi_old_m.get_value(q), phi_old_extr_m.get_value(q)));
1004 *   const auto& avg_p_old = 0.5*(phi_old_press_p.get_value(q) + phi_old_press_m.get_value(q));
1005 *  
1006 *   phi_p.submit_value((a21/Re*avg_grad_u_old - a21*avg_tensor_product_u_n)*n_plus - avg_p_old*n_plus, q);
1007 *   phi_m.submit_value(-(a21/Re*avg_grad_u_old - a21*avg_tensor_product_u_n)*n_plus + avg_p_old*n_plus, q);
1008 *   }
1009 *   phi_p.integrate_scatter(EvaluationFlags::values, dst);
1010 *   phi_m.integrate_scatter(EvaluationFlags::values, dst);
1011 *   }
1012 *   }
1013 *   else {
1014 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1016 *   phi_m(data, false, 0),
1017 *   phi_old_p(data, true, 0),
1018 *   phi_old_m(data, false, 0),
1019 *   phi_int_p(data, true, 0),
1020 *   phi_int_m(data, false, 0);
1021 *   FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number> phi_old_press_p(data, true, 1),
1022 *   phi_old_press_m(data, false, 1);
1023 *  
1024 *   /*--- We loop over the faces in the range ---*/
1025 *   for(unsigned int face = face_range.first; face < face_range.second; ++ face) {
1026 *   phi_old_p.reinit(face);
1027 *   phi_old_p.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1028 *   phi_old_m.reinit(face);
1029 *   phi_old_m.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1030 *   phi_int_p.reinit(face);
1031 *   phi_int_p.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
1032 *   phi_int_m.reinit(face);
1033 *   phi_int_m.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
1034 *   phi_old_press_p.reinit(face);
1035 *   phi_old_press_p.gather_evaluate(src[2], EvaluationFlags::values);
1036 *   phi_old_press_m.reinit(face);
1037 *   phi_old_press_m.gather_evaluate(src[2], EvaluationFlags::values);
1038 *   phi_p.reinit(face);
1039 *   phi_m.reinit(face);
1040 *  
1041 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1042 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1043 *   const auto& n_plus = phi_p.get_normal_vector(q);
1044 *  
1045 *   const auto& avg_grad_u_old = 0.5*(phi_old_p.get_gradient(q) + phi_old_m.get_gradient(q));
1046 *   const auto& avg_grad_u_int = 0.5*(phi_int_p.get_gradient(q) + phi_int_m.get_gradient(q));
1047 *   const auto& avg_tensor_product_u_n = 0.5*(outer_product(phi_old_p.get_value(q), phi_old_p.get_value(q)) +
1048 *   outer_product(phi_old_m.get_value(q), phi_old_m.get_value(q)));
1049 *   const auto& avg_tensor_product_u_n_gamma = 0.5*(outer_product(phi_int_p.get_value(q), phi_int_p.get_value(q)) +
1050 *   outer_product(phi_int_m.get_value(q), phi_int_m.get_value(q)));
1051 *   const auto& avg_p_old = 0.5*(phi_old_press_p.get_value(q) + phi_old_press_m.get_value(q));
1052 *  
1053 *   phi_p.submit_value((a31/Re*avg_grad_u_old + a32/Re*avg_grad_u_int -
1054 *   a31*avg_tensor_product_u_n - a32*avg_tensor_product_u_n_gamma)*n_plus - avg_p_old*n_plus, q);
1055 *   phi_m.submit_value(-(a31/Re*avg_grad_u_old + a32/Re*avg_grad_u_int -
1056 *   a31*avg_tensor_product_u_n - a32*avg_tensor_product_u_n_gamma)*n_plus + avg_p_old*n_plus, q);
1057 *   }
1058 *   phi_p.integrate_scatter(EvaluationFlags::values, dst);
1059 *   phi_m.integrate_scatter(EvaluationFlags::values, dst);
1060 *   }
1061 *   }
1062 *   }
1063 *  
1064 *  
1065 * @endcode
1066 *
1067 * The followinf function assembles rhs boundary term for the velocity
1068 *
1069
1070 *
1071 *
1072 * @code
1073 *   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>
1074 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1075 *   assemble_rhs_boundary_term_velocity(const MatrixFree<dim, Number>& data,
1076 *   Vec& dst,
1077 *   const std::vector<Vec>& src,
1078 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1079 *   if(TR_BDF2_stage == 1) {
1080 *   /*--- We first start by declaring the suitable instances to read already available quantities. Clearly on the boundary
1081 *   the second argument has to be true. ---*/
1083 *   phi_old(data, true, 0),
1084 *   phi_old_extr(data, true, 0);
1086 *  
1087 *   /*--- We loop over the faces in the range ---*/
1088 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1089 *   phi_old.reinit(face);
1090 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1091 *   phi_old_extr.reinit(face);
1092 *   phi_old_extr.gather_evaluate(src[1], EvaluationFlags::values);
1093 *   phi_old_press.reinit(face);
1094 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
1095 *   phi.reinit(face);
1096 *  
1097 *   const auto boundary_id = data.get_boundary_id(face); /*--- Get the id in order to impose the proper boundary condition ---*/
1098 *   const auto coef_jump = (boundary_id == 1) ? 0.0 : C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1099 *   const double aux_coeff = (boundary_id == 1) ? 0.0 : 1.0;
1100 *  
1101 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1102 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1103 *   const auto& n_plus = phi.get_normal_vector(q);
1104 *  
1105 *   const auto& grad_u_old = phi_old.get_gradient(q);
1106 *   const auto& tensor_product_u_n = outer_product(phi_old.get_value(q), phi_old_extr.get_value(q));
1107 *   const auto& p_old = phi_old_press.get_value(q);
1108 *   const auto& point_vectorized = phi.quadrature_point(q);
1109 *   auto u_int_m = Tensor<1, dim, VectorizedArray<Number>>();
1110 *   if(boundary_id == 0) {
1111 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1112 *   Point<dim> point; /*--- The point returned by the 'quadrature_point' function is not an instance of Point
1113 *   and so it is not ready to be directly used. We need to pay attention to the
1114 *   vectorization ---*/
1115 *   for(unsigned int d = 0; d < dim; ++d)
1116 *   point[d] = point_vectorized[d][v];
1117 *   for(unsigned int d = 0; d < dim; ++d)
1118 *   u_int_m[d][v] = vel_boundary_inflow.value(point, d);
1119 *   }
1120 *   }
1121 *   const auto tensor_product_u_int_m = outer_product(u_int_m, phi_old_extr.get_value(q));
1122 *   const auto lambda = (boundary_id == 1) ? 0.0 : std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
1123 *  
1124 *   phi.submit_value((a21/Re*grad_u_old - a21*tensor_product_u_n)*n_plus - p_old*n_plus +
1125 *   a22/Re*2.0*coef_jump*u_int_m -
1126 *   aux_coeff*a22*tensor_product_u_int_m*n_plus + a22*lambda*u_int_m, q);
1127 *   phi.submit_normal_derivative(-aux_coeff*theta_v*a22/Re*u_int_m, q); /*--- This is equivalent to multiply to the gradient
1128 *   with outer product and use 'submit_gradient' ---*/
1129 *   }
1130 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1131 *   }
1132 *   }
1133 *   else {
1134 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1136 *   phi_old(data, true, 0),
1137 *   phi_int(data, true, 0),
1138 *   phi_int_extr(data, true, 0);
1140 *  
1141 *   /*--- We loop over the faces in the range ---*/
1142 *   for(unsigned int face = face_range.first; face < face_range.second; ++ face) {
1143 *   phi_old.reinit(face);
1144 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1145 *   phi_int.reinit(face);
1146 *   phi_int.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
1147 *   phi_old_press.reinit(face);
1148 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
1149 *   phi_int_extr.reinit(face);
1150 *   phi_int_extr.gather_evaluate(src[3], EvaluationFlags::values);
1151 *   phi.reinit(face);
1152 *  
1153 *   const auto boundary_id = data.get_boundary_id(face);
1154 *   const auto coef_jump = (boundary_id == 1) ? 0.0 : C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1155 *   const double aux_coeff = (boundary_id == 1) ? 0.0 : 1.0;
1156 *  
1157 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1158 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1159 *   const auto& n_plus = phi.get_normal_vector(q);
1160 *  
1161 *   const auto& grad_u_old = phi_old.get_gradient(q);
1162 *   const auto& grad_u_int = phi_int.get_gradient(q);
1163 *   const auto& tensor_product_u_n = outer_product(phi_old.get_value(q), phi_old.get_value(q));
1164 *   const auto& tensor_product_u_n_gamma = outer_product(phi_int.get_value(q), phi_int.get_value(q));
1165 *   const auto& p_old = phi_old_press.get_value(q);
1166 *   const auto& point_vectorized = phi.quadrature_point(q);
1167 *   auto u_m = Tensor<1, dim, VectorizedArray<Number>>();
1168 *   if(boundary_id == 0) {
1169 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1170 *   Point<dim> point;
1171 *   for(unsigned int d = 0; d < dim; ++d)
1172 *   point[d] = point_vectorized[d][v];
1173 *   for(unsigned int d = 0; d < dim; ++d)
1174 *   u_m[d][v] = vel_boundary_inflow.value(point, d);
1175 *   }
1176 *   }
1177 *   const auto tensor_product_u_m = outer_product(u_m, phi_int_extr.get_value(q));
1178 *   const auto lambda = (boundary_id == 1) ? 0.0 : std::abs(scalar_product(phi_int_extr.get_value(q), n_plus));
1179 *  
1180 *   phi.submit_value((a31/Re*grad_u_old + a32/Re*grad_u_int -
1181 *   a31*tensor_product_u_n - a32*tensor_product_u_n_gamma)*n_plus - p_old*n_plus +
1182 *   a33/Re*2.0*coef_jump*u_m -
1183 *   aux_coeff*a33*tensor_product_u_m*n_plus + a33*lambda*u_m, q);
1184 *   phi.submit_normal_derivative(-aux_coeff*theta_v*a33/Re*u_m, q);
1185 *   }
1186 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1187 *   }
1188 *   }
1189 *   }
1190 *  
1191 *  
1192 * @endcode
1193 *
1194 * Put together all the previous steps for velocity. This is done automatically by the loop function of 'MatrixFree' class
1195 *
1196
1197 *
1198 *
1199 * @code
1200 *   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>
1201 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1202 *   vmult_rhs_velocity(Vec& dst, const std::vector<Vec>& src) const {
1203 *   for(auto& vec : src)
1204 *   vec.update_ghost_values();
1205 *  
1206 *   this->data->loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_velocity,
1207 *   &NavierStokesProjectionOperator::assemble_rhs_face_term_velocity,
1208 *   &NavierStokesProjectionOperator::assemble_rhs_boundary_term_velocity,
1209 *   this, dst, src, true,
1212 *   }
1213 *  
1214 *  
1215 * @endcode
1216 *
1217 * Now we focus on computing the rhs for the projection step for the pressure with the same ratio.
1218 * The following function assembles rhs cell term for the pressure
1219 *
1220
1221 *
1222 *
1223 * @code
1224 *   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>
1225 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1226 *   assemble_rhs_cell_term_pressure(const MatrixFree<dim, Number>& data,
1227 *   Vec& dst,
1228 *   const std::vector<Vec>& src,
1229 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1230 *   /*--- We first start by declaring the suitable instances to read already available quantities.
1231 *   The third parameter specifies that we want to use the second quadrature formula stored. ---*/
1233 *   phi_old(data, 1, 1);
1235 *  
1236 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0e6*gamma*dt*gamma*dt : 1.0e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
1237 *  
1238 *   const double coeff_2 = (TR_BDF2_stage == 1) ? gamma*dt : (1.0 - gamma)*dt;
1239 *  
1240 *   /*--- We loop over cells in the range ---*/
1241 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1242 *   phi_proj.reinit(cell);
1243 *   phi_proj.gather_evaluate(src[0], EvaluationFlags::values);
1244 *   phi_old.reinit(cell);
1245 *   phi_old.gather_evaluate(src[1], EvaluationFlags::values);
1246 *   phi.reinit(cell);
1247 *  
1248 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1249 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1250 *   const auto& u_star_star = phi_proj.get_value(q);
1251 *   const auto& p_old = phi_old.get_value(q);
1252 *  
1253 *   phi.submit_value(1.0/coeff*p_old, q);
1254 *   phi.submit_gradient(1.0/coeff_2*u_star_star, q);
1255 *   }
1256 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1257 *   }
1258 *   }
1259 *  
1260 *  
1261 * @endcode
1262 *
1263 * The following function assembles rhs face term for the pressure
1264 *
1265
1266 *
1267 *
1268 * @code
1269 *   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>
1270 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1271 *   assemble_rhs_face_term_pressure(const MatrixFree<dim, Number>& data,
1272 *   Vec& dst,
1273 *   const std::vector<Vec>& src,
1274 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1275 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1277 *   phi_m(data, false, 1, 1);
1279 *   phi_proj_m(data, false, 0, 1);
1280 *  
1281 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0/(gamma*dt) : 1.0/((1.0 - gamma)*dt);
1282 *  
1283 *   /*--- We loop over faces in the range ---*/
1284 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1285 *   phi_proj_p.reinit(face);
1286 *   phi_proj_p.gather_evaluate(src[0], EvaluationFlags::values);
1287 *   phi_proj_m.reinit(face);
1288 *   phi_proj_m.gather_evaluate(src[0], EvaluationFlags::values);
1289 *   phi_p.reinit(face);
1290 *   phi_m.reinit(face);
1291 *  
1292 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1293 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1294 *   const auto& n_plus = phi_p.get_normal_vector(q);
1295 *   const auto& avg_u_star_star = 0.5*(phi_proj_p.get_value(q) + phi_proj_m.get_value(q));
1296 *  
1297 *   phi_p.submit_value(-coeff*scalar_product(avg_u_star_star, n_plus), q);
1298 *   phi_m.submit_value(coeff*scalar_product(avg_u_star_star, n_plus), q);
1299 *   }
1300 *   phi_p.integrate_scatter(EvaluationFlags::values, dst);
1301 *   phi_m.integrate_scatter(EvaluationFlags::values, dst);
1302 *   }
1303 *   }
1304 *  
1305 *  
1306 * @endcode
1307 *
1308 * The following function assembles rhs boundary term for the pressure
1309 *
1310
1311 *
1312 *
1313 * @code
1314 *   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>
1315 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1316 *   assemble_rhs_boundary_term_pressure(const MatrixFree<dim, Number>& data,
1317 *   Vec& dst,
1318 *   const std::vector<Vec>& src,
1319 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1320 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1323 *  
1324 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0/(gamma*dt) : 1.0/((1.0 - gamma)*dt);
1325 *  
1326 *   /*--- We loop over faces in the range ---*/
1327 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1328 *   phi_proj.reinit(face);
1329 *   phi_proj.gather_evaluate(src[0], EvaluationFlags::values);
1330 *   phi.reinit(face);
1331 *  
1332 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1333 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1334 *   const auto& n_plus = phi.get_normal_vector(q);
1335 *  
1336 *   phi.submit_value(-coeff*scalar_product(phi_proj.get_value(q), n_plus), q);
1337 *   }
1338 *   phi.integrate_scatter(EvaluationFlags::values, dst);
1339 *   }
1340 *   }
1341 *  
1342 *  
1343 * @endcode
1344 *
1345 * Put together all the previous steps for pressure
1346 *
1347
1348 *
1349 *
1350 * @code
1351 *   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>
1352 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1353 *   vmult_rhs_pressure(Vec& dst, const std::vector<Vec>& src) const {
1354 *   for(auto& vec : src)
1355 *   vec.update_ghost_values();
1356 *  
1357 *   this->data->loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_pressure,
1358 *   &NavierStokesProjectionOperator::assemble_rhs_face_term_pressure,
1359 *   &NavierStokesProjectionOperator::assemble_rhs_boundary_term_pressure,
1360 *   this, dst, src, true,
1363 *   }
1364 *  
1365 *  
1366 * @endcode
1367 *
1368 * Now we need to build the 'matrices', i.e. the bilinear forms. We start by
1369 * assembling the cell term for the velocity
1370 *
1371
1372 *
1373 *
1374 * @code
1375 *   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>
1376 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1377 *   assemble_cell_term_velocity(const MatrixFree<dim, Number>& data,
1378 *   Vec& dst,
1379 *   const Vec& src,
1380 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1381 *   if(TR_BDF2_stage == 1) {
1382 *   /*--- We first start by declaring the suitable instances to read already available quantities. Moreover 'phi' in
1383 *   this case serves for a bilinear form and so it will not used only to submit but also to read the src ---*/
1385 *   phi_old_extr(data, 0);
1386 *  
1387 *   /*--- We loop over all cells in the range ---*/
1388 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1389 *   phi.reinit(cell);
1390 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1391 *   phi_old_extr.reinit(cell);
1392 *   phi_old_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1393 *  
1394 *   /*--- Now we loop over all quadrature points ---*/
1395 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1396 *   const auto& u_int = phi.get_value(q);
1397 *   const auto& grad_u_int = phi.get_gradient(q);
1398 *   const auto& u_n_gamma_ov_2 = phi_old_extr.get_value(q);
1399 *   const auto& tensor_product_u_int = outer_product(u_int, u_n_gamma_ov_2);
1400 *  
1401 *   phi.submit_value(1.0/(gamma*dt)*u_int, q);
1402 *   phi.submit_gradient(-a22*tensor_product_u_int + a22/Re*grad_u_int, q);
1403 *   }
1404 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1405 *   }
1406 *   }
1407 *   else {
1408 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1410 *   phi_int_extr(data, 0);
1411 *  
1412 *   /*--- We loop over all cells in the range ---*/
1413 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1414 *   phi.reinit(cell);
1415 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1416 *   phi_int_extr.reinit(cell);
1417 *   phi_int_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1418 *  
1419 *   /*--- Now we loop over all quadrature points ---*/
1420 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1421 *   const auto& u_curr = phi.get_value(q);
1422 *   const auto& grad_u_curr = phi.get_gradient(q);
1423 *   const auto& u_n1_int = phi_int_extr.get_value(q);
1424 *   const auto& tensor_product_u_curr = outer_product(u_curr, u_n1_int);
1425 *  
1426 *   phi.submit_value(1.0/((1.0 - gamma)*dt)*u_curr, q);
1427 *   phi.submit_gradient(-a33*tensor_product_u_curr + a33/Re*grad_u_curr, q);
1428 *   }
1429 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1430 *   }
1431 *   }
1432 *   }
1433 *  
1434 *  
1435 * @endcode
1436 *
1437 * The following function assembles face term for the velocity
1438 *
1439
1440 *
1441 *
1442 * @code
1443 *   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>
1444 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1445 *   assemble_face_term_velocity(const MatrixFree<dim, Number>& data,
1446 *   Vec& dst,
1447 *   const Vec& src,
1448 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1449 *   if(TR_BDF2_stage == 1) {
1450 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1452 *   phi_m(data, false, 0),
1453 *   phi_old_extr_p(data, true, 0),
1454 *   phi_old_extr_m(data, false, 0);
1455 *  
1456 *   /*--- We loop over all faces in the range ---*/
1457 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1458 *   phi_p.reinit(face);
1459 *   phi_p.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1460 *   phi_m.reinit(face);
1461 *   phi_m.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1462 *   phi_old_extr_p.reinit(face);
1463 *   phi_old_extr_p.gather_evaluate(u_extr, EvaluationFlags::values);
1464 *   phi_old_extr_m.reinit(face);
1465 *   phi_old_extr_m.gather_evaluate(u_extr, EvaluationFlags::values);
1466 *  
1467 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
1468 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
1469 *  
1470 *   /*--- Now we loop over all quadrature points ---*/
1471 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1472 *   const auto& n_plus = phi_p.get_normal_vector(q);
1473 *  
1474 *   const auto& avg_grad_u_int = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
1475 *   const auto& jump_u_int = phi_p.get_value(q) - phi_m.get_value(q);
1476 *   const auto& avg_tensor_product_u_int = 0.5*(outer_product(phi_p.get_value(q), phi_old_extr_p.get_value(q)) +
1477 *   outer_product(phi_m.get_value(q), phi_old_extr_m.get_value(q)));
1478 *   const auto lambda = std::max(std::abs(scalar_product(phi_old_extr_p.get_value(q), n_plus)),
1479 *   std::abs(scalar_product(phi_old_extr_m.get_value(q), n_plus)));
1480 *  
1481 *   phi_p.submit_value(a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) +
1482 *   a22*avg_tensor_product_u_int*n_plus + 0.5*a22*lambda*jump_u_int, q);
1483 *   phi_m.submit_value(-a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) -
1484 *   a22*avg_tensor_product_u_int*n_plus - 0.5*a22*lambda*jump_u_int, q);
1485 *   phi_p.submit_normal_derivative(-theta_v*a22/Re*0.5*jump_u_int, q);
1486 *   phi_m.submit_normal_derivative(-theta_v*a22/Re*0.5*jump_u_int, q);
1487 *   }
1488 *   phi_p.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1489 *   phi_m.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1490 *   }
1491 *   }
1492 *   else {
1493 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1495 *   phi_m(data, false, 0),
1496 *   phi_extr_p(data, true, 0),
1497 *   phi_extr_m(data, false, 0);
1498 *  
1499 *   /*--- We loop over all faces in the range ---*/
1500 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1501 *   phi_p.reinit(face);
1502 *   phi_p.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1503 *   phi_m.reinit(face);
1504 *   phi_m.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1505 *   phi_extr_p.reinit(face);
1506 *   phi_extr_p.gather_evaluate(u_extr, EvaluationFlags::values);
1507 *   phi_extr_m.reinit(face);
1508 *   phi_extr_m.gather_evaluate(u_extr, EvaluationFlags::values);
1509 *  
1510 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
1511 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
1512 *  
1513 *   /*--- Now we loop over all quadrature points ---*/
1514 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1515 *   const auto& n_plus = phi_p.get_normal_vector(q);
1516 *  
1517 *   const auto& avg_grad_u = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
1518 *   const auto& jump_u = phi_p.get_value(q) - phi_m.get_value(q);
1519 *   const auto& avg_tensor_product_u = 0.5*(outer_product(phi_p.get_value(q), phi_extr_p.get_value(q)) +
1520 *   outer_product(phi_m.get_value(q), phi_extr_m.get_value(q)));
1521 *   const auto lambda = std::max(std::abs(scalar_product(phi_extr_p.get_value(q), n_plus)),
1522 *   std::abs(scalar_product(phi_extr_m.get_value(q), n_plus)));
1523 *  
1524 *   phi_p.submit_value(a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) +
1525 *   a33*avg_tensor_product_u*n_plus + 0.5*a33*lambda*jump_u, q);
1526 *   phi_m.submit_value(-a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) -
1527 *   a33*avg_tensor_product_u*n_plus - 0.5*a33*lambda*jump_u, q);
1528 *   phi_p.submit_normal_derivative(-theta_v*a33/Re*0.5*jump_u, q);
1529 *   phi_m.submit_normal_derivative(-theta_v*a33/Re*0.5*jump_u, q);
1530 *   }
1531 *   phi_p.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1532 *   phi_m.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1533 *   }
1534 *   }
1535 *   }
1536 *  
1537 *  
1538 * @endcode
1539 *
1540 * The following function assembles boundary term for the velocity
1541 *
1542
1543 *
1544 *
1545 * @code
1546 *   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>
1547 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1548 *   assemble_boundary_term_velocity(const MatrixFree<dim, Number>& data,
1549 *   Vec& dst,
1550 *   const Vec& src,
1551 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1552 *   if(TR_BDF2_stage == 1) {
1553 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1555 *   phi_old_extr(data, true, 0);
1556 *  
1557 *   /*--- We loop over all faces in the range ---*/
1558 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1559 *   phi.reinit(face);
1560 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1561 *   phi_old_extr.reinit(face);
1562 *   phi_old_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1563 *  
1564 *   const auto boundary_id = data.get_boundary_id(face);
1565 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1566 *  
1567 *   /*--- The application of the mirror principle is not so trivial because we have a Dirichlet condition
1568 *   on a single component for the outflow; so we distinguish the two cases ---*/
1569 *   if(boundary_id != 1) {
1570 *   const double coef_trasp = 0.0;
1571 *  
1572 *   /*--- Now we loop over all quadrature points ---*/
1573 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1574 *   const auto& n_plus = phi.get_normal_vector(q);
1575 *   const auto& grad_u_int = phi.get_gradient(q);
1576 *   const auto& u_int = phi.get_value(q);
1577 *   const auto& tensor_product_u_int = outer_product(phi.get_value(q), phi_old_extr.get_value(q));
1578 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
1579 *  
1580 *   phi.submit_value(a22/Re*(-grad_u_int*n_plus + 2.0*coef_jump*u_int) +
1581 *   a22*coef_trasp*tensor_product_u_int*n_plus + a22*lambda*u_int, q);
1582 *   phi.submit_normal_derivative(-theta_v*a22/Re*u_int, q);
1583 *   }
1584 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1585 *   }
1586 *   else {
1587 *   /*--- Now we loop over all quadrature points ---*/
1588 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1589 *   const auto& n_plus = phi.get_normal_vector(q);
1590 *   const auto& grad_u_int = phi.get_gradient(q);
1591 *   const auto& u_int = phi.get_value(q);
1592 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
1593 *  
1594 *   const auto& point_vectorized = phi.quadrature_point(q);
1595 *   auto u_int_m = u_int;
1596 *   auto grad_u_int_m = grad_u_int;
1597 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1598 *   Point<dim> point;
1599 *   for(unsigned int d = 0; d < dim; ++d)
1600 *   point[d] = point_vectorized[d][v];
1601 *  
1602 *   u_int_m[1][v] = -u_int_m[1][v];
1603 *  
1604 *   grad_u_int_m[0][0][v] = -grad_u_int_m[0][0][v];
1605 *   grad_u_int_m[0][1][v] = -grad_u_int_m[0][1][v];
1606 *   }
1607 *  
1608 *   phi.submit_value(a22/Re*(-(0.5*(grad_u_int + grad_u_int_m))*n_plus + coef_jump*(u_int - u_int_m)) +
1609 *   a22*outer_product(0.5*(u_int + u_int_m), phi_old_extr.get_value(q))*n_plus +
1610 *   a22*0.5*lambda*(u_int - u_int_m), q);
1611 *   phi.submit_normal_derivative(-theta_v*a22/Re*(u_int - u_int_m), q);
1612 *   }
1613 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1614 *   }
1615 *   }
1616 *   }
1617 *   else {
1618 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1620 *   phi_extr(data, true, 0);
1621 *  
1622 *   /*--- We loop over all faces in the range ---*/
1623 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1624 *   phi.reinit(face);
1625 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1626 *   phi_extr.reinit(face);
1627 *   phi_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1628 *  
1629 *   const auto boundary_id = data.get_boundary_id(face);
1630 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1631 *  
1632 *   if(boundary_id != 1) {
1633 *   const double coef_trasp = 0.0;
1634 *  
1635 *   /*--- Now we loop over all quadrature points ---*/
1636 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1637 *   const auto& n_plus = phi.get_normal_vector(q);
1638 *   const auto& grad_u = phi.get_gradient(q);
1639 *   const auto& u = phi.get_value(q);
1640 *   const auto& tensor_product_u = outer_product(phi.get_value(q), phi_extr.get_value(q));
1641 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
1642 *  
1643 *   phi.submit_value(a33/Re*(-grad_u*n_plus + 2.0*coef_jump*u) +
1644 *   a33*coef_trasp*tensor_product_u*n_plus + a33*lambda*u, q);
1645 *   phi.submit_normal_derivative(-theta_v*a33/Re*u, q);
1646 *   }
1647 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1648 *   }
1649 *   else {
1650 *   /*--- Now we loop over all quadrature points ---*/
1651 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1652 *   const auto& n_plus = phi.get_normal_vector(q);
1653 *   const auto& grad_u = phi.get_gradient(q);
1654 *   const auto& u = phi.get_value(q);
1655 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
1656 *  
1657 *   const auto& point_vectorized = phi.quadrature_point(q);
1658 *   auto u_m = u;
1659 *   auto grad_u_m = grad_u;
1660 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1661 *   Point<dim> point;
1662 *   for(unsigned int d = 0; d < dim; ++d)
1663 *   point[d] = point_vectorized[d][v];
1664 *  
1665 *   u_m[1][v] = -u_m[1][v];
1666 *  
1667 *   grad_u_m[0][0][v] = -grad_u_m[0][0][v];
1668 *   grad_u_m[0][1][v] = -grad_u_m[0][1][v];
1669 *   }
1670 *  
1671 *   phi.submit_value(a33/Re*(-(0.5*(grad_u + grad_u_m))*n_plus + coef_jump*(u - u_m)) +
1672 *   a33*outer_product(0.5*(u + u_m), phi_extr.get_value(q))*n_plus + a33*0.5*lambda*(u - u_m), q);
1673 *   phi.submit_normal_derivative(-theta_v*a33/Re*(u - u_m), q);
1674 *   }
1675 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1676 *   }
1677 *   }
1678 *   }
1679 *   }
1680 *  
1681 *  
1682 * @endcode
1683 *
1684 * Next, we focus on 'matrices' to compute the pressure. We first assemble cell term for the pressure
1685 *
1686
1687 *
1688 *
1689 * @code
1690 *   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>
1691 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1692 *   assemble_cell_term_pressure(const MatrixFree<dim, Number>& data,
1693 *   Vec& dst,
1694 *   const Vec& src,
1695 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1696 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1698 *  
1699 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0e6*gamma*dt*gamma*dt : 1.0e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
1700 *  
1701 *   /*--- Loop over all cells in the range ---*/
1702 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1703 *   phi.reinit(cell);
1704 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1705 *  
1706 *   /*--- Now we loop over all quadrature points ---*/
1707 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1708 *   phi.submit_gradient(phi.get_gradient(q), q);
1709 *   phi.submit_value(1.0/coeff*phi.get_value(q), q);
1710 *   }
1711 *  
1712 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1713 *   }
1714 *   }
1715 *  
1716 *  
1717 * @endcode
1718 *
1719 * The following function assembles face term for the pressure
1720 *
1721
1722 *
1723 *
1724 * @code
1725 *   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>
1726 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1727 *   assemble_face_term_pressure(const MatrixFree<dim, Number>& data,
1728 *   Vec& dst,
1729 *   const Vec& src,
1730 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1731 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1733 *   phi_m(data, false, 1, 1);
1734 *  
1735 *   /*--- Loop over all faces in the range ---*/
1736 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1737 *   phi_p.reinit(face);
1738 *   phi_p.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1739 *   phi_m.reinit(face);
1740 *   phi_m.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1741 *  
1742 *   const auto coef_jump = C_p*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
1743 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
1744 *  
1745 *   /*--- Loop over quadrature points ---*/
1746 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1747 *   const auto& n_plus = phi_p.get_normal_vector(q);
1748 *  
1749 *   const auto& avg_grad_pres = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
1750 *   const auto& jump_pres = phi_p.get_value(q) - phi_m.get_value(q);
1751 *  
1752 *   phi_p.submit_value(-scalar_product(avg_grad_pres, n_plus) + coef_jump*jump_pres, q);
1753 *   phi_m.submit_value(scalar_product(avg_grad_pres, n_plus) - coef_jump*jump_pres, q);
1754 *   phi_p.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
1755 *   phi_m.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
1756 *   }
1757 *   phi_p.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1758 *   phi_m.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1759 *   }
1760 *   }
1761 *  
1762 *  
1763 * @endcode
1764 *
1765 * The following function assembles boundary term for the pressure
1766 *
1767
1768 *
1769 *
1770 * @code
1771 *   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>
1772 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1773 *   assemble_boundary_term_pressure(const MatrixFree<dim, Number>& data,
1774 *   Vec& dst,
1775 *   const Vec& src,
1776 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1778 *  
1779 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1780 *   const auto boundary_id = data.get_boundary_id(face);
1781 *  
1782 *   if(boundary_id == 1) {
1783 *   phi.reinit(face);
1784 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1785 *  
1786 *   const auto coef_jump = C_p*std::abs((phi.get_normal_vector(0)*phi.inverse_jacobian(0))[dim - 1]);
1787 *  
1788 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1789 *   const auto& n_plus = phi.get_normal_vector(q);
1790 *  
1791 *   const auto& grad_pres = phi.get_gradient(q);
1792 *   const auto& pres = phi.get_value(q);
1793 *  
1794 *   phi.submit_value(-scalar_product(grad_pres, n_plus) + coef_jump*pres , q);
1795 *   phi.submit_normal_derivative(-theta_p*pres, q);
1796 *   }
1797 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1798 *   }
1799 *   }
1800 *   }
1801 *  
1802 *  
1803 * @endcode
1804 *
1805 * Before coding the 'apply_add' function, which is the one that will perform the loop, we focus on
1806 * the linear system that arises to project the gradient of the pressure into the velocity space.
1807 * The following function assembles rhs cell term for the projection of gradient of pressure. Since no
1808 * integration by parts is performed, only a cell term contribution is present.
1809 *
1810
1811 *
1812 *
1813 * @code
1814 *   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>
1815 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1816 *   assemble_rhs_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
1817 *   Vec& dst,
1818 *   const Vec& src,
1819 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1820 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1823 *  
1824 *   /*--- Loop over all cells in the range ---*/
1825 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1826 *   phi_pres.reinit(cell);
1827 *   phi_pres.gather_evaluate(src, EvaluationFlags::gradients);
1828 *   phi.reinit(cell);
1829 *  
1830 *   /*--- Loop over quadrature points ---*/
1831 *   for(unsigned int q = 0; q < phi.n_q_points; ++q)
1832 *   phi.submit_value(phi_pres.get_gradient(q), q);
1833 *  
1834 *   phi.integrate_scatter(EvaluationFlags::values, dst);
1835 *   }
1836 *   }
1837 *  
1838 *  
1839 * @endcode
1840 *
1841 * Put together all the previous steps for projection of pressure gradient. Here we loop only over cells
1842 *
1843
1844 *
1845 *
1846 * @code
1847 *   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>
1848 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1849 *   vmult_grad_p_projection(Vec& dst, const Vec& src) const {
1850 *   this->data->cell_loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_projection_grad_p,
1851 *   this, dst, src, true);
1852 *   }
1853 *  
1854 *  
1855 * @endcode
1856 *
1857 * Assemble now cell term for the projection of gradient of pressure. This is nothing but a mass matrix
1858 *
1859
1860 *
1861 *
1862 * @code
1863 *   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>
1864 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1865 *   assemble_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
1866 *   Vec& dst,
1867 *   const Vec& src,
1868 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1870 *  
1871 *   /*--- Loop over all cells in the range ---*/
1872 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1873 *   phi.reinit(cell);
1874 *   phi.gather_evaluate(src, EvaluationFlags::values);
1875 *  
1876 *   /*--- Loop over quadrature points ---*/
1877 *   for(unsigned int q = 0; q < phi.n_q_points; ++q)
1878 *   phi.submit_value(phi.get_value(q), q);
1879 *  
1880 *   phi.integrate_scatter(EvaluationFlags::values, dst);
1881 *   }
1882 *   }
1883 *  
1884 *  
1885 * @endcode
1886 *
1887 * Put together all previous steps. This is the overridden function that effectively performs the
1888 * matrix-vector multiplication.
1889 *
1890
1891 *
1892 *
1893 * @code
1894 *   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>
1895 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1896 *   apply_add(Vec& dst, const Vec& src) const {
1897 *   if(NS_stage == 1) {
1898 *   this->data->loop(&NavierStokesProjectionOperator::assemble_cell_term_velocity,
1899 *   &NavierStokesProjectionOperator::assemble_face_term_velocity,
1900 *   &NavierStokesProjectionOperator::assemble_boundary_term_velocity,
1901 *   this, dst, src, false,
1904 *   }
1905 *   else if(NS_stage == 2) {
1906 *   this->data->loop(&NavierStokesProjectionOperator::assemble_cell_term_pressure,
1907 *   &NavierStokesProjectionOperator::assemble_face_term_pressure,
1908 *   &NavierStokesProjectionOperator::assemble_boundary_term_pressure,
1909 *   this, dst, src, false,
1912 *   }
1913 *   else if(NS_stage == 3) {
1914 *   this->data->cell_loop(&NavierStokesProjectionOperator::assemble_cell_term_projection_grad_p,
1915 *   this, dst, src, false); /*--- Since we have only a cell term contribution, we use cell_loop ---*/
1916 *   }
1917 *   else
1918 *   Assert(false, ExcNotImplemented());
1919 *   }
1920 *  
1921 *  
1922 * @endcode
1923 *
1924 * Finally, we focus on computing the diagonal for preconditioners and we start by assembling
1925 * the diagonal cell term for the velocity. Since we do not have access to the entries of the matrix,
1926 * 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.
1927 * This is why 'src' will result as unused.
1928 *
1929
1930 *
1931 *
1932 * @code
1933 *   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>
1934 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1935 *   assemble_diagonal_cell_term_velocity(const MatrixFree<dim, Number>& data,
1936 *   Vec& dst,
1937 *   const unsigned int& ,
1938 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1939 *   if(TR_BDF2_stage == 1) {
1941 *   phi_old_extr(data, 0);
1942 *  
1944 *   /*--- Build a vector of ones to be tested (here we will see the velocity as a whole vector, since
1945 *   dof_handler_velocity is vectorial and so the dof values are vectors). ---*/
1947 *   for(unsigned int d = 0; d < dim; ++d)
1948 *   tmp[d] = make_vectorized_array<Number>(1.0);
1949 *  
1950 *   /*--- Loop over cells in the range ---*/
1951 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1952 *   phi_old_extr.reinit(cell);
1953 *   phi_old_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1954 *   phi.reinit(cell);
1955 *  
1956 *   /*--- Loop over dofs ---*/
1957 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
1958 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
1959 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j); /*--- Set all dofs to zero ---*/
1960 *   phi.submit_dof_value(tmp, i); /*--- Set dof i equal to one ---*/
1962 *  
1963 *   /*--- Loop over quadrature points ---*/
1964 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1965 *   const auto& u_int = phi.get_value(q);
1966 *   const auto& grad_u_int = phi.get_gradient(q);
1967 *   const auto& u_n_gamma_ov_2 = phi_old_extr.get_value(q);
1968 *   const auto& tensor_product_u_int = outer_product(u_int, u_n_gamma_ov_2);
1969 *  
1970 *   phi.submit_value(1.0/(gamma*dt)*u_int, q);
1971 *   phi.submit_gradient(-a22*tensor_product_u_int + a22/Re*grad_u_int, q);
1972 *   }
1974 *   diagonal[i] = phi.get_dof_value(i);
1975 *   }
1976 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
1977 *   phi.submit_dof_value(diagonal[i], i);
1978 *   phi.distribute_local_to_global(dst);
1979 *   }
1980 *   }
1981 *   else {
1983 *   phi_int_extr(data, 0);
1984 *  
1987 *   for(unsigned int d = 0; d < dim; ++d)
1988 *   tmp[d] = make_vectorized_array<Number>(1.0);
1989 *  
1990 *   /*--- Loop over cells in the range ---*/
1991 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1992 *   phi_int_extr.reinit(cell);
1993 *   phi_int_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1994 *   phi.reinit(cell);
1995 *  
1996 *   /*--- Loop over dofs ---*/
1997 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
1998 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
1999 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2000 *   phi.submit_dof_value(tmp, i);
2002 *  
2003 *   /*--- Loop over quadrature points ---*/
2004 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2005 *   const auto& u_curr = phi.get_value(q);
2006 *   const auto& grad_u_curr = phi.get_gradient(q);
2007 *   const auto& u_n1_int = phi_int_extr.get_value(q);
2008 *   const auto& tensor_product_u_curr = outer_product(u_curr, u_n1_int);
2009 *  
2010 *   phi.submit_value(1.0/((1.0 - gamma)*dt)*u_curr, q);
2011 *   phi.submit_gradient(-a33*tensor_product_u_curr + a33/Re*grad_u_curr, q);
2012 *   }
2014 *   diagonal[i] = phi.get_dof_value(i);
2015 *   }
2016 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2017 *   phi.submit_dof_value(diagonal[i], i);
2018 *   phi.distribute_local_to_global(dst);
2019 *   }
2020 *   }
2021 *   }
2022 *  
2023 *  
2024 * @endcode
2025 *
2026 * The following function assembles diagonal face term for the velocity
2027 *
2028
2029 *
2030 *
2031 * @code
2032 *   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>
2033 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2034 *   assemble_diagonal_face_term_velocity(const MatrixFree<dim, Number>& data,
2035 *   Vec& dst,
2036 *   const unsigned int& ,
2037 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2038 *   if(TR_BDF2_stage == 1) {
2040 *   phi_m(data, false, 0),
2041 *   phi_old_extr_p(data, true, 0),
2042 *   phi_old_extr_m(data, false, 0);
2043 *  
2044 *   AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component); /*--- We just assert for safety that dimension match,
2045 *   in the sense that we have selected the proper
2046 *   space ---*/
2047 *   AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal_p(phi_p.dofs_per_component),
2048 *   diagonal_m(phi_m.dofs_per_component);
2049 *  
2051 *   for(unsigned int d = 0; d < dim; ++d)
2052 *   tmp[d] = make_vectorized_array<Number>(1.0); /*--- We build the usual vector of ones that we will use as dof value ---*/
2053 *  
2054 *   /*--- Now we loop over faces ---*/
2055 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2056 *   phi_old_extr_p.reinit(face);
2057 *   phi_old_extr_p.gather_evaluate(u_extr, EvaluationFlags::values);
2058 *   phi_old_extr_m.reinit(face);
2059 *   phi_old_extr_m.gather_evaluate(u_extr, EvaluationFlags::values);
2060 *   phi_p.reinit(face);
2061 *   phi_m.reinit(face);
2062 *  
2063 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
2064 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
2065 *  
2066 *   /*--- Loop over dofs. We will set all equal to zero apart from the current one ---*/
2067 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2068 *   for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
2069 *   phi_p.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2070 *   phi_m.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2071 *   }
2072 *   phi_p.submit_dof_value(tmp, i);
2074 *   phi_m.submit_dof_value(tmp, i);
2076 *  
2077 *   /*--- Loop over quadrature points to compute the integral ---*/
2078 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
2079 *   const auto& n_plus = phi_p.get_normal_vector(q);
2080 *   const auto& avg_grad_u_int = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
2081 *   const auto& jump_u_int = phi_p.get_value(q) - phi_m.get_value(q);
2082 *   const auto& avg_tensor_product_u_int = 0.5*(outer_product(phi_p.get_value(q), phi_old_extr_p.get_value(q)) +
2083 *   outer_product(phi_m.get_value(q), phi_old_extr_m.get_value(q)));
2084 *   const auto lambda = std::max(std::abs(scalar_product(phi_old_extr_p.get_value(q), n_plus)),
2085 *   std::abs(scalar_product(phi_old_extr_m.get_value(q), n_plus)));
2086 *  
2087 *   phi_p.submit_value(a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) +
2088 *   a22*avg_tensor_product_u_int*n_plus + 0.5*a22*lambda*jump_u_int , q);
2089 *   phi_m.submit_value(-a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) -
2090 *   a22*avg_tensor_product_u_int*n_plus - 0.5*a22*lambda*jump_u_int, q);
2091 *   phi_p.submit_normal_derivative(-theta_v*0.5*a22/Re*jump_u_int, q);
2092 *   phi_m.submit_normal_derivative(-theta_v*0.5*a22/Re*jump_u_int, q);
2093 *   }
2095 *   diagonal_p[i] = phi_p.get_dof_value(i);
2097 *   diagonal_m[i] = phi_m.get_dof_value(i);
2098 *   }
2099 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2100 *   phi_p.submit_dof_value(diagonal_p[i], i);
2101 *   phi_m.submit_dof_value(diagonal_m[i], i);
2102 *   }
2103 *   phi_p.distribute_local_to_global(dst);
2104 *   phi_m.distribute_local_to_global(dst);
2105 *   }
2106 *   }
2107 *   else {
2109 *   phi_m(data, false, 0),
2110 *   phi_extr_p(data, true, 0),
2111 *   phi_extr_m(data, false, 0);
2112 *  
2113 *   AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component);
2114 *   AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal_p(phi_p.dofs_per_component),
2115 *   diagonal_m(phi_m.dofs_per_component);
2117 *   for(unsigned int d = 0; d < dim; ++d)
2118 *   tmp[d] = make_vectorized_array<Number>(1.0);
2119 *  
2120 *   /*--- Now we loop over faces ---*/
2121 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2122 *   phi_extr_p.reinit(face);
2123 *   phi_extr_p.gather_evaluate(u_extr, EvaluationFlags::values);
2124 *   phi_extr_m.reinit(face);
2125 *   phi_extr_m.gather_evaluate(u_extr, EvaluationFlags::values);
2126 *   phi_p.reinit(face);
2127 *   phi_m.reinit(face);
2128 *  
2129 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
2130 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
2131 *  
2132 *   /*--- Loop over dofs. We will set all equal to zero apart from the current one ---*/
2133 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2134 *   for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
2135 *   phi_p.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2136 *   phi_m.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2137 *   }
2138 *   phi_p.submit_dof_value(tmp, i);
2140 *   phi_m.submit_dof_value(tmp, i);
2142 *  
2143 *   /*--- Loop over quadrature points to compute the integral ---*/
2144 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
2145 *   const auto& n_plus = phi_p.get_normal_vector(q);
2146 *   const auto& avg_grad_u = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
2147 *   const auto& jump_u = phi_p.get_value(q) - phi_m.get_value(q);
2148 *   const auto& avg_tensor_product_u = 0.5*(outer_product(phi_p.get_value(q), phi_extr_p.get_value(q)) +
2149 *   outer_product(phi_m.get_value(q), phi_extr_m.get_value(q)));
2150 *   const auto lambda = std::max(std::abs(scalar_product(phi_extr_p.get_value(q), n_plus)),
2151 *   std::abs(scalar_product(phi_extr_m.get_value(q), n_plus)));
2152 *  
2153 *   phi_p.submit_value(a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) +
2154 *   a33*avg_tensor_product_u*n_plus + 0.5*a33*lambda*jump_u, q);
2155 *   phi_m.submit_value(-a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) -
2156 *   a33*avg_tensor_product_u*n_plus - 0.5*a33*lambda*jump_u, q);
2157 *   phi_p.submit_normal_derivative(-theta_v*0.5*a33/Re*jump_u, q);
2158 *   phi_m.submit_normal_derivative(-theta_v*0.5*a33/Re*jump_u, q);
2159 *   }
2161 *   diagonal_p[i] = phi_p.get_dof_value(i);
2163 *   diagonal_m[i] = phi_m.get_dof_value(i);
2164 *   }
2165 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2166 *   phi_p.submit_dof_value(diagonal_p[i], i);
2167 *   phi_m.submit_dof_value(diagonal_m[i], i);
2168 *   }
2169 *   phi_p.distribute_local_to_global(dst);
2170 *   phi_m.distribute_local_to_global(dst);
2171 *   }
2172 *   }
2173 *   }
2174 *  
2175 *  
2176 * @endcode
2177 *
2178 * The following function assembles boundary term for the velocity
2179 *
2180
2181 *
2182 *
2183 * @code
2184 *   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>
2185 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2186 *   assemble_diagonal_boundary_term_velocity(const MatrixFree<dim, Number>& data,
2187 *   Vec& dst,
2188 *   const unsigned int& ,
2189 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2190 *   if(TR_BDF2_stage == 1) {
2192 *   phi_old_extr(data, true, 0);
2193 *  
2196 *   for(unsigned int d = 0; d < dim; ++d)
2197 *   tmp[d] = make_vectorized_array<Number>(1.0);
2198 *  
2199 *   /*--- Loop over all faces in the range ---*/
2200 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2201 *   phi_old_extr.reinit(face);
2202 *   phi_old_extr.gather_evaluate(u_extr, EvaluationFlags::values);
2203 *   phi.reinit(face);
2204 *  
2205 *   const auto boundary_id = data.get_boundary_id(face);
2206 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
2207 *  
2208 *   if(boundary_id != 1) {
2209 *   const double coef_trasp = 0.0;
2210 *  
2211 *   /*--- Loop over all dofs ---*/
2212 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2213 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2214 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2215 *   phi.submit_dof_value(tmp, i);
2217 *  
2218 *   /*--- Loop over quadrature points to compute the integral ---*/
2219 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2220 *   const auto& n_plus = phi.get_normal_vector(q);
2221 *   const auto& grad_u_int = phi.get_gradient(q);
2222 *   const auto& u_int = phi.get_value(q);
2223 *   const auto& tensor_product_u_int = outer_product(phi.get_value(q), phi_old_extr.get_value(q));
2224 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
2225 *  
2226 *   phi.submit_value(a22/Re*(-grad_u_int*n_plus + 2.0*coef_jump*u_int) +
2227 *   a22*coef_trasp*tensor_product_u_int*n_plus + a22*lambda*u_int, q);
2228 *   phi.submit_normal_derivative(-theta_v*a22/Re*u_int, q);
2229 *   }
2231 *   diagonal[i] = phi.get_dof_value(i);
2232 *   }
2233 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2234 *   phi.submit_dof_value(diagonal[i], i);
2235 *   phi.distribute_local_to_global(dst);
2236 *   }
2237 *   else {
2238 *   /*--- Loop over all dofs ---*/
2239 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2240 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2241 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2242 *   phi.submit_dof_value(tmp, i);
2244 *  
2245 *   /*--- Loop over quadrature points to compute the integral ---*/
2246 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2247 *   const auto& n_plus = phi.get_normal_vector(q);
2248 *   const auto& grad_u_int = phi.get_gradient(q);
2249 *   const auto& u_int = phi.get_value(q);
2250 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
2251 *  
2252 *   const auto& point_vectorized = phi.quadrature_point(q);
2253 *   auto u_int_m = u_int;
2254 *   auto grad_u_int_m = grad_u_int;
2255 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
2256 *   Point<dim> point;
2257 *   for(unsigned int d = 0; d < dim; ++d)
2258 *   point[d] = point_vectorized[d][v];
2259 *  
2260 *   u_int_m[1][v] = -u_int_m[1][v];
2261 *  
2262 *   grad_u_int_m[0][0][v] = -grad_u_int_m[0][0][v];
2263 *   grad_u_int_m[0][1][v] = -grad_u_int_m[0][1][v];
2264 *   }
2265 *  
2266 *   phi.submit_value(a22/Re*(-(0.5*(grad_u_int + grad_u_int_m))*n_plus + coef_jump*(u_int - u_int_m)) +
2267 *   a22*outer_product(0.5*(u_int + u_int_m), phi_old_extr.get_value(q))*n_plus +
2268 *   a22*0.5*lambda*(u_int - u_int_m), q);
2269 *   phi.submit_normal_derivative(-theta_v*a22/Re*(u_int - u_int_m), q);
2270 *   }
2272 *   diagonal[i] = phi.get_dof_value(i);
2273 *   }
2274 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2275 *   phi.submit_dof_value(diagonal[i], i);
2276 *   phi.distribute_local_to_global(dst);
2277 *   }
2278 *   }
2279 *   }
2280 *   else {
2282 *   phi_extr(data, true, 0);
2283 *  
2286 *   for(unsigned int d = 0; d < dim; ++d)
2287 *   tmp[d] = make_vectorized_array<Number>(1.0);
2288 *  
2289 *   /*--- Loop over all faces in the range ---*/
2290 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2291 *   phi_extr.reinit(face);
2292 *   phi_extr.gather_evaluate(u_extr, EvaluationFlags::values);
2293 *   phi.reinit(face);
2294 *  
2295 *   const auto boundary_id = data.get_boundary_id(face);
2296 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
2297 *  
2298 *   if(boundary_id != 1) {
2299 *   const double coef_trasp = 0.0;
2300 *  
2301 *   /*--- Loop over all dofs ---*/
2302 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2303 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2304 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2305 *   phi.submit_dof_value(tmp, i);
2307 *  
2308 *   /*--- Loop over quadrature points to compute the integral ---*/
2309 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2310 *   const auto& n_plus = phi.get_normal_vector(q);
2311 *   const auto& grad_u = phi.get_gradient(q);
2312 *   const auto& u = phi.get_value(q);
2313 *   const auto& tensor_product_u = outer_product(phi.get_value(q), phi_extr.get_value(q));
2314 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
2315 *  
2316 *   phi.submit_value(a33/Re*(-grad_u*n_plus + 2.0*coef_jump*u) +
2317 *   a33*coef_trasp*tensor_product_u*n_plus + a33*lambda*u, q);
2318 *   phi.submit_normal_derivative(-theta_v*a33/Re*u, q);
2319 *   }
2321 *   diagonal[i] = phi.get_dof_value(i);
2322 *   }
2323 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2324 *   phi.submit_dof_value(diagonal[i], i);
2325 *   phi.distribute_local_to_global(dst);
2326 *   }
2327 *   else {
2328 *   /*--- Loop over all dofs ---*/
2329 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2330 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2331 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2332 *   phi.submit_dof_value(tmp, i);
2334 *  
2335 *   /*--- Loop over quadrature points to compute the integral ---*/
2336 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2337 *   const auto& n_plus = phi.get_normal_vector(q);
2338 *   const auto& grad_u = phi.get_gradient(q);
2339 *   const auto& u = phi.get_value(q);
2340 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
2341 *  
2342 *   const auto& point_vectorized = phi.quadrature_point(q);
2343 *   auto u_m = u;
2344 *   auto grad_u_m = grad_u;
2345 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
2346 *   Point<dim> point;
2347 *   for(unsigned int d = 0; d < dim; ++d)
2348 *   point[d] = point_vectorized[d][v];
2349 *  
2350 *   u_m[1][v] = -u_m[1][v];
2351 *  
2352 *   grad_u_m[0][0][v] = -grad_u_m[0][0][v];
2353 *   grad_u_m[0][1][v] = -grad_u_m[0][1][v];
2354 *   }
2355 *  
2356 *   phi.submit_value(a33/Re*(-(0.5*(grad_u + grad_u_m))*n_plus + coef_jump*(u - u_m)) +
2357 *   a33*outer_product(0.5*(u + u_m), phi_extr.get_value(q))*n_plus +
2358 *   a33*0.5*lambda*(u - u_m), q);
2359 *   phi.submit_normal_derivative(-theta_v*a33/Re*(u - u_m), q);
2360 *   }
2362 *   diagonal[i] = phi.get_dof_value(i);
2363 *   }
2364 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2365 *   phi.submit_dof_value(diagonal[i], i);
2366 *   phi.distribute_local_to_global(dst);
2367 *   }
2368 *   }
2369 *   }
2370 *   }
2371 *  
2372 *  
2373 * @endcode
2374 *
2375 * Now we consider the pressure related bilinear forms. We first assemble diagonal cell term for the pressure
2376 *
2377
2378 *
2379 *
2380 * @code
2381 *   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>
2382 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2383 *   assemble_diagonal_cell_term_pressure(const MatrixFree<dim, Number>& data,
2384 *   Vec& dst,
2385 *   const unsigned int& ,
2386 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
2388 *  
2389 *   AlignedVector<VectorizedArray<Number>> diagonal(phi.dofs_per_component); /*--- Here we are using dofs_per_component but
2390 *   it coincides with dofs_per_cell since it is
2391 *   scalar finite element space ---*/
2392 *  
2393 *   const double coeff = (TR_BDF2_stage == 1) ? 1e6*gamma*dt*gamma*dt : 1e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
2394 *  
2395 *   /*--- Loop over all cells in the range ---*/
2396 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
2397 *   phi.reinit(cell);
2398 *  
2399 *   /*--- Loop over all dofs ---*/
2400 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2401 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2402 *   phi.submit_dof_value(VectorizedArray<Number>(), j); /*--- We set all dofs to zero ---*/
2403 *   phi.submit_dof_value(make_vectorized_array<Number>(1.0), i); /*--- Now we set the current one to 1; since it is scalar,
2404 *   we can directly use 'make_vectorized_array' without
2405 *   relying on 'Tensor' ---*/
2407 *  
2408 *   /*--- Loop over quadrature points ---*/
2409 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2410 *   phi.submit_value(1.0/coeff*phi.get_value(q), q);
2411 *   phi.submit_gradient(phi.get_gradient(q), q);
2412 *   }
2414 *   diagonal[i] = phi.get_dof_value(i);
2415 *   }
2416 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2417 *   phi.submit_dof_value(diagonal[i], i);
2418 *  
2419 *   phi.distribute_local_to_global(dst);
2420 *   }
2421 *   }
2422 *  
2423 *  
2424 * @endcode
2425 *
2426 * The following function assembles diagonal face term for the pressure
2427 *
2428
2429 *
2430 *
2431 * @code
2432 *   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>
2433 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2434 *   assemble_diagonal_face_term_pressure(const MatrixFree<dim, Number>& data,
2435 *   Vec& dst,
2436 *   const unsigned int& ,
2437 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2439 *   phi_m(data, false, 1, 1);
2440 *  
2441 *   AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component);
2442 *   AlignedVector<VectorizedArray<Number>> diagonal_p(phi_p.dofs_per_component),
2443 *   diagonal_m(phi_m.dofs_per_component); /*--- Again, we just assert for safety that dimension
2444 *   match, in the sense that we have selected
2445 *   the proper space ---*/
2446 *  
2447 *   /*--- Loop over all faces ---*/
2448 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2449 *   phi_p.reinit(face);
2450 *   phi_m.reinit(face);
2451 *  
2452 *   const auto coef_jump = C_p*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
2453 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
2454 *  
2455 *   /*--- Loop over all dofs ---*/
2456 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2457 *   for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
2458 *   phi_p.submit_dof_value(VectorizedArray<Number>(), j);
2459 *   phi_m.submit_dof_value(VectorizedArray<Number>(), j);
2460 *   }
2461 *   phi_p.submit_dof_value(make_vectorized_array<Number>(1.0), i);
2462 *   phi_m.submit_dof_value(make_vectorized_array<Number>(1.0), i);
2465 *  
2466 *   /*--- Loop over all quadrature points to compute the integral ---*/
2467 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
2468 *   const auto& n_plus = phi_p.get_normal_vector(q);
2469 *  
2470 *   const auto& avg_grad_pres = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
2471 *   const auto& jump_pres = phi_p.get_value(q) - phi_m.get_value(q);
2472 *  
2473 *   phi_p.submit_value(-scalar_product(avg_grad_pres, n_plus) + coef_jump*jump_pres, q);
2474 *   phi_m.submit_value(scalar_product(avg_grad_pres, n_plus) - coef_jump*jump_pres, q);
2475 *   phi_p.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
2476 *   phi_m.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
2477 *   }
2479 *   diagonal_p[i] = phi_p.get_dof_value(i);
2481 *   diagonal_m[i] = phi_m.get_dof_value(i);
2482 *   }
2483 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2484 *   phi_p.submit_dof_value(diagonal_p[i], i);
2485 *   phi_m.submit_dof_value(diagonal_m[i], i);
2486 *   }
2487 *   phi_p.distribute_local_to_global(dst);
2488 *   phi_m.distribute_local_to_global(dst);
2489 *   }
2490 *   }
2491 *  
2492 *  
2493 * @endcode
2494 *
2495 * Eventually, we assemble diagonal boundary term for the pressure
2496 *
2497
2498 *
2499 *
2500 * @code
2501 *   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>
2502 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2503 *   assemble_diagonal_boundary_term_pressure(const MatrixFree<dim, Number>& data,
2504 *   Vec& dst,
2505 *   const unsigned int& ,
2506 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2508 *  
2509 *   AlignedVector<VectorizedArray<Number>> diagonal(phi.dofs_per_component);
2510 *  
2511 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2512 *   const auto boundary_id = data.get_boundary_id(face);
2513 *  
2514 *   if(boundary_id == 1) {
2515 *   phi.reinit(face);
2516 *  
2517 *   const auto coef_jump = C_p*std::abs((phi.get_normal_vector(0)*phi.inverse_jacobian(0))[dim - 1]);
2518 *  
2519 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2520 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2521 *   phi.submit_dof_value(VectorizedArray<Number>(), j);
2522 *   phi.submit_dof_value(make_vectorized_array<Number>(1.0), i);
2524 *  
2525 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2526 *   const auto& n_plus = phi.get_normal_vector(q);
2527 *  
2528 *   const auto& grad_pres = phi.get_gradient(q);
2529 *   const auto& pres = phi.get_value(q);
2530 *  
2531 *   phi.submit_value(-scalar_product(grad_pres, n_plus) + 2.0*coef_jump*pres , q);
2532 *   phi.submit_normal_derivative(-theta_p*pres, q);
2533 *   }
2535 *   diagonal[i] = phi.get_dof_value(i);
2536 *   }
2537 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2538 *   phi.submit_dof_value(diagonal[i], i);
2539 *   phi.distribute_local_to_global(dst);
2540 *   }
2541 *   }
2542 *   }
2543 *  
2544 *  
2545 * @endcode
2546 *
2547 * Put together all previous steps. We create a dummy auxliary vector that serves for the src input argument in
2548 * the previous functions that as we have seen before is unused. Then everything is done by the 'loop' function
2549 * and it is saved in the field 'inverse_diagonal_entries' already present in the base class. Anyway since there is
2550 * only one field, we need to resize properly depending on whether we are considering the velocity or the pressure.
2551 *
2552
2553 *
2554 *
2555 * @code
2556 *   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>
2557 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2558 *   compute_diagonal() {
2559 *   Assert(NS_stage == 1 || NS_stage == 2, ExcInternalError());
2560 *  
2561 *   this->inverse_diagonal_entries.reset(new DiagonalMatrix<Vec>());
2562 *   auto& inverse_diagonal = this->inverse_diagonal_entries->get_vector();
2563 *  
2564 *   if(NS_stage == 1) {
2565 *   ::MatrixFreeTools::compute_diagonal<dim, Number, VectorizedArray<Number>>
2566 *   (*(this->data),
2567 *   inverse_diagonal,
2568 *   [&](const auto& data, auto& dst, const auto& src, const auto& cell_range) {
2569 *   (this->assemble_diagonal_cell_term_velocity)(data, dst, src, cell_range);
2570 *   },
2571 *   [&](const auto& data, auto& dst, const auto& src, const auto& face_range) {
2572 *   (this->assemble_diagonal_face_term_velocity)(data, dst, src, face_range);
2573 *   },
2574 *   [&](const auto& data, auto& dst, const auto& src, const auto& boundary_range) {
2575 *   (this->assemble_diagonal_boundary_term_velocity)(data, dst, src, boundary_range);
2576 *   },
2577 *   0);
2578 *   }
2579 *   else if(NS_stage == 2) {
2580 *   ::MatrixFreeTools::compute_diagonal<dim, Number, VectorizedArray<Number>>
2581 *   (*(this->data),
2582 *   inverse_diagonal,
2583 *   [&](const auto& data, auto& dst, const auto& src, const auto& cell_range) {
2584 *   (this->assemble_diagonal_cell_term_pressure)(data, dst, src, cell_range);
2585 *   },
2586 *   [&](const auto& data, auto& dst, const auto& src, const auto& face_range) {
2587 *   (this->assemble_diagonal_face_term_pressure)(data, dst, src, face_range);
2588 *   },
2589 *   [&](const auto& data, auto& dst, const auto& src, const auto& boundary_range) {
2590 *   (this->assemble_diagonal_boundary_term_pressure)(data, dst, src, boundary_range);
2591 *   },
2592 *   1);
2593 *   }
2594 *  
2595 *   for(unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i) {
2596 *   Assert(inverse_diagonal.local_element(i) != 0.0,
2597 *   ExcMessage("No diagonal entry in a definite operator should be zero"));
2598 *   inverse_diagonal.local_element(i) = 1.0/inverse_diagonal.local_element(i);
2599 *   }
2600 *   }
2601 *  
2602 *  
2603 * @endcode
2604 *
2605 *
2606 * <a name="navier_stokes_TRBDF2_DG.cc-"></a>
2607 * @sect{The <code>NavierStokesProjection</code> class}
2608 *
2609
2610 *
2611 * Now we are ready for the main class of the program. It implements the calls to the various steps
2612 * of the projection method for Navier-Stokes equations.
2613 *
2614
2615 *
2616 *
2617 * @code
2618 *   template<int dim>
2619 *   class NavierStokesProjection {
2620 *   public:
2621 *   NavierStokesProjection(RunTimeParameters::Data_Storage& data);
2622 *  
2623 *   void run(const bool verbose = false, const unsigned int output_interval = 10);
2624 *  
2625 *   protected:
2626 *   const double t_0;
2627 *   const double T;
2628 *   const double gamma; //--- TR-BDF2 parameter
2629 *   unsigned int TR_BDF2_stage; //--- Flag to check at which current stage of TR-BDF2 are
2630 *   const double Re;
2631 *   double dt;
2632 *  
2633 *   EquationData::Velocity<dim> vel_init;
2634 *   EquationData::Pressure<dim> pres_init; /*--- Instance of 'Velocity' and 'Pressure' classes to initialize. ---*/
2635 *  
2637 *  
2638 *   /*--- Finite Element spaces ---*/
2639 *   FESystem<dim> fe_velocity;
2640 *   FESystem<dim> fe_pressure;
2641 *  
2642 *   /*--- Handler for dofs ---*/
2643 *   DoFHandler<dim> dof_handler_velocity;
2644 *   DoFHandler<dim> dof_handler_pressure;
2645 *  
2646 *   /*--- Quadrature formulas for velocity and pressure, respectively ---*/
2647 *   QGauss<dim> quadrature_pressure;
2648 *   QGauss<dim> quadrature_velocity;
2649 *  
2650 *   /*--- Now we define all the vectors for the solution. We start from the pressure
2651 *   with p^n, p^(n+gamma) and a vector for rhs ---*/
2655 *  
2656 *   /*--- Next, we move to the velocity, with u^n, u^(n-1), u^(n+gamma/2),
2657 *   u^(n+gamma) and other two auxiliary vectors as well as the rhs ---*/
2666 *  
2667 *   Vector<double> Linfty_error_per_cell_vel;
2668 *  
2669 *   DeclException2(ExcInvalidTimeStep,
2670 *   double,
2671 *   double,
2672 *   << " The time step " << arg1 << " is out of range."
2673 *   << std::endl
2674 *   << " The permitted range is (0," << arg2 << "]");
2675 *  
2676 *   void create_triangulation(const unsigned int n_refines);
2677 *  
2678 *   void setup_dofs();
2679 *  
2680 *   void initialize();
2681 *  
2682 *   void interpolate_velocity();
2683 *  
2684 *   void diffusion_step();
2685 *  
2686 *   void projection_step();
2687 *  
2688 *   void project_grad(const unsigned int flag);
2689 *  
2690 *   double get_maximal_velocity();
2691 *  
2692 *   double get_maximal_difference_velocity();
2693 *  
2694 *   void output_results(const unsigned int step);
2695 *  
2696 *   void refine_mesh();
2697 *  
2698 *   void interpolate_max_res(const unsigned int level);
2699 *  
2700 *   void save_max_res();
2701 *  
2702 *   private:
2703 *   void compute_lift_and_drag();
2704 *  
2705 *   /*--- Technical member to handle the various steps ---*/
2706 *   std::shared_ptr<MatrixFree<dim, double>> matrix_free_storage;
2707 *  
2708 *   /*--- Now we need an instance of the class implemented before with the weak form ---*/
2709 *   NavierStokesProjectionOperator<dim, EquationData::degree_p, EquationData::degree_p + 1,
2710 *   EquationData::degree_p + 1, EquationData::degree_p + 2,
2711 *   LinearAlgebra::distributed::Vector<double>> navier_stokes_matrix;
2712 *  
2713 *   /*--- This is an instance for geometric multigrid preconditioner ---*/
2714 *   MGLevelObject<NavierStokesProjectionOperator<dim, EquationData::degree_p, EquationData::degree_p + 1,
2715 *   EquationData::degree_p + 1, EquationData::degree_p + 2,
2717 *  
2718 *   /*--- Here we define two 'AffineConstraints' instance, one for each finite element space.
2719 *   This is just a technical issue, due to MatrixFree requirements. In general
2720 *   this class is used to impose boundary conditions (or any kind of constraints), but in this case, since
2721 *   we are using a weak imposition of bcs, everything is already in the weak forms and so these instances
2722 *   will be default constructed ---*/
2723 *   AffineConstraints<double> constraints_velocity,
2724 *   constraints_pressure;
2725 *  
2726 *   /*--- Now a bunch of variables handled by 'ParamHandler' introduced at the beginning of the code ---*/
2727 *   unsigned int max_its;
2728 *   double eps;
2729 *  
2730 *   unsigned int max_loc_refinements;
2731 *   unsigned int min_loc_refinements;
2732 *   unsigned int refinement_iterations;
2733 *  
2734 *   std::string saving_dir;
2735 *  
2736 *   /*--- Finally, some output related streams ---*/
2737 *   ConditionalOStream pcout;
2738 *  
2739 *   std::ofstream time_out;
2740 *   ConditionalOStream ptime_out;
2741 *   TimerOutput time_table;
2742 *  
2743 *   std::ofstream output_n_dofs_velocity;
2744 *   std::ofstream output_n_dofs_pressure;
2745 *  
2746 *   std::ofstream output_lift;
2747 *   std::ofstream output_drag;
2748 *   };
2749 *  
2750 *  
2751 * @endcode
2752 *
2753 * In the constructor, we just read all the data from the
2754 * <code>Data_Storage</code> object that is passed as an argument, verify that
2755 * the data we read are reasonable and, finally, create the triangulation and
2756 * load the initial data.
2757 *
2758
2759 *
2760 *
2761 * @code
2762 *   template<int dim>
2763 *   NavierStokesProjection<dim>::NavierStokesProjection(RunTimeParameters::Data_Storage& data):
2764 *   t_0(data.initial_time),
2765 *   T(data.final_time),
2766 *   gamma(2.0 - std::sqrt(2.0)), //--- Save also in the NavierStokes class the TR-BDF2 parameter value
2767 *   TR_BDF2_stage(1), //--- Initialize the flag for the TR_BDF2 stage
2768 *   Re(data.Reynolds),
2769 *   dt(data.dt),
2770 *   vel_init(data.initial_time),
2771 *   pres_init(data.initial_time),
2772 *   triangulation(MPI_COMM_WORLD, parallel::distributed::Triangulation<dim>::limit_level_difference_at_vertices,
2773 *   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
2774 *   fe_velocity(FE_DGQ<dim>(EquationData::degree_p + 1), dim),
2775 *   fe_pressure(FE_DGQ<dim>(EquationData::degree_p), 1),
2776 *   dof_handler_velocity(triangulation),
2777 *   dof_handler_pressure(triangulation),
2778 *   quadrature_pressure(EquationData::degree_p + 1),
2779 *   quadrature_velocity(EquationData::degree_p + 2),
2780 *   navier_stokes_matrix(data),
2781 *   max_its(data.max_iterations),
2782 *   eps(data.eps),
2783 *   max_loc_refinements(data.max_loc_refinements),
2784 *   min_loc_refinements(data.min_loc_refinements),
2785 *   refinement_iterations(data.refinement_iterations),
2786 *   saving_dir(data.dir),
2787 *   pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
2788 *   time_out("./" + data.dir + "/time_analysis_" +
2789 *   Utilities::int_to_string(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)) + "proc.dat"),
2790 *   ptime_out(time_out, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
2791 *   time_table(ptime_out, TimerOutput::summary, TimerOutput::cpu_and_wall_times),
2792 *   output_n_dofs_velocity("./" + data.dir + "/n_dofs_velocity.dat", std::ofstream::out),
2793 *   output_n_dofs_pressure("./" + data.dir + "/n_dofs_pressure.dat", std::ofstream::out),
2794 *   output_lift("./" + data.dir + "/lift.dat", std::ofstream::out),
2795 *   output_drag("./" + data.dir + "/drag.dat", std::ofstream::out) {
2796 *   if(EquationData::degree_p < 1) {
2797 *   pcout
2798 *   << " WARNING: The chosen pair of finite element spaces is not stable."
2799 *   << std::endl
2800 *   << " The obtained results will be nonsense" << std::endl;
2801 *   }
2802 *  
2803 *   AssertThrow(!((dt <= 0.0) || (dt > 0.5*T)), ExcInvalidTimeStep(dt, 0.5*T));
2804 *  
2805 *   matrix_free_storage = std::make_shared<MatrixFree<dim, double>>();
2806 *  
2807 *   create_triangulation(data.n_refines);
2808 *   setup_dofs();
2809 *   initialize();
2810 *   }
2811 *  
2812 *  
2813 * @endcode
2814 *
2815 * The method that creates the triangulation and refines it the needed number
2816 * of times.
2817 *
2818
2819 *
2820 *
2821 * @code
2822 *   template<int dim>
2823 *   void NavierStokesProjection<dim>::create_triangulation(const unsigned int n_refines) {
2824 *   TimerOutput::Scope t(time_table, "Create triangulation");
2825 *  
2826 *   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);
2827 *   /*--- We strongly advice to check the documentation to verify the meaning of all input parameters. ---*/
2828 *  
2829 *   pcout << "Number of refines = " << n_refines << std::endl;
2830 *   triangulation.refine_global(n_refines);
2831 *   }
2832 *  
2833 *  
2834 * @endcode
2835 *
2836 * After creating the triangulation, it creates the mesh dependent
2837 * data, i.e. it distributes degrees of freedom, and
2838 * initializes the vectors that we will use.
2839 *
2840
2841 *
2842 *
2843 * @code
2844 *   template<int dim>
2845 *   void NavierStokesProjection<dim>::setup_dofs() {
2846 *   pcout << "Number of active cells: " << triangulation.n_global_active_cells() << std::endl;
2847 *   pcout << "Number of levels: " << triangulation.n_global_levels() << std::endl;
2848 *  
2849 *   /*--- Distribute dofs and prepare for multigrid ---*/
2850 *   dof_handler_velocity.distribute_dofs(fe_velocity);
2851 *   dof_handler_pressure.distribute_dofs(fe_pressure);
2852 *  
2853 *   pcout << "dim (X_h) = " << dof_handler_velocity.n_dofs()
2854 *   << std::endl
2855 *   << "dim (M_h) = " << dof_handler_pressure.n_dofs()
2856 *   << std::endl
2857 *   << "Re = " << Re << std::endl
2858 *   << std::endl;
2859 *  
2860 *   if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
2861 *   output_n_dofs_velocity << dof_handler_velocity.n_dofs() << std::endl;
2862 *   output_n_dofs_pressure << dof_handler_pressure.n_dofs() << std::endl;
2863 *   }
2864 *  
2865 *   typename MatrixFree<dim, double>::AdditionalData additional_data;
2866 *   additional_data.mapping_update_flags = (update_gradients | update_JxW_values |
2868 *   additional_data.mapping_update_flags_inner_faces = (update_gradients | update_JxW_values | update_quadrature_points |
2870 *   additional_data.mapping_update_flags_boundary_faces = (update_gradients | update_JxW_values | update_quadrature_points |
2872 *   additional_data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
2873 *  
2874 *   std::vector<const DoFHandler<dim>*> dof_handlers; /*--- Vector of dof_handlers to feed the 'MatrixFree'. Here the order
2875 *   counts and enters into the game as parameter of FEEvaluation and
2876 *   FEFaceEvaluation in the previous class ---*/
2877 *   dof_handlers.push_back(&dof_handler_velocity);
2878 *   dof_handlers.push_back(&dof_handler_pressure);
2879 *  
2880 *   constraints_velocity.clear();
2881 *   constraints_velocity.close();
2882 *   constraints_pressure.clear();
2883 *   constraints_pressure.close();
2884 *   std::vector<const AffineConstraints<double>*> constraints;
2885 *   constraints.push_back(&constraints_velocity);
2886 *   constraints.push_back(&constraints_pressure);
2887 *  
2888 *   std::vector<QGauss<1>> quadratures; /*--- We cannot directly use 'quadrature_velocity' and 'quadrature_pressure',
2889 *   because the 'MatrixFree' structure wants a quadrature formula for 1D
2890 *   (this is way the template parameter of the previous class was called 'n_q_points_1d_p'
2891 *   and 'n_q_points_1d_v' and the reason of '1' as QGauss template parameter). ---*/
2892 *   quadratures.push_back(QGauss<1>(EquationData::degree_p + 2));
2893 *   quadratures.push_back(QGauss<1>(EquationData::degree_p + 1));
2894 *  
2895 *   /*--- Initialize the matrix-free structure and size properly the vectors. Here again the
2896 *   second input argument of the 'initialize_dof_vector' method depends on the order of 'dof_handlers' ---*/
2897 *   matrix_free_storage->reinit(MappingQ1<dim>(),dof_handlers, constraints, quadratures, additional_data);
2898 *   matrix_free_storage->initialize_dof_vector(u_star, 0);
2899 *   matrix_free_storage->initialize_dof_vector(rhs_u, 0);
2900 *   matrix_free_storage->initialize_dof_vector(u_n, 0);
2901 *   matrix_free_storage->initialize_dof_vector(u_extr, 0);
2902 *   matrix_free_storage->initialize_dof_vector(u_n_minus_1, 0);
2903 *   matrix_free_storage->initialize_dof_vector(u_n_gamma, 0);
2904 *   matrix_free_storage->initialize_dof_vector(u_tmp, 0);
2905 *   matrix_free_storage->initialize_dof_vector(grad_pres_int, 0);
2906 *  
2907 *   matrix_free_storage->initialize_dof_vector(pres_int, 1);
2908 *   matrix_free_storage->initialize_dof_vector(pres_n, 1);
2909 *   matrix_free_storage->initialize_dof_vector(rhs_p, 1);
2910 *  
2911 *   /*--- Initialize the multigrid structure. We dedicate ad hoc 'dof_handlers_mg' and 'constraints_mg' because
2912 *   we use float as type. Moreover we can initialize already with the index of the finite element of the pressure;
2913 *   anyway we need by requirement to declare also structures for the velocity for coherence (basically because
2914 *   the index of finite element space has to be the same, so the pressure has to be the second).---*/
2915 *   mg_matrices.clear_elements();
2916 *   dof_handler_velocity.distribute_mg_dofs();
2917 *   dof_handler_pressure.distribute_mg_dofs();
2918 *  
2919 *   const unsigned int nlevels = triangulation.n_global_levels();
2920 *   mg_matrices.resize(0, nlevels - 1);
2921 *   for(unsigned int level = 0; level < nlevels; ++level) {
2922 *   typename MatrixFree<dim, float>::AdditionalData additional_data_mg;
2924 *   additional_data_mg.mapping_update_flags = (update_gradients | update_JxW_values);
2925 *   additional_data_mg.mapping_update_flags_inner_faces = (update_gradients | update_JxW_values);
2926 *   additional_data_mg.mapping_update_flags_boundary_faces = (update_gradients | update_JxW_values);
2927 *   additional_data_mg.mg_level = level;
2928 *  
2929 *   std::vector<const DoFHandler<dim>*> dof_handlers_mg;
2930 *   dof_handlers_mg.push_back(&dof_handler_velocity);
2931 *   dof_handlers_mg.push_back(&dof_handler_pressure);
2932 *   std::vector<const AffineConstraints<float>*> constraints_mg;
2933 *   AffineConstraints<float> constraints_velocity_mg;
2934 *   constraints_velocity_mg.clear();
2935 *   constraints_velocity_mg.close();
2936 *   constraints_mg.push_back(&constraints_velocity_mg);
2937 *   AffineConstraints<float> constraints_pressure_mg;
2938 *   constraints_pressure_mg.clear();
2939 *   constraints_pressure_mg.close();
2940 *   constraints_mg.push_back(&constraints_pressure_mg);
2941 *  
2942 *   std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(new MatrixFree<dim, float>());
2943 *   mg_mf_storage_level->reinit(MappingQ1<dim>(),dof_handlers_mg, constraints_mg, quadratures, additional_data_mg);
2944 *   const std::vector<unsigned int> tmp = {1};
2945 *   mg_matrices[level].initialize(mg_mf_storage_level, tmp, tmp);
2946 *   mg_matrices[level].set_dt(dt);
2947 *   mg_matrices[level].set_NS_stage(2);
2948 *   }
2949 *  
2950 *   Linfty_error_per_cell_vel.reinit(triangulation.n_active_cells());
2951 *   }
2952 *  
2953 *  
2954 * @endcode
2955 *
2956 * This method loads the initial data. It simply uses the class <code>Pressure</code> instance for the pressure
2957 * and the class <code>Velocity</code> instance for the velocity.
2958 *
2959
2960 *
2961 *
2962 * @code
2963 *   template<int dim>
2964 *   void NavierStokesProjection<dim>::initialize() {
2965 *   TimerOutput::Scope t(time_table, "Initialize pressure and velocity");
2966 *  
2967 *   VectorTools::interpolate(dof_handler_pressure, pres_init, pres_n);
2968 *  
2969 *   VectorTools::interpolate(dof_handler_velocity, vel_init, u_n_minus_1);
2970 *   VectorTools::interpolate(dof_handler_velocity, vel_init, u_n);
2971 *   }
2972 *  
2973 *  
2974 * @endcode
2975 *
2976 * This function computes the extrapolated velocity to be used in the momentum predictor
2977 *
2978
2979 *
2980 *
2981 * @code
2982 *   template<int dim>
2983 *   void NavierStokesProjection<dim>::interpolate_velocity() {
2984 *   TimerOutput::Scope t(time_table, "Interpolate velocity");
2985 *  
2986 * @endcode
2987 *
2988 * --- TR-BDF2 first step
2989 *
2990 * @code
2991 *   if(TR_BDF2_stage == 1) {
2992 *   u_extr.equ(1.0 + gamma/(2.0*(1.0 - gamma)), u_n);
2993 *   u_tmp.equ(gamma/(2.0*(1.0 - gamma)), u_n_minus_1);
2994 *   u_extr -= u_tmp;
2995 *   }
2996 * @endcode
2997 *
2998 * --- TR-BDF2 second step
2999 *
3000 * @code
3001 *   else {
3002 *   u_extr.equ(1.0 + (1.0 - gamma)/gamma, u_n_gamma);
3003 *   u_tmp.equ((1.0 - gamma)/gamma, u_n);
3004 *   u_extr -= u_tmp;
3005 *   }
3006 *   }
3007 *  
3008 *  
3009 * @endcode
3010 *
3011 * We are finally ready to solve the diffusion step.
3012 *
3013
3014 *
3015 *
3016 * @code
3017 *   template<int dim>
3018 *   void NavierStokesProjection<dim>::diffusion_step() {
3019 *   TimerOutput::Scope t(time_table, "Diffusion step");
3020 *  
3021 *   /*--- We first specify that we want to deal with velocity dof_handler (index 0, since it is the first one
3022 *   in the 'dof_handlers' vector) ---*/
3023 *   const std::vector<unsigned int> tmp = {0};
3024 *   navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
3025 *  
3026 *   /*--- Next, we specify at we are at stage 1, namely the diffusion step ---*/
3027 *   navier_stokes_matrix.set_NS_stage(1);
3028 *  
3029 *   /*--- Now, we compute the right-hand side and we set the convective velocity. The necessity of 'set_u_extr' is
3030 *   that this quantity is required in the bilinear forms and we can't use a vector of src like on the right-hand side,
3031 *   so it has to be available ---*/
3032 *   if(TR_BDF2_stage == 1) {
3033 *   navier_stokes_matrix.vmult_rhs_velocity(rhs_u, {u_n, u_extr, pres_n});
3034 *   navier_stokes_matrix.set_u_extr(u_extr);
3035 *   u_star = u_extr;
3036 *   }
3037 *   else {
3038 *   navier_stokes_matrix.vmult_rhs_velocity(rhs_u, {u_n, u_n_gamma, pres_int, u_extr});
3039 *   navier_stokes_matrix.set_u_extr(u_extr);
3040 *   u_star = u_extr;
3041 *   }
3042 *  
3043 *   /*--- Build the linear solver; in this case we specify the maximum number of iterations and residual ---*/
3044 *   SolverControl solver_control(max_its, eps*rhs_u.l2_norm());
3046 *  
3047 *   /*--- Build a Jacobi preconditioner and solve ---*/
3048 *   PreconditionJacobi<NavierStokesProjectionOperator<dim,
3049 *   EquationData::degree_p,
3050 *   EquationData::degree_p + 1,
3051 *   EquationData::degree_p + 1,
3052 *   EquationData::degree_p + 2,
3053 *   LinearAlgebra::distributed::Vector<double>>> preconditioner;
3054 *   navier_stokes_matrix.compute_diagonal();
3055 *   preconditioner.initialize(navier_stokes_matrix);
3056 *  
3057 *   gmres.solve(navier_stokes_matrix, u_star, rhs_u, preconditioner);
3058 *   }
3059 *  
3060 *  
3061 * @endcode
3062 *
3063 * Next, we solve the projection step.
3064 *
3065
3066 *
3067 *
3068 * @code
3069 *   template<int dim>
3070 *   void NavierStokesProjection<dim>::projection_step() {
3071 *   TimerOutput::Scope t(time_table, "Projection step pressure");
3072 *  
3073 *   /*--- We start in the same way of 'diffusion_step': we first reinitialize with the index of FE space,
3074 *   we specify that this is the second stage and we compute the right-hand side ---*/
3075 *   const std::vector<unsigned int> tmp = {1};
3076 *   navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
3077 *  
3078 *   navier_stokes_matrix.set_NS_stage(2);
3079 *  
3080 *   if(TR_BDF2_stage == 1)
3081 *   navier_stokes_matrix.vmult_rhs_pressure(rhs_p, {u_star, pres_n});
3082 *   else
3083 *   navier_stokes_matrix.vmult_rhs_pressure(rhs_p, {u_star, pres_int});
3084 *  
3085 *   /*--- Build the linear solver (Conjugate Gradient in this case) ---*/
3086 *   SolverControl solver_control(max_its, eps*rhs_p.l2_norm());
3088 *  
3089 *   /*--- Build the preconditioner (as in @ref step_37 "step-37") ---*/
3090 *   MGTransferMatrixFree<dim, float> mg_transfer;
3091 *   mg_transfer.build(dof_handler_pressure);
3092 *  
3093 *   using SmootherType = PreconditionChebyshev<NavierStokesProjectionOperator<dim,
3094 *   EquationData::degree_p,
3095 *   EquationData::degree_p + 1,
3096 *   EquationData::degree_p + 1,
3097 *   EquationData::degree_p + 2,
3102 *   smoother_data.resize(0, triangulation.n_global_levels() - 1);
3103 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level) {
3104 *   if(level > 0) {
3105 *   smoother_data[level].smoothing_range = 15.0;
3106 *   smoother_data[level].degree = 3;
3107 *   smoother_data[level].eig_cg_n_iterations = 10;
3108 *   }
3109 *   else {
3110 *   smoother_data[0].smoothing_range = 2e-2;
3111 *   smoother_data[0].degree = numbers::invalid_unsigned_int;
3112 *   smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
3113 *   }
3114 *   mg_matrices[level].compute_diagonal();
3115 *   smoother_data[level].preconditioner = mg_matrices[level].get_matrix_diagonal_inverse();
3116 *   }
3117 *   mg_smoother.initialize(mg_matrices, smoother_data);
3118 *  
3120 *   SolverCG<LinearAlgebra::distributed::Vector<float>> cg_mg(solver_control);
3123 *   NavierStokesProjectionOperator<dim,
3124 *   EquationData::degree_p,
3125 *   EquationData::degree_p + 1,
3126 *   EquationData::degree_p + 1,
3127 *   EquationData::degree_p + 2,
3129 *   PreconditionIdentity> mg_coarse(cg_mg, mg_matrices[0], identity);
3130 *  
3131 *   mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(mg_matrices);
3132 *  
3133 *   Multigrid<LinearAlgebra::distributed::Vector<float>> mg(mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
3134 *  
3135 *   PreconditionMG<dim,
3137 *   MGTransferMatrixFree<dim, float>> preconditioner(dof_handler_pressure, mg, mg_transfer);
3138 *  
3139 *   /*--- Solve the linear system ---*/
3140 *   if(TR_BDF2_stage == 1) {
3141 *   pres_int = pres_n;
3142 *   cg.solve(navier_stokes_matrix, pres_int, rhs_p, preconditioner);
3143 *   }
3144 *   else {
3145 *   pres_n = pres_int;
3146 *   cg.solve(navier_stokes_matrix, pres_n, rhs_p, preconditioner);
3147 *   }
3148 *   }
3149 *  
3150 *  
3151 * @endcode
3152 *
3153 * This implements the projection step for the gradient of pressure
3154 *
3155
3156 *
3157 *
3158 * @code
3159 *   template<int dim>
3160 *   void NavierStokesProjection<dim>::project_grad(const unsigned int flag) {
3161 *   TimerOutput::Scope t(time_table, "Gradient of pressure projection");
3162 *  
3163 *   /*--- The input parameter flag is used just to specify where we want to save the result ---*/
3164 *   AssertIndexRange(flag, 3);
3165 *   Assert(flag > 0, ExcInternalError());
3166 *  
3167 *   /*--- We need to select the dof handler related to the velocity since the result lives there ---*/
3168 *   const std::vector<unsigned int> tmp = {0};
3169 *   navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
3170 *  
3171 *   if(flag == 1)
3172 *   navier_stokes_matrix.vmult_grad_p_projection(rhs_u, pres_n);
3173 *   else if(flag == 2)
3174 *   navier_stokes_matrix.vmult_grad_p_projection(rhs_u, pres_int);
3175 *  
3176 *   /*--- We conventionally decide that the this corresponds to third stage ---*/
3177 *   navier_stokes_matrix.set_NS_stage(3);
3178 *  
3179 *   /*--- Solve the system ---*/
3180 *   SolverControl solver_control(max_its, 1e-12*rhs_u.l2_norm());
3182 *   cg.solve(navier_stokes_matrix, u_tmp, rhs_u, PreconditionIdentity());
3183 *   }
3184 *  
3185 *  
3186 * @endcode
3187 *
3188 * The following function is used in determining the maximal velocity
3189 * in order to compute the Courant number.
3190 *
3191
3192 *
3193 *
3194 * @code
3195 *   template<int dim>
3196 *   double NavierStokesProjection<dim>::get_maximal_velocity() {
3197 *   return u_n.linfty_norm();
3198 *   }
3199 *  
3200 *  
3201 * @endcode
3202 *
3203 * The following function is used in determining the maximal nodal difference
3204 * between old and current velocity value in order to see if we have reached steady-state.
3205 *
3206
3207 *
3208 *
3209 * @code
3210 *   template<int dim>
3211 *   double NavierStokesProjection<dim>::get_maximal_difference_velocity() {
3212 *   u_tmp = u_n;
3213 *   u_tmp -= u_n_minus_1;
3214 *  
3215 *   return u_tmp.linfty_norm();
3216 *   }
3217 *  
3218 *  
3219 * @endcode
3220 *
3221 * This method plots the current solution. The main difficulty is that we want
3222 * to create a single output file that contains the data for all velocity
3223 * components and the pressure. On the other hand, velocities and the pressure
3224 * live on separate DoFHandler objects, so we need to pay attention when we use
3225 * 'add_data_vector' to select the proper space.
3226 *
3227
3228 *
3229 *
3230 * @code
3231 *   template<int dim>
3232 *   void NavierStokesProjection<dim>::output_results(const unsigned int step) {
3233 *   TimerOutput::Scope t(time_table, "Output results");
3234 *  
3235 *   DataOut<dim> data_out;
3236 *  
3237 *   std::vector<std::string> velocity_names(dim, "v");
3238 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3239 *   component_interpretation_velocity(dim, DataComponentInterpretation::component_is_part_of_vector);
3240 *   u_n.update_ghost_values();
3241 *   data_out.add_data_vector(dof_handler_velocity, u_n, velocity_names, component_interpretation_velocity);
3242 *   pres_n.update_ghost_values();
3243 *   data_out.add_data_vector(dof_handler_pressure, pres_n, "p", {DataComponentInterpretation::component_is_scalar});
3244 *  
3245 *   std::vector<std::string> velocity_names_old(dim, "v_old");
3246 *   u_n_minus_1.update_ghost_values();
3247 *   data_out.add_data_vector(dof_handler_velocity, u_n_minus_1, velocity_names_old, component_interpretation_velocity);
3248 *  
3249 *   /*--- Here we rely on the postprocessor we have built ---*/
3250 *   PostprocessorVorticity<dim> postprocessor;
3251 *   data_out.add_data_vector(dof_handler_velocity, u_n, postprocessor);
3252 *  
3253 *   data_out.build_patches(MappingQ1<dim>(), 1, DataOut<dim>::curved_inner_cells);
3254 *  
3255 *   const std::string output = "./" + saving_dir + "/solution-" + Utilities::int_to_string(step, 5) + ".vtu";
3256 *   data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
3257 *   }
3258 *  
3259 *  
3260 * @endcode
3261 *
3262 *
3263 * <a name="navier_stokes_TRBDF2_DG.cc-"></a>
3264 * @sect{<code>NavierStokesProjection::compute_lift_and_drag</code>}
3265 *
3266
3267 *
3268 * This routine computes the lift and the drag forces in a non-dimensional framework
3269 * (so basically for the classical coefficients, it is necessary to multiply by a factor 2).
3270 *
3271
3272 *
3273 *
3274 * @code
3275 *   template<int dim>
3276 *   void NavierStokesProjection<dim>::compute_lift_and_drag() {
3277 *   QGauss<dim - 1> face_quadrature_formula(EquationData::degree_p + 2);
3278 *   const int n_q_points = face_quadrature_formula.size();
3279 *  
3280 *   std::vector<double> pressure_values(n_q_points);
3281 *   std::vector<std::vector<Tensor<1, dim>>> velocity_gradients(n_q_points, std::vector<Tensor<1, dim>>(dim));
3282 *  
3283 *   Tensor<1, dim> normal_vector;
3284 *   Tensor<2, dim> fluid_stress;
3285 *   Tensor<2, dim> fluid_pressure;
3286 *   Tensor<1, dim> forces;
3287 *  
3288 *   /*--- We need to compute the integral over the cylinder boundary, so we need to use 'FEFaceValues' instances.
3289 *   For the velocity we need the gradients, for the pressure the values. ---*/
3290 *   FEFaceValues<dim> fe_face_values_velocity(fe_velocity, face_quadrature_formula,
3293 *   FEFaceValues<dim> fe_face_values_pressure(fe_pressure, face_quadrature_formula, update_values);
3294 *  
3295 *   double local_drag = 0.0;
3296 *   double local_lift = 0.0;
3297 *  
3298 *   /*--- We need to perform a unique loop because the whole stress tensor takes into account contributions of
3299 *   velocity and pressure obviously. However, the two dof_handlers are different, so we need to create an ad-hoc
3300 *   iterator for the pressure that we update manually. It is guaranteed that the cells are visited in the same order
3301 *   (see the documentation) ---*/
3302 *   auto tmp_cell = dof_handler_pressure.begin_active();
3303 *   for(const auto& cell : dof_handler_velocity.active_cell_iterators()) {
3304 *   if(cell->is_locally_owned()) {
3305 *   for(unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
3306 *   if(cell->face(face)->at_boundary() && cell->face(face)->boundary_id() == 4) {
3307 *   fe_face_values_velocity.reinit(cell, face);
3308 *   fe_face_values_pressure.reinit(tmp_cell, face);
3309 *  
3310 *   fe_face_values_velocity.get_function_gradients(u_n, velocity_gradients); /*--- velocity gradients ---*/
3311 *   fe_face_values_pressure.get_function_values(pres_n, pressure_values); /*--- pressure values ---*/
3312 *  
3313 *   for(int q = 0; q < n_q_points; q++) {
3314 *   normal_vector = -fe_face_values_velocity.normal_vector(q);
3315 *  
3316 *   for(unsigned int d = 0; d < dim; ++ d) {
3317 *   fluid_pressure[d][d] = pressure_values[q];
3318 *   for(unsigned int k = 0; k < dim; ++k)
3319 *   fluid_stress[d][k] = 1.0/Re*velocity_gradients[q][d][k];
3320 *   }
3321 *   fluid_stress = fluid_stress - fluid_pressure;
3322 *  
3323 *   forces = fluid_stress*normal_vector*fe_face_values_velocity.JxW(q);
3324 *  
3325 *   local_drag += forces[0];
3326 *   local_lift += forces[1];
3327 *   }
3328 *   }
3329 *   }
3330 *   }
3331 *   ++tmp_cell;
3332 *   }
3333 *  
3334 *   /*--- At the end, each processor has computed the contribution to the boundary cells it owns and, therefore,
3335 *   we need to sum up all the contributions. ---*/
3336 *   const double lift = Utilities::MPI::sum(local_lift, MPI_COMM_WORLD);
3337 *   const double drag = Utilities::MPI::sum(local_drag, MPI_COMM_WORLD);
3338 *   if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
3339 *   output_lift << lift << std::endl;
3340 *   output_drag << drag << std::endl;
3341 *   }
3342 *   }
3343 *  
3344 *  
3345 * @endcode
3346 *
3347 *
3348 * <a name="navier_stokes_TRBDF2_DG.cc-"></a>
3349 * @sect{ <code>NavierStokesProjection::refine_mesh</code>}
3350 *
3351
3352 *
3353 * After finding a good initial guess on the coarse mesh, we hope to
3354 * decrease the error through refining the mesh. We also need to transfer the current solution to the
3355 * next mesh using the SolutionTransfer class.
3356 *
3357
3358 *
3359 *
3360 * @code
3361 *   template <int dim>
3362 *   void NavierStokesProjection<dim>::refine_mesh() {
3363 *   TimerOutput::Scope t(time_table, "Refine mesh");
3364 *  
3365 *   /*--- We first create a proper vector for computing estimator ---*/
3366 *   const IndexSet locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler_velocity);
3368 *   tmp_velocity.reinit(dof_handler_velocity.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
3369 *   tmp_velocity = u_n;
3370 *   tmp_velocity.update_ghost_values();
3371 *  
3372 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
3373 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
3374 *  
3375 *   /*--- This is basically the indicator per cell computation (see @ref step_50 "step-50"). Since it is not so complciated
3376 *   we implement it through a lambda expression ---*/
3377 *   const auto cell_worker = [&](const Iterator& cell,
3378 *   ScratchData<dim>& scratch_data,
3379 *   CopyData& copy_data) {
3380 *   FEValues<dim>& fe_values = scratch_data.fe_values; /*--- Here we finally use the 'FEValues' inside ScratchData ---*/
3381 *   fe_values.reinit(cell);
3382 *  
3383 *   /*--- Compute the gradients for all quadrature points ---*/
3384 *   std::vector<std::vector<Tensor<1, dim>>> gradients(fe_values.n_quadrature_points, std::vector<Tensor<1, dim>>(dim));
3385 *   fe_values.get_function_gradients(tmp_velocity, gradients);
3386 *   copy_data.cell_index = cell->active_cell_index();
3387 *   double vorticity_norm_square = 0.0;
3388 *   /*--- Loop over quadrature points and evaluate the integral multiplying the vorticty
3389 *   by the weights and the determinant of the Jacobian (which are included in 'JxW') ---*/
3390 *   for(unsigned k = 0; k < fe_values.n_quadrature_points; ++k) {
3391 *   const double vorticity = gradients[k][1][0] - gradients[k][0][1];
3392 *   vorticity_norm_square += vorticity*vorticity*fe_values.JxW(k);
3393 *   }
3394 *   copy_data.value = cell->diameter()*cell->diameter()*vorticity_norm_square;
3395 *   };
3396 *  
3398 *  
3399 *   const auto copier = [&](const CopyData &copy_data) {
3400 *   if(copy_data.cell_index != numbers::invalid_unsigned_int)
3401 *   estimated_error_per_cell[copy_data.cell_index] += copy_data.value;
3402 *   };
3403 *  
3404 *   /*--- Now everything is 'automagically' handled by 'mesh_loop' ---*/
3405 *   ScratchData<dim> scratch_data(fe_velocity, EquationData::degree_p + 2, cell_flags);
3406 *   CopyData copy_data;
3407 *   MeshWorker::mesh_loop(dof_handler_velocity.begin_active(),
3408 *   dof_handler_velocity.end(),
3409 *   cell_worker,
3410 *   copier,
3411 *   scratch_data,
3412 *   copy_data,
3414 *  
3415 *   /*--- Refine grid. In case the refinement level is above a certain value (or the coarsening level is below)
3416 *   we clear the flags. ---*/
3418 *   for(const auto& cell: triangulation.active_cell_iterators()) {
3419 *   if(cell->refine_flag_set() && static_cast<unsigned int>(cell->level()) == max_loc_refinements)
3420 *   cell->clear_refine_flag();
3421 *   if(cell->coarsen_flag_set() && static_cast<unsigned int>(cell->level()) == min_loc_refinements)
3422 *   cell->clear_coarsen_flag();
3423 *   }
3425 *  
3426 *   /*--- Now we prepare the object for transferring, basically saving the old quantities using SolutionTransfer.
3427 *   Since the 'prepare_for_coarsening_and_refinement' method can be called only once, but we have two vectors
3428 *   for dof_handler_velocity, we need to put them in an auxiliary vector. ---*/
3429 *   std::vector<const LinearAlgebra::distributed::Vector<double>*> velocities;
3430 *   velocities.push_back(&u_n);
3431 *   velocities.push_back(&u_n_minus_1);
3433 *   solution_transfer_velocity(dof_handler_velocity);
3434 *   solution_transfer_velocity.prepare_for_coarsening_and_refinement(velocities);
3436 *   solution_transfer_pressure(dof_handler_pressure);
3437 *   solution_transfer_pressure.prepare_for_coarsening_and_refinement(pres_n);
3438 *  
3439 *   triangulation.execute_coarsening_and_refinement(); /*--- Effectively perform the remeshing ---*/
3440 *  
3441 *   /*--- First DoFHandler objects are set up within the new grid ----*/
3442 *   setup_dofs();
3443 *  
3444 *   /*--- Interpolate current solutions to new mesh. This is done using auxliary vectors just for safety,
3445 *   but the new u_n or pres_n could be used. Again, the only point is that the function 'interpolate'
3446 *   can be called once and so the vectors related to 'dof_handler_velocity' have to collected in an auxiliary vector. ---*/
3447 *   LinearAlgebra::distributed::Vector<double> transfer_velocity,
3448 *   transfer_velocity_minus_1,
3449 *   transfer_pressure;
3450 *   transfer_velocity.reinit(u_n);
3451 *   transfer_velocity.zero_out_ghost_values();
3452 *   transfer_velocity_minus_1.reinit(u_n_minus_1);
3453 *   transfer_velocity_minus_1.zero_out_ghost_values();
3454 *   transfer_pressure.reinit(pres_n);
3455 *   transfer_pressure.zero_out_ghost_values();
3456 *  
3457 *   std::vector<LinearAlgebra::distributed::Vector<double>*> transfer_velocities;
3458 *   transfer_velocities.push_back(&transfer_velocity);
3459 *   transfer_velocities.push_back(&transfer_velocity_minus_1);
3460 *   solution_transfer_velocity.interpolate(transfer_velocities);
3461 *   transfer_velocity.update_ghost_values();
3462 *   transfer_velocity_minus_1.update_ghost_values();
3463 *   solution_transfer_pressure.interpolate(transfer_pressure);
3464 *   transfer_pressure.update_ghost_values();
3465 *  
3466 *   u_n = transfer_velocity;
3467 *   u_n_minus_1 = transfer_velocity_minus_1;
3468 *   pres_n = transfer_pressure;
3469 *   }
3470 *  
3471 *  
3472 * @endcode
3473 *
3474 * Interpolate the locally refined solution to a mesh with maximal resolution
3475 * and transfer velocity and pressure.
3476 *
3477
3478 *
3479 *
3480 * @code
3481 *   template<int dim>
3482 *   void NavierStokesProjection<dim>::interpolate_max_res(const unsigned int level) {
3484 *   solution_transfer_velocity(dof_handler_velocity);
3485 *   std::vector<const LinearAlgebra::distributed::Vector<double>*> velocities;
3486 *   velocities.push_back(&u_n);
3487 *   velocities.push_back(&u_n_minus_1);
3488 *   solution_transfer_velocity.prepare_for_coarsening_and_refinement(velocities);
3489 *  
3491 *   solution_transfer_pressure(dof_handler_pressure);
3492 *   solution_transfer_pressure.prepare_for_coarsening_and_refinement(pres_n);
3493 *  
3494 *   for(const auto& cell: triangulation.active_cell_iterators_on_level(level)) {
3495 *   if(cell->is_locally_owned())
3496 *   cell->set_refine_flag();
3497 *   }
3499 *  
3500 *   setup_dofs();
3501 *  
3502 *   LinearAlgebra::distributed::Vector<double> transfer_velocity, transfer_velocity_minus_1,
3503 *   transfer_pressure;
3504 *  
3505 *   transfer_velocity.reinit(u_n);
3506 *   transfer_velocity.zero_out_ghost_values();
3507 *   transfer_velocity_minus_1.reinit(u_n_minus_1);
3508 *   transfer_velocity_minus_1.zero_out_ghost_values();
3509 *  
3510 *   transfer_pressure.reinit(pres_n);
3511 *   transfer_pressure.zero_out_ghost_values();
3512 *  
3513 *   std::vector<LinearAlgebra::distributed::Vector<double>*> transfer_velocities;
3514 *  
3515 *   transfer_velocities.push_back(&transfer_velocity);
3516 *   transfer_velocities.push_back(&transfer_velocity_minus_1);
3517 *   solution_transfer_velocity.interpolate(transfer_velocities);
3518 *   transfer_velocity.update_ghost_values();
3519 *   transfer_velocity_minus_1.update_ghost_values();
3520 *  
3521 *   solution_transfer_pressure.interpolate(transfer_pressure);
3522 *   transfer_pressure.update_ghost_values();
3523 *  
3524 *   u_n = transfer_velocity;
3525 *   u_n_minus_1 = transfer_velocity_minus_1;
3526 *   pres_n = transfer_pressure;
3527 *   }
3528 *  
3529 *  
3530 * @endcode
3531 *
3532 * Save maximum resolution to a mesh adapted.
3533 *
3534
3535 *
3536 *
3537 * @code
3538 *   template<int dim>
3539 *   void NavierStokesProjection<dim>::save_max_res() {
3540 *   parallel::distributed::Triangulation<dim> triangulation_tmp(MPI_COMM_WORLD);
3541 *   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);
3542 *   triangulation_tmp.refine_global(triangulation.n_global_levels() - 1);
3543 *  
3544 *   DoFHandler<dim> dof_handler_velocity_tmp(triangulation_tmp);
3545 *   DoFHandler<dim> dof_handler_pressure_tmp(triangulation_tmp);
3546 *   dof_handler_velocity_tmp.distribute_dofs(fe_velocity);
3547 *   dof_handler_pressure_tmp.distribute_dofs(fe_pressure);
3548 *  
3550 *   pres_n_tmp;
3551 *   u_n_tmp.reinit(dof_handler_velocity_tmp.n_dofs());
3552 *   pres_n_tmp.reinit(dof_handler_pressure_tmp.n_dofs());
3553 *  
3554 *   DataOut<dim> data_out;
3555 *   std::vector<std::string> velocity_names(dim, "v");
3556 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3557 *   component_interpretation_velocity(dim, DataComponentInterpretation::component_is_part_of_vector);
3558 *   VectorTools::interpolate_to_different_mesh(dof_handler_velocity, u_n, dof_handler_velocity_tmp, u_n_tmp);
3559 *   u_n_tmp.update_ghost_values();
3560 *   data_out.add_data_vector(dof_handler_velocity_tmp, u_n_tmp, velocity_names, component_interpretation_velocity);
3561 *   VectorTools::interpolate_to_different_mesh(dof_handler_pressure, pres_n, dof_handler_pressure_tmp, pres_n_tmp);
3562 *   pres_n_tmp.update_ghost_values();
3563 *   data_out.add_data_vector(dof_handler_pressure_tmp, pres_n_tmp, "p", {DataComponentInterpretation::component_is_scalar});
3564 *   PostprocessorVorticity<dim> postprocessor;
3565 *   data_out.add_data_vector(dof_handler_velocity_tmp, u_n_tmp, postprocessor);
3566 *  
3567 *   data_out.build_patches(MappingQ1<dim>(), 1, DataOut<dim>::curved_inner_cells);
3568 *   const std::string output = "./" + saving_dir + "/solution_max_res_end.vtu";
3569 *   data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
3570 *   }
3571 *  
3572 *  
3573 * @endcode
3574 *
3575 *
3576 * <a name="navier_stokes_TRBDF2_DG.cc-"></a>
3577 * @sect{ <code>NavierStokesProjection::run</code> }
3578 *
3579
3580 *
3581 * This is the time marching function, which starting at <code>t_0</code>
3582 * advances in time using the projection method with time step <code>dt</code>
3583 * until <code>T</code>.
3584 *
3585
3586 *
3587 * Its second parameter, <code>verbose</code> indicates whether the function
3588 * should output information what it is doing at any given moment:
3589 * we use the ConditionalOStream class to do that for us.
3590 *
3591
3592 *
3593 *
3594 * @code
3595 *   template<int dim>
3596 *   void NavierStokesProjection<dim>::run(const bool verbose, const unsigned int output_interval) {
3597 *   ConditionalOStream verbose_cout(std::cout, verbose && Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0);
3598 *  
3599 *   output_results(1);
3600 *   double time = t_0 + dt;
3601 *   unsigned int n = 1;
3602 *   while(std::abs(T - time) > 1e-10) {
3603 *   time += dt;
3604 *   n++;
3605 *   pcout << "Step = " << n << " Time = " << time << std::endl;
3606 *  
3607 *   /*--- First stage of TR-BDF2 and we start by setting the proper flag ---*/
3608 *   TR_BDF2_stage = 1;
3609 *   navier_stokes_matrix.set_TR_BDF2_stage(TR_BDF2_stage);
3610 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
3611 *   mg_matrices[level].set_TR_BDF2_stage(TR_BDF2_stage);
3612 *  
3613 *   verbose_cout << " Interpolating the velocity stage 1" << std::endl;
3614 *   interpolate_velocity();
3615 *  
3616 *   verbose_cout << " Diffusion Step stage 1 " << std::endl;
3617 *   diffusion_step();
3618 *  
3619 *   verbose_cout << " Projection Step stage 1" << std::endl;
3620 *   project_grad(1);
3621 *   u_tmp.equ(gamma*dt, u_tmp);
3622 *   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 ---*/
3623 *   projection_step();
3624 *  
3625 *   verbose_cout << " Updating the Velocity stage 1" << std::endl;
3626 *   u_n_gamma.equ(1.0, u_star);
3627 *   project_grad(2);
3628 *   grad_pres_int.equ(1.0, u_tmp); /*--- We save grad(pres_int), because we will need it soon ---*/
3629 *   u_tmp.equ(-gamma*dt, u_tmp);
3630 *   u_n_gamma += u_tmp; /*--- u_n_gamma = u_star - gamma*dt*grad(pres_int) ---*/
3631 *   u_n_minus_1 = u_n;
3632 *  
3633 *   /*--- Second stage of TR-BDF2 ---*/
3634 *   TR_BDF2_stage = 2;
3635 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
3636 *   mg_matrices[level].set_TR_BDF2_stage(TR_BDF2_stage);
3637 *   navier_stokes_matrix.set_TR_BDF2_stage(TR_BDF2_stage);
3638 *  
3639 *   verbose_cout << " Interpolating the velocity stage 2" << std::endl;
3640 *   interpolate_velocity();
3641 *  
3642 *   verbose_cout << " Diffusion Step stage 2 " << std::endl;
3643 *   diffusion_step();
3644 *  
3645 *   verbose_cout << " Projection Step stage 2" << std::endl;
3646 *   u_tmp.equ((1.0 - gamma)*dt, grad_pres_int);
3647 *   u_star += u_tmp; /*--- In the rhs of the projection step we need u_star + (1 - gamma)*dt*grad(pres_int) ---*/
3648 *   projection_step();
3649 *  
3650 *   verbose_cout << " Updating the Velocity stage 2" << std::endl;
3651 *   u_n.equ(1.0, u_star);
3652 *   project_grad(1);
3653 *   u_tmp.equ((gamma - 1.0)*dt, u_tmp);
3654 *   u_n += u_tmp; /*--- u_n = u_star - (1 - gamma)*dt*grad(pres_n) ---*/
3655 *  
3656 *   const double max_vel = get_maximal_velocity();
3657 *   pcout<< "Maximal velocity = " << max_vel << std::endl;
3658 *   /*--- The Courant number is computed taking into account the polynomial degree for the velocity ---*/
3659 *   pcout << "CFL = " << dt*max_vel*(EquationData::degree_p + 1)*
3661 *   compute_lift_and_drag();
3662 *   if(n % output_interval == 0) {
3663 *   verbose_cout << "Plotting Solution final" << std::endl;
3664 *   output_results(n);
3665 *   }
3666 *   /*--- In case dt is not a multiple of T, we reduce dt in order to end up at T ---*/
3667 *   if(T - time < dt && T - time > 1e-10) {
3668 *   dt = T - time;
3669 *   navier_stokes_matrix.set_dt(dt);
3670 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
3671 *   mg_matrices[level].set_dt(dt);
3672 *   }
3673 *   /*--- Perform the refinement if desired ---*/
3674 *   if(refinement_iterations > 0 && n % refinement_iterations == 0) {
3675 *   verbose_cout << "Refining mesh" << std::endl;
3676 *   refine_mesh();
3677 *   }
3678 *   }
3679 *   if(n % output_interval != 0) {
3680 *   verbose_cout << "Plotting Solution final" << std::endl;
3681 *   output_results(n);
3682 *   }
3683 *   if(refinement_iterations > 0) {
3684 *   for(unsigned int lev = 0; lev < triangulation.n_global_levels() - 1; ++ lev)
3685 *   interpolate_max_res(lev);
3686 *   save_max_res();
3687 *   }
3688 *   }
3689 *  
3690 *   } // namespace NS_TRBDF2
3691 *  
3692 *  
3693 * @endcode
3694 *
3695 *
3696 * <a name="navier_stokes_TRBDF2_DG.cc-"></a>
3697 * @sect{ The main function }
3698 *
3699
3700 *
3701 * The main function looks very much like in all the other tutorial programs. We first initialize MPI,
3702 * we initialize the class 'NavierStokesProjection' with the dimension as template parameter and then
3703 * let the method 'run' do the job.
3704 *
3705
3706 *
3707 *
3708 * @code
3709 *   int main(int argc, char *argv[]) {
3710 *   try {
3711 *   using namespace NS_TRBDF2;
3712 *  
3713 *   RunTimeParameters::Data_Storage data;
3714 *   data.read_data("parameter-file.prm");
3715 *  
3716 *   Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, -1);
3717 *  
3718 *   const auto& curr_rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
3719 *   deallog.depth_console(data.verbose && curr_rank == 0 ? 2 : 0);
3720 *  
3721 *   NavierStokesProjection<2> test(data);
3722 *   test.run(data.verbose, data.output_interval);
3723 *  
3724 *   if(curr_rank == 0)
3725 *   std::cout << "----------------------------------------------------"
3726 *   << std::endl
3727 *   << "Apparently everything went fine!" << std::endl
3728 *   << "Don't forget to brush your teeth :-)" << std::endl
3729 *   << std::endl;
3730 *  
3731 *   return 0;
3732 *   }
3733 *   catch(std::exception &exc) {
3734 *   std::cerr << std::endl
3735 *   << std::endl
3736 *   << "----------------------------------------------------"
3737 *   << std::endl;
3738 *   std::cerr << "Exception on processing: " << std::endl
3739 *   << exc.what() << std::endl
3740 *   << "Aborting!" << std::endl
3741 *   << "----------------------------------------------------"
3742 *   << std::endl;
3743 *   return 1;
3744 *   }
3745 *   catch(...) {
3746 *   std::cerr << std::endl
3747 *   << std::endl
3748 *   << "----------------------------------------------------"
3749 *   << std::endl;
3750 *   std::cerr << "Unknown exception!" << std::endl
3751 *   << "Aborting!" << std::endl
3752 *   << "----------------------------------------------------"
3753 *   << std::endl;
3754 *   return 1;
3755 *   }
3756 *  
3757 *   }
3758 * @endcode
3759
3760
3761<a name="ann-runtime_parameters.h"></a>
3762<h1>Annotated version of runtime_parameters.h</h1>
3763 *
3764 *
3765 *
3766 *
3767 * @code
3768 *   /* -----------------------------------------------------------------------------
3769 *   *
3770 *   * SPDX-License-Identifier: LGPL-2.1-or-later
3771 *   * Copyright (C) 2022 by Giuseppe Orlando
3772 *   *
3773 *   * This file is part of the deal.II code gallery.
3774 *   *
3775 *   * -----------------------------------------------------------------------------
3776 *   */
3777 *  
3778 * @endcode
3779 *
3780 * We start by including all the necessary deal.II header files
3781 *
3782
3783 *
3784 *
3785 * @code
3786 *   #include <deal.II/base/parameter_handler.h>
3787 *  
3788 * @endcode
3789 *
3790 *
3791 * <a name="runtime_parameters.h-"></a>
3792 * @sect{Run time parameters}
3793 *
3794
3795 *
3796 * Since our method has several parameters that can be fine-tuned we put them
3797 * into an external file, so that they can be determined at run-time.
3798 *
3799
3800 *
3801 *
3802 * @code
3803 *   namespace RunTimeParameters {
3804 *   using namespace dealii;
3805 *  
3806 *   class Data_Storage {
3807 *   public:
3808 *   Data_Storage();
3809 *  
3810 *   void read_data(const std::string& filename);
3811 *  
3812 *   double initial_time;
3813 *   double final_time;
3814 *  
3815 *   double Reynolds;
3816 *   double dt;
3817 *  
3818 *   unsigned int n_refines; /*--- Number of refinements ---*/
3819 *   unsigned int max_loc_refinements; /*--- Number of maximum local refinements allowed ---*/
3820 *   unsigned int min_loc_refinements; /*--- Number of minimum local refinements allowed
3821 *   once reached that level ---*/
3822 *  
3823 *   /*--- Parameters related to the linear solver ---*/
3824 *   unsigned int max_iterations;
3825 *   double eps;
3826 *  
3827 *   bool verbose;
3828 *   unsigned int output_interval;
3829 *  
3830 *   std::string dir; /*--- Auxiliary string variable for output storage ---*/
3831 *  
3832 *   unsigned int refinement_iterations; /*--- Auxiliary variable about how many steps perform remeshing ---*/
3833 *  
3834 *   protected:
3835 *   ParameterHandler prm;
3836 *   };
3837 *  
3838 * @endcode
3839 *
3840 * In the constructor of this class we declare all the parameters in suitable (but arbitrary) subsections.
3841 *
3842
3843 *
3844 *
3845 * @code
3846 *   Data_Storage::Data_Storage(): initial_time(0.0),
3847 *   final_time(1.0),
3848 *   Reynolds(1.0),
3849 *   dt(5e-4),
3850 *   n_refines(0),
3851 *   max_loc_refinements(0),
3852 *   min_loc_refinements(0),
3853 *   max_iterations(1000),
3854 *   eps(1e-12),
3855 *   verbose(true),
3856 *   output_interval(15),
3857 *   refinement_iterations(0) {
3858 *   prm.enter_subsection("Physical data");
3859 *   {
3860 *   prm.declare_entry("initial_time",
3861 *   "0.0",
3862 *   Patterns::Double(0.0),
3863 *   " The initial time of the simulation. ");
3864 *   prm.declare_entry("final_time",
3865 *   "1.0",
3866 *   Patterns::Double(0.0),
3867 *   " The final time of the simulation. ");
3868 *   prm.declare_entry("Reynolds",
3869 *   "1.0",
3870 *   Patterns::Double(0.0),
3871 *   " The Reynolds number. ");
3872 *   }
3873 *   prm.leave_subsection();
3874 *  
3875 *   prm.enter_subsection("Time step data");
3876 *   {
3877 *   prm.declare_entry("dt",
3878 *   "5e-4",
3879 *   Patterns::Double(0.0),
3880 *   " The time step size. ");
3881 *   }
3882 *   prm.leave_subsection();
3883 *  
3884 *   prm.enter_subsection("Space discretization");
3885 *   {
3886 *   prm.declare_entry("n_of_refines",
3887 *   "100",
3888 *   Patterns::Integer(0, 1500),
3889 *   " The number of cells we want on each direction of the mesh. ");
3890 *   prm.declare_entry("max_loc_refinements",
3891 *   "4",
3892 *   Patterns::Integer(0, 10),
3893 *   " The number of maximum local refinements. ");
3894 *   prm.declare_entry("min_loc_refinements",
3895 *   "2",
3896 *   Patterns::Integer(0, 10),
3897 *   " The number of minimum local refinements. ");
3898 *   }
3899 *   prm.leave_subsection();
3900 *  
3901 *   prm.enter_subsection("Data solve");
3902 *   {
3903 *   prm.declare_entry("max_iterations",
3904 *   "1000",
3905 *   Patterns::Integer(1, 30000),
3906 *   " The maximal number of iterations linear solvers must make. ");
3907 *   prm.declare_entry("eps",
3908 *   "1e-12",
3909 *   Patterns::Double(0.0),
3910 *   " The stopping criterion. ");
3911 *   }
3912 *   prm.leave_subsection();
3913 *  
3914 *   prm.declare_entry("refinement_iterations",
3915 *   "0",
3916 *   Patterns::Integer(0),
3917 *   " This number indicates how often we need to "
3918 *   "refine the mesh");
3919 *  
3920 *   prm.declare_entry("saving directory", "SimTest");
3921 *  
3922 *   prm.declare_entry("verbose",
3923 *   "true",
3924 *   Patterns::Bool(),
3925 *   " This indicates whether the output of the solution "
3926 *   "process should be verbose. ");
3927 *  
3928 *   prm.declare_entry("output_interval",
3929 *   "1",
3930 *   Patterns::Integer(1),
3931 *   " This indicates between how many time steps we print "
3932 *   "the solution. ");
3933 *   }
3934 *  
3935 * @endcode
3936 *
3937 * We need now a routine to read all declared parameters in the constructor
3938 *
3939
3940 *
3941 *
3942 * @code
3943 *   void Data_Storage::read_data(const std::string& filename) {
3944 *   std::ifstream file(filename);
3945 *   AssertThrow(file, ExcFileNotOpen(filename));
3946 *  
3947 *   prm.parse_input(file);
3948 *  
3949 *   prm.enter_subsection("Physical data");
3950 *   {
3951 *   initial_time = prm.get_double("initial_time");
3952 *   final_time = prm.get_double("final_time");
3953 *   Reynolds = prm.get_double("Reynolds");
3954 *   }
3955 *   prm.leave_subsection();
3956 *  
3957 *   prm.enter_subsection("Time step data");
3958 *   {
3959 *   dt = prm.get_double("dt");
3960 *   }
3961 *   prm.leave_subsection();
3962 *  
3963 *   prm.enter_subsection("Space discretization");
3964 *   {
3965 *   n_refines = prm.get_integer("n_of_refines");
3966 *   max_loc_refinements = prm.get_integer("max_loc_refinements");
3967 *   min_loc_refinements = prm.get_integer("min_loc_refinements");
3968 *   }
3969 *   prm.leave_subsection();
3970 *  
3971 *   prm.enter_subsection("Data solve");
3972 *   {
3973 *   max_iterations = prm.get_integer("max_iterations");
3974 *   eps = prm.get_double("eps");
3975 *   }
3976 *   prm.leave_subsection();
3977 *  
3978 *   dir = prm.get("saving directory");
3979 *  
3980 *   refinement_iterations = prm.get_integer("refinement_iterations");
3981 *  
3982 *   verbose = prm.get_bool("verbose");
3983 *  
3984 *   output_interval = prm.get_integer("output_interval");
3985 *   }
3986 *  
3987 *   } // namespace RunTimeParameters
3988 * @endcode
3989
3990
3991*/
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:349
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
void build(const DoFHandler< dim > &dof_handler, const std::vector< std::shared_ptr< const Utilities::MPI::Partitioner > > &external_partitioners=std::vector< std::shared_ptr< const Utilities::MPI::Partitioner > >())
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void clear()
unsigned int mg_level
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
Definition point.h:111
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
virtual types::global_cell_index n_global_active_cells() const override
Definition tria_base.cc:151
virtual unsigned int n_global_levels() const override
Definition tria_base.cc:141
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
virtual bool prepare_coarsening_and_refinement() override
Definition tria.cc:2805
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
unsigned int cell_index
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:281
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
UpdateFlags
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
LogStream deallog
Definition logstream.cc:36
const Event initial
Definition event.cc:64
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > &center=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &cell_operation, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >()), const bool project_to_boundary_first=false)
void interpolate_to_different_mesh(const DoFHandler< dim, spacedim > &dof_handler_1, const VectorType &u1, const DoFHandler< dim, spacedim > &dof_handler_2, VectorType &u2)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
int(& functions)(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
Definition mg.h:81
static const unsigned int invalid_unsigned_int
Definition types.h:220
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int boundary_id
Definition types.h:144
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector<::Vector< double > > solution_values
std::vector< std::vector< Tensor< 1, spacedim > > > solution_gradients
TasksParallelScheme tasks_parallel_scheme
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
std_cxx20::type_identity< T > identity
VectorizedArray< Number, width > make_vectorized_array(const Number &u)