Reference documentation for deal.II version 9.2.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\}}\)
step-63.h
Go to the documentation of this file.
1 
523  * ParameterHandler prm;
524  *
525  * prm.declare_entry("Epsilon",
526  * "0.005",
527  * Patterns::Double(0),
528  * "Diffusion parameter");
529  *
530  * prm.declare_entry("Fe degree",
531  * "1",
532  * Patterns::Integer(1),
533  * "Finite Element degree");
534  * prm.declare_entry("Smoother type",
535  * "block SOR",
536  * Patterns::Selection("SOR|Jacobi|block SOR|block Jacobi"),
537  * "Select smoother: SOR|Jacobi|block SOR|block Jacobi");
538  * prm.declare_entry("Smoothing steps",
539  * "2",
540  * Patterns::Integer(1),
541  * "Number of smoothing steps");
542  * prm.declare_entry(
543  * "DoF renumbering",
544  * "downstream",
545  * Patterns::Selection("none|downstream|upstream|random"),
546  * "Select DoF renumbering: none|downstream|upstream|random");
547  * prm.declare_entry("With streamline diffusion",
548  * "true",
549  * Patterns::Bool(),
550  * "Enable streamline diffusion stabilization: true|false");
551  * prm.declare_entry("Output",
552  * "true",
553  * Patterns::Bool(),
554  * "Generate graphical output: true|false");
555  *
556  * /* ...and then try to read their values from the input file: */
557  * if (prm_filename.empty())
558  * {
559  * prm.print_parameters(std::cout, ParameterHandler::Text);
560  * AssertThrow(
561  * false, ExcMessage("Please pass a .prm file as the first argument!"));
562  * }
563  *
564  * prm.parse_input(prm_filename);
565  *
566  * epsilon = prm.get_double("Epsilon");
567  * fe_degree = prm.get_integer("Fe degree");
568  * smoother_type = prm.get("Smoother type");
569  * smoothing_steps = prm.get_integer("Smoothing steps");
570  *
571  * const std::string renumbering = prm.get("DoF renumbering");
572  * if (renumbering == "none")
573  * dof_renumbering = DoFRenumberingStrategy::none;
574  * else if (renumbering == "downstream")
575  * dof_renumbering = DoFRenumberingStrategy::downstream;
576  * else if (renumbering == "upstream")
577  * dof_renumbering = DoFRenumberingStrategy::upstream;
578  * else if (renumbering == "random")
579  * dof_renumbering = DoFRenumberingStrategy::random;
580  * else
581  * AssertThrow(false,
582  * ExcMessage("The <DoF renumbering> parameter has "
583  * "an invalid value."));
584  *
585  * with_streamline_diffusion = prm.get_bool("With streamline diffusion");
586  * output = prm.get_bool("Output");
587  * }
588  *
589  *
590  * @endcode
591  *
592  *
593  * <a name="Cellpermutations"></a>
594  * <h3>Cell permutations</h3>
595  *
596 
597  *
598  * The ordering in which cells and degrees of freedom are traversed
599  * will play a role in the speed of convergence for multiplicative
600  * methods. Here we define functions which return a specific ordering
601  * of cells to be used by the block smoothers.
602  *
603 
604  *
605  * For each type of cell ordering, we define a function for the
606  * active mesh and one for a level mesh (i.e., for the cells at one
607  * level of a multigrid hierarchy). While the only reordering
608  * necessary for solving the system will be on the level meshes, we
609  * include the active reordering for visualization purposes in
610  * output_results().
611  *
612 
613  *
614  * For the two downstream ordering functions, we first create an
615  * array with all of the relevant cells that we then sort in
616  * downstream direction using a "comparator" object. The output of
617  * the functions is then simply an array of the indices of the cells
618  * in the just computed order.
619  *
620  * @code
621  * template <int dim>
622  * std::vector<unsigned int>
623  * create_downstream_cell_ordering(const DoFHandler<dim> &dof_handler,
624  * const Tensor<1, dim> direction,
625  * const unsigned int level)
626  * {
627  * std::vector<typename DoFHandler<dim>::level_cell_iterator> ordered_cells;
628  * ordered_cells.reserve(dof_handler.get_triangulation().n_cells(level));
629  * for (const auto &cell : dof_handler.cell_iterators_on_level(level))
630  * ordered_cells.push_back(cell);
631  *
632  * const DoFRenumbering::
633  * CompareDownstream<typename DoFHandler<dim>::level_cell_iterator, dim>
634  * comparator(direction);
635  * std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
636  *
637  * std::vector<unsigned> ordered_indices;
638  * ordered_indices.reserve(dof_handler.get_triangulation().n_cells(level));
639  *
640  * for (const auto &cell : ordered_cells)
641  * ordered_indices.push_back(cell->index());
642  *
643  * return ordered_indices;
644  * }
645  *
646  *
647  *
648  * template <int dim>
649  * std::vector<unsigned int>
650  * create_downstream_cell_ordering(const DoFHandler<dim> &dof_handler,
651  * const Tensor<1, dim> direction)
652  * {
653  * std::vector<typename DoFHandler<dim>::active_cell_iterator> ordered_cells;
654  * ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
655  * for (const auto &cell : dof_handler.active_cell_iterators())
656  * ordered_cells.push_back(cell);
657  *
658  * const DoFRenumbering::
659  * CompareDownstream<typename DoFHandler<dim>::active_cell_iterator, dim>
660  * comparator(direction);
661  * std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
662  *
663  * std::vector<unsigned int> ordered_indices;
664  * ordered_indices.reserve(dof_handler.get_triangulation().n_active_cells());
665  *
666  * for (const auto &cell : ordered_cells)
667  * ordered_indices.push_back(cell->index());
668  *
669  * return ordered_indices;
670  * }
671  *
672  *
673  * @endcode
674  *
675  * The functions that produce a random ordering are similar in
676  * spirit in that they first put information about all cells into an
677  * array. But then, instead of sorting them, they shuffle the
678  * elements randomly using the facilities C++ offers to generate
679  * random numbers. The way this is done is by iterating over all
680  * elements of the array, drawing a random number for another
681  * element before that, and then exchanging these elements. The
682  * result is a random shuffle of the elements of the array.
683  *
684  * @code
685  * template <int dim>
686  * std::vector<unsigned int>
687  * create_random_cell_ordering(const DoFHandler<dim> &dof_handler,
688  * const unsigned int level)
689  * {
690  * std::vector<unsigned int> ordered_cells;
691  * ordered_cells.reserve(dof_handler.get_triangulation().n_cells(level));
692  * for (const auto &cell : dof_handler.cell_iterators_on_level(level))
693  * ordered_cells.push_back(cell->index());
694  *
695  * std::mt19937 random_number_generator;
696  * std::shuffle(ordered_cells.begin(),
697  * ordered_cells.end(),
698  * random_number_generator);
699  *
700  * return ordered_cells;
701  * }
702  *
703  *
704  *
705  * template <int dim>
706  * std::vector<unsigned int>
707  * create_random_cell_ordering(const DoFHandler<dim> &dof_handler)
708  * {
709  * std::vector<unsigned int> ordered_cells;
710  * ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
711  * for (const auto &cell : dof_handler.active_cell_iterators())
712  * ordered_cells.push_back(cell->index());
713  *
714  * std::mt19937 random_number_generator;
715  * std::shuffle(ordered_cells.begin(),
716  * ordered_cells.end(),
717  * random_number_generator);
718  *
719  * return ordered_cells;
720  * }
721  *
722  *
723  * @endcode
724  *
725  *
726  * <a name="Righthandsideandboundaryvalues"></a>
727  * <h3>Right-hand side and boundary values</h3>
728  *
729 
730  *
731  * The problem solved in this tutorial is an adaptation of Ex. 3.1.3 found
732  * on pg. 118 of <a
733  * href="https://global.oup.com/academic/product/finite-elements-and-fast-iterative-solvers-9780199678808">
734  * Finite Elements and Fast Iterative Solvers: with Applications in
735  * Incompressible Fluid Dynamics by Elman, Silvester, and Wathen</a>. The
736  * main difference being that we add a hole in the center of our domain with
737  * zero Dirichlet boundary conditions.
738  *
739 
740  *
741  * For a complete description, we need classes that implement the
742  * zero right-hand side first (we could of course have just used
744  *
745  * @code
746  * template <int dim>
747  * class RightHandSide : public Function<dim>
748  * {
749  * public:
750  * virtual double value(const Point<dim> & p,
751  * const unsigned int component = 0) const override;
752  *
753  * virtual void value_list(const std::vector<Point<dim>> &points,
754  * std::vector<double> & values,
755  * const unsigned int component = 0) const override;
756  * };
757  *
758  *
759  *
760  * template <int dim>
761  * double RightHandSide<dim>::value(const Point<dim> &,
762  * const unsigned int component) const
763  * {
764  * Assert(component == 0, ExcIndexRange(component, 0, 1));
765  * (void)component;
766  *
767  * return 0.0;
768  * }
769  *
770  *
771  *
772  * template <int dim>
773  * void RightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
774  * std::vector<double> & values,
775  * const unsigned int component) const
776  * {
777  * Assert(values.size() == points.size(),
778  * ExcDimensionMismatch(values.size(), points.size()));
779  *
780  * for (unsigned int i = 0; i < points.size(); ++i)
781  * values[i] = RightHandSide<dim>::value(points[i], component);
782  * }
783  *
784  *
785  * @endcode
786  *
787  * We also have Dirichlet boundary conditions. On a connected portion of the
788  * outer, square boundary we set the value to 1, and we set the value to 0
789  * everywhere else (including the inner, circular boundary):
790  *
791  * @code
792  * template <int dim>
793  * class BoundaryValues : public Function<dim>
794  * {
795  * public:
796  * virtual double value(const Point<dim> & p,
797  * const unsigned int component = 0) const override;
798  *
799  * virtual void value_list(const std::vector<Point<dim>> &points,
800  * std::vector<double> & values,
801  * const unsigned int component = 0) const override;
802  * };
803  *
804  *
805  *
806  * template <int dim>
807  * double BoundaryValues<dim>::value(const Point<dim> & p,
808  * const unsigned int component) const
809  * {
810  * Assert(component == 0, ExcIndexRange(component, 0, 1));
811  * (void)component;
812  *
813  * @endcode
814  *
815  * Set boundary to 1 if @f$x=1@f$, or if @f$x>0.5@f$ and @f$y=-1@f$.
816  *
817  * @code
818  * if (std::fabs(p[0] - 1) < 1e-8 ||
819  * (std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5))
820  * {
821  * return 1.0;
822  * }
823  * else
824  * {
825  * return 0.0;
826  * }
827  * }
828  *
829  *
830  *
831  * template <int dim>
832  * void BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
833  * std::vector<double> & values,
834  * const unsigned int component) const
835  * {
836  * Assert(values.size() == points.size(),
837  * ExcDimensionMismatch(values.size(), points.size()));
838  *
839  * for (unsigned int i = 0; i < points.size(); ++i)
840  * values[i] = BoundaryValues<dim>::value(points[i], component);
841  * }
842  *
843  *
844  *
845  * @endcode
846  *
847  *
848  * <a name="Streamlinediffusionimplementation"></a>
849  * <h3>Streamline diffusion implementation</h3>
850  *
851 
852  *
853  * The streamline diffusion method has a stabilization constant that
854  * we need to be able to compute. The choice of how this parameter
855  * is computed is taken from <a
856  * href="https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27">On
857  * Discontinuity-Capturing Methods for Convection-Diffusion
858  * Equations by Volker John and Petr Knobloch</a>.
859  *
860  * @code
861  * template <int dim>
862  * double compute_stabilization_delta(const double hk,
863  * const double eps,
864  * const Tensor<1, dim> dir,
865  * const double pk)
866  * {
867  * const double Peclet = dir.norm() * hk / (2.0 * eps * pk);
868  * const double coth =
869  * (1.0 + std::exp(-2.0 * Peclet)) / (1.0 - std::exp(-2.0 * Peclet));
870  *
871  * return hk / (2.0 * dir.norm() * pk) * (coth - 1.0 / Peclet);
872  * }
873  *
874  *
875  * @endcode
876  *
877  *
878  * <a name="codeAdvectionProlemcodeclass"></a>
879  * <h3><code>AdvectionProlem</code> class</h3>
880  *
881 
882  *
883  * This is the main class of the program, and should look very similar to
884  * @ref step_16 "step-16". The major difference is that, since we are defining our multigrid
885  * smoother at runtime, we choose to define a function `create_smoother()` and
886  * a class object `mg_smoother` which is a `std::unique_ptr` to a smoother
887  * that is derived from MGSmoother. Note that for smoothers derived from
888  * RelaxationBlock, we must include a `smoother_data` object for each level.
889  * This will contain information about the cell ordering and the method of
890  * inverting cell matrices.
891  *
892 
893  *
894  *
895  * @code
896  * template <int dim>
897  * class AdvectionProblem
898  * {
899  * public:
900  * AdvectionProblem(const Settings &settings);
901  * void run();
902  *
903  * private:
904  * void setup_system();
905  *
906  * template <class IteratorType>
907  * void assemble_cell(const IteratorType &cell,
908  * ScratchData<dim> & scratch_data,
909  * CopyData & copy_data);
910  * void assemble_system_and_multigrid();
911  *
912  * void setup_smoother();
913  *
914  * void solve();
915  * void refine_grid();
916  * void output_results(const unsigned int cycle) const;
917  *
919  * DoFHandler<dim> dof_handler;
920  *
921  * const FE_Q<dim> fe;
922  * const MappingQ<dim> mapping;
923  *
924  * AffineConstraints<double> constraints;
925  *
926  * SparsityPattern sparsity_pattern;
927  * SparseMatrix<double> system_matrix;
928  *
929  * Vector<double> solution;
930  * Vector<double> system_rhs;
931  *
932  * MGLevelObject<SparsityPattern> mg_sparsity_patterns;
933  * MGLevelObject<SparsityPattern> mg_interface_sparsity_patterns;
934  *
935  * MGLevelObject<SparseMatrix<double>> mg_matrices;
936  * MGLevelObject<SparseMatrix<double>> mg_interface_in;
937  * MGLevelObject<SparseMatrix<double>> mg_interface_out;
938  *
939  * mg::Matrix<Vector<double>> mg_matrix;
940  * mg::Matrix<Vector<double>> mg_interface_matrix_in;
941  * mg::Matrix<Vector<double>> mg_interface_matrix_out;
942  *
943  * std::unique_ptr<MGSmoother<Vector<double>>> mg_smoother;
944  *
945  * using SmootherType =
947  * using SmootherAdditionalDataType = SmootherType::AdditionalData;
949  *
950  * MGConstrainedDoFs mg_constrained_dofs;
951  *
952  * Tensor<1, dim> advection_direction;
953  *
954  * const Settings settings;
955  * };
956  *
957  *
958  *
959  * template <int dim>
960  * AdvectionProblem<dim>::AdvectionProblem(const Settings &settings)
962  * , dof_handler(triangulation)
963  * , fe(settings.fe_degree)
964  * , mapping(settings.fe_degree)
965  * , settings(settings)
966  * {
967  * advection_direction[0] = -std::sin(numbers::PI / 6.0);
968  * if (dim >= 2)
969  * advection_direction[1] = std::cos(numbers::PI / 6.0);
970  * if (dim >= 3)
971  * AssertThrow(false, ExcNotImplemented());
972  * }
973  *
974  *
975  * @endcode
976  *
977  *
978  * <a name="codeAdvectionProblemsetup_systemcode"></a>
979  * <h4><code>AdvectionProblem::setup_system()</code></h4>
980  *
981 
982  *
983  * Here we first set up the DoFHandler, AffineConstraints, and
984  * SparsityPattern objects for both active and multigrid level meshes.
985  *
986 
987  *
988  * We could renumber the active DoFs with the DoFRenumbering class,
989  * but the smoothers only act on multigrid levels and as such, this
990  * would not matter for the computations. Instead, we will renumber the
991  * DoFs on each multigrid level below.
992  *
993  * @code
994  * template <int dim>
995  * void AdvectionProblem<dim>::setup_system()
996  * {
997  * const unsigned int n_levels = triangulation.n_levels();
998  *
999  * dof_handler.distribute_dofs(fe);
1000  *
1001  * solution.reinit(dof_handler.n_dofs());
1002  * system_rhs.reinit(dof_handler.n_dofs());
1003  *
1004  * constraints.clear();
1005  * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1006  *
1008  * mapping, dof_handler, 0, BoundaryValues<dim>(), constraints);
1010  * mapping, dof_handler, 1, BoundaryValues<dim>(), constraints);
1011  * constraints.close();
1012  *
1013  * DynamicSparsityPattern dsp(dof_handler.n_dofs());
1014  * DoFTools::make_sparsity_pattern(dof_handler,
1015  * dsp,
1016  * constraints,
1017  * /*keep_constrained_dofs = */ false);
1018  *
1019  * sparsity_pattern.copy_from(dsp);
1020  * system_matrix.reinit(sparsity_pattern);
1021  *
1022  * dof_handler.distribute_mg_dofs();
1023  *
1024  * @endcode
1025  *
1026  * Having enumerated the global degrees of freedom as well as (in
1027  * the last line above) the level degrees of freedom, let us
1028  * renumber the level degrees of freedom to get a better smoother
1029  * as explained in the introduction. The first block below
1030  * renumbers DoFs on each level in downstream or upstream
1031  * direction if needed. This is only necessary for point smoothers
1032  * (SOR and Jacobi) as the block smoothers operate on cells (see
1033  * `create_smoother()`). The blocks below then also implement
1034  * random numbering.
1035  *
1036  * @code
1037  * if (settings.smoother_type == "SOR" || settings.smoother_type == "Jacobi")
1038  * {
1039  * if (settings.dof_renumbering ==
1041  * settings.dof_renumbering ==
1042  * Settings::DoFRenumberingStrategy::upstream)
1043  * {
1044  * const Tensor<1, dim> direction =
1045  * (settings.dof_renumbering ==
1046  * Settings::DoFRenumberingStrategy::upstream ?
1047  * -1.0 :
1048  * 1.0) *
1049  * advection_direction;
1050  *
1051  * for (unsigned int level = 0; level < n_levels; ++level)
1052  * DoFRenumbering::downstream(dof_handler,
1053  * level,
1054  * direction,
1055  * /*dof_wise_renumbering = */ true);
1056  * }
1057  * else if (settings.dof_renumbering ==
1059  * {
1060  * for (unsigned int level = 0; level < n_levels; ++level)
1061  * DoFRenumbering::random(dof_handler, level);
1062  * }
1063  * else
1064  * Assert(false, ExcNotImplemented());
1065  * }
1066  *
1067  * @endcode
1068  *
1069  * The rest of the function just sets up data structures. The last
1070  * lines of the code below is unlike the other GMG tutorials, as
1071  * it sets up both the interface in and out matrices. We need this
1072  * since our problem is non-symmetric.
1073  *
1074  * @code
1075  * mg_constrained_dofs.clear();
1076  * mg_constrained_dofs.initialize(dof_handler);
1077  *
1078  * mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, {0, 1});
1079  *
1080  * mg_matrices.resize(0, n_levels - 1);
1081  * mg_matrices.clear_elements();
1082  * mg_interface_in.resize(0, n_levels - 1);
1083  * mg_interface_in.clear_elements();
1084  * mg_interface_out.resize(0, n_levels - 1);
1085  * mg_interface_out.clear_elements();
1086  * mg_sparsity_patterns.resize(0, n_levels - 1);
1087  * mg_interface_sparsity_patterns.resize(0, n_levels - 1);
1088  *
1089  * for (unsigned int level = 0; level < n_levels; ++level)
1090  * {
1091  * {
1092  * DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
1093  * dof_handler.n_dofs(level));
1094  * MGTools::make_sparsity_pattern(dof_handler, dsp, level);
1095  * mg_sparsity_patterns[level].copy_from(dsp);
1096  * mg_matrices[level].reinit(mg_sparsity_patterns[level]);
1097  * }
1098  * {
1099  * DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
1100  * dof_handler.n_dofs(level));
1102  * mg_constrained_dofs,
1103  * dsp,
1104  * level);
1105  * mg_interface_sparsity_patterns[level].copy_from(dsp);
1106  *
1107  * mg_interface_in[level].reinit(mg_interface_sparsity_patterns[level]);
1108  * mg_interface_out[level].reinit(mg_interface_sparsity_patterns[level]);
1109  * }
1110  * }
1111  * }
1112  *
1113  *
1114  * @endcode
1115  *
1116  *
1117  * <a name="codeAdvectionProblemassemble_cellcode"></a>
1118  * <h4><code>AdvectionProblem::assemble_cell()</code></h4>
1119  *
1120 
1121  *
1122  * Here we define the assembly of the linear system on each cell to
1123  * be used by the mesh_loop() function below. This one function
1124  * assembles the cell matrix for either an active or a level cell
1125  * (whatever it is passed as its first argument), and only assembles
1126  * a right-hand side if called with an active cell.
1127  *
1128 
1129  *
1130  *
1131  * @code
1132  * template <int dim>
1133  * template <class IteratorType>
1134  * void AdvectionProblem<dim>::assemble_cell(const IteratorType &cell,
1135  * ScratchData<dim> & scratch_data,
1136  * CopyData & copy_data)
1137  * {
1138  * copy_data.level = cell->level();
1139  *
1140  * const unsigned int dofs_per_cell =
1141  * scratch_data.fe_values.get_fe().dofs_per_cell;
1142  * copy_data.dofs_per_cell = dofs_per_cell;
1143  * copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1144  *
1145  * const unsigned int n_q_points =
1146  * scratch_data.fe_values.get_quadrature().size();
1147  *
1148  * if (cell->is_level_cell() == false)
1149  * copy_data.cell_rhs.reinit(dofs_per_cell);
1150  *
1151  * copy_data.local_dof_indices.resize(dofs_per_cell);
1152  * cell->get_active_or_mg_dof_indices(copy_data.local_dof_indices);
1153  *
1154  * scratch_data.fe_values.reinit(cell);
1155  *
1156  * RightHandSide<dim> right_hand_side;
1157  * std::vector<double> rhs_values(n_q_points);
1158  *
1159  * right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(),
1160  * rhs_values);
1161  *
1162  * @endcode
1163  *
1164  * If we are using streamline diffusion we must add its contribution
1165  * to both the cell matrix and the cell right-hand side. If we are not
1166  * using streamline diffusion, setting @f$\delta=0@f$ negates this contribution
1167  * below and we are left with the standard, Galerkin finite element
1168  * assembly.
1169  *
1170  * @code
1171  * const double delta = (settings.with_streamline_diffusion ?
1172  * compute_stabilization_delta(cell->diameter(),
1173  * settings.epsilon,
1174  * advection_direction,
1175  * settings.fe_degree) :
1176  * 0.0);
1177  *
1178  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1179  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1180  * {
1181  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1182  * {
1183  * @endcode
1184  *
1185  * The assembly of the local matrix has two parts. First
1186  * the Galerkin contribution:
1187  *
1188  * @code
1189  * copy_data.cell_matrix(i, j) +=
1190  * (settings.epsilon *
1191  * scratch_data.fe_values.shape_grad(i, q_point) *
1192  * scratch_data.fe_values.shape_grad(j, q_point) *
1193  * scratch_data.fe_values.JxW(q_point)) +
1194  * (scratch_data.fe_values.shape_value(i, q_point) *
1195  * (advection_direction *
1196  * scratch_data.fe_values.shape_grad(j, q_point)) *
1197  * scratch_data.fe_values.JxW(q_point))
1198  * @endcode
1199  *
1200  * and then the streamline diffusion contribution:
1201  *
1202  * @code
1203  * + delta *
1204  * (advection_direction *
1205  * scratch_data.fe_values.shape_grad(j, q_point)) *
1206  * (advection_direction *
1207  * scratch_data.fe_values.shape_grad(i, q_point)) *
1208  * scratch_data.fe_values.JxW(q_point) -
1209  * delta * settings.epsilon *
1210  * trace(scratch_data.fe_values.shape_hessian(j, q_point)) *
1211  * (advection_direction *
1212  * scratch_data.fe_values.shape_grad(i, q_point)) *
1213  * scratch_data.fe_values.JxW(q_point);
1214  * }
1215  * if (cell->is_level_cell() == false)
1216  * {
1217  * @endcode
1218  *
1219  * The same applies to the right hand side. First the
1220  * Galerkin contribution:
1221  *
1222  * @code
1223  * copy_data.cell_rhs(i) +=
1224  * scratch_data.fe_values.shape_value(i, q_point) *
1225  * rhs_values[q_point] * scratch_data.fe_values.JxW(q_point)
1226  * @endcode
1227  *
1228  * and then the streamline diffusion contribution:
1229  *
1230  * @code
1231  * + delta * rhs_values[q_point] * advection_direction *
1232  * scratch_data.fe_values.shape_grad(i, q_point) *
1233  * scratch_data.fe_values.JxW(q_point);
1234  * }
1235  * }
1236  * }
1237  *
1238  *
1239  * @endcode
1240  *
1241  *
1242  * <a name="codeAdvectionProblemassemble_system_and_multigridcode"></a>
1243  * <h4><code>AdvectionProblem::assemble_system_and_multigrid()</code></h4>
1244  *
1245 
1246  *
1247  * Here we employ MeshWorker::mesh_loop() to go over cells and assemble the
1248  * system_matrix, system_rhs, and all mg_matrices for us.
1249  *
1250 
1251  *
1252  *
1253  * @code
1254  * template <int dim>
1255  * void AdvectionProblem<dim>::assemble_system_and_multigrid()
1256  * {
1257  * const auto cell_worker_active =
1258  * [&](const decltype(dof_handler.begin_active()) &cell,
1259  * ScratchData<dim> & scratch_data,
1260  * CopyData & copy_data) {
1261  * this->assemble_cell(cell, scratch_data, copy_data);
1262  * };
1263  *
1264  * const auto copier_active = [&](const CopyData &copy_data) {
1265  * constraints.distribute_local_to_global(copy_data.cell_matrix,
1266  * copy_data.cell_rhs,
1267  * copy_data.local_dof_indices,
1268  * system_matrix,
1269  * system_rhs);
1270  * };
1271  *
1272  *
1273  * MeshWorker::mesh_loop(dof_handler.begin_active(),
1274  * dof_handler.end(),
1275  * cell_worker_active,
1276  * copier_active,
1277  * ScratchData<dim>(fe, fe.degree + 1),
1278  * CopyData(),
1280  *
1281  * @endcode
1282  *
1283  * Unlike the constraints for the active level, we choose to create
1284  * constraint objects for each multigrid level local to this function
1285  * since they are never needed elsewhere in the program.
1286  *
1287  * @code
1288  * std::vector<AffineConstraints<double>> boundary_constraints(
1289  * triangulation.n_global_levels());
1290  * for (unsigned int level = 0; level < triangulation.n_global_levels();
1291  * ++level)
1292  * {
1293  * IndexSet locally_owned_level_dof_indices;
1295  * dof_handler, level, locally_owned_level_dof_indices);
1296  * boundary_constraints[level].reinit(locally_owned_level_dof_indices);
1297  * boundary_constraints[level].add_lines(
1298  * mg_constrained_dofs.get_refinement_edge_indices(level));
1299  * boundary_constraints[level].add_lines(
1300  * mg_constrained_dofs.get_boundary_indices(level));
1301  * boundary_constraints[level].close();
1302  * }
1303  *
1304  * const auto cell_worker_mg =
1305  * [&](const decltype(dof_handler.begin_mg()) &cell,
1306  * ScratchData<dim> & scratch_data,
1307  * CopyData & copy_data) {
1308  * this->assemble_cell(cell, scratch_data, copy_data);
1309  * };
1310  *
1311  * const auto copier_mg = [&](const CopyData &copy_data) {
1312  * boundary_constraints[copy_data.level].distribute_local_to_global(
1313  * copy_data.cell_matrix,
1314  * copy_data.local_dof_indices,
1315  * mg_matrices[copy_data.level]);
1316  *
1317  * @endcode
1318  *
1319  * If @f$(i,j)@f$ is an `interface_out` dof pair, then @f$(j,i)@f$ is an
1320  * `interface_in` dof pair. Note: For `interface_in`, we load
1321  * the transpose of the interface entries, i.e., the entry for
1322  * dof pair @f$(j,i)@f$ is stored in `interface_in(i,j)`. This is an
1323  * optimization for the symmetric case which allows only one
1324  * matrix to be used when setting the edge_matrices in
1325  * solve(). Here, however, since our problem is non-symmetric,
1326  * we must store both `interface_in` and `interface_out`
1327  * matrices.
1328  *
1329  * @code
1330  * for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i)
1331  * for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j)
1332  * if (mg_constrained_dofs.is_interface_matrix_entry(
1333  * copy_data.level,
1334  * copy_data.local_dof_indices[i],
1335  * copy_data.local_dof_indices[j]))
1336  * {
1337  * mg_interface_out[copy_data.level].add(
1338  * copy_data.local_dof_indices[i],
1339  * copy_data.local_dof_indices[j],
1340  * copy_data.cell_matrix(i, j));
1341  * mg_interface_in[copy_data.level].add(
1342  * copy_data.local_dof_indices[i],
1343  * copy_data.local_dof_indices[j],
1344  * copy_data.cell_matrix(j, i));
1345  * }
1346  * };
1347  *
1348  * MeshWorker::mesh_loop(dof_handler.begin_mg(),
1349  * dof_handler.end_mg(),
1350  * cell_worker_mg,
1351  * copier_mg,
1352  * ScratchData<dim>(fe, fe.degree + 1),
1353  * CopyData(),
1355  * }
1356  *
1357  *
1358  * @endcode
1359  *
1360  *
1361  * <a name="codeAdvectionProblemsetup_smoothercode"></a>
1362  * <h4><code>AdvectionProblem::setup_smoother()</code></h4>
1363  *
1364 
1365  *
1366  * Next, we set up the smoother based on the settings in the `.prm` file. The
1367  * two options that are of significance is the number of pre- and
1368  * post-smoothing steps on each level of the multigrid v-cycle and the
1369  * relaxation parameter.
1370  *
1371 
1372  *
1373  * Since multiplicative methods tend to be more powerful than additive method,
1374  * fewer smoothing steps are required to see convergence independent of mesh
1375  * size. The same holds for block smoothers over point smoothers. This is
1376  * reflected in the choice for the number of smoothing steps for each type of
1377  * smoother below.
1378  *
1379 
1380  *
1381  * The relaxation parameter for point smoothers is chosen based on trial and
1382  * error, and reflects values necessary to keep the iteration counts in
1383  * the GMRES solve constant (or as close as possible) as we refine the mesh.
1384  * The two values given for both "Jacobi" and "SOR" in the `.prm` files are
1385  * for degree 1 and degree 3 finite elements. If the user wants to change to
1386  * another degree, they may need to adjust these numbers. For block smoothers,
1387  * this parameter has a more straightforward interpretation, namely that for
1388  * additive methods in 2D, a DoF can have a repeated contribution from up to 4
1389  * cells, therefore we must relax these methods by 0.25 to compensate. This is
1390  * not an issue for multiplicative methods as each cell's inverse application
1391  * carries new information to all its DoFs.
1392  *
1393 
1394  *
1395  * Finally, as mentioned above, the point smoothers only operate on DoFs, and
1396  * the block smoothers on cells, so only the block smoothers need to be given
1397  * information regarding cell orderings. DoF ordering for point smoothers has
1398  * already been taken care of in `setup_system()`.
1399  *
1400 
1401  *
1402  *
1403  * @code
1404  * template <int dim>
1405  * void AdvectionProblem<dim>::setup_smoother()
1406  * {
1407  * if (settings.smoother_type == "SOR")
1408  * {
1409  * using Smoother = PreconditionSOR<SparseMatrix<double>>;
1410  *
1411  * auto smoother =
1412  * std_cxx14::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1413  * Smoother,
1414  * Vector<double>>>();
1415  * smoother->initialize(mg_matrices,
1416  * Smoother::AdditionalData(fe.degree == 1 ? 1.0 :
1417  * 0.62));
1418  * smoother->set_steps(settings.smoothing_steps);
1419  * mg_smoother = std::move(smoother);
1420  * }
1421  * else if (settings.smoother_type == "Jacobi")
1422  * {
1423  * using Smoother = PreconditionJacobi<SparseMatrix<double>>;
1424  * auto smoother =
1425  * std_cxx14::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1426  * Smoother,
1427  * Vector<double>>>();
1428  * smoother->initialize(mg_matrices,
1429  * Smoother::AdditionalData(fe.degree == 1 ? 0.6667 :
1430  * 0.47));
1431  * smoother->set_steps(settings.smoothing_steps);
1432  * mg_smoother = std::move(smoother);
1433  * }
1434  * else if (settings.smoother_type == "block SOR" ||
1435  * settings.smoother_type == "block Jacobi")
1436  * {
1437  * smoother_data.resize(0, triangulation.n_levels() - 1);
1438  *
1439  * for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
1440  * {
1441  * DoFTools::make_cell_patches(smoother_data[level].block_list,
1442  * dof_handler,
1443  * level);
1444  *
1445  * smoother_data[level].relaxation =
1446  * (settings.smoother_type == "block SOR" ? 1.0 : 0.25);
1447  * smoother_data[level].inversion = PreconditionBlockBase<double>::svd;
1448  *
1449  * std::vector<unsigned int> ordered_indices;
1450  * switch (settings.dof_renumbering)
1451  * {
1452  * case Settings::DoFRenumberingStrategy::downstream:
1453  * ordered_indices =
1454  * create_downstream_cell_ordering(dof_handler,
1455  * advection_direction,
1456  * level);
1457  * break;
1458  *
1459  * case Settings::DoFRenumberingStrategy::upstream:
1460  * ordered_indices =
1461  * create_downstream_cell_ordering(dof_handler,
1462  * -1.0 * advection_direction,
1463  * level);
1464  * break;
1465  *
1466  * case Settings::DoFRenumberingStrategy::random:
1467  * ordered_indices =
1468  * create_random_cell_ordering(dof_handler, level);
1469  * break;
1470  *
1471  * case Settings::DoFRenumberingStrategy::none:
1472  * break;
1473  *
1474  * default:
1475  * AssertThrow(false, ExcNotImplemented());
1476  * break;
1477  * }
1478  *
1479  * smoother_data[level].order =
1480  * std::vector<std::vector<unsigned int>>(1, ordered_indices);
1481  * }
1482  *
1483  * if (settings.smoother_type == "block SOR")
1484  * {
1485  * auto smoother = std_cxx14::make_unique<MGSmootherPrecondition<
1486  * SparseMatrix<double>,
1487  * RelaxationBlockSOR<SparseMatrix<double>, double, Vector<double>>,
1488  * Vector<double>>>();
1489  * smoother->initialize(mg_matrices, smoother_data);
1490  * smoother->set_steps(settings.smoothing_steps);
1491  * mg_smoother = std::move(smoother);
1492  * }
1493  * else if (settings.smoother_type == "block Jacobi")
1494  * {
1495  * auto smoother = std_cxx14::make_unique<
1496  * MGSmootherPrecondition<SparseMatrix<double>,
1497  * RelaxationBlockJacobi<SparseMatrix<double>,
1498  * double,
1499  * Vector<double>>,
1500  * Vector<double>>>();
1501  * smoother->initialize(mg_matrices, smoother_data);
1502  * smoother->set_steps(settings.smoothing_steps);
1503  * mg_smoother = std::move(smoother);
1504  * }
1505  * }
1506  * else
1507  * AssertThrow(false, ExcNotImplemented());
1508  * }
1509  *
1510  *
1511  * @endcode
1512  *
1513  *
1514  * <a name="codeAdvectionProblemsolvecode"></a>
1515  * <h4><code>AdvectionProblem::solve()</code></h4>
1516  *
1517 
1518  *
1519  * Before we can solve the system, we must first set up the multigrid
1520  * preconditioner. This requires the setup of the transfer between levels,
1521  * the coarse matrix solver, and the smoother. This setup follows almost
1522  * identically to @ref step_16 "step-16", the main difference being the various smoothers
1523  * defined above and the fact that we need different interface edge matrices
1524  * for in and out since our problem is non-symmetric. (In reality, for this
1525  * tutorial these interface matrices are empty since we are only using global
1526  * refinement, and thus have no refinement edges. However, we have still
1527  * included both here since if one made the simple switch to an adaptively
1528  * refined method, the program would still run correctly.)
1529  *
1530 
1531  *
1532  * The last thing to note is that since our problem is non-symmetric, we must
1533  * use an appropriate Krylov subspace method. We choose here to
1534  * use GMRES since it offers the guarantee of residual reduction in each
1535  * iteration. The major disavantage of GMRES is that, for each iteration,
1536  * the number of stored temporary vectors increases by one, and one also needs
1537  * to compute a scalar product with all previously stored vectors. This is
1538  * rather expensive. This requirement is relaxed by using the restarted GMRES
1539  * method which puts a cap on the number of vectors we are required to store
1540  * at any one time (here we restart after 50 temporary vectors, or 48
1541  * iterations). This then has the disadvantage that we lose information we
1542  * have gathered throughout the iteration and therefore we could see slower
1543  * convergence. As a consequence, where to restart is a question of balancing
1544  * memory consumption, CPU effort, and convergence speed.
1545  * However, the goal of this tutorial is to have very low
1546  * iteration counts by using a powerful GMG preconditioner, so we have picked
1547  * the restart length such that all of the results shown below converge prior
1548  * to restart happening, and thus we have a standard GMRES method. If the user
1549  * is interested, another suitable method offered in deal.II would be
1550  * BiCGStab.
1551  *
1552 
1553  *
1554  *
1555  * @code
1556  * template <int dim>
1557  * void AdvectionProblem<dim>::solve()
1558  * {
1559  * const unsigned int max_iters = 200;
1560  * const double solve_tolerance = 1e-8 * system_rhs.l2_norm();
1561  * SolverControl solver_control(max_iters, solve_tolerance, true, true);
1562  * solver_control.enable_history_data();
1563  *
1564  * using Transfer = MGTransferPrebuilt<Vector<double>>;
1565  * Transfer mg_transfer(mg_constrained_dofs);
1566  * mg_transfer.build(dof_handler);
1567  *
1568  * FullMatrix<double> coarse_matrix;
1569  * coarse_matrix.copy_from(mg_matrices[0]);
1570  * MGCoarseGridHouseholder<double, Vector<double>> coarse_grid_solver;
1571  * coarse_grid_solver.initialize(coarse_matrix);
1572  *
1573  * setup_smoother();
1574  *
1575  * mg_matrix.initialize(mg_matrices);
1576  * mg_interface_matrix_in.initialize(mg_interface_in);
1577  * mg_interface_matrix_out.initialize(mg_interface_out);
1578  *
1579  * Multigrid<Vector<double>> mg(
1580  * mg_matrix, coarse_grid_solver, mg_transfer, *mg_smoother, *mg_smoother);
1581  * mg.set_edge_matrices(mg_interface_matrix_out, mg_interface_matrix_in);
1582  *
1583  * PreconditionMG<dim, Vector<double>, Transfer> preconditioner(dof_handler,
1584  * mg,
1585  * mg_transfer);
1586  *
1587  * std::cout << " Solving with GMRES to tol " << solve_tolerance << "..."
1588  * << std::endl;
1589  * SolverGMRES<Vector<double>> solver(
1590  * solver_control, SolverGMRES<Vector<double>>::AdditionalData(50, true));
1591  *
1592  * Timer time;
1593  * time.start();
1594  * solver.solve(system_matrix, solution, system_rhs, preconditioner);
1595  * time.stop();
1596  *
1597  * std::cout << " converged in " << solver_control.last_step()
1598  * << " iterations"
1599  * << " in " << time.last_wall_time() << " seconds " << std::endl;
1600  *
1601  * constraints.distribute(solution);
1602  *
1603  * mg_smoother.release();
1604  * }
1605  *
1606  *
1607  * @endcode
1608  *
1609  *
1610  * <a name="codeAdvectionProblemoutput_resultscode"></a>
1611  * <h4><code>AdvectionProblem::output_results()</code></h4>
1612  *
1613 
1614  *
1615  * The final function of interest generates graphical output.
1616  * Here we output the solution and cell ordering in a .vtu format.
1617  *
1618 
1619  *
1620  * At the top of the function, we generate an index for each cell to
1621  * visualize the ordering used by the smoothers. Note that we do
1622  * this only for the active cells instead of the levels, where the
1623  * smoothers are actually used. For the point smoothers we renumber
1624  * DoFs instead of cells, so this is only an approximation of what
1625  * happens in reality. Finally, the random ordering is not the
1626  * random ordering we actually use (see `create_smoother()` for that).
1627  *
1628 
1629  *
1630  * The (integer) ordering of cells is then copied into a (floating
1631  * point) vector for graphical output.
1632  *
1633  * @code
1634  * template <int dim>
1635  * void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1636  * {
1637  * const unsigned int n_active_cells = triangulation.n_active_cells();
1638  * Vector<double> cell_indices(n_active_cells);
1639  * {
1640  * std::vector<unsigned int> ordered_indices;
1641  * switch (settings.dof_renumbering)
1642  * {
1643  * case Settings::DoFRenumberingStrategy::downstream:
1644  * ordered_indices =
1645  * create_downstream_cell_ordering(dof_handler, advection_direction);
1646  * break;
1647  *
1648  * case Settings::DoFRenumberingStrategy::upstream:
1649  * ordered_indices =
1650  * create_downstream_cell_ordering(dof_handler,
1651  * -1.0 * advection_direction);
1652  * break;
1653  *
1654  * case Settings::DoFRenumberingStrategy::random:
1655  * ordered_indices = create_random_cell_ordering(dof_handler);
1656  * break;
1657  *
1658  * case Settings::DoFRenumberingStrategy::none:
1659  * ordered_indices.resize(n_active_cells);
1660  * for (unsigned int i = 0; i < n_active_cells; ++i)
1661  * ordered_indices[i] = i;
1662  * break;
1663  *
1664  * default:
1665  * AssertThrow(false, ExcNotImplemented());
1666  * break;
1667  * }
1668  *
1669  * for (unsigned int i = 0; i < n_active_cells; ++i)
1670  * cell_indices(ordered_indices[i]) = static_cast<double>(i);
1671  * }
1672  *
1673  * @endcode
1674  *
1675  * The remainder of the function is then straightforward, given
1676  * previous tutorial programs:
1677  *
1678  * @code
1679  * DataOut<dim> data_out;
1680  * data_out.attach_dof_handler(dof_handler);
1681  * data_out.add_data_vector(solution, "solution");
1682  * data_out.add_data_vector(cell_indices, "cell_index");
1683  * data_out.build_patches();
1684  *
1685  * const std::string filename =
1686  * "solution-" + Utilities::int_to_string(cycle) + ".vtu";
1687  * std::ofstream output(filename.c_str());
1688  * data_out.write_vtu(output);
1689  * }
1690  *
1691  *
1692  * @endcode
1693  *
1694  *
1695  * <a name="codeAdvectionProblemruncode"></a>
1696  * <h4><code>AdvectionProblem::run()</code></h4>
1697  *
1698 
1699  *
1700  * As in most tutorials, this function creates/refines the mesh and calls
1701  * the various functions defined above to set up, assemble, solve, and output
1702  * the results.
1703  *
1704 
1705  *
1706  * In cycle zero, we generate the mesh for the on the square
1707  * <code>[-1,1]^dim</code> with a hole of radius 3/10 units centered
1708  * at the origin. For objects with `manifold_id` equal to one
1709  * (namely, the faces adjacent to the hole), we assign a spherical
1710  * manifold.
1711  *
1712 
1713  *
1714  *
1715  * @code
1716  * template <int dim>
1717  * void AdvectionProblem<dim>::run()
1718  * {
1719  * for (unsigned int cycle = 0; cycle < (settings.fe_degree == 1 ? 7 : 5);
1720  * ++cycle)
1721  * {
1722  * std::cout << " Cycle " << cycle << ':' << std::endl;
1723  *
1724  * if (cycle == 0)
1725  * {
1726  * GridGenerator::hyper_cube_with_cylindrical_hole(triangulation,
1727  * 0.3,
1728  * 1.0);
1729  *
1730  * const SphericalManifold<dim> manifold_description(Point<dim>(0, 0));
1731  * triangulation.set_manifold(1, manifold_description);
1732  * }
1733  *
1734  * triangulation.refine_global();
1735  *
1736  * setup_system();
1737  *
1738  * std::cout << " Number of active cells: "
1739  * << triangulation.n_active_cells() << " ("
1740  * << triangulation.n_levels() << " levels)" << std::endl;
1741  * std::cout << " Number of degrees of freedom: "
1742  * << dof_handler.n_dofs() << std::endl;
1743  *
1744  * assemble_system_and_multigrid();
1745  *
1746  * solve();
1747  *
1748  * if (settings.output)
1749  * output_results(cycle);
1750  *
1751  * std::cout << std::endl;
1752  * }
1753  * }
1754  * } // namespace Step63
1755  *
1756  *
1757  * @endcode
1758  *
1759  *
1760  * <a name="Thecodemaincodefunction"></a>
1761  * <h3>The <code>main</code> function</h3>
1762  *
1763 
1764  *
1765  * Finally, the main function is like most tutorials. The only
1766  * interesting bit is that we require the user to pass a `.prm` file
1767  * as a sole command line argument. If no parameter file is given, the
1768  * program will output the contents of a sample parameter file with
1769  * all default values to the screen that the user can then copy and
1770  * paste into their own `.prm` file.
1771  *
1772 
1773  *
1774  *
1775  * @code
1776  * int main(int argc, char *argv[])
1777  * {
1778  * try
1779  * {
1780  * Step63::Settings settings;
1781  * settings.get_parameters((argc > 1) ? (argv[1]) : "");
1782  *
1783  * Step63::AdvectionProblem<2> advection_problem_2d(settings);
1784  * advection_problem_2d.run();
1785  * }
1786  * catch (std::exception &exc)
1787  * {
1788  * std::cerr << std::endl
1789  * << std::endl
1790  * << "----------------------------------------------------"
1791  * << std::endl;
1792  * std::cerr << "Exception on processing: " << std::endl
1793  * << exc.what() << std::endl
1794  * << "Aborting!" << std::endl
1795  * << "----------------------------------------------------"
1796  * << std::endl;
1797  * return 1;
1798  * }
1799  * catch (...)
1800  * {
1801  * std::cerr << std::endl
1802  * << std::endl
1803  * << "----------------------------------------------------"
1804  * << std::endl;
1805  * std::cerr << "Unknown exception!" << std::endl
1806  * << "Aborting!" << std::endl
1807  * << "----------------------------------------------------"
1808  * << std::endl;
1809  * return 1;
1810  * }
1811  *
1812  * return 0;
1813  * }
1814  * @endcode
1815 <a name="Results"></a><h1>Results</h1>
1816 
1817 
1818 <a name="GMRESIterationNumbers"></a><h3> GMRES Iteration Numbers </h3>
1819 
1820 
1821 The major advantage for GMG is that it is an @f$\mathcal{O}(n)@f$ method,
1822 that is, the complexity of the problem increases linearly with the
1823 problem size. To show then that the linear solver presented in this
1824 tutorial is in fact @f$\mathcal{O}(n)@f$, all one needs to do is show that
1825 the iteration counts for the GMRES solve stay roughly constant as we
1826 refine the mesh.
1827 
1828 Each of the following tables gives the GMRES iteration counts to reduce the
1829 initial residual by a factor of @f$10^8@f$. We selected a sufficient number of smoothing steps
1830 (based on the method) to get iteration numbers independent of mesh size. As
1831 can be seen from the tables below, the method is indeed @f$\mathcal{O}(n)@f$.
1832 
1833 <a name="DoFCellRenumbering"></a><h4> DoF/Cell Renumbering </h4>
1834 
1835 
1836 The point-wise smoothers ("Jacobi" and "SOR") get applied in the order the
1837 DoFs are numbered on each level. We can influence this using the
1838 DoFRenumbering namespace. The block smoothers are applied based on the
1839 ordering we set in `setup_smoother()`. We can visualize this numbering. The
1840 following pictures show the cell numbering of the active cells in downstream,
1841 random, and upstream numbering (left to right):
1842 
1843 <img src="https://www.dealii.org/images/steps/developer/step-63-cell-order.png" alt="">
1844 
1845 Let us start with the additive smoothers. The following table shows
1846 the number of iterations necessary to obtain convergence from GMRES:
1847 
1848 <table align="center" class="doxtable">
1849 <tr>
1850  <th></th>
1851  <th></th>
1852  <th colspan="1">@f$Q_1@f$</th>
1853  <th colspan="7">Smoother (smoothing steps)</th>
1854 </tr>
1855 <tr>
1856  <th></th>
1857  <th></th>
1858  <th></th>
1859  <th colspan="3">Jacobi (6)</th>
1860  <th></th>
1861  <th colspan="3">Block Jacobi (3)</th>
1862 </tr>
1863 <tr>
1864  <th></th>
1865  <th></th>
1866  <th></th>
1867  <th colspan="3">Renumbering Strategy</th>
1868  <th></th>
1869  <th colspan="3">Renumbering Strategy</th>
1870 </tr>
1871 <tr>
1872  <th>Cells</th>
1873  <th></th>
1874  <th>DoFs</th>
1875  <th>Downstream</th>
1876  <th>Random</th>
1877  <th>Upstream</th>
1878  <th></th>
1879  <th>Downstream</th>
1880  <th>Random</th>
1881  <th>Upstream</th>
1882 </tr>
1883 <tr>
1884  <th>32</th>
1885  <th></th>
1886  <th>48</th>
1887  <td>3</th>
1888  <td>3</th>
1889  <td>3</th>
1890  <th></th>
1891  <td>3</th>
1892  <td>3</th>
1893  <td>3</th>
1894 </tr>
1895 <tr>
1896  <th>128</th>
1897  <th></th>
1898  <th>160</th>
1899  <td>6</th>
1900  <td>6</th>
1901  <td>6</th>
1902  <th></th>
1903  <td>6</th>
1904  <td>6</th>
1905  <td>6</th>
1906 </tr>
1907 <tr>
1908  <th>512</th>
1909  <th></th>
1910  <th>576</th>
1911  <td>11</th>
1912  <td>11</th>
1913  <td>11</th>
1914  <th></th>
1915  <td>9</th>
1916  <td>9</th>
1917  <td>9</th>
1918 </tr>
1919 <tr>
1920  <th>2048</th>
1921  <th></th>
1922  <th>2176</th>
1923  <td>15</th>
1924  <td>15</th>
1925  <td>15</th>
1926  <th></th>
1927  <td>13</th>
1928  <td>13</th>
1929  <td>13</th>
1930 </tr>
1931 <tr>
1932  <th>8192</th>
1933  <th></th>
1934  <th>8448</th>
1935  <td>18</th>
1936  <td>18</th>
1937  <td>18</th>
1938  <th></th>
1939  <td>15</th>
1940  <td>15</th>
1941  <td>15</th>
1942 </tr>
1943 <tr>
1944  <th>32768</th>
1945  <th></th>
1946  <th>33280</th>
1947  <td>20</th>
1948  <td>20</th>
1949  <td>20</th>
1950  <th></th>
1951  <td>16</th>
1952  <td>16</th>
1953  <td>16</th>
1954 </tr>
1955 <tr>
1956  <th>131072</th>
1957  <th></th>
1958  <th>132096</th>
1959  <td>20</th>
1960  <td>20</th>
1961  <td>20</th>
1962  <th></th>
1963  <td>16</th>
1964  <td>16</th>
1965  <td>16</th>
1966 </tr>
1967 </table>
1968 
1969 We see that renumbering the
1970 DoFs/cells has no effect on convergence speed. This is because these
1971 smoothers compute operations on each DoF (point-smoother) or cell
1972 (block-smoother) independently and add up the results. Since we can
1973 define these smoothers as an application of a sum of matrices, and
1974 matrix addition is commutative, the order at which we sum the
1975 different components will not affect the end result.
1976 
1977 On the other hand, the situation is different for multiplicative smoothers:
1978 
1979 <table align="center" class="doxtable">
1980 <tr>
1981  <th></th>
1982  <th></th>
1983  <th colspan="1">@f$Q_1@f$</th>
1984  <th colspan="7">Smoother (smoothing steps)</th>
1985 </tr>
1986 <tr>
1987  <th></th>
1988  <th></th>
1989  <th></th>
1990  <th colspan="3">SOR (3)</th>
1991  <th></th>
1992  <th colspan="3">Block SOR (1)</th>
1993 </tr>
1994 <tr>
1995  <th></th>
1996  <th></th>
1997  <th></th>
1998  <th colspan="3">Renumbering Strategy</th>
1999  <th></th>
2000  <th colspan="3">Renumbering Strategy</th>
2001 </tr>
2002 <tr>
2003  <th>Cells</th>
2004  <th></th>
2005  <th>DoFs</th>
2006  <th>Downstream</th>
2007  <th>Random</th>
2008  <th>Upstream</th>
2009  <th></th>
2010  <th>Downstream</th>
2011  <th>Random</th>
2012  <th>Upstream</th>
2013 </tr>
2014 <tr>
2015  <th>32</th>
2016  <th></th>
2017  <th>48</th>
2018  <td>2</th>
2019  <td>2</th>
2020  <td>3</th>
2021  <th></th>
2022  <td>2</th>
2023  <td>2</th>
2024  <td>3</th>
2025 </tr>
2026 <tr>
2027  <th>128</th>
2028  <th></th>
2029  <th>160</th>
2030  <td>5</th>
2031  <td>5</th>
2032  <td>7</th>
2033  <th></th>
2034  <td>5</th>
2035  <td>5</th>
2036  <td>7</th>
2037 </tr>
2038 <tr>
2039  <th>512</th>
2040  <th></th>
2041  <th>576</th>
2042  <td>7</th>
2043  <td>9</th>
2044  <td>11</th>
2045  <th></th>
2046  <td>7</th>
2047  <td>7</th>
2048  <td>12</th>
2049 </tr>
2050 <tr>
2051  <th>2048</th>
2052  <th></th>
2053  <th>2176</th>
2054  <td>10</th>
2055  <td>12</th>
2056  <td>15</th>
2057  <th></th>
2058  <td>8</th>
2059  <td>10</th>
2060  <td>17</th>
2061 </tr>
2062 <tr>
2063  <th>8192</th>
2064  <th></th>
2065  <th>8448</th>
2066  <td>11</th>
2067  <td>15</th>
2068  <td>19</th>
2069  <th></th>
2070  <td>10</th>
2071  <td>11</th>
2072  <td>20</th>
2073 </tr>
2074 <tr>
2075  <th>32768</th>
2076  <th></th>
2077  <th>33280</th>
2078  <td>12</th>
2079  <td>16</th>
2080  <td>20</th>
2081  <th></th>
2082  <td>10</th>
2083  <td>12</th>
2084  <td>21</th>
2085 </tr>
2086 <tr>
2087  <th>131072</th>
2088  <th></th>
2089  <th>132096</th>
2090  <td>12</th>
2091  <td>16</th>
2092  <td>19</th>
2093  <th></th>
2094  <td>11</th>
2095  <td>12</th>
2096  <td>21</th>
2097 </tr>
2098 </table>
2099 
2100 Here, we can speed up
2101 convergence by renumbering the DoFs/cells in the advection direction,
2102 and similarly, we can slow down convergence if we do the renumbering
2103 in the opposite direction. This is because advection-dominated
2104 problems have a directional flow of information (in the advection
2105 direction) which, given the right renumbering of DoFs/cells,
2106 multiplicative methods are able to capture.
2107 
2108 This feature of multiplicative methods is, however, dependent on the
2109 value of @f$\varepsilon@f$. As we increase @f$\varepsilon@f$ and the problem
2110 becomes more diffusion-dominated, we have a more uniform propagation
2111 of information over the mesh and there is a diminished advantage for
2112 renumbering in the advection direction. On the opposite end, in the
2113 extreme case of @f$\varepsilon=0@f$ (advection-only), we have a 1st-order
2114 PDE and multiplicative methods with the right renumbering become
2115 effective solvers: A correct downstream numbering may lead to methods
2116 that require only a single iteration because information can be
2117 propagated from the inflow boundary downstream, with no information
2118 transport in the opposite direction. (Note, however, that in the case
2119 of @f$\varepsilon=0@f$, special care must be taken for the boundary
2120 conditions in this case).
2121 
2122 
2123 <a name="Pointvsblocksmoothers"></a><h4> %Point vs. block smoothers </h4>
2124 
2125 
2126 We will limit the results to runs using the downstream
2127 renumbering. Here is a cross comparison of all four smoothers for both
2128 @f$Q_1@f$ and @f$Q_3@f$ elements:
2129 
2130 <table align="center" class="doxtable">
2131 <tr>
2132  <th></th>
2133  <td></th>
2134  <th colspan="1">@f$Q_1@f$</th>
2135  <th colspan="4">Smoother (smoothing steps)</th>
2136  <th></th>
2137  <th colspan="1">@f$Q_3@f$</th>
2138  <th colspan="4">Smoother (smoothing steps)</th>
2139 </tr>
2140 <tr>
2141  <th colspan="1">Cells</th>
2142  <td></th>
2143  <th colspan="1">DoFs</th>
2144  <th colspan="1">Jacobi (6)</th>
2145  <th colspan="1">Block Jacobi (3)</th>
2146  <th colspan="1">SOR (3)</th>
2147  <th colspan="1">Block SOR (1)</th>
2148  <th></th>
2149  <th colspan="1">DoFs</th>
2150  <th colspan="1">Jacobi (6)</th>
2151  <th colspan="1">Block Jacobi (3)</th>
2152  <th colspan="1">SOR (3)</th>
2153  <th colspan="1">Block SOR (1)</th>
2154 </tr>
2155 <tr>
2156  <th>32</th>
2157  <td></th>
2158  <th>48</th>
2159  <td>3</th>
2160  <td>3</th>
2161  <td>2</th>
2162  <td>2</th>
2163  <td></th>
2164  <th>336</th>
2165  <td>15</th>
2166  <td>14</th>
2167  <td>15</th>
2168  <td>6</th>
2169 </tr>
2170 <tr>
2171  <th>128</th>
2172  <td></th>
2173  <th>160</th>
2174  <td>6</th>
2175  <td>6</th>
2176  <td>5</th>
2177  <td>5</th>
2178  <td></th>
2179  <th>1248</th>
2180  <td>23</th>
2181  <td>18</th>
2182  <td>21</th>
2183  <td>9</th>
2184 </tr>
2185 <tr>
2186  <th>512</th>
2187  <td></th>
2188  <th>576</th>
2189  <td>11</th>
2190  <td>9</th>
2191  <td>7</th>
2192  <td>7</th>
2193  <td></th>
2194  <th>4800</th>
2195  <td>29</th>
2196  <td>21</th>
2197  <td>28</th>
2198  <td>9</th>
2199 </tr>
2200 <tr>
2201  <th>2048</th>
2202  <td></th>
2203  <th>2176</th>
2204  <td>15</th>
2205  <td>13</th>
2206  <td>10</th>
2207  <td>8</th>
2208  <td></th>
2209  <th>18816</th>
2210  <td>33</th>
2211  <td>22</th>
2212  <td>32</th>
2213  <td>9</th>
2214 </tr>
2215 <tr>
2216  <th>8192</th>
2217  <td></th>
2218  <th>8448</th>
2219  <td>18</th>
2220  <td>15</th>
2221  <td>11</th>
2222  <td>10</th>
2223  <td></th>
2224  <th>74496</th>
2225  <td>35</th>
2226  <td>22</th>
2227  <td>34</th>
2228  <td>10</th>
2229 </tr>
2230 <tr>
2231  <th>32768</th>
2232  <td></th>
2233  <th>33280</th>
2234  <td>20</th>
2235  <td>16</th>
2236  <td>12</th>
2237  <td>10</th>
2238  <td></th>
2239 </tr>
2240 <tr>
2241  <th>131072</th>
2242  <td></th>
2243  <th>132096</th>
2244  <td>20</th>
2245  <td>16</th>
2246  <td>12</th>
2247  <td>11</th>
2248  <td></th>
2249 </tr>
2250 </table>
2251 
2252 We see that for @f$Q_1@f$, both multiplicative smoothers require a smaller
2253 combination of smoothing steps and iteration counts than either
2254 additive smoother. However, when we increase the degree to a @f$Q_3@f$
2255 element, there is a clear advantage for the block smoothers in terms
2256 of the number of smoothing steps and iterations required to
2257 solve. Specifically, the block SOR smoother gives constant iteration
2258 counts over the degree, and the block Jacobi smoother only sees about
2259 a 38% increase in iterations compared to 75% and 183% for Jacobi and
2260 SOR respectively.
2261 
2262 <a name="Cost"></a><h3> Cost </h3>
2263 
2264 
2265 Iteration counts do not tell the full story in the optimality of a one
2266 smoother over another. Obviously we must examine the cost of an
2267 iteration. Block smoothers here are at a disadvantage as they are
2268 having to construct and invert a cell matrix for each cell. Here is a
2269 comparison of solve times for a @f$Q_3@f$ element with 74,496 DoFs:
2270 
2271 <table align="center" class="doxtable">
2272 <tr>
2273  <th colspan="1">@f$Q_3@f$</th>
2274  <th colspan="4">Smoother (smoothing steps)</th>
2275 </tr>
2276 <tr>
2277  <th colspan="1">DoFs</th>
2278  <th colspan="1">Jacobi (6)</th>
2279  <th colspan="1">Block Jacobi (3)</th>
2280  <th colspan="1">SOR (3)</th>
2281  <th colspan="1">Block SOR (1)</th>
2282 </tr>
2283 <tr>
2284  <th>74496</th>
2285  <td>0.68s</th>
2286  <td>5.82s</th>
2287  <td>1.18s</th>
2288  <td>1.02s</th>
2289 </tr>
2290 </table>
2291 
2292 The smoother that requires the most iterations (Jacobi) actually takes
2293 the shortest time (roughly 2/3 the time of the next fastest
2294 method). This is because all that is required to apply a Jacobi
2295 smoothing step is multiplication by a diagonal matrix which is very
2296 cheap. On the other hand, while SOR requires over 3x more iterations
2297 (each with 3x more smoothing steps) than block SOR, the times are
2298 roughly equivalent, implying that a smoothing step of block SOR is
2299 roughly 9x slower than a smoothing step of SOR. Lastly, block Jacobi
2300 is almost 6x more expensive than block SOR, which intuitively makes
2301 sense from the fact that 1 step of each method has the same cost
2302 (inverting the cell matrices and either adding or multiply them
2303 together), and block Jacobi has 3 times the number of smoothing steps per
2304 iteration with 2 times the iterations.
2305 
2306 
2307 <a name="Additionalpoints"></a><h3> Additional points </h3>
2308 
2309 
2310 There are a few more important points to mention:
2311 
2312 <ol>
2313 <li> For a mesh distributed in parallel, multiplicative methods cannot
2314 be executed over the entire domain. This is because they operate one
2315 cell at a time, and downstream cells can only be handled once upstream
2316 cells have already been done. This is fine on a single processor: The
2317 processor just goes through the list of cells one after the
2318 other. However, in parallel, it would imply that some processors are
2319 idle because upstream processors have not finished doing the work on
2320 cells upstream from the ones owned by the current processor. Once the
2321 upstream processors are done, the downstream ones can start, but by
2322 that time the upstream processors have no work left. In other words,
2323 most of the time during these smoother steps, most processors are in
2324 fact idle. This is not how one obtains good parallel scalability!
2325 
2326 One can use a hybrid method where
2327 a multiplicative smoother is applied on each subdomain, but as you
2328 increase the number of subdomains, the method approaches the behavior
2329 of an additive method. This is a major disadvantage to these methods.
2330 </li>
2331 
2332 <li> Current research into block smoothers suggest that soon we will be
2333 able to compute the inverse of the cell matrices much cheaper than
2334 what is currently being done inside deal.II. This research is based on
2335 the fast diagonalization method (dating back to the 1960s) and has
2336 been used in the spectral community for around 20 years (see, e.g., <a
2337 href="https://doi.org/10.1007/s10915-004-4787-3"> Hybrid
2338 Multigrid/Schwarz Algorithms for the Spectral Element Method by Lottes
2339 and Fischer</a>). There are currently efforts to generalize these
2340 methods to DG and make them more robust. Also, it seems that one
2341 should be able to take advantage of matrix-free implementations and
2342 the fact that, in the interior of the domain, cell matrices tend to
2343 look very similar, allowing fewer matrix inverse computations.
2344 </li>
2345 </ol>
2346 
2347 Combining 1. and 2. gives a good reason for expecting that a method
2348 like block Jacobi could become very powerful in the future, even
2349 though currently for these examples it is quite slow.
2350 
2351 
2352 <a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
2353 
2354 
2355 <a name="ConstantiterationsforQsub5sub"></a><h4> Constant iterations for Q<sub>5</sub> </h4>
2356 
2357 
2358 Change the number of smoothing steps and the smoother relaxation
2359 parameter (set in <code>Smoother::AdditionalData()</code> inside
2360 <code>create_smoother()</code>, only necessary for point smoothers) so
2361 that we maintain a constant number of iterations for a @f$Q_5@f$ element.
2362 
2363 <a name="Effectivenessofrenumberingforchangingepsilon"></a><h4> Effectiveness of renumbering for changing epsilon </h4>
2364 
2365 
2366 Increase/decrease the parameter "Epsilon" in the `.prm` files of the
2367 multiplicative methods and observe for which values renumbering no
2368 longer influences convergence speed.
2369 
2370 <a name="Meshadaptivity"></a><h4> Mesh adaptivity </h4>
2371 
2372 
2373 The code is set up to work correctly with an adaptively refined mesh (the
2374 interface matrices are created and set). Devise a suitable refinement
2375 criterium or try the KellyErrorEstimator class.
2376 
2377 
2378 
2379  *
2380  *
2381 <a name="PlainProg"></a>
2382 <h1> The plain program</h1>
2383 @include "step-63.cc"
2384 */
SymmetricTensor::trace
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2705
Physics::Elasticity::Kinematics::epsilon
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
std::exp
inline ::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5368
ParameterHandler::get
std::string get(const std::string &entry_string) const
Definition: parameter_handler.cc:975
IndexSet
Definition: index_set.h:74
MappingQ
Definition: mapping_manifold.h:33
ParameterHandler::get_double
double get_double(const std::string &entry_name) const
Definition: parameter_handler.cc:1056
ParameterHandler::Text
@ Text
Definition: parameter_handler.h:894
MGSmoother
Definition: mg_smoother.h:50
DoFRenumbering::random
void random(DoFHandlerType &dof_handler)
Definition: dof_renumbering.cc:2102
FE_Q
Definition: fe_q.h:554
Differentiation::SD::coth
Expression coth(const Expression &x)
Definition: symengine_math.cc:217
MeshWorker::assemble_own_cells
@ assemble_own_cells
Definition: assemble_flags.h:50
ParameterHandler::get_bool
bool get_bool(const std::string &entry_name) const
Definition: parameter_handler.cc:1101
MGConstrainedDoFs
Definition: mg_constrained_dofs.h:45
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Triangulation< dim >
ParameterHandler::declare_entry
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
Definition: parameter_handler.cc:784
DoFRenumbering::downstream
void downstream(DoFHandlerType &dof_handler, const Tensor< 1, DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering=false)
Definition: dof_renumbering.cc:1735
SparseMatrix< double >
MeshWorker::mesh_loop
void mesh_loop(const CellIteratorType &begin, const typename identity< CellIteratorType >::type &end, const typename identity< std::function< void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type &cell_worker, const typename identity< std::function< void(const CopyData &)>>::type &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const typename identity< std::function< void(const CellIteratorBaseType &, const unsigned int, ScratchData &, CopyData &)>>::type &boundary_worker=std::function< void(const CellIteratorBaseType &, const unsigned int, ScratchData &, CopyData &)>(), const typename identity< std::function< void(const CellIteratorBaseType &, const unsigned int, const unsigned int, const CellIteratorBaseType &, const unsigned int, const unsigned int, ScratchData &, CopyData &)>>::type &face_worker=std::function< void(const CellIteratorBaseType &, const unsigned int, const unsigned int, const CellIteratorBaseType &, const unsigned int, const unsigned int, ScratchData &, CopyData &)>(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: mesh_loop.h:237
Patterns::Bool
Definition: patterns.h:984
Function::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
MGTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:565
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
std::cos
inline ::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5324
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition: symengine_math.cc:273
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
DataOutBase::none
@ none
Definition: data_out_base.h:1557
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
MGTools::make_interface_sparsity_pattern
void make_interface_sparsity_pattern(const DoFHandlerType &dof_handler, const MGConstrainedDoFs &mg_constrained_dofs, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:1023
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
ParameterHandler::get_integer
long int get_integer(const std::string &entry_string) const
Definition: parameter_handler.cc:1013
Tensor< 1, dim >
mg::Matrix
Definition: mg_matrix.h:47
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
ParameterHandler::print_parameters
std::ostream & print_parameters(std::ostream &out, const OutputStyle style) const
Definition: parameter_handler.cc:1238
SparsityPattern
Definition: sparsity_pattern.h:865
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
numbers
Definition: numbers.h:207
DoFRenumbering
Definition: dof_renumbering.h:345
value
static const bool value
Definition: dof_tools_constraints.cc:433
AffineConstraints< double >
DataOutBase::eps
@ eps
Definition: data_out_base.h:1582
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
MGLevelObject< SparsityPattern >
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
Functions::ZeroFunction
Definition: function.h:513
ParameterHandler
Definition: parameter_handler.h:845
DoFTools::extract_locally_relevant_level_dofs
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1215
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
Patterns::Selection
Definition: patterns.h:383
TriangulationDescription::Settings
Settings
Definition: tria_description.h:253
first
Point< 2 > first
Definition: grid_out.cc:4352
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector< double >
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
center
Point< 3 > center
Definition: data_out_base.cc:169
ParameterHandler::parse_input
virtual void parse_input(std::istream &input, const std::string &filename="input file", const std::string &last_line="", const bool skip_undefined=false)
Definition: parameter_handler.cc:399
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
RelaxationBlock
Definition: relaxation_block.h:59
Patterns::Double
Definition: patterns.h:293
int
std::sin
inline ::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5297
Patterns::Integer
Definition: patterns.h:190