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