Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+00:00
\(\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
advection_reaction_estimator.h
Go to the documentation of this file.
1
214 *  
215 *   #ifndef INCLUDE_DG_UPWIND_H_
216 *   #define INCLUDE_DG_UPWIND_H_
217 *  
218 *  
219 *   #include <deal.II/base/function.h>
220 *   #include <deal.II/base/quadrature_lib.h>
221 *  
222 *   #include <deal.II/dofs/dof_handler.h>
223 *   #include <deal.II/dofs/dof_tools.h>
224 *  
225 *   #include <deal.II/fe/fe_dgq.h>
226 *   #include <deal.II/fe/fe_q.h>
227 *   #include <deal.II/fe/fe_values.h>
228 *   #include <deal.II/fe/mapping_q1.h>
229 *  
230 *   #include <deal.II/grid/grid_generator.h>
231 *   #include <deal.II/grid/grid_out.h>
232 *   #include <deal.II/grid/grid_refinement.h>
233 *   #include <deal.II/grid/tria.h>
234 *  
235 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
236 *   #include <deal.II/lac/full_matrix.h>
237 *   #include <deal.II/lac/sparse_matrix.h>
238 *   #include <deal.II/lac/vector.h>
239 *  
240 *   #include <deal.II/numerics/data_out.h>
241 *   #include <deal.II/numerics/vector_tools.h>
242 * @endcode
243 *
244 * This header is needed for FEInterfaceValues to compute integrals on
245 * interfaces:
246 *
247 * @code
248 *   #include <deal.II/fe/fe_interface_values.h>
249 * @endcode
250 *
251 * Solver
252 *
253 * @code
254 *   #include <deal.II/lac/precondition_block.h>
255 *   #include <deal.II/lac/solver_richardson.h>
256 *   #include <deal.II/lac/sparse_direct.h>
257 * @endcode
258 *
259 * We are going to use gradients as refinement indicator.
260 *
261 * @code
262 *   #include <deal.II/numerics/derivative_approximation.h>
263 * @endcode
264 *
265 * Using using the mesh_loop from the MeshWorker framework
266 *
267 * @code
268 *   #include <deal.II/base/convergence_table.h>
269 *  
270 *   #include <deal.II/meshworker/mesh_loop.h>
271 *  
272 * @endcode
273 *
274 * To enable parameter handling
275 *
276 * @code
277 *   #include <deal.II/base/function_parser.h>
278 *   #include <deal.II/base/parameter_acceptor.h>
279 *   #include <deal.II/base/parameter_handler.h>
280 *   #include <deal.II/base/parsed_convergence_table.h>
281 *   #include <deal.II/base/symbolic_function.h>
282 *  
283 *   #include <deal.II/meshworker/copy_data.h>
284 *   #include <deal.II/meshworker/mesh_loop.h>
285 *   #include <deal.II/meshworker/scratch_data.h>
286 *  
287 *   #include <fstream>
288 *   #include <iostream>
289 *   using namespace dealii;
290 *  
291 * @endcode
292 *
293 * This is a struct used only for throwing an exception when theta parameter is
294 * not okay.
295 *
296 * @code
297 *   struct theta_exc
298 *   {
299 *   std::string message;
300 *   theta_exc(std::string &&s)
301 *   : message{std::move(s)} {};
302 *   const char *
303 *   what() const
304 *   {
305 *   return message.c_str();
306 *   }
307 *   };
308 *  
309 *  
310 *   template <int dim>
311 *   class AdvectionReaction : ParameterAcceptor
312 *   {
313 *   public:
314 *   AdvectionReaction();
315 *   void
316 *   initialize_params(const std::string &filename);
317 *   void
318 *   run();
319 *  
320 *   private:
321 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
322 *   void
323 *   parse_string(const std::string &parameters);
324 *   void
325 *   setup_system();
326 *   void
327 *   assemble_system();
328 *   void
329 *   solve();
330 *   void
331 *   refine_grid();
332 *   void
333 *   output_results(const unsigned int cycle) const;
334 *   void
335 *   compute_error();
336 *   double
337 *   compute_energy_norm();
338 *   void
339 *   compute_local_projection_and_estimate();
340 *  
342 *   const MappingQ1<dim> mapping;
343 *  
344 * @endcode
345 *
346 * Furthermore we want to use DG elements.
347 *
348 * @code
349 *   std::unique_ptr<FE_DGQ<dim>> fe;
350 *   DoFHandler<dim> dof_handler;
351 *  
352 *   SparsityPattern sparsity_pattern;
353 *   SparseMatrix<double> system_matrix;
354 *  
355 *   Vector<double> solution;
356 *   Vector<double> right_hand_side;
357 *   Vector<double> energy_norm_square_per_cell;
358 *   Vector<double> error_indicator_per_cell;
359 *  
360 * @endcode
361 *
362 * So far we declared the usual objects. Hereafter we declare
363 * `FunctionParser<dim>` objects
364 *
365 * @code
366 *   FunctionParser<dim> exact_solution;
367 *   FunctionParser<dim> boundary_conditions;
368 *   FunctionParser<dim> rhs;
369 *   FunctionParser<dim> advection_coeff;
370 *  
371 *   unsigned int fe_degree = 1;
372 *  
373 * @endcode
374 *
375 * and then we define default values that will be parsed from the following
376 * strings
377 *
378 * @code
379 *   std::string exact_solution_expression =
380 *   "tanh(100*(x+y-0.5))"; // internal layer solution
381 *   std::string rhs_expression =
382 *   "-200*tanh(100*x + 100*y - 50.0)^2 + tanh(100*x + 100*y - 50.0) + 200";
383 *   std::string advection_coefficient_expression = "1.0";
384 *   std::string boundary_conditions_expression = "tanh(100*x + 100*y - 50.0)";
385 *   std::string refinement = "residual";
386 *   std::string output_filename = "DG_advection_reaction_estimator";
387 *   std::map<std::string, double> constants;
388 *   ParsedConvergenceTable error_table;
389 *  
390 *   bool use_direct_solver = true;
391 *   unsigned int n_refinement_cycles = 8;
392 *   unsigned int n_global_refinements = 3;
393 *   double theta = 0.5; // default is 0.5 so that I have classical upwind flux
394 *   };
395 *  
396 *   #endif /* INCLUDE_DG_UPWIND_H_ */
397 * @endcode
398
399
400<a name="ann-main.cc"></a>
401<h1>Annotated version of main.cc</h1>
402 *
403 *
404 *
405 *
406 * @code
407 *   /* -----------------------------------------------------------------------------
408 *   *
409 *   * SPDX-License-Identifier: LGPL-2.1-or-later
410 *   * Copyright (C) 2015 by Marco Feder
411 *   *
412 *   * This file is part of the deal.II code gallery.
413 *   *
414 *   * -----------------------------------------------------------------------------
415 *   */
416 *  
417 *   #include "include/DG_advection_reaction.h"
418 *  
419 *   int
420 *   main(int argc, char **argv)
421 *   {
422 *   try
423 *   {
424 *   std::string par_name = "";
425 *   if (argc > 1)
426 *   {
427 *   par_name = argv[1];
428 *   }
429 *   else
430 *   {
431 *   par_name = "parameters.prm";
432 *   }
433 *   deallog.depth_console(2);
434 *   AdvectionReaction<2> problem;
435 *   problem.initialize_params(par_name);
436 *   problem.run();
437 *   }
438 *   catch (std::exception &exc)
439 *   {
440 *   std::cerr << std::endl
441 *   << std::endl
442 *   << "----------------------------------------------------"
443 *   << std::endl;
444 *   std::cerr << "Exception on processing: " << std::endl
445 *   << exc.what() << std::endl
446 *   << "Aborting!" << std::endl
447 *   << "----------------------------------------------------"
448 *   << std::endl;
449 *   return 1;
450 *   }
451 *   catch (const theta_exc &theta_range)
452 *   {
453 *   std::cerr << std::endl
454 *   << std::endl
455 *   << "----------------------------------------------------"
456 *   << std::endl;
457 *   std::cerr << "Exception on processing: " << std::endl
458 *   << theta_range.what() << std::endl
459 *   << "Aborting!" << std::endl
460 *   << "----------------------------------------------------"
461 *   << std::endl;
462 *   return 1;
463 *   }
464 *   catch (...)
465 *   {
466 *   std::cerr << std::endl
467 *   << std::endl
468 *   << "----------------------------------------------------"
469 *   << std::endl;
470 *   std::cerr << "Unknown exception!" << std::endl
471 *   << "Aborting!" << std::endl
472 *   << "----------------------------------------------------"
473 *   << std::endl;
474 *   return 1;
475 *   }
476 *  
477 *   return 0;
478 *   }
479 * @endcode
480
481
482<a name="ann-source/DG_advection_reaction.cc"></a>
483<h1>Annotated version of source/DG_advection_reaction.cc</h1>
484 *
485 *
486 *
487 *
488 * @code
489 *   /* -----------------------------------------------------------------------------
490 *   *
491 *   * SPDX-License-Identifier: LGPL-2.1-or-later
492 *   * Copyright (C) 2015 by Marco Feder
493 *   *
494 *   * This file is part of the deal.II code gallery.
495 *   *
496 *   * -----------------------------------------------------------------------------
497 *   */
498 *  
499 *   #include "../include/DG_advection_reaction.h"
500 *  
501 * @endcode
502 *
503 * Compute and returns the wind field b
504 *
505 * @code
506 *   template <int dim>
507 *   Tensor<1, dim>
508 *   beta(const Point<dim> &p)
509 *   {
510 *   Assert(dim > 1, ExcNotImplemented());
511 *   (void)p; // suppress warnings from p
512 *   Tensor<1, dim> wind_field;
513 *   wind_field[0] = 1.0;
514 *   wind_field[1] = 1.0;
515 *  
516 *   return wind_field;
517 *   }
518 *  
519 * @endcode
520 *
521 *
522 * <a name="source/DG_advection_reaction.cc-TheScratchDataandCopyDataclasses"></a>
523 * <h3>The ScratchData and CopyData classes</h3>
524 *
525
526 *
527 * The following objects are the scratch and copy objects we use in the call
528 * to MeshWorker::mesh_loop(). The new object is the FEInterfaceValues object,
529 * that works similar to FEValues or FEFacesValues, except that it acts on
530 * an interface between two cells and allows us to assemble the interface
531 * terms in our weak form.
532 *
533 * @code
534 *   template <int dim>
535 *   struct ScratchData
536 *   {
537 *   ScratchData(const Mapping<dim> &mapping,
538 *   const FiniteElement<dim> &fe,
539 *   const Quadrature<dim> &quadrature,
540 *   const Quadrature<dim - 1> &quadrature_face,
541 *   const UpdateFlags update_flags = update_values |
542 *   update_gradients |
545 *   const UpdateFlags interface_update_flags =
548 *   : fe_values(mapping, fe, quadrature, update_flags)
549 *   , fe_interface_values(mapping, fe, quadrature_face, interface_update_flags)
550 *   {}
551 *  
552 *   ScratchData(const ScratchData<dim> &scratch_data)
553 *   : fe_values(scratch_data.fe_values.get_mapping(),
554 *   scratch_data.fe_values.get_fe(),
555 *   scratch_data.fe_values.get_quadrature(),
556 *   scratch_data.fe_values.get_update_flags())
557 *   , fe_interface_values(scratch_data.fe_interface_values.get_mapping(),
558 *   scratch_data.fe_interface_values.get_fe(),
559 *   scratch_data.fe_interface_values.get_quadrature(),
560 *   scratch_data.fe_interface_values.get_update_flags())
561 *   {}
562 *  
563 *   FEValues<dim> fe_values;
564 *   FEInterfaceValues<dim> fe_interface_values;
565 *   };
566 *  
567 *  
568 *  
569 *   struct CopyDataFace
570 *   {
572 *   std::vector<types::global_dof_index> joint_dof_indices;
573 *   std::array<double, 2> values;
574 *   std::array<unsigned int, 2> cell_indices;
575 *   };
576 *  
577 *  
578 *  
579 *   struct CopyData
580 *   {
582 *   Vector<double> cell_rhs;
583 *   std::vector<types::global_dof_index> local_dof_indices;
584 *   std::vector<CopyDataFace> face_data;
585 *  
586 *   double value;
587 *   double value_estimator;
588 *   unsigned int cell_index;
589 *  
590 *   FullMatrix<double> cell_mass_matrix;
591 *   Vector<double> cell_mass_rhs;
592 *  
593 *   template <class Iterator>
594 *   void
595 *   reinit(const Iterator &cell, unsigned int dofs_per_cell)
596 *   {
597 *   cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
598 *   cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
599 *  
600 *   cell_rhs.reinit(dofs_per_cell);
601 *   cell_mass_rhs.reinit(dofs_per_cell);
602 *  
603 *   local_dof_indices.resize(dofs_per_cell);
604 *   cell->get_dof_indices(local_dof_indices);
605 *   }
606 *   };
607 *  
608 *  
609 *  
610 * @endcode
611 *
612 *
613 * <a name="source/DG_advection_reaction.cc-Auxiliaryfunction"></a>
614 * <h3>Auxiliary function</h3>
615 * This auxiliary function is taken from @ref step_74 "step-74" and it's used to
616 * compute the jump of the finite element function @f$u_h@f$ on a face.
617 *
618 * @code
619 *   template <int dim>
620 *   void
621 *   get_function_jump(const FEInterfaceValues<dim> &fe_iv,
622 *   const Vector<double> &solution,
623 *   std::vector<double> &jump)
624 *   {
625 *   const unsigned int n_q = fe_iv.n_quadrature_points;
626 *   std::array<std::vector<double>, 2> face_values;
627 *   jump.resize(n_q);
628 *   for (unsigned int i = 0; i < 2; ++i)
629 *   {
630 *   face_values[i].resize(n_q);
631 *   fe_iv.get_fe_face_values(i).get_function_values(solution, face_values[i]);
632 *   }
633 *   for (unsigned int q = 0; q < n_q; ++q)
634 *   jump[q] = face_values[0][q] - face_values[1][q];
635 *   }
636 *  
637 *  
638 *  
639 *   template <int dim>
640 *   AdvectionReaction<dim>::AdvectionReaction()
641 *   : mapping()
642 *   , dof_handler(triangulation)
643 *   {
644 *   Assert(dim > 1, ExcMessage("Not implemented in 1D."));
645 *   add_parameter("Finite element degree", fe_degree);
646 *   add_parameter("Problem constants", constants);
647 *   add_parameter("Output filename", output_filename);
648 *   add_parameter("Use direct solver", use_direct_solver);
649 *   add_parameter("Number of refinement cycles", n_refinement_cycles);
650 *   add_parameter("Number of global refinement", n_global_refinements);
651 *   add_parameter("Refinement", refinement);
652 *   add_parameter("Exact solution expression", exact_solution_expression);
653 *   add_parameter("Boundary conditions expression",
654 *   boundary_conditions_expression);
655 *   add_parameter("Theta", theta);
656 *   add_parameter("Advection coefficient expression",
657 *   advection_coefficient_expression);
658 *   add_parameter("Right hand side expression", rhs_expression);
659 *  
660 *   this->prm.enter_subsection("Error table");
661 *   error_table.add_parameters(this->prm);
662 *   this->prm.leave_subsection();
663 *   }
664 *  
665 *  
666 *  
667 *   template <int dim>
668 *   void
669 *   AdvectionReaction<dim>::initialize_params(const std::string &filename)
670 *   {
671 *   ParameterAcceptor::initialize(filename,
672 *   "last_used_parameters.prm",
673 *   ParameterHandler::Short);
674 *   if (theta < 0.0 || theta > 10.0 || std::abs(theta) < 1e-12)
675 *   {
676 *   throw(
677 *   theta_exc("Theta parameter is not in a suitable range: see paper by "
678 *   "Brezzi, Marini, Suli for an extended discussion"));
679 *   }
680 *   }
681 *  
682 *  
683 *  
684 *   template <int dim>
685 *   void
686 *   AdvectionReaction<dim>::parse_string(const std::string &parameters)
687 *   {
688 *   ParameterAcceptor::prm.parse_input_from_string(parameters);
689 *   ParameterAcceptor::parse_all_parameters();
690 *   }
691 *  
692 *  
693 *  
694 *   template <int dim>
695 *   void
696 *   AdvectionReaction<dim>::setup_system()
697 *   {
698 * @endcode
699 *
700 * first need to distribute the DoFs.
701 *
702 * @code
703 *   if (!fe)
704 *   {
705 *   fe = std::make_unique<FE_DGQ<dim>>(fe_degree);
706 *   const auto vars = dim == 2 ? "x,y" : "x,y,z";
707 *   exact_solution.initialize(vars, exact_solution_expression, constants);
708 *   rhs.initialize(vars, rhs_expression, constants);
709 *   advection_coeff.initialize(vars,
710 *   advection_coefficient_expression,
711 *   constants);
712 *   boundary_conditions.initialize(vars,
713 *   boundary_conditions_expression,
714 *   constants);
715 *   }
716 *   dof_handler.distribute_dofs(*fe);
717 *  
718 * @endcode
719 *
720 * To build the sparsity pattern for DG discretizations, we can call the
721 * function analogue to DoFTools::make_sparsity_pattern, which is called
722 * DoFTools::make_flux_sparsity_pattern:
723 *
724 * @code
725 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
726 *   DoFTools::make_flux_sparsity_pattern(dof_handler,
727 *   dsp); // DG sparsity pattern generator
728 *   sparsity_pattern.copy_from(dsp);
729 *  
730 * @endcode
731 *
732 * Finally, we set up the structure of all components of the linear system.
733 *
734 * @code
735 *   system_matrix.reinit(sparsity_pattern);
736 *   solution.reinit(dof_handler.n_dofs());
737 *   right_hand_side.reinit(dof_handler.n_dofs());
738 *   }
739 *  
740 *  
741 *  
742 * @endcode
743 *
744 * in the call to MeshWorker::mesh_loop() we only need to specify what should
745 * happen on
746 * each cell, each boundary face, and each interior face. These three tasks
747 * are handled by the lambda functions inside the function below.
748 *
749
750 *
751 *
752 * @code
753 *   template <int dim>
754 *   void
755 *   AdvectionReaction<dim>::assemble_system()
756 *   {
757 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
758 *  
759 *   const QGauss<dim> quadrature = fe->tensor_degree() + 1;
760 *   const QGauss<dim - 1> quadrature_face = fe->tensor_degree() + 1;
761 *  
762 * @endcode
763 *
764 * This is the function that will be executed for each cell.
765 *
766 * @code
767 *   const auto cell_worker = [&](const Iterator &cell,
768 *   ScratchData<dim> &scratch_data,
769 *   CopyData &copy_data) {
770 *   FEValues<dim> fe_values_continuous(*fe,
771 *   quadrature,
772 *   update_values | update_gradients |
773 *   update_quadrature_points |
774 *   update_JxW_values);
775 *  
776 *   const unsigned int n_dofs =
777 *   scratch_data.fe_values.get_fe().n_dofs_per_cell();
778 *   copy_data.reinit(cell, n_dofs);
779 *   scratch_data.fe_values.reinit(cell);
780 *  
781 *   const auto &q_points = scratch_data.fe_values.get_quadrature_points();
782 *  
783 *   const FEValues<dim> &fe_v = scratch_data.fe_values;
784 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
785 *  
786 *   for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
787 *   {
788 *   auto beta_q = beta(q_points[point]);
789 *   for (unsigned int i = 0; i < n_dofs; ++i)
790 *   {
791 *   for (unsigned int j = 0; j < n_dofs; ++j)
792 *   {
793 *   copy_data.cell_matrix(i, j) +=
794 *   (-beta_q // -\beta
795 *   * fe_v.shape_grad(i, point) // \nabla \phi_i
796 *   * fe_v.shape_value(j, point) // \phi_j
797 *   + advection_coeff.value(q_points[point]) * // gamma
798 *   fe_v.shape_value(i, point) // phi_i
799 *   * fe_v.shape_value(j, point) // phi_j
800 *   ) *
801 *   JxW[point]; // dx
802 *   }
803 *   copy_data.cell_rhs(i) += rhs.value(q_points[point]) // f(x_q)
804 *   * fe_v.shape_value(i, point) // phi_i(x_q)
805 *   * JxW[point]; // dx
806 *   }
807 *   }
808 *   };
809 *  
810 * @endcode
811 *
812 * This is the function called for boundary faces and consists of a normal
813 * integration using FEFaceValues. New is the logic to decide if the term
814 * goes into the system matrix (outflow) or the right-hand side (inflow).
815 *
816 * @code
817 *   const auto boundary_worker = [&](const Iterator &cell,
818 *   const unsigned int &face_no,
819 *   ScratchData<dim> &scratch_data,
820 *   CopyData &copy_data) {
821 *   scratch_data.fe_interface_values.reinit(cell, face_no);
822 *   const FEFaceValuesBase<dim> &fe_face =
823 *   scratch_data.fe_interface_values.get_fe_face_values(0);
824 *  
825 *   const auto &q_points = fe_face.get_quadrature_points();
826 *  
827 *   const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
828 *   const std::vector<double> &JxW = fe_face.get_JxW_values();
829 *   const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
830 *  
831 *   std::vector<double> g(q_points.size());
832 *   exact_solution.value_list(q_points, g);
833 *  
834 *   for (unsigned int point = 0; point < q_points.size(); ++point)
835 *   {
836 *   const double beta_dot_n = beta(q_points[point]) * normals[point];
837 *  
838 *   if (beta_dot_n > 0)
839 *   {
840 *   for (unsigned int i = 0; i < n_facet_dofs; ++i)
841 *   for (unsigned int j = 0; j < n_facet_dofs; ++j)
842 *   copy_data.cell_matrix(i, j) +=
843 *   fe_face.shape_value(i,
844 *   point) // \phi_i
845 *   * fe_face.shape_value(j, point) // \phi_j
846 *   * beta_dot_n // \beta . n
847 *   * JxW[point]; // dx
848 *   }
849 *   else
850 *   for (unsigned int i = 0; i < n_facet_dofs; ++i)
851 *   copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
852 *   * g[point] // g*/
853 *   * beta_dot_n // \beta . n
854 *   * JxW[point]; // dx
855 *   }
856 *   };
857 *  
858 * @endcode
859 *
860 * This is the function called on interior faces. The arguments specify
861 * cells, face and subface indices (for adaptive refinement). We just pass
862 * them along to the reinit() function of FEInterfaceValues.
863 *
864 * @code
865 *   const auto face_worker = [&](const Iterator &cell,
866 *   const unsigned int &f,
867 *   const unsigned int &sf,
868 *   const Iterator &ncell,
869 *   const unsigned int &nf,
870 *   const unsigned int &nsf,
871 *   ScratchData<dim> &scratch_data,
872 *   CopyData &copy_data) {
873 *   FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
874 *   fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
875 *   const auto &q_points = fe_iv.get_quadrature_points();
876 *  
877 *   copy_data.face_data.emplace_back();
878 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
879 *  
880 *   const unsigned int n_dofs = fe_iv.n_current_interface_dofs();
881 *   copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
882 *  
883 *   copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
884 *  
885 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
886 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
887 *  
888 *   for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
889 *   {
890 *   const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint];
891 *   for (unsigned int i = 0; i < n_dofs; ++i)
892 *   {
893 *   for (unsigned int j = 0; j < n_dofs; ++j)
894 *   {
895 *   copy_data_face.cell_matrix(i, j) +=
896 *   (beta(q_points[qpoint]) * normals[qpoint] *
897 *   fe_iv.average_of_shape_values(j, qpoint) *
898 *   fe_iv.jump_in_shape_values(i, qpoint) +
899 *   theta * std::abs(beta_dot_n) *
900 *   fe_iv.jump_in_shape_values(j, qpoint) *
901 *   fe_iv.jump_in_shape_values(i, qpoint)) *
902 *   JxW[qpoint];
903 *   }
904 *   }
905 *   }
906 *   };
907 *  
908 * @endcode
909 *
910 * The following lambda function will handle copying the data from the
911 * cell and face assembly into the global matrix and right-hand side.
912 *
913
914 *
915 * While we would not need an AffineConstraints object, because there are
916 * no hanging node constraints in DG discretizations, we use an empty
917 * object here as this allows us to use its `copy_local_to_global`
918 * functionality.
919 *
920 * @code
921 *   const AffineConstraints<double> constraints;
922 *  
923 *   const auto copier = [&](const CopyData &c) {
924 *   constraints.distribute_local_to_global(c.cell_matrix,
925 *   c.cell_rhs,
926 *   c.local_dof_indices,
927 *   system_matrix,
928 *   right_hand_side);
929 *  
930 *   for (auto &cdf : c.face_data)
931 *   {
932 *   constraints.distribute_local_to_global(cdf.cell_matrix,
933 *   cdf.joint_dof_indices,
934 *   system_matrix);
935 *   }
936 *   };
937 *  
938 *   ScratchData<dim> scratch_data(mapping, *fe, quadrature, quadrature_face);
939 *   CopyData copy_data;
940 *  
941 * @endcode
942 *
943 * Here, we finally handle the assembly. We pass in ScratchData and
944 * CopyData objects, the lambda functions from above, an specify that we
945 * want to assemble interior faces once.
946 *
947 * @code
948 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
949 *   dof_handler.end(),
950 *   cell_worker,
951 *   copier,
952 *   scratch_data,
953 *   copy_data,
954 *   MeshWorker::assemble_own_cells |
955 *   MeshWorker::assemble_boundary_faces |
956 *   MeshWorker::assemble_own_interior_faces_once,
957 *   boundary_worker,
958 *   face_worker);
959 *   }
960 *  
961 *  
962 *  
963 *   template <int dim>
964 *   void
965 *   AdvectionReaction<dim>::solve()
966 *   {
967 *   if (use_direct_solver)
968 *   {
969 *   SparseDirectUMFPACK system_matrix_inverse;
970 *   system_matrix_inverse.initialize(system_matrix);
971 *   system_matrix_inverse.vmult(solution, right_hand_side);
972 *   }
973 *   else
974 *   {
975 * @endcode
976 *
977 * Here we have a classic iterative solver, as done in many tutorials:
978 *
979 * @code
980 *   SolverControl solver_control(1000, 1e-15);
981 *   SolverRichardson<Vector<double>> solver(solver_control);
982 *   PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
983 *   preconditioner.initialize(system_matrix, fe->n_dofs_per_cell());
984 *   solver.solve(system_matrix, solution, right_hand_side, preconditioner);
985 *   std::cout << " Solver converged in " << solver_control.last_step()
986 *   << " iterations." << std::endl;
987 *   }
988 *   }
989 *  
990 *  
991 *  
992 * @endcode
993 *
994 *
995 * <a name="source/DG_advection_reaction.cc-Meshrefinement"></a>
996 * <h3>Mesh refinement</h3>
997 * We refine the grid according the proposed estimator or with an approximation
998 * to the gradient of the solution. The first option is the default one (you can
999 * see it in the header file)
1000 *
1001 * @code
1002 *   template <int dim>
1003 *   void
1004 *   AdvectionReaction<dim>::refine_grid()
1005 *   {
1006 *   if (refinement == "residual")
1007 *   {
1008 * @endcode
1009 *
1010 * If the `refinement` string is `"residual"`, then we first compute the
1011 * local projection
1012 *
1013 * @code
1014 *   compute_local_projection_and_estimate();
1015 * @endcode
1016 *
1017 * We then set the refinement fraction and as usual we execute the
1018 * refinement.
1019 *
1020 * @code
1021 *   const double refinement_fraction = 0.6;
1022 *   GridRefinement::refine_and_coarsen_fixed_fraction(
1023 *   triangulation, error_indicator_per_cell, refinement_fraction, 0.0);
1024 *   triangulation.execute_coarsening_and_refinement();
1025 *   }
1026 *   else if (refinement == "gradient")
1027 *   {
1028 *   Vector<float> gradient_indicator(triangulation.n_active_cells());
1029 *  
1030 * @endcode
1031 *
1032 * Now the approximate gradients are computed
1033 *
1034 * @code
1035 *   DerivativeApproximation::approximate_gradient(mapping,
1036 *   dof_handler,
1037 *   solution,
1038 *   gradient_indicator);
1039 *  
1040 * @endcode
1041 *
1042 * and they are cell-wise scaled by the factor @f$h^{1+d/2}@f$
1043 *
1044 * @code
1045 *   unsigned int cell_no = 0;
1046 *   for (const auto &cell : dof_handler.active_cell_iterators())
1047 *   gradient_indicator(cell_no++) *=
1048 *   std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
1049 *  
1050 * @endcode
1051 *
1052 * Finally they serve as refinement indicator.
1053 *
1054 * @code
1055 *   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
1056 *   gradient_indicator,
1057 *   0.25,
1058 *   0.0);
1059 *  
1060 *   triangulation.execute_coarsening_and_refinement();
1061 *   std::cout << gradient_indicator.l2_norm() << '\n';
1062 *   }
1063 *   else if (refinement == "global")
1064 *   {
1065 *   triangulation.refine_global(
1066 *   1); // just for testing on uniformly refined meshes
1067 *   }
1068 *   else
1069 *   {
1070 *   Assert(false, ExcInternalError());
1071 *   }
1072 *   }
1073 *  
1074 *  
1075 *  
1076 * @endcode
1077 *
1078 * The output of this program consists of a vtk file of the adaptively
1079 * refined grids and the numerical solutions.
1080 *
1081 * @code
1082 *   template <int dim>
1083 *   void
1084 *   AdvectionReaction<dim>::output_results(const unsigned int cycle) const
1085 *   {
1086 *   const std::string filename = "solution-" + std::to_string(cycle) + ".vtk";
1087 *   std::cout << " Writing solution to <" << filename << ">" << std::endl;
1088 *   std::ofstream output(filename);
1089 *  
1090 *   DataOut<dim> data_out;
1091 *   data_out.attach_dof_handler(dof_handler);
1092 *   data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
1093 *   data_out.build_patches(mapping);
1094 *   data_out.write_vtk(output);
1095 *   }
1096 *  
1097 *   template <int dim>
1098 *   void
1099 *   AdvectionReaction<dim>::compute_error()
1100 *   {
1101 *   error_table.error_from_exact(
1102 *   mapping,
1103 *   dof_handler,
1104 *   solution,
1105 *   exact_solution); // be careful: a FD approximation of the gradient is used
1106 * @endcode
1107 *
1108 * to compute the H^1 norm if Solution<dim> doesn't
1109 * implements the Gradient function
1110 *
1111 * @code
1112 *   }
1113 *  
1114 *  
1115 *  
1116 * @endcode
1117 *
1118 *
1119 * <a name="source/DG_advection_reaction.cc-Computetheenergynorm"></a>
1120 * <h3>Compute the energy norm</h3>
1121 * The energy norm is defined as @f$ |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 +
1122 * \sum_{F \in \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2
1123 * \Bigr)^{\frac{1}{2}}@f$ Notice that in the current case we have @f$c_f = \frac{|b
1124 * \cdot n|}{2}@f$ Like in the assembly, all the contributions are handled
1125 * separately by using ScratchData and CopyData objects.
1126 *
1127 * @code
1128 *   template <int dim>
1129 *   double
1130 *   AdvectionReaction<dim>::compute_energy_norm()
1131 *   {
1132 *   energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
1133 *  
1134 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
1135 *  
1136 * @endcode
1137 *
1138 * We start off by adding cell contributions
1139 *
1140 * @code
1141 *   const auto cell_worker = [&](const Iterator &cell,
1142 *   ScratchData<dim> &scratch_data,
1143 *   CopyData &copy_data) {
1144 *   const unsigned int n_dofs =
1145 *   scratch_data.fe_values.get_fe().n_dofs_per_cell();
1146 *   copy_data.reinit(cell, n_dofs);
1147 *   scratch_data.fe_values.reinit(cell);
1148 *  
1149 *   copy_data.cell_index = cell->active_cell_index();
1150 *  
1151 *   const auto &q_points = scratch_data.fe_values.get_quadrature_points();
1152 *   const FEValues<dim> &fe_v = scratch_data.fe_values;
1153 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
1154 *  
1155 *   double error_square_norm{0.0};
1156 *   std::vector<double> sol_u(fe_v.n_quadrature_points);
1157 *   fe_v.get_function_values(solution, sol_u);
1158 *  
1159 *   for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
1160 *   {
1161 *   const double diff =
1162 *   (sol_u[point] - exact_solution.value(q_points[point]));
1163 *   error_square_norm += diff * diff * JxW[point];
1164 *   }
1165 *   copy_data.value = error_square_norm;
1166 *   };
1167 *  
1168 * @endcode
1169 *
1170 * Here we add contributions coming from the internal faces
1171 *
1172 * @code
1173 *   const auto face_worker = [&](const Iterator &cell,
1174 *   const unsigned int &f,
1175 *   const unsigned int &sf,
1176 *   const Iterator &ncell,
1177 *   const unsigned int &nf,
1178 *   const unsigned int &nsf,
1179 *   ScratchData<dim> &scratch_data,
1180 *   CopyData &copy_data) {
1181 *   FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
1182 *   fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
1183 *  
1184 *   copy_data.face_data.emplace_back();
1185 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
1186 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
1187 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
1188 *  
1189 *   const auto &q_points = fe_iv.get_quadrature_points();
1190 *   const unsigned n_q_points = q_points.size();
1191 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
1192 *   std::vector<double> g(n_q_points);
1193 *  
1194 *   std::vector<double> jump(n_q_points);
1195 *   get_function_jump(fe_iv, solution, jump);
1196 *  
1197 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
1198 *  
1199 *   double error_jump_square{0.0};
1200 *   for (unsigned int point = 0; point < n_q_points; ++point)
1201 *   {
1202 *   const double beta_dot_n =
1203 *   theta * std::abs(beta(q_points[point]) * normals[point]);
1204 *   error_jump_square +=
1205 *   beta_dot_n * jump[point] * jump[point] * JxW[point];
1206 *   }
1207 *  
1208 *   copy_data.value = error_jump_square;
1209 *   };
1210 *  
1211 * @endcode
1212 *
1213 * Finally, we add the boundary contributions
1214 *
1215 * @code
1216 *   const auto boundary_worker = [&](const Iterator &cell,
1217 *   const unsigned int &face_no,
1218 *   ScratchData<dim> &scratch_data,
1219 *   CopyData &copy_data) {
1220 *   scratch_data.fe_interface_values.reinit(cell, face_no);
1221 *   const FEFaceValuesBase<dim> &fe_fv =
1222 *   scratch_data.fe_interface_values.get_fe_face_values(0);
1223 *   const auto &q_points = fe_fv.get_quadrature_points();
1224 *   const unsigned n_q_points = q_points.size();
1225 *   const std::vector<double> &JxW = fe_fv.get_JxW_values();
1226 *  
1227 *   std::vector<double> g(n_q_points);
1228 *  
1229 *   std::vector<double> sol_u(n_q_points);
1230 *   fe_fv.get_function_values(solution, sol_u);
1231 *  
1232 *   const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
1233 *  
1234 *   double difference_norm_square = 0.;
1235 *   for (unsigned int point = 0; point < q_points.size(); ++point)
1236 *   {
1237 *   const double beta_dot_n =
1238 *   theta * std::abs(beta(q_points[point]) * normals[point]);
1239 *   const double diff =
1240 *   (boundary_conditions.value(q_points[point]) - sol_u[point]);
1241 *   difference_norm_square += beta_dot_n * diff * diff * JxW[point];
1242 *   }
1243 *   copy_data.value = difference_norm_square;
1244 *   };
1245 *  
1246 *   const auto copier = [&](const auto &copy_data) {
1247 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
1248 *   {
1249 *   energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
1250 *   }
1251 *   for (auto &cdf : copy_data.face_data)
1252 *   for (unsigned int j = 0; j < 2; ++j)
1253 *   energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
1254 *   };
1255 *  
1256 *   ScratchData<dim> scratch_data(mapping,
1257 *   *fe,
1258 *   QGauss<dim>{fe->tensor_degree() + 1},
1259 *   QGauss<dim - 1>{fe->tensor_degree() + 1});
1260 *  
1261 *   CopyData copy_data;
1262 *  
1263 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1264 *   dof_handler.end(),
1265 *   cell_worker,
1266 *   copier,
1267 *   scratch_data,
1268 *   copy_data,
1272 *   boundary_worker,
1273 *   face_worker);
1274 *  
1275 *   const double energy_error = std::sqrt(energy_norm_square_per_cell.l1_norm());
1276 *   return energy_error;
1277 *   }
1278 *  
1279 *  
1280 *  
1281 * @endcode
1282 *
1283 *
1284 * <a name="source/DG_advection_reaction.cc-Computingtheestimator"></a>
1285 * <h3>Computing the estimator</h3>
1286 * In the estimator, we have to compute the term @f$||f- c u_h - \Pi(f- c
1287 * u_h)||_{T}^{2}@f$ over a generic cell @f$T@f$. To achieve this, we first need to
1288 * compute the projection involving the finite element function @f$u_h@f$. Using the
1289 * definition of orthogonal projection, we're required to solve cellwise
1290 * @f$(v_h,f-c u_h)_T = (v_h,\Pi)_T \qquad \forall v_h \in V_h@f$ for @f$\Pi@f$, which
1291 * means that we have to build a mass-matrix on each cell. Once we have the
1292 * projection, which is a finite element function, we can add its contribution
1293 * in the <code>cell_worker</code> lambda. As done in @ref step_74 "step-74", the square of the
1294 * error indicator is computed.
1295 *
1296
1297 *
1298 *
1299 * @code
1300 *   template <int dim>
1301 *   void
1302 *   AdvectionReaction<dim>::compute_local_projection_and_estimate()
1303 *   {
1304 * @endcode
1305 *
1306 * Compute the term @f$||f-c u_h - \Pi(f- cu_h)||_T^2@f$
1307 *
1308 * @code
1309 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
1310 *   error_indicator_per_cell.reinit(triangulation.n_active_cells());
1311 *  
1312 *   const auto cell_worker = [&](const Iterator &cell,
1313 *   ScratchData<dim> &scratch_data,
1314 *   CopyData &copy_data) {
1315 *   const unsigned int n_dofs =
1316 *   scratch_data.fe_values.get_fe().n_dofs_per_cell();
1317 *  
1318 *   copy_data.reinit(cell, n_dofs);
1319 *   scratch_data.fe_values.reinit(cell);
1320 *   copy_data.cell_index = cell->active_cell_index();
1321 *  
1322 *   const auto &q_points = scratch_data.fe_values.get_quadrature_points();
1323 *   const unsigned n_q_points = q_points.size();
1324 *  
1325 *   const FEValues<dim> &fe_v = scratch_data.fe_values;
1326 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
1327 *  
1328 *   std::vector<double> sol_u_at_quadrature_points(fe_v.n_quadrature_points);
1329 *   fe_v.get_function_values(solution, sol_u_at_quadrature_points);
1330 *  
1331 * @endcode
1332 *
1333 * Compute local L^2 projection of @f$f- c u_h@f$ over the local finite element
1334 * space
1335 *
1336 * @code
1337 *   for (unsigned int point = 0; point < n_q_points; ++point)
1338 *   {
1339 *   for (unsigned int i = 0; i < n_dofs; ++i)
1340 *   {
1341 *   for (unsigned int j = 0; j < n_dofs; ++j)
1342 *   {
1343 *   copy_data.cell_mass_matrix(i, j) +=
1344 *   fe_v.shape_value(i, point) * // phi_i(x_q)
1345 *   fe_v.shape_value(j, point) * // phi_j(x_q)
1346 *   JxW[point]; // dx(x_q)
1347 *   }
1348 *   copy_data.cell_mass_rhs(i) +=
1349 *   (rhs.value(q_points[point]) * // f(x_q)
1350 *   fe_v.shape_value(i, point) // phi_i(x_q)
1351 *   - advection_coeff.value(q_points[point]) *
1352 *   fe_v.shape_value(i, point) * // c*phi_i(x_q)
1353 *   sol_u_at_quadrature_points[point]) * // u_h(x_q)
1354 *   JxW[point]; // dx
1355 *   }
1356 *   }
1357 *  
1358 *   FullMatrix<double> inverse(fe_v.n_quadrature_points,
1359 *   fe_v.n_quadrature_points);
1360 *   inverse.invert(copy_data.cell_mass_matrix);
1361 *   Vector<double> proj(fe_v.n_quadrature_points); // projection of (f-c*U_h) on
1362 * @endcode
1363 *
1364 * the local fe_space
1365 *
1366 * @code
1367 *   inverse.vmult(proj, copy_data.cell_mass_rhs); // M^{-1}*rhs = proj
1368 *  
1369 *   double square_norm_over_cell = 0.0;
1370 *   for (unsigned int point = 0; point < n_q_points; ++point)
1371 *   {
1372 *   const double diff = rhs.value(q_points[point]) -
1373 *   sol_u_at_quadrature_points[point] - proj[point];
1374 *   square_norm_over_cell += diff * diff * JxW[point];
1375 *   }
1376 *   copy_data.value_estimator = square_norm_over_cell;
1377 *   };
1378 *  
1379 * @endcode
1380 *
1381 * Finally we have the boundary term with @f$||\beta (g-u_h^+)||^2@f$
1382 *
1383 * @code
1384 *   const auto boundary_worker = [&](const Iterator &cell,
1385 *   const unsigned int &face_no,
1386 *   ScratchData<dim> &scratch_data,
1387 *   CopyData &copy_data) {
1388 *   scratch_data.fe_interface_values.reinit(cell, face_no);
1389 *   const FEFaceValuesBase<dim> &fe_fv =
1390 *   scratch_data.fe_interface_values.get_fe_face_values(0);
1391 *   const auto &q_points = fe_fv.get_quadrature_points();
1392 *   const unsigned n_q_points = q_points.size();
1393 *   const std::vector<double> &JxW = fe_fv.get_JxW_values();
1394 *  
1395 *   std::vector<double> g(n_q_points);
1396 *   exact_solution.value_list(q_points, g);
1397 *  
1398 *   std::vector<double> sol_u(n_q_points);
1399 *   fe_fv.get_function_values(solution, sol_u);
1400 *  
1401 *   const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
1402 *  
1403 *   double square_norm_over_bdary_face = 0.;
1404 *   for (unsigned int point = 0; point < q_points.size(); ++point)
1405 *   {
1406 *   const double beta_dot_n = beta(q_points[point]) * normals[point];
1407 *  
1408 *   if (beta_dot_n < 0) //\partial_{-T} \cap \partial_{- \Omega}
1409 *   {
1410 *   const double diff =
1411 *   std::abs(beta_dot_n) * (g[point] - sol_u[point]);
1412 *   square_norm_over_bdary_face += diff * diff * JxW[point];
1413 *   }
1414 *   }
1415 *   copy_data.value_estimator += square_norm_over_bdary_face;
1416 *   };
1417 *  
1418 * @endcode
1419 *
1420 * Then compute the interior face terms with @f$|| \sqrt{b \cdot n}[u_h]||^2@f$
1421 *
1422 * @code
1423 *   const auto face_worker = [&](const Iterator &cell,
1424 *   const unsigned int &f,
1425 *   const unsigned int &sf,
1426 *   const Iterator &ncell,
1427 *   const unsigned int &nf,
1428 *   const unsigned int &nsf,
1429 *   ScratchData<dim> &scratch_data,
1430 *   CopyData &copy_data) {
1431 *   FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
1432 *   fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
1433 *  
1434 *   copy_data.face_data.emplace_back();
1435 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
1436 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
1437 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
1438 *  
1439 *   const auto &q_points = fe_iv.get_quadrature_points();
1440 *   const unsigned n_q_points = q_points.size();
1441 *  
1442 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
1443 *   std::vector<double> g(n_q_points);
1444 *  
1445 *   std::vector<double> jump(n_q_points);
1446 *   get_function_jump(fe_iv, solution, jump);
1447 *  
1448 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
1449 *  
1450 *   double error_jump_square{0.0};
1451 *   for (unsigned int point = 0; point < n_q_points; ++point)
1452 *   {
1453 *   const double beta_dot_n = beta(q_points[point]) * normals[point];
1454 *   if (beta_dot_n < 0)
1455 *   {
1456 *   error_jump_square +=
1457 *   std::abs(beta_dot_n) * jump[point] * jump[point] * JxW[point];
1458 *   }
1459 *   }
1460 *  
1461 *   copy_data_face.values[0] = error_jump_square;
1462 *   copy_data_face.values[1] = copy_data_face.values[0];
1463 *   };
1464 *  
1465 *   ScratchData<dim> scratch_data(mapping,
1466 *   *fe,
1467 *   QGauss<dim>{fe->tensor_degree() + 1},
1468 *   QGauss<dim - 1>{fe->tensor_degree() + 1});
1469 *  
1470 *   const auto copier = [&](const auto &copy_data) {
1471 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
1472 *   {
1473 *   error_indicator_per_cell[copy_data.cell_index] +=
1474 *   copy_data.value_estimator;
1475 *   }
1476 *   for (auto &cdf : copy_data.face_data)
1477 *   {
1478 *   for (unsigned int j = 0; j < 2; ++j)
1479 *   {
1480 *   error_indicator_per_cell[cdf.cell_indices[j]] += cdf.values[j];
1481 *   }
1482 *   }
1483 *   };
1484 *  
1485 * @endcode
1486 *
1487 * Here, we finally handle the assembly of the Mass matrix (M)_{ij} = (\phi_j,
1488 * \phi_i)_T. We pass in ScratchData and CopyData objects
1489 *
1490 * @code
1491 *   CopyData copy_data;
1492 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1493 *   dof_handler.end(),
1494 *   cell_worker,
1495 *   copier,
1496 *   scratch_data,
1497 *   copy_data,
1498 *   MeshWorker::assemble_own_cells |
1499 *   MeshWorker::assemble_boundary_faces |
1500 *   MeshWorker::assemble_own_interior_faces_once,
1501 *   boundary_worker,
1502 *   face_worker);
1503 *   }
1504 *  
1505 *  
1506 *  
1507 * @endcode
1508 *
1509 * Usual <code>run</code> function, which runs over several refinement cycles
1510 *
1511 * @code
1512 *   template <int dim>
1513 *   void
1514 *   AdvectionReaction<dim>::run()
1515 *   {
1516 *   std::vector<double> energy_errors;
1517 *   std::vector<int> dofs_hist;
1518 *   for (unsigned int cycle = 0; cycle < n_refinement_cycles; ++cycle)
1519 *   {
1520 *   std::cout << "Cycle " << cycle << std::endl;
1521 *  
1522 *   if (cycle == 0)
1523 *   {
1524 *   GridGenerator::hyper_cube(triangulation);
1525 *   triangulation.refine_global(n_global_refinements);
1526 *   }
1527 *   else
1528 *   {
1529 *   refine_grid();
1530 *   }
1531 *   std::cout << " Number of active cells: "
1532 *   << triangulation.n_active_cells() << std::endl;
1533 *   std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1534 *   << std::endl;
1535 *  
1536 *   setup_system();
1537 *   assemble_system();
1538 *   solve();
1539 *   compute_error();
1540 *   output_results(cycle);
1541 *  
1542 *   energy_errors.emplace_back(compute_energy_norm());
1543 *   dofs_hist.emplace_back(triangulation.n_active_cells());
1544 *   }
1545 *   error_table.output_table(std::cout);
1546 *  
1547 *   for (unsigned int i = 0; i < n_refinement_cycles; ++i)
1548 *   std::cout << "Cycle " << i << "\t" << energy_errors[i] << std::endl;
1549 *   }
1550 * @endcode
1551 *
1552 * Explicit instantiation
1553 *
1554 * @code
1555 *   template class AdvectionReaction<1>;
1556 *   template class AdvectionReaction<2>;
1557 *   template class AdvectionReaction<3>;
1558 * @endcode
1559
1560
1561*/
void reinit(const CellIteratorType &cell, const unsigned int face_no, const unsigned int sub_face_no, const CellNeighborIteratorType &cell_neighbor, const unsigned int face_no_neighbor, const unsigned int sub_face_no_neighbor, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
const std::vector< double > & get_JxW_values() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:349
Abstract base class for mapping classes.
Definition mapping.h:318
The ParsedConvergenceTable class.
Definition point.h:111
Point< 2 > first
Definition grid_out.cc:4623
unsigned int cell_index
#define Assert(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
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
spacedim & mapping
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
const std::vector< bool > & used
const Triangulation< dim, spacedim > & tria
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
@ assemble_boundary_faces
@ assemble_own_interior_faces_once
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
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)
void copy(const T *begin, const T *end, U *dest)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static const unsigned int invalid_unsigned_int
Definition types.h:220
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation