Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-6.h
Go to the documentation of this file.
1
627 2)
628 *   , dof_handler(triangulation)
629 *   {}
630 *  
631 *  
632 *  
633 * @endcode
634 *
635 *
636 * <a name="step_6-Step6setup_system"></a>
637 * <h4>Step6::setup_system</h4>
638 *
639
640 *
641 * The next function sets up all the variables that describe the linear
642 * finite element problem, such as the DoFHandler, matrices, and
643 * vectors. The difference to what we did in @ref step_5 "step-5" is only that we now also
644 * have to take care of hanging node constraints. These constraints are
645 * handled almost exclusively by the library, i.e. you only need to know
646 * that they exist and how to get them, but you do not have to know how they
647 * are formed or what exactly is done with them.
648 *
649
650 *
651 * At the beginning of the function, you find all the things that are the same
652 * as in @ref step_5 "step-5": setting up the degrees of freedom (this time we have
653 * quadratic elements, but there is no difference from a user code perspective
654 * to the linear -- or any other degree, for that matter -- case), generating
655 * the sparsity pattern, and initializing the solution and right hand side
656 * vectors. Note that the sparsity pattern will have significantly more
657 * entries per row now, since there are now 9 degrees of freedom per cell
658 * (rather than only four), that can couple with each other.
659 *
660 * @code
661 *   template <int dim>
662 *   void Step6<dim>::setup_system()
663 *   {
664 *   dof_handler.distribute_dofs(fe);
665 *  
666 *   solution.reinit(dof_handler.n_dofs());
667 *   system_rhs.reinit(dof_handler.n_dofs());
668 *  
669 * @endcode
670 *
671 * We may now populate the AffineConstraints object with the hanging node
672 * constraints. Since we will call this function in a loop we first clear
673 * the current set of constraints from the last system and then compute new
674 * ones:
675 *
676 * @code
677 *   constraints.clear();
678 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
679 *  
680 *  
681 * @endcode
682 *
683 * Now we are ready to interpolate the boundary values with indicator 0 (the
684 * whole boundary) and store the resulting constraints in our
685 * <code>constraints</code> object. Note that we do not to apply the
686 * boundary conditions after assembly, like we did in earlier steps: instead
687 * we put all constraints on our function space in the AffineConstraints
688 * object. We can add constraints to the AffineConstraints object in either
689 * order: if two constraints conflict then the constraint matrix either abort
690 * or throw an exception via the Assert macro.
691 *
692 * @code
694 *   types::boundary_id(0),
695 *   Functions::ZeroFunction<dim>(),
696 *   constraints);
697 *  
698 * @endcode
699 *
700 * After all constraints have been added, they need to be sorted and
701 * rearranged to perform some actions more efficiently. This postprocessing
702 * is done using the <code>close()</code> function, after which no further
703 * constraints may be added any more:
704 *
705 * @code
706 *   constraints.close();
707 *  
708 * @endcode
709 *
710 * Now we first build our compressed sparsity pattern like we did in the
711 * previous examples. Nevertheless, we do not copy it to the final sparsity
712 * pattern immediately. Note that we call a variant of
713 * make_sparsity_pattern that takes the AffineConstraints object as the third
714 * argument. We are letting the routine know that we will never write into
715 * the locations given by <code>constraints</code> by setting the argument
716 * <code>keep_constrained_dofs</code> to false (in other words, that we will
717 * never write into entries of the matrix that correspond to constrained
718 * degrees of freedom). If we were to condense the
719 * constraints after assembling, we would have to pass <code>true</code>
720 * instead because then we would first write into these locations only to
721 * later set them to zero again during condensation.
722 *
723 * @code
724 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
725 *   DoFTools::make_sparsity_pattern(dof_handler,
726 *   dsp,
727 *   constraints,
728 *   /*keep_constrained_dofs = */ false);
729 *  
730 * @endcode
731 *
732 * Now all non-zero entries of the matrix are known (i.e. those from
733 * regularly assembling the matrix and those that were introduced by
734 * eliminating constraints). We may copy our intermediate object to the
735 * sparsity pattern:
736 *
737 * @code
738 *   sparsity_pattern.copy_from(dsp);
739 *  
740 * @endcode
741 *
742 * We may now, finally, initialize the sparse matrix:
743 *
744 * @code
745 *   system_matrix.reinit(sparsity_pattern);
746 *   }
747 *  
748 *  
749 * @endcode
750 *
751 *
752 * <a name="step_6-Step6assemble_system"></a>
753 * <h4>Step6::assemble_system</h4>
754 *
755
756 *
757 * Next, we have to assemble the matrix. However, to copy the local matrix and
758 * vector on each cell into the global system, we are no longer using a
759 * hand-written loop. Instead, we use
760 * AffineConstraints::distribute_local_to_global() that internally executes
761 * this loop while performing Gaussian elimination on rows and columns
762 * corresponding to constrained degrees on freedom.
763 *
764
765 *
766 * The rest of the code that forms the local contributions remains
767 * unchanged. It is worth noting, however, that under the hood several things
768 * are different than before. First, the variable <code>dofs_per_cell</code>
769 * and return value of <code>quadrature_formula.size()</code> now are 9 each,
770 * where they were 4 before. Introducing such variables as abbreviations is a
771 * good strategy to make code work with different elements without having to
772 * change too much code. Secondly, the <code>fe_values</code> object of course
773 * needs to do other things as well, since the shape functions are now
774 * quadratic, rather than linear, in each coordinate variable. Again, however,
775 * this is something that is completely handled by the library.
776 *
777 * @code
778 *   template <int dim>
779 *   void Step6<dim>::assemble_system()
780 *   {
781 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
782 *  
783 *   FEValues<dim> fe_values(fe,
784 *   quadrature_formula,
787 *  
788 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
789 *  
790 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
791 *   Vector<double> cell_rhs(dofs_per_cell);
792 *  
793 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
794 *  
795 *   for (const auto &cell : dof_handler.active_cell_iterators())
796 *   {
797 *   fe_values.reinit(cell);
798 *  
799 *   cell_matrix = 0;
800 *   cell_rhs = 0;
801 *  
802 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
803 *   {
804 *   const double current_coefficient =
805 *   coefficient(fe_values.quadrature_point(q_index));
806 *   for (const unsigned int i : fe_values.dof_indices())
807 *   {
808 *   for (const unsigned int j : fe_values.dof_indices())
809 *   cell_matrix(i, j) +=
810 *   (current_coefficient * // a(x_q)
811 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
812 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
813 *   fe_values.JxW(q_index)); // dx
814 *  
815 *   cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
816 *   1.0 * // f(x)
817 *   fe_values.JxW(q_index)); // dx
818 *   }
819 *   }
820 *  
821 * @endcode
822 *
823 * Finally, transfer the contributions from @p cell_matrix and
824 * @p cell_rhs into the global objects.
825 *
826 * @code
827 *   cell->get_dof_indices(local_dof_indices);
828 *   constraints.distribute_local_to_global(
829 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
830 *   }
831 * @endcode
832 *
833 * Now we are done assembling the linear system. The constraint matrix took
834 * care of applying the boundary conditions and also eliminated hanging node
835 * constraints. The constrained nodes are still in the linear system (there
836 * is a nonzero entry, chosen in a way that the matrix is well conditioned,
837 * on the diagonal of the matrix and all other entries for this line are set
838 * to zero) but the computed values are invalid (i.e., the corresponding
839 * entries in <code>system_rhs</code> are currently meaningless). We compute
840 * the correct values for these nodes at the end of the <code>solve</code>
841 * function.
842 *
843 * @code
844 *   }
845 *  
846 *  
847 * @endcode
848 *
849 *
850 * <a name="step_6-Step6solve"></a>
851 * <h4>Step6::solve</h4>
852 *
853
854 *
855 * We continue with gradual improvements. The function that solves the linear
856 * system again uses the SSOR preconditioner, and is again unchanged except
857 * that we have to incorporate hanging node constraints. As mentioned above,
858 * the degrees of freedom from the AffineConstraints object corresponding to
859 * hanging node constraints and boundary values have been removed from the
860 * linear system by giving the rows and columns of the matrix a special
861 * treatment. This way, the values for these degrees of freedom have wrong,
862 * but well-defined values after solving the linear system. What we then have
863 * to do is to use the constraints to assign to them the values that they
864 * should have. This process, called <code>distributing</code> constraints,
865 * computes the values of constrained nodes from the values of the
866 * unconstrained ones, and requires only a single additional function call
867 * that you find at the end of this function:
868 *
869
870 *
871 *
872 * @code
873 *   template <int dim>
874 *   void Step6<dim>::solve()
875 *   {
876 *   SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
877 *   SolverCG<Vector<double>> solver(solver_control);
878 *  
879 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
880 *   preconditioner.initialize(system_matrix, 1.2);
881 *  
882 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
883 *  
884 *   constraints.distribute(solution);
885 *   }
886 *  
887 *  
888 * @endcode
889 *
890 *
891 * <a name="step_6-Step6refine_grid"></a>
892 * <h4>Step6::refine_grid</h4>
893 *
894
895 *
896 * We use a sophisticated error estimation scheme to refine the mesh instead
897 * of global refinement. We will use the KellyErrorEstimator class which
898 * implements an error estimator for the Laplace equation; it can in principle
899 * handle variable coefficients, but we will not use these advanced features,
900 * but rather use its most simple form since we are not interested in
901 * quantitative results but only in a quick way to generate locally refined
902 * grids.
903 *
904
905 *
906 * Although the error estimator derived by Kelly et al. was originally
907 * developed for the Laplace equation, we have found that it is also well
908 * suited to quickly generate locally refined grids for a wide class of
909 * problems. This error estimator uses the solution gradient's jump at
910 * cell faces (which is a measure for the second derivatives) and
911 * scales it by the size of the cell. It is therefore a measure for the local
912 * smoothness of the solution at the place of each cell and it is thus
913 * understandable that it yields reasonable grids also for hyperbolic
914 * transport problems or the wave equation as well, although these grids are
915 * certainly suboptimal compared to approaches specially tailored to the
916 * problem. This error estimator may therefore be understood as a quick way to
917 * test an adaptive program.
918 *
919
920 *
921 * The way the estimator works is to take a <code>DoFHandler</code> object
922 * describing the degrees of freedom and a vector of values for each degree of
923 * freedom as input and compute a single indicator value for each active cell
924 * of the triangulation (i.e. one value for each of the active cells). To do
925 * so, it needs two additional pieces of information: a face quadrature formula,
926 * i.e., a quadrature formula on <code>dim-1</code> dimensional objects. We use
927 * a 3-point Gauss rule again, a choice that is consistent and appropriate with
928 * the bi-quadratic finite element shape functions in this program.
929 * (What constitutes a suitable quadrature rule here of course depends on
930 * knowledge of the way the error estimator evaluates the solution field. As
931 * said above, the jump of the gradient is integrated over each face, which
932 * would be a quadratic function on each face for the quadratic elements in
933 * use in this example. In fact, however, it is the square of the jump of the
934 * gradient, as explained in the documentation of that class, and that is a
935 * quartic function, for which a 3 point Gauss formula is sufficient since it
936 * integrates polynomials up to order 5 exactly.)
937 *
938
939 *
940 * Secondly, the function wants a list of boundary indicators for those
941 * boundaries where we have imposed Neumann values of the kind
942 * @f$\partial_n u(\mathbf x) = h(\mathbf x)@f$, along with a function @f$h(\mathbf
943 * x)@f$ for each such boundary. This information is represented by a map from
944 * boundary indicators to function objects describing the Neumann boundary
945 * values. In the present example program, we do not use Neumann boundary
946 * values, so this map is empty, and in fact constructed using the default
947 * constructor of the map in the place where the function call expects the
948 * respective function argument.
949 *
950
951 *
952 * The output is a vector of values for all active cells. While it may
953 * make sense to compute the <b>value</b> of a solution degree of freedom
954 * very accurately, it is usually not necessary to compute the <b>error
955 * indicator</b> corresponding to the solution on a cell particularly
956 * accurately. We therefore typically use a vector of floats instead of a vector
957 * of doubles to represent error indicators.
958 *
959 * @code
960 *   template <int dim>
961 *   void Step6<dim>::refine_grid()
962 *   {
963 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
964 *  
965 *   KellyErrorEstimator<dim>::estimate(dof_handler,
966 *   QGauss<dim - 1>(fe.degree + 1),
967 *   {},
968 *   solution,
969 *   estimated_error_per_cell);
970 *  
971 * @endcode
972 *
973 * The above function returned one error indicator value for each cell in
974 * the <code>estimated_error_per_cell</code> array. Refinement is now done
975 * as follows: refine those 30 per cent of the cells with the highest error
976 * values, and coarsen the 3 per cent of cells with the lowest values.
977 *
978
979 *
980 * One can easily verify that if the second number were zero, this would
981 * approximately result in a doubling of cells in each step in two space
982 * dimensions, since for each of the 30 per cent of cells, four new would be
983 * replaced, while the remaining 70 per cent of cells remain untouched. In
984 * practice, some more cells are usually produced since it is disallowed
985 * that a cell is refined twice while the neighbor cell is not refined; in
986 * that case, the neighbor cell would be refined as well.
987 *
988
989 *
990 * In many applications, the number of cells to be coarsened would be set to
991 * something larger than only three per cent. A non-zero value is useful
992 * especially if for some reason the initial (coarse) grid is already rather
993 * refined. In that case, it might be necessary to refine it in some
994 * regions, while coarsening in some other regions is useful. In our case
995 * here, the initial grid is very coarse, so coarsening is only necessary in
996 * a few regions where over-refinement may have taken place. Thus a small,
997 * non-zero value is appropriate here.
998 *
999
1000 *
1001 * The following function now takes these refinement indicators and flags
1002 * some cells of the triangulation for refinement or coarsening using the
1003 * method described above. It is from a class that implements several
1004 * different algorithms to refine a triangulation based on cell-wise error
1005 * indicators.
1006 *
1007 * @code
1008 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1009 *   estimated_error_per_cell,
1010 *   0.3,
1011 *   0.03);
1012 *  
1013 * @endcode
1014 *
1015 * After the previous function has exited, some cells are flagged for
1016 * refinement, and some other for coarsening. The refinement or coarsening
1017 * itself is not performed by now, however, since there are cases where
1018 * further modifications of these flags is useful. Here, we don't want to do
1019 * any such thing, so we can tell the triangulation to perform the actions
1020 * for which the cells are flagged:
1021 *
1022 * @code
1024 *   }
1025 *  
1026 *  
1027 * @endcode
1028 *
1029 *
1030 * <a name="step_6-Step6output_results"></a>
1031 * <h4>Step6::output_results</h4>
1032 *
1033
1034 *
1035 * At the end of computations on each grid, and just before we continue the
1036 * next cycle with mesh refinement, we want to output the results from this
1037 * cycle.
1038 *
1039
1040 *
1041 * We have already seen in @ref step_1 "step-1" how this can be achieved for the
1042 * mesh itself. Here, we change a few things:
1043 * <ol>
1044 * <li>We use two different formats: gnuplot and VTU.</li>
1045 * <li>We embed the cycle number in the output file name.</li>
1046 * <li>For gnuplot output, we set up a GridOutFlags::Gnuplot object to
1047 * provide a few extra visualization arguments so that edges appear
1048 * curved. This is explained in further detail in @ref step_10 "step-10".</li>
1049 * </ol>
1050 *
1051 * @code
1052 *   template <int dim>
1053 *   void Step6<dim>::output_results(const unsigned int cycle) const
1054 *   {
1055 *   {
1056 *   GridOut grid_out;
1057 *   std::ofstream output("grid-" + std::to_string(cycle) + ".gnuplot");
1058 *   GridOutFlags::Gnuplot gnuplot_flags(false, 5);
1059 *   grid_out.set_flags(gnuplot_flags);
1060 *   MappingQ<dim> mapping(3);
1061 *   grid_out.write_gnuplot(triangulation, output, &mapping);
1062 *   }
1063 *  
1064 *   {
1065 *   DataOut<dim> data_out;
1066 *   data_out.attach_dof_handler(dof_handler);
1067 *   data_out.add_data_vector(solution, "solution");
1068 *   data_out.build_patches();
1069 *  
1070 *   std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
1071 *   data_out.write_vtu(output);
1072 *   }
1073 *   }
1074 *  
1075 *  
1076 * @endcode
1077 *
1078 *
1079 * <a name="step_6-Step6run"></a>
1080 * <h4>Step6::run</h4>
1081 *
1082
1083 *
1084 * The final function before <code>main()</code> is again the main driver of
1085 * the class, <code>run()</code>. It is similar to the one of @ref step_5 "step-5", except
1086 * that we generate a file in the program again instead of reading it from
1087 * disk, in that we adaptively instead of globally refine the mesh, and that
1088 * we output the solution on the final mesh in the present function.
1089 *
1090
1091 *
1092 * The first block in the main loop of the function deals with mesh generation.
1093 * If this is the first cycle of the program, instead of reading the grid from
1094 * a file on disk as in the previous example, we now again create it using a
1095 * library function. The domain is again a circle with center at the origin and
1096 * a radius of one (these are the two hidden arguments to the function, which
1097 * have default values).
1098 *
1099
1100 *
1101 * You will notice by looking at the coarse grid that it is of inferior
1102 * quality than the one which we read from the file in the previous example:
1103 * the cells are less equally formed. However, using the library function this
1104 * program works in any space dimension, which was not the case before.
1105 *
1106
1107 *
1108 * In case we find that this is not the first cycle, we want to refine the
1109 * grid. Unlike the global refinement employed in the last example program, we
1110 * now use the adaptive procedure described above.
1111 *
1112
1113 *
1114 * The rest of the loop looks as before:
1115 *
1116 * @code
1117 *   template <int dim>
1118 *   void Step6<dim>::run()
1119 *   {
1120 *   for (unsigned int cycle = 0; cycle < 8; ++cycle)
1121 *   {
1122 *   std::cout << "Cycle " << cycle << ':' << std::endl;
1123 *  
1124 *   if (cycle == 0)
1125 *   {
1128 *   }
1129 *   else
1130 *   refine_grid();
1131 *  
1132 *  
1133 *   std::cout << " Number of active cells: "
1134 *   << triangulation.n_active_cells() << std::endl;
1135 *  
1136 *   setup_system();
1137 *  
1138 *   std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1139 *   << std::endl;
1140 *  
1141 *   assemble_system();
1142 *   solve();
1143 *   output_results(cycle);
1144 *   }
1145 *   }
1146 *  
1147 *  
1148 * @endcode
1149 *
1150 *
1151 * <a name="step_6-Thecodemaincodefunction"></a>
1152 * <h3>The <code>main</code> function</h3>
1153 *
1154
1155 *
1156 * The main function is unaltered in its functionality from the previous
1157 * example, but we have taken a step of additional caution. Sometimes,
1158 * something goes wrong (such as insufficient disk space upon writing an
1159 * output file, not enough memory when trying to allocate a vector or a
1160 * matrix, or if we can't read from or write to a file for whatever reason),
1161 * and in these cases the library will throw exceptions. Since these are
1162 * run-time problems, not programming errors that can be fixed once and for
1163 * all, this kind of exceptions is not switched off in optimized mode, in
1164 * contrast to the <code>Assert</code> macro which we have used to test
1165 * against programming errors. If uncaught, these exceptions propagate the
1166 * call tree up to the <code>main</code> function, and if they are not caught
1167 * there either, the program is aborted. In many cases, like if there is not
1168 * enough memory or disk space, we can't do anything but we can at least print
1169 * some text trying to explain the reason why the program failed. A way to do
1170 * so is shown in the following. It is certainly useful to write any larger
1171 * program in this way, and you can do so by more or less copying this
1172 * function except for the <code>try</code> block that actually encodes the
1173 * functionality particular to the present application.
1174 *
1175 * @code
1176 *   int main()
1177 *   {
1178 * @endcode
1179 *
1180 * The general idea behind the layout of this function is as follows: let's
1181 * try to run the program as we did before...
1182 *
1183 * @code
1184 *   try
1185 *   {
1186 *   Step6<2> laplace_problem_2d;
1187 *   laplace_problem_2d.run();
1188 *   }
1189 * @endcode
1190 *
1191 * ...and if this should fail, try to gather as much information as
1192 * possible. Specifically, if the exception that was thrown is an object of
1193 * a class that is derived from the C++ standard class
1194 * <code>exception</code>, then we can use the <code>what</code> member
1195 * function to get a string which describes the reason why the exception was
1196 * thrown.
1197 *
1198
1199 *
1200 * The deal.II exception classes are all derived from the standard class,
1201 * and in particular, the <code>exc.what()</code> function will return
1202 * approximately the same string as would be generated if the exception was
1203 * thrown using the <code>Assert</code> macro. You have seen the output of
1204 * such an exception in the previous example, and you then know that it
1205 * contains the file and line number of where the exception occurred, and
1206 * some other information. This is also what the following statements would
1207 * print.
1208 *
1209
1210 *
1211 * Apart from this, there isn't much that we can do except exiting the
1212 * program with an error code (this is what the <code>return 1;</code>
1213 * does):
1214 *
1215 * @code
1216 *   catch (std::exception &exc)
1217 *   {
1218 *   std::cerr << std::endl
1219 *   << std::endl
1220 *   << "----------------------------------------------------"
1221 *   << std::endl;
1222 *   std::cerr << "Exception on processing: " << std::endl
1223 *   << exc.what() << std::endl
1224 *   << "Aborting!" << std::endl
1225 *   << "----------------------------------------------------"
1226 *   << std::endl;
1227 *  
1228 *   return 1;
1229 *   }
1230 * @endcode
1231 *
1232 * If the exception that was thrown somewhere was not an object of a class
1233 * derived from the standard <code>exception</code> class, then we can't do
1234 * anything at all. We then simply print an error message and exit.
1235 *
1236 * @code
1237 *   catch (...)
1238 *   {
1239 *   std::cerr << std::endl
1240 *   << std::endl
1241 *   << "----------------------------------------------------"
1242 *   << std::endl;
1243 *   std::cerr << "Unknown exception!" << std::endl
1244 *   << "Aborting!" << std::endl
1245 *   << "----------------------------------------------------"
1246 *   << std::endl;
1247 *   return 1;
1248 *   }
1249 *  
1250 * @endcode
1251 *
1252 * If we got to this point, there was no exception which propagated up to
1253 * the main function (there may have been exceptions, but they were caught
1254 * somewhere in the program or the library). Therefore, the program
1255 * performed as was expected and we can return without error.
1256 *
1257 * @code
1258 *   return 0;
1259 *   }
1260 * @endcode
1261<a name="step_6-Results"></a><h1>Results</h1>
1262
1263
1264
1265The output of the program looks as follows:
1266@code
1267Cycle 0:
1268 Number of active cells: 20
1269 Number of degrees of freedom: 89
1270Cycle 1:
1271 Number of active cells: 38
1272 Number of degrees of freedom: 183
1273Cycle 2:
1274 Number of active cells: 71
1275 Number of degrees of freedom: 325
1276Cycle 3:
1277 Number of active cells: 146
1278 Number of degrees of freedom: 689
1279Cycle 4:
1280 Number of active cells: 284
1281 Number of degrees of freedom: 1373
1282Cycle 5:
1283 Number of active cells: 599
1284 Number of degrees of freedom: 2769
1285Cycle 6:
1286 Number of active cells: 1241
1287 Number of degrees of freedom: 5833
1288Cycle 7:
1289 Number of active cells: 2507
1290 Number of degrees of freedom: 11916
1291@endcode
1292
1293
1294
1295As intended, the number of cells roughly doubles in each cycle. The
1296number of degrees is slightly more than four times the number of
1297cells; one would expect a factor of exactly four in two spatial
1298dimensions on an infinite grid (since the spacing between the degrees
1299of freedom is half the cell width: one additional degree of freedom on
1300each edge and one in the middle of each cell), but it is larger than
1301that factor due to the finite size of the mesh and due to additional
1302degrees of freedom which are introduced by hanging nodes and local
1303refinement.
1304
1305
1306
1307The program outputs the solution and mesh in each cycle of the
1308refinement loop. The solution looks as follows:
1309
1310<img src="https://www.dealii.org/images/steps/developer/step-6.solution.9.2.png" alt="">
1311
1312It is interesting to follow how the program arrives at the final mesh:
1313
1314<div class="twocolumn" style="width: 80%">
1315 <div>
1316 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_0.svg"
1317 alt="Initial grid: the five-cell circle grid with one global refinement."
1318 width="300" height="300">
1319 </div>
1320 <div>
1321 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_1.svg"
1322 alt="First grid: the five-cell circle grid with two global refinements."
1323 width="300" height="300">
1324 </div>
1325 <div>
1326 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_2.svg"
1327 alt="Second grid: the five-cell circle grid with one adaptive refinement."
1328 width="300" height="300">
1329 </div>
1330 <div>
1331 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_3.svg"
1332 alt="Third grid: the five-cell circle grid with two adaptive
1333 refinements, showing clustering around the inner circle."
1334 width="300" height="300">
1335 </div>
1336 <div>
1337 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_4.svg"
1338 alt="Fourth grid: the five-cell circle grid with three adaptive
1339 refinements, showing clustering around the inner circle."
1340 width="300" height="300">
1341 </div>
1342 <div>
1343 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_5.svg"
1344 alt="Fifth grid: the five-cell circle grid with four adaptive
1345 refinements, showing clustering around the inner circle."
1346 width="300" height="300">
1347 </div>
1348 <div>
1349 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_6.svg"
1350 alt="Sixth grid: the five-cell circle grid with five adaptive
1351 refinements, showing clustering around the inner circle."
1352 width="300" height="300">
1353 </div>
1354 <div>
1355 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_7.svg"
1356 alt="Last grid: the five-cell circle grid with six adaptive
1357 refinements, showing that most cells are clustered around the inner circle."
1358 width="300" height="300">
1359 </div>
1360</div>
1361
1362
1363It is clearly visible that the region where the solution has a kink,
1364i.e. the circle at radial distance 0.5 from the center, is
1365refined most. Furthermore, the central region where the solution is
1366very smooth and almost flat, is almost not refined at all, but this
1367results from the fact that we did not take into account that the
1368coefficient is large there. The region outside is refined rather
1369arbitrarily, since the second derivative is constant there and refinement
1370is therefore mostly based on the size of the cells and their deviation
1371from the optimal square.
1372
1373
1374
1375<a name="step-6-extensions"></a>
1376<a name="step_6-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1377
1378
1379<a name="step_6-Solversandpreconditioners"></a><h4>Solvers and preconditioners</h4>
1380
1381
1382One thing that is always worth playing around with if one solves
1383problems of appreciable size (much bigger than the one we have here)
1384is to try different solvers or preconditioners. In the current case,
1385the linear system is symmetric and positive definite, which makes the
1386Conjugate Gradient (CG) algorithm pretty much the canonical choice for solving. However,
1387the SSOR preconditioner we use in the <code>solve()</code> function is
1388up for grabs.
1389
1390In deal.II, it is relatively simple to change the preconditioner.
1391Several simple preconditioner choices are accessible
1392by changing the following two lines
1393@code
1394 PreconditionSSOR<SparseMatrix<double>> preconditioner;
1395 preconditioner.initialize(system_matrix, 1.2);
1396@endcode
1397of code in the program. For example, switching this into
1398@code
1399 PreconditionSSOR<SparseMatrix<double>> preconditioner;
1400 preconditioner.initialize(system_matrix, 1.0);
1401@endcode
1402allows us to try out a different relaxation parameter for SSOR (1.0 instead
1403of the 1.2 in the original program). Similarly, by using
1404@code
1405 PreconditionJacobi<SparseMatrix<double>> preconditioner;
1406 preconditioner.initialize(system_matrix);
1407@endcode
1408we can use Jacobi as a preconditioner. And by using
1409@code
1410 SparseILU<double> preconditioner;
1411 preconditioner.initialize(system_matrix);
1412@endcode
1413we can use a simple incomplete LU decomposition without any thresholding or
1414strengthening of the diagonal (to use this preconditioner, you have to also
1415add the header file <code>deal.II/lac/sparse_ilu.h</code> to the include list
1416at the top of the file).
1417
1418Using these various different preconditioners, we can compare the
1419number of CG iterations needed (available through the
1420<code>solver_control.last_step()</code> call, see
1421@ref step_4 "step-4") as well as CPU time needed (using the Timer class,
1422discussed, for example, in @ref step_28 "step-28") and get the
1423following results (left: iterations; right: CPU time):
1424
1425<table width="60%" align="center">
1426 <tr>
1427 <td align="center">
1428 <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_iterations.png" alt="">
1429 </td>
1430 <td align="center">
1431 <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_time.png" alt="">
1432 </td>
1433 </tr>
1434</table>
1435
1436As we can see, all preconditioners behave pretty much the same on this
1437simple problem, with the number of iterations growing like @f${\cal
1438O}(N^{1/2})@f$ and because each iteration requires around @f${\cal
1439O}(N)@f$ operations the total CPU time grows like @f${\cal
1440O}(N^{3/2})@f$ (for the few smallest meshes, the CPU time is so small
1441that it doesn't record). Note that even though it is the simplest
1442method, Jacobi is the fastest for this problem.
1443
1444The situation changes slightly when the finite element is not a
1445bi-quadratic one (i.e., polynomial degree two) as selected in the
1446constructor of this program, but a bi-linear one (polynomial degree one).
1447If one makes this change, the results are as follows:
1448
1449<table width="60%" align="center">
1450 <tr>
1451 <td align="center">
1452 <img src="https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_iterations.png" alt="">
1453 </td>
1454 <td align="center">
1455 <img src="https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_time.png" alt="">
1456 </td>
1457 </tr>
1458</table>
1459
1460In other words, while the increase in iterations and CPU time is as
1461before, Jacobi is now the method that requires the most iterations; it
1462is still the fastest one, however, owing to the simplicity of the
1463operations it has to perform. This is not to say that Jacobi
1464is actually a good preconditioner -- for problems of appreciable size, it is
1465definitely not, and other methods will be substantially better -- but really
1466only that it is fast because its implementation is so simple that it can
1467compensate for a larger number of iterations.
1468
1469The message to take away from this is not that simplicity in
1470preconditioners is always best. While this may be true for the current
1471problem, it definitely is not once we move to more complicated
1472problems (elasticity or Stokes, for examples @ref step_8 "step-8" or
1473@ref step_22 "step-22"). Secondly, all of these preconditioners still
1474lead to an increase in the number of iterations as the number @f$N@f$ of
1475degrees of freedom grows, for example @f${\cal O}(N^\alpha)@f$; this, in
1476turn, leads to a total growth in effort as @f${\cal O}(N^{1+\alpha})@f$
1477since each iteration takes @f${\cal O}(N)@f$ work. This behavior is
1478undesirable: we would really like to solve linear systems with @f$N@f$
1479unknowns in a total of @f${\cal O}(N)@f$ work; there is a class
1480of preconditioners that can achieve this, namely geometric (@ref step_16 "step-16",
1481@ref step_37 "step-37", @ref step_39 "step-39")
1482or algebraic multigrid (@ref step_31 "step-31", @ref step_40 "step-40", and several others)
1483preconditioners. They are, however, significantly more complex than
1484the preconditioners outlined above, and so we will leave their use
1485to these later tutorial programs. The point to make, however, is
1486that "real" finite element programs do not use the preconditioners
1487we mention above: These are simply shown for expository purposes.
1488
1489Finally, the last message to take
1490home is that when the data shown above was generated (in 2018), linear
1491systems with 100,000 unknowns are
1492easily solved on a desktop or laptop machine in about a second, making
1493the solution of relatively simple 2d problems even to very high
1494accuracy not that big a task as it used to be in the
1495past. At the same time, the situation for 3d problems continues to be
1496quite different: A uniform 2d mesh with 100,000 unknowns corresponds to
1497a grid with about @f$300 \times 300@f$ nodes; the corresponding 3d mesh has
1498@f$300 \times 300 \times 300@f$ nodes and 30 million unknowns. Because
1499finite element matrices in 3d have many more nonzero entries than in 2d,
1500solving these linear systems will not only take 300 times as much CPU
1501time, but substantially longer. In other words, achieving the
1502same resolution in 3d *is* quite a large problem, and solving it
1503within a reasonable amount of time *will* require much more work to
1504implement better linear solvers. As mentioned above, multigrid methods
1505and matrix-free methods (see, for example, @ref step_37 "step-37"), along with
1506parallelization (@ref step_40 "step-40") will be necessary, but are then also able
1507to comfortably solve such linear systems.
1508
1509
1510
1511<a name="step_6-Abettermesh"></a><h4>A better mesh</h4>
1512
1513
1514If you look at the meshes above, you will see even though the domain is the
1515unit disk, and the jump in the coefficient lies along a circle, the cells
1516that make up the mesh do not track this geometry well. The reason, already hinted
1517at in @ref step_1 "step-1", is that in the absence of other information,
1518the Triangulation class only sees a bunch of
1519coarse grid cells but has, of course, no real idea what kind of geometry they
1520might represent when looked at together. For this reason, we need to tell
1521the Triangulation what to do when a cell is refined: where should the new
1522vertices at the edge midpoints and the cell midpoint be located so that the
1523child cells better represent the desired geometry than the parent cell.
1524
1525To visualize what the triangulation actually knows about the geometry,
1526it is not enough to just output the location of vertices and draw a
1527straight line for each edge; instead, we have to output both interior
1528and boundary lines as multiple segments so that they look
1529curved. We can do this by making one change to the gnuplot part of
1530<code>output_results</code>:
1531@code
1532{
1533 GridOut grid_out;
1534 std::ofstream output("grid-" + std::to_string(cycle) + ".gnuplot");
1535 GridOutFlags::Gnuplot gnuplot_flags(false, 5, /*curved_interior_cells*/true);
1536 grid_out.set_flags(gnuplot_flags);
1537 MappingQ<dim> mapping(3);
1538 grid_out.write_gnuplot(triangulation, output, &mapping);
1539}
1540@endcode
1541
1542In the code above, we already do this for faces that sit at the boundary: this
1543happens automatically since we use GridGenerator::hyper_ball, which attaches a
1544SphericalManifold to the boundary of the domain. To make the mesh
1545<i>interior</i> also track a circular domain, we need to work a bit harder,
1546though. First, recall that our coarse mesh consists of a central square
1547cell and four cells around it. Now first consider what would happen if we
1548also attached the SphericalManifold object not only to the four exterior faces
1549but also the four cells at the perimeter as well as all of their faces. We can
1550do this by adding the following snippet (testing that the center of a cell is
1551larger than a small multiple, say one tenth, of the cell diameter away from
1552center of the mesh only fails for the central square of the mesh):
1553@code
1555// after GridGenerator::hyper_ball is called the Triangulation has
1556// a SphericalManifold with id 0. We can use it again on the interior.
1557const Point<dim> mesh_center;
1558for (const auto &cell : triangulation.active_cell_iterators())
1559 if (mesh_center.distance (cell->center()) > cell->diameter()/10)
1560 cell->set_all_manifold_ids(0);
1561
1563@endcode
1564
1565After a few global refinement steps, this would lead to a mesh of the following
1566kind:
1567
1568
1569 <div class="onecolumn" style="width: 80%">
1570 <div>
1571 <img src="https://www.dealii.org/images/steps/developer/step_6_bad_grid_4.svg"
1572 alt="Grid where some central cells are nearly triangular."
1573 width="300" height="300">
1574 </div>
1575 </div>
1576
1577This is not a good mesh: the central cell has been refined in such a way that
1578the children located in the four corners of the original central cell
1579<i>degenerate</i>: they all tend towards triangles as mesh refinement
1580continues. This means that the Jacobian matrix of the transformation from
1581reference cell to actual cell degenerates for these cells, and because
1582all error estimates for finite element solutions contain the norm of the
1583inverse of the Jacobian matrix, you will get very large errors on these
1584cells and, in the limit as mesh refinement, a loss of convergence order because
1585the cells in these corners become worse and worse under mesh refinement.
1586
1587So we need something smarter. To this end, consider the following solution
1588originally developed by Konstantin Ladutenko. We will use the following code:
1589@code
1591
1592const Point<dim> mesh_center;
1593const double core_radius = 1.0/5.0,
1594 inner_radius = 1.0/3.0;
1595
1596// Step 1: Shrink the inner cell
1597//
1598// We cannot get a circle out of the inner cell because of
1599// the degeneration problem mentioned above. Rather, shrink
1600// the inner cell to a core radius of 1/5 that stays
1601// sufficiently far away from the place where the
1602// coefficient will have a discontinuity and where we want
1603// to have cell interfaces that actually lie on a circle.
1604// We do this shrinking by just scaling the location of each
1605// of the vertices, given that the center of the circle is
1606// simply the origin of the coordinate system.
1607for (const auto &cell : triangulation.active_cell_iterators())
1608 if (mesh_center.distance(cell->center()) < 1e-5)
1609 {
1610 for (const auto v : cell->vertex_indices())
1611 cell->vertex(v) *= core_radius/mesh_center.distance(cell->vertex(v));
1612 }
1613
1614// Step 2: Refine all cells except the central one
1615for (const auto &cell : triangulation.active_cell_iterators())
1616 if (mesh_center.distance(cell->center()) >= 1e-5)
1617 cell->set_refine_flag();
1619
1620// Step 3: Resize the inner children of the outer cells
1621//
1622// The previous step replaced each of the four outer cells
1623// by its four children, but the radial distance at which we
1624// have intersected is not what we want to later refinement
1625// steps. Consequently, move the vertices that were just
1626// created in radial direction to a place where we need
1627// them.
1628for (const auto &cell : triangulation.active_cell_iterators())
1629 for (const auto v : cell->vertex_indices())
1630 {
1631 const double dist = mesh_center.distance(cell->vertex(v));
1632 if (dist > core_radius*1.0001 && dist < 0.9999)
1633 cell->vertex(v) *= inner_radius/dist;
1634 }
1635
1636// Step 4: Apply curved manifold description
1637//
1638// As discussed above, we can not expect to subdivide the
1639// inner four cells (or their faces) onto concentric rings,
1640// but we can do so for all other cells that are located
1641// outside the inner radius. To this end, we loop over all
1642// cells and determine whether it is in this zone. If it
1643// isn't, then we set the manifold description of the cell
1644// and all of its bounding faces to the one that describes
1645// the spherical manifold already introduced above and that
1646// will be used for all further mesh refinement.
1647for (const auto &cell : triangulation.active_cell_iterators())
1648 {
1649 bool is_in_inner_circle = false;
1650 for (const auto v : cell->vertex_indices())
1651 if (mesh_center.distance(cell->vertex(v)) < inner_radius)
1652 {
1653 is_in_inner_circle = true;
1654 break;
1655 }
1656
1657 if (is_in_inner_circle == false)
1658 // The Triangulation already has a SphericalManifold with
1659 // manifold id 0 (see the documentation of
1660 // GridGenerator::hyper_ball) so we just attach it to the outer
1661 // ring here:
1662 cell->set_all_manifold_ids(0);
1663 }
1664@endcode
1665
1666This code then generates the following, much better sequence of meshes:
1667
1668<div class="twocolumn" style="width: 80%">
1669 <div>
1670 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_0_ladutenko.svg"
1671 alt="Initial grid: the Ladutenko grid with one global refinement."
1672 width="300" height="300">
1673 </div>
1674 <div>
1675 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_1_ladutenko.svg"
1676 alt="First adaptively refined Ladutenko grid."
1677 width="300" height="300">
1678 </div>
1679 <div>
1680 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_2_ladutenko.svg"
1681 alt="Second adaptively refined Ladutenko grid."
1682 width="300" height="300">
1683 </div>
1684 <div>
1685 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_3_ladutenko.svg"
1686 alt="Third adaptively refined Ladutenko grid."
1687 width="300" height="300">
1688 </div>
1689 <div>
1690 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_4_ladutenko.svg"
1691 alt="Fourth adaptively refined Ladutenko grid. The cells are clustered
1692 along the inner circle."
1693 width="300" height="300">
1694 </div>
1695 <div>
1696 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_5_ladutenko.svg"
1697 alt="Fifth adaptively refined Ladutenko grid: the cells are clustered
1698 along the inner circle."
1699 width="300" height="300">
1700 </div>
1701</div>
1702
1703Creating good meshes, and in particular making them fit the geometry you
1704want, is a complex topic in itself. You can find much more on this in
1705@ref step_49 "step-49", @ref step_53 "step-53", and @ref step_54 "step-54", among other tutorial programs that cover
1706the issue. @ref step_65 "step-65" shows another, less manual way to achieve a mesh
1707well fit to the problem here.
1708Information on curved domains can also be found in the
1709documentation topic on @ref manifold "Manifold descriptions".
1710
1711Why does it make sense to choose a mesh that tracks the internal
1712interface? There are a number of reasons, but the most essential one
1713comes down to what we actually integrate in our bilinear
1714form. Conceptually, we want to integrate the term @f$A_{ij}^K=\int_K
1715a(\mathbf x) \nabla \varphi_i(\mathbf x) \nabla \varphi_j(\mathbf x) ; dx@f$ as the
1716contribution of cell @f$K@f$ to the matrix entry @f$A_{ij}@f$. We can not
1717compute it exactly and have to resort to quadrature. We know that
1718quadrature is accurate if the integrand is smooth. That is because
1719quadrature in essence computes a polynomial approximation to the
1720integrand that coincides with the integrand in the quadrature points,
1721and then computes the volume under this polynomial as an approximation
1722to the volume under the original integrand. This polynomial
1723interpolant is accurate if the integrand is smooth on a cell, but it
1724is usually rather inaccurate if the integrand is discontinuous on a
1725cell.
1726
1727Consequently, it is worthwhile to align cells in such a way that the
1728interfaces across which the coefficient is discontinuous are aligned
1729with cell interfaces. This way, the coefficient is constant on each
1730cell, following which the integrand will be smooth, and its polynomial
1731approximation and the quadrature approximation of the integral will
1732both be accurate. Note that such an alignment is common in many
1733practical cases, so deal.II provides a number of functions (such as
1734@ref GlossMaterialId "material_id") to help manage such a scenario.
1735Refer to @ref step_28 "step-28" and @ref step_46 "step-46" for examples of how material ids can be
1736applied.
1737
1738Finally, let us consider the case of a coefficient that has a smooth
1739and non-uniform distribution in space. We can repeat once again all of
1740the above discussion on the representation of such a function with the
1741quadrature. So, to simulate it accurately there are a few readily
1742available options: you could reduce the cell size, increase the order
1743of the polynomial used in the quadrature formula, select a more
1744appropriate quadrature formula, or perform a combination of these
1745steps. The key is that providing the best fit of the coefficient's
1746spatial dependence with the quadrature polynomial will lead to a more
1747accurate finite element solution of the PDE.
1748
1749As a final note: The discussion in the previous paragraphs shows, we here
1750have a very concrete way of stating what we think of a good mesh -- it should
1751be aligned with the jump in the coefficient. But one could also have asked
1752this kind of question in a more general setting: Given some equation with
1753a smooth solution and smooth coefficients, can we say what a good mesh
1754would look like? This is a question for which the answer is easier to state
1755in intuitive terms than mathematically: A good mesh has cells that all,
1756by and large, look like squares (or cubes, in 3d). A bad mesh would contain
1757cells that are very elongated in some directions or, more generally, for which
1758there are cells that have both short and long edges. There are many ways
1759in which one could assign a numerical quality index to each cell that measures
1760whether the cell is "good" or "bad"; some of these are often chosen because
1761they are cheap and easy to compute, whereas others are based on what enters
1762into proofs of convergence. An example of the former would be the ratio of
1763the longest to the shortest edge of a cell: In the ideal case, that ratio
1764would be one; bad cells have values much larger than one. An example of the
1765latter kind would consider the gradient (the "Jacobian") of the mapping
1766from the reference cell @f$\hat K=[0,1]^d@f$ to the real cell @f$K@f$; this
1767gradient is a matrix, and a quantity that enters into error estimates
1768is the maximum over all points on the reference cell of the ratio of the
1769largest to the smallest eigenvalue of this matrix. It is again not difficult
1770to see that this ratio is constant if the cell @f$K@f$ is an affine image of
1771@f$\hat K@f$, and that it is one for squares and cubes.
1772
1773In practice, it might be interesting to visualize such quality measures.
1774The function GridTools::compute_aspect_ratio_of_cells() provides one
1775way to get this kind of information. Even better, visualization tools
1776such as VisIt often allow you to visualize this sort of information
1777for a variety of measures from within the visualization software
1778itself; in the case of VisIt, just add a "pseudo-color" plot and select
1779one of the mesh quality measures instead of the solution field.
1780
1781
1782<a name="step_6-Playingwiththeregularityofthesolution"></a><h4>Playing with the regularity of the solution</h4>
1783
1784
1785From a mathematical perspective, solutions of the Laplace equation
1786@f[
1787 -\Delta u = f
1788@f]
1789on smoothly bounded, convex domains are known to be smooth themselves. The exact degree
1790of smoothness, i.e., the function space in which the solution lives, depends
1791on how smooth exactly the boundary of the domain is, and how smooth the right
1792hand side is. Some regularity of the solution may be lost at the boundary, but
1793one generally has that the solution is twice more differentiable in
1794compact subsets of the domain than the right hand side.
1795If, in particular, the right hand side satisfies @f$f\in C^\infty(\Omega)@f$, then
1796@f$u \in C^\infty(\Omega_i)@f$ where @f$\Omega_i@f$ is any compact subset of @f$\Omega@f$
1797(@f$\Omega@f$ is an open domain, so a compact subset needs to keep a positive distance
1798from @f$\partial\Omega@f$).
1799
1800The situation we chose for the current example is different, however: we look
1801at an equation with a non-constant coefficient @f$a(\mathbf x)@f$:
1802@f[
1803 -\nabla \cdot (a \nabla u) = f.
1804@f]
1805Here, if @f$a@f$ is not smooth, then the solution will not be smooth either,
1806regardless of @f$f@f$. In particular, we expect that wherever @f$a@f$ is discontinuous
1807along a line (or along a plane in 3d),
1808the solution will have a kink. This is easy to see: if for example @f$f@f$
1809is continuous, then @f$f=-\nabla \cdot (a \nabla u)@f$ needs to be
1810continuous. This means that @f$a \nabla u@f$ must be continuously differentiable
1811(not have a kink). Consequently, if @f$a@f$ has a discontinuity, then @f$\nabla u@f$
1812must have an opposite discontinuity so that the two exactly cancel and their
1813product yields a function without a discontinuity. But for @f$\nabla u@f$ to have
1814a discontinuity, @f$u@f$ must have a kink. This is of course exactly what is
1815happening in the current example, and easy to observe in the pictures of the
1816solution.
1817
1818In general, if the coefficient @f$a(\mathbf x)@f$ is discontinuous along a line in 2d,
1819or a plane in 3d, then the solution may have a kink, but the gradient of the
1820solution will not go to infinity. That means, that the solution is at least
1821still in the <a href="https://en.wikipedia.org/wiki/Sobolev_space">Sobolev space</a>
1822@f$W^{1,\infty}@f$ (i.e., roughly speaking, in the
1823space of functions whose derivatives are bounded). On the other hand,
1824we know that in the most
1825extreme cases -- i.e., where the domain has reentrant corners, the
1826right hand side only satisfies @f$f\in H^{-1}@f$, or the coefficient @f$a@f$ is only in
1827@f$L^\infty@f$ -- all we can expect is that @f$u\in H^1@f$ (i.e., the
1828<a
1829href="https://en.wikipedia.org/wiki/Sobolev_space#Sobolev_spaces_with_integer_k">Sobolev
1830space</a> of functions whose derivative is square integrable), a much larger space than
1831@f$W^{1,\infty}@f$. It is not very difficult to create cases where
1832the solution is in a space @f$H^{1+s}@f$ where we can get @f$s@f$ to become as small
1833as we want. Such cases are often used to test adaptive finite element
1834methods because the mesh will have to resolve the singularity that causes
1835the solution to not be in @f$W^{1,\infty}@f$ any more.
1836
1837The typical example one uses for this is called the <i>Kellogg problem</i>
1838(referring to @cite Kel74), which in the commonly used form has a coefficient
1839@f$a(\mathbf x)@f$ that has different values in the four quadrants of the plane
1840(or eight different values in the octants of @f${\mathbb R}^3@f$). The exact degree
1841of regularity (the @f$s@f$ in the index of the Sobolev space above) depends on the
1842values of @f$a(\mathbf x)@f$ coming together at the origin, and by choosing the
1843jumps large enough, the regularity of the solution can be made as close as
1844desired to @f$H^1@f$.
1845
1846To implement something like this, one could replace the coefficient
1847function by the following (shown here only for the 2d case):
1848@code
1849template <int dim>
1850double coefficient (const Point<dim> &p)
1851{
1852 if ((p[0] < 0) && (p[1] < 0)) // lower left quadrant
1853 return 1;
1854 else if ((p[0] >= 0) && (p[1] < 0)) // lower right quadrant
1855 return 10;
1856 else if ((p[0] < 0) && (p[1] >= 0)) // upper left quadrant
1857 return 100;
1858 else if ((p[0] >= 0) && (p[1] >= 0)) // upper right quadrant
1859 return 1000;
1860 else
1861 DEAL_II_ASSERT_UNREACHABLE();
1862}
1863@endcode
1864(Adding the `DEAL_II_ASSERT_UNREACHABLE();` call at the end ensures that either an exception
1865is thrown or that the program aborts if we ever get to that point
1866-- which of course we shouldn't,
1867but this is a good way to insure yourself: We all make mistakes by
1868sometimes not thinking of all cases, for example by checking
1869for <code>p[0]</code> to be less than and greater than zero,
1870rather than greater-or-equal to zero, and thereby forgetting
1871some cases that would otherwise lead to bugs that are awkward
1872to find.
1873
1874By playing with such cases where four or more sectors come
1875together and on which the coefficient has different values, one can
1876construct cases where the solution has singularities at the
1877origin. One can also see how the meshes are refined in such cases.
1878 *
1879 *
1880<a name="step_6-PlainProg"></a>
1881<h1> The plain program</h1>
1882@include "step-6.cc"
1883*/
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
void set_flags(const GridOutFlags::DX &flags)
Definition grid_out.cc:471
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:4608
Definition point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
Point< 3 > center
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int vertex_indices[2]
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
double volume(const Triangulation< dim, spacedim > &tria)
double diameter(const Triangulation< dim, spacedim > &tria)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
void free(T *&pointer)
Definition cuda.h:96
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< 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={})
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void abort(const ExceptionBase &exc) noexcept
void copy(const T *begin, const T *end, U *dest)
int(& functions)(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
STL namespace.
Definition types.h:32
unsigned int boundary_id
Definition types.h:144
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation