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