717 * Now all non-zero entries of the
matrix are known (i.e. those from
718 * regularly assembling the matrix and those that were introduced by
719 * eliminating constraints). We may
copy our intermediate
object to the
723 * sparsity_pattern.copy_from(dsp);
727 * We may now,
finally, initialize the sparse
matrix:
730 * system_matrix.reinit(sparsity_pattern);
737 * <a name=
"Step6assemble_system"></a>
738 * <h4>Step6::assemble_system</h4>
743 * vector on each cell into the global system, we are no longer
using a
744 * hand-written
loop. Instead, we use
746 * this loop while performing Gaussian elimination on rows and columns
747 * corresponding to constrained degrees on freedom.
751 * The rest of the code that forms the local contributions remains
752 * unchanged. It is worth noting, however, that under the hood several things
753 * are different than before. First, the variable <code>dofs_per_cell</code>
754 * and return value of <code>quadrature_formula.size()</code> now are 9 each,
755 * where they were 4 before. Introducing such variables as abbreviations is a
756 * good strategy to make code work with different elements without having to
757 * change too much code. Secondly, the <code>fe_values</code>
object of course
758 * needs to do other things as well, since the shape functions are now
759 * quadratic, rather than linear, in each coordinate variable. Again, however,
760 * this is something that is completely handled by the library.
764 *
void Step6<dim>::assemble_system()
766 *
const QGauss<dim> quadrature_formula(fe.degree + 1);
769 * quadrature_formula,
773 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
778 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
780 *
for (
const auto &cell : dof_handler.active_cell_iterators())
785 * fe_values.reinit(cell);
787 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
789 * const double current_coefficient =
790 * coefficient(fe_values.quadrature_point(q_index));
791 *
for (
const unsigned int i : fe_values.dof_indices())
793 * for (const unsigned
int j : fe_values.dof_indices())
795 * (current_coefficient *
796 * fe_values.shape_grad(i, q_index) *
797 * fe_values.shape_grad(j, q_index) *
798 * fe_values.JxW(q_index));
800 * cell_rhs(i) += (1.0 *
801 * fe_values.shape_value(i, q_index) *
802 * fe_values.JxW(q_index));
808 * Finally, transfer the contributions from @p
cell_matrix and
809 * @p cell_rhs into the global objects.
812 * cell->get_dof_indices(local_dof_indices);
813 * constraints.distribute_local_to_global(
814 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
818 * Now we are done assembling the linear system. The constraint
matrix took
819 * care of applying the boundary conditions and also eliminated hanging node
820 * constraints. The constrained nodes are still in the linear system (there
821 * is a nonzero entry, chosen in a way that the matrix is well conditioned,
822 * on the diagonal of the matrix and all other entries
for this line are set
823 * to zero) but the computed
values are
invalid (i.e., the corresponding
824 * entries in <code>system_rhs</code> are currently meaningless). We compute
825 * the correct
values for these nodes at the
end of the <code>solve</code>
835 * <a name=
"Step6solve"></a>
836 * <h4>Step6::solve</h4>
840 * We
continue with gradual improvements. The function that solves the linear
841 * system again uses the SSOR preconditioner, and is again unchanged except
842 * that we have to incorporate hanging node constraints. As mentioned above,
844 * hanging node constraints and boundary
values have been removed from the
845 * linear system by giving the rows and columns of the
matrix a special
846 * treatment. This way, the
values for these degrees of freedom have wrong,
847 * but well-defined
values after solving the linear system. What we then have
848 * to
do is to use the constraints to assign to them the
values that they
849 * should have. This process, called <code>distributing</code> constraints,
850 * computes the
values of constrained nodes from the
values of the
852 * that you find at the
end of
this function:
859 *
void Step6<dim>::solve()
865 * preconditioner.
initialize(system_matrix, 1.2);
867 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
869 * constraints.distribute(solution);
876 * <a name=
"Step6refine_grid"></a>
877 * <h4>Step6::refine_grid</h4>
881 * We use a sophisticated error estimation scheme to
refine the mesh instead
883 * implements an error estimator
for the Laplace equation; it can in principle
884 * handle variable coefficients, but we will not use these advanced features,
885 * but rather use its most simple form since we are not interested in
886 * quantitative results but only in a quick way to generate locally refined
891 * Although the error estimator derived by Kelly et al. was originally
892 * developed
for the Laplace equation, we have found that it is also well
893 * suited to quickly generate locally refined grids
for a wide
class of
894 * problems. This error estimator uses the solution
gradient's jump at
895 * cell faces (which is a measure for the second derivatives) and
896 * scales it by the size of the cell. It is therefore a measure for the local
897 * smoothness of the solution at the place of each cell and it is thus
898 * understandable that it yields reasonable grids also for hyperbolic
899 * transport problems or the wave equation as well, although these grids are
900 * certainly suboptimal compared to approaches specially tailored to the
901 * problem. This error estimator may therefore be understood as a quick way to
902 * test an adaptive program.
906 * The way the estimator works is to take a <code>DoFHandler</code> object
907 * describing the degrees of freedom and a vector of values for each degree of
908 * freedom as input and compute a single indicator value for each active cell
909 * of the triangulation (i.e. one value for each of the active cells). To do
910 * so, it needs two additional pieces of information: a face quadrature formula,
911 * i.e., a quadrature formula on <code>dim-1</code> dimensional objects. We use
912 * a 3-point Gauss rule again, a choice that is consistent and appropriate with
913 * the bi-quadratic finite element shape functions in this program.
914 * (What constitutes a suitable quadrature rule here of course depends on
915 * knowledge of the way the error estimator evaluates the solution field. As
916 * said above, the jump of the gradient is integrated over each face, which
917 * would be a quadratic function on each face for the quadratic elements in
918 * use in this example. In fact, however, it is the square of the jump of the
919 * gradient, as explained in the documentation of that class, and that is a
920 * quartic function, for which a 3 point Gauss formula is sufficient since it
921 * integrates polynomials up to order 5 exactly.)
925 * Secondly, the function wants a list of boundary indicators for those
926 * boundaries where we have imposed Neumann values of the kind
927 * @f$\partial_n u(\mathbf x) = h(\mathbf x)@f$, along with a function @f$h(\mathbf
928 * x)@f$ for each such boundary. This information is represented by a map from
929 * boundary indicators to function objects describing the Neumann boundary
930 * values. In the present example program, we do not use Neumann boundary
931 * values, so this map is empty, and in fact constructed using the default
932 * constructor of the map in the place where the function call expects the
933 * respective function argument.
937 * The output is a vector of values for all active cells. While it may
938 * make sense to compute the <b>value</b> of a solution degree of freedom
939 * very accurately, it is usually not necessary to compute the <b>error
940 * indicator</b> corresponding to the solution on a cell particularly
941 * accurately. We therefore typically use a vector of floats instead of a vector
942 * of doubles to represent error indicators.
946 * void Step6<dim>::refine_grid()
948 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
950 * KellyErrorEstimator<dim>::estimate(dof_handler,
951 * QGauss<dim - 1>(fe.degree + 1),
954 * estimated_error_per_cell);
958 * The above function returned one error indicator value for each cell in
959 * the <code>estimated_error_per_cell</code> array. Refinement is now done
960 * as follows: refine those 30 per cent of the cells with the highest error
961 * values, and coarsen the 3 per cent of cells with the lowest values.
965 * One can easily verify that if the second number were zero, this would
966 * approximately result in a doubling of cells in each step in two space
967 * dimensions, since for each of the 30 per cent of cells, four new would be
968 * replaced, while the remaining 70 per cent of cells remain untouched. In
969 * practice, some more cells are usually produced since it is disallowed
970 * that a cell is refined twice while the neighbor cell is not refined; in
971 * that case, the neighbor cell would be refined as well.
975 * In many applications, the number of cells to be coarsened would be set to
976 * something larger than only three per cent. A non-zero value is useful
977 * especially if for some reason the initial (coarse) grid is already rather
978 * refined. In that case, it might be necessary to refine it in some
979 * regions, while coarsening in some other regions is useful. In our case
980 * here, the initial grid is very coarse, so coarsening is only necessary in
981 * a few regions where over-refinement may have taken place. Thus a small,
982 * non-zero value is appropriate here.
986 * The following function now takes these refinement indicators and flags
987 * some cells of the triangulation for refinement or coarsening using the
988 * method described above. It is from a class that implements several
989 * different algorithms to refine a triangulation based on cell-wise error
993 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
994 * estimated_error_per_cell,
1000 * After the previous function has exited, some cells are flagged for
1001 * refinement, and some other for coarsening. The refinement or coarsening
1002 * itself is not performed by now, however, since there are cases where
1003 * further modifications of these flags is useful. Here, we don't want to
do
1004 * any such thing, so we can tell the
triangulation to perform the actions
1005 *
for which the cells are flagged:
1015 * <a name=
"Step6output_results"></a>
1016 * <h4>Step6::output_results</h4>
1020 * At the
end of computations on each grid, and just before we
continue the
1021 * next cycle with mesh refinement, we want to output the results from
this
1026 * We have already seen in @ref step_1
"step-1" how
this can be achieved
for the
1027 * mesh itself. Here, we change a few things:
1029 * <li>We use two different formats:
gnuplot and VTU.</li>
1030 * <li>We embed the cycle number in the output file name.</li>
1032 * provide a few extra visualization arguments so that edges appear
1033 * curved. This is explained in further detail in @ref step_10
"step-10".</li>
1037 *
template <
int dim>
1038 *
void Step6<dim>::output_results(
const unsigned int cycle)
const
1042 * std::ofstream output(
"grid-" + std::to_string(cycle) +
".gnuplot");
1044 * grid_out.set_flags(gnuplot_flags);
1052 * data_out.add_data_vector(solution,
"solution");
1053 * data_out.build_patches();
1055 * std::ofstream output(
"solution-" + std::to_string(cycle) +
".vtu");
1056 * data_out.write_vtu(output);
1064 * <a name=
"Step6run"></a>
1065 * <h4>Step6::run</h4>
1069 * The
final function before <code>main()</code> is again the main driver of
1070 * the
class, <code>
run()</code>. It is similar to the one of @ref step_5
"step-5", except
1071 * that we generate a file in the program again instead of reading it from
1072 * disk, in that we adaptively instead of globally
refine the mesh, and that
1073 * we output the solution on the
final mesh in the present function.
1077 * The
first block in the main
loop of the function deals with mesh generation.
1078 * If
this is the
first cycle of the program, instead of reading the grid from
1079 * a file on disk as in the previous example, we now again create it
using a
1080 * library function. The domain is again a circle with
center at the origin and
1081 * a radius of one (these are the two hidden arguments to the function, which
1082 * have
default values).
1086 * You will notice by looking at the coarse grid that it is of inferior
1087 * quality than the one which we read from the file in the previous example:
1088 * the cells are less equally formed. However,
using the library function
this
1089 * program works in any space dimension, which was not the
case before.
1093 * In
case we find that
this is not the
first cycle, we want to
refine the
1094 * grid. Unlike the global refinement employed in the last example program, we
1095 * now use the adaptive procedure described above.
1099 * The rest of the
loop looks as before:
1102 *
template <
int dim>
1103 *
void Step6<dim>::run()
1105 *
for (
unsigned int cycle = 0; cycle < 8; ++cycle)
1107 * std::cout <<
"Cycle " << cycle <<
':' << std::endl;
1118 * std::cout <<
" Number of active cells: "
1123 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1126 * assemble_system();
1128 * output_results(cycle);
1136 * <a name=
"Thecodemaincodefunction"></a>
1137 * <h3>The <code>main</code> function</h3>
1141 * The main function is unaltered in its functionality from the previous
1142 * example, but we have taken a step of additional caution. Sometimes,
1143 * something goes wrong (such as insufficient disk space upon writing an
1144 * output file, not enough memory when trying to allocate a vector or a
1145 * matrix, or
if we can
't read from or write to a file for whatever reason),
1146 * and in these cases the library will throw exceptions. Since these are
1147 * run-time problems, not programming errors that can be fixed once and for
1148 * all, this kind of exceptions is not switched off in optimized mode, in
1149 * contrast to the <code>Assert</code> macro which we have used to test
1150 * against programming errors. If uncaught, these exceptions propagate the
1151 * call tree up to the <code>main</code> function, and if they are not caught
1152 * there either, the program is aborted. In many cases, like if there is not
1153 * enough memory or disk space, we can't
do anything but we can at least print
1154 * some text trying to explain the reason why the program failed. A way to
do
1155 * so is shown in the following. It is certainly useful to write any larger
1156 * program in
this way, and you can
do so by more or less copying
this
1157 * function except
for the <code>try</code> block that actually encodes the
1158 * functionality particular to the present application.
1165 * The
general idea behind the layout of
this function is as follows: let
's
1166 * try to run the program as we did before...
1171 * Step6<2> laplace_problem_2d;
1172 * laplace_problem_2d.run();
1176 * ...and if this should fail, try to gather as much information as
1177 * possible. Specifically, if the exception that was thrown is an object of
1178 * a class that is derived from the C++ standard class
1179 * <code>exception</code>, then we can use the <code>what</code> member
1180 * function to get a string which describes the reason why the exception was
1185 * The deal.II exception classes are all derived from the standard class,
1186 * and in particular, the <code>exc.what()</code> function will return
1187 * approximately the same string as would be generated if the exception was
1188 * thrown using the <code>Assert</code> macro. You have seen the output of
1189 * such an exception in the previous example, and you then know that it
1190 * contains the file and line number of where the exception occurred, and
1191 * some other information. This is also what the following statements would
1196 * Apart from this, there isn't much that we can
do except exiting the
1197 * program with an error code (
this is what the <code>
return 1;</code>
1201 * catch (
std::exception &exc)
1205 * <<
"----------------------------------------------------"
1207 * std::cerr <<
"Exception on processing: " << std::endl
1208 * << exc.what() << std::endl
1209 * <<
"Aborting!" << std::endl
1210 * <<
"----------------------------------------------------"
1217 * If the exception that was thrown somewhere was not an
object of a
class
1218 * derived from the standard <code>exception</code>
class, then we can
't do
1219 * anything at all. We then simply print an error message and exit.
1224 * std::cerr << std::endl
1226 * << "----------------------------------------------------"
1228 * std::cerr << "Unknown exception!" << std::endl
1229 * << "Aborting!" << std::endl
1230 * << "----------------------------------------------------"
1237 * If we got to this point, there was no exception which propagated up to
1238 * the main function (there may have been exceptions, but they were caught
1239 * somewhere in the program or the library). Therefore, the program
1240 * performed as was expected and we can return without error.
1246<a name="Results"></a><h1>Results</h1>
1250The output of the program looks as follows:
1253 Number of active cells: 20
1254 Number of degrees of freedom: 89
1256 Number of active cells: 44
1257 Number of degrees of freedom: 209
1259 Number of active cells: 92
1260 Number of degrees of freedom: 449
1262 Number of active cells: 200
1263 Number of degrees of freedom: 921
1265 Number of active cells: 440
1266 Number of degrees of freedom: 2017
1268 Number of active cells: 956
1269 Number of degrees of freedom: 4425
1271 Number of active cells: 1916
1272 Number of degrees of freedom: 8993
1274 Number of active cells: 3860
1275 Number of degrees of freedom: 18353
1280As intended, the number of cells roughly doubles in each cycle. The
1281number of degrees is slightly more than four times the number of
1282cells; one would expect a factor of exactly four in two spatial
1283dimensions on an infinite grid (since the spacing between the degrees
1284of freedom is half the cell width: one additional degree of freedom on
1285each edge and one in the middle of each cell), but it is larger than
1286that factor due to the finite size of the mesh and due to additional
1287degrees of freedom which are introduced by hanging nodes and local
1292The program outputs the solution and mesh in each cycle of the
1293refinement loop. The solution looks as follows:
1295<img src="https://www.dealii.org/images/steps/developer/step-6.solution.9.2.png" alt="">
1297It is interesting to follow how the program arrives at the final mesh:
1299<div class="twocolumn" style="width: 80%">
1301 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_0.svg"
1302 alt="Initial grid: the five-cell circle grid with one global refinement."
1303 width="300" height="300">
1306 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_1.svg"
1307 alt="First grid: the five-cell circle grid with two global refinements."
1308 width="300" height="300">
1311 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_2.svg"
1312 alt="Second grid: the five-cell circle grid with one adaptive refinement."
1313 width="300" height="300">
1316 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_3.svg"
1317 alt="Third grid: the five-cell circle grid with two adaptive
1318 refinements, showing clustering around the inner circle."
1319 width="300" height="300">
1322 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_4.svg"
1323 alt="Fourth grid: the five-cell circle grid with three adaptive
1324 refinements, showing clustering around the inner circle."
1325 width="300" height="300">
1328 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_5.svg"
1329 alt="Fifth grid: the five-cell circle grid with four adaptive
1330 refinements, showing clustering around the inner circle."
1331 width="300" height="300">
1334 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_6.svg"
1335 alt="Sixth grid: the five-cell circle grid with five adaptive
1336 refinements, showing clustering around the inner circle."
1337 width="300" height="300">
1340 <img src="https://www.dealii.org/images/steps/developer/step_6_grid_7.svg"
1341 alt="Last grid: the five-cell circle grid with six adaptive
1342 refinements, showing that most cells are clustered around the inner circle."
1343 width="300" height="300">
1348It is clearly visible that the region where the solution has a kink,
1349i.e. the circle at radial distance 0.5 from the center, is
1350refined most. Furthermore, the central region where the solution is
1351very smooth and almost flat, is almost not refined at all, but this
1352results from the fact that we did not take into account that the
1353coefficient is large there. The region outside is refined rather
1354arbitrarily, since the second derivative is constant there and refinement
1355is therefore mostly based on the size of the cells and their deviation
1356from the optimal square.
1360<a name="extensions"></a>
1361<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1364<a name="Solversandpreconditioners"></a><h4>Solvers and preconditioners</h4>
1367One thing that is always worth playing around with if one solves
1368problems of appreciable size (much bigger than the one we have here)
1369is to try different solvers or preconditioners. In the current case,
1370the linear system is symmetric and positive definite, which makes the
1371Conjugate Gradient (CG) algorithm pretty much the canonical choice for solving. However,
1372the SSOR preconditioner we use in the <code>solve()</code> function is
1375In deal.II, it is relatively simple to change the preconditioner.
1376Several simple preconditioner choices are accessible
1377by changing the following two lines
1379 PreconditionSSOR<SparseMatrix<double>> preconditioner;
1380 preconditioner.initialize(system_matrix, 1.2);
1382of code in the program. For example, switching this into
1384 PreconditionSSOR<SparseMatrix<double>> preconditioner;
1385 preconditioner.initialize(system_matrix, 1.0);
1387allows us to try out a different relaxation parameter for SSOR (1.0 instead
1388of the 1.2 in the original program). Similarly, by using
1390 PreconditionJacobi<SparseMatrix<double>> preconditioner;
1391 preconditioner.initialize(system_matrix);
1393we can use Jacobi as a preconditioner. And by using
1395 SparseILU<double> preconditioner;
1396 preconditioner.initialize(system_matrix);
1398we can use a simple incomplete LU decomposition without any thresholding or
1399strengthening of the diagonal (to use this preconditioner, you have to also
1400add the header file <code>deal.II/lac/sparse_ilu.h</code> to the include list
1401at the top of the file).
1403Using these various different preconditioners, we can compare the
1404number of CG iterations needed (available through the
1405<code>solver_control.last_step()</code> call, see
1406@ref step_4 "step-4") as well as CPU time needed (using the Timer class,
1407discussed, for example, in @ref step_28 "step-28") and get the
1408following results (left: iterations; right: CPU time):
1410<table width="60%" align="center">
1413 <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_iterations.png" alt="">
1416 <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_time.png" alt="">
1421As we can see, all preconditioners behave pretty much the same on this
1422simple problem, with the number of iterations growing like @f${\cal
1423O}(N^{1/2})@f$ and because each iteration requires around @f${\cal
1424O}(N)@f$ operations the total CPU time grows like @f${\cal
1425O}(N^{3/2})@f$ (for the few smallest meshes, the CPU time is so small
1426that it doesn't record). Note that even though it is the simplest
1427method, Jacobi is the fastest
for this problem.
1429The situation changes slightly when the finite element is not a
1430bi-quadratic one (i.e., polynomial degree two) as selected in the
1431constructor of
this program, but a bi-linear one (polynomial degree one).
1432If one makes
this change, the results are as follows:
1434<table width=
"60%" align=
"center">
1437 <img src=
"https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_iterations.png" alt=
"">
1440 <img src=
"https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_time.png" alt=
"">
1445In other words,
while the increase in iterations and CPU time is as
1446before, Jacobi is now the method that
requires the most iterations; it
1447is still the fastest one, however, owing to the simplicity of the
1448operations it has to perform. This is not to say that Jacobi
1449is actually a good preconditioner --
for problems of appreciable size, it is
1450definitely not, and other methods will be substantially better -- but really
1451only that it is fast because its implementation is so simple that it can
1452compensate
for a larger number of iterations.
1454The message to take away from
this is not that simplicity in
1455preconditioners is
always best. While
this may be
true for the current
1456problem, it definitely is not once we move to more complicated
1457problems (elasticity or Stokes,
for examples @ref step_8
"step-8" or
1458@ref step_22
"step-22"). Secondly, all of these preconditioners still
1459lead to an increase in the number of iterations as the number @f$N@f$ of
1460degrees of freedom grows,
for example @f${\cal O}(N^\alpha)@f$;
this, in
1461turn, leads to a total growth in effort as @f${\cal O}(N^{1+\alpha})@f$
1462since each iteration takes @f${\cal O}(N)@f$ work. This behavior is
1463undesirable: we would really like to solve linear systems with @f$N@f$
1464unknowns in a total of @f${\cal O}(N)@f$ work; there is a
class
1465of preconditioners that can achieve
this, namely geometric (@ref step_16
"step-16",
1466@ref step_37
"step-37", @ref step_39
"step-39")
1467or algebraic multigrid (@ref step_31
"step-31", @ref step_40
"step-40", and several others)
1468preconditioners. They are, however, significantly more complex than
1469the preconditioners outlined above, and so we will leave their use
1470to these later tutorial programs. The
point to make, however, is
1471that
"real" finite element programs
do not use the preconditioners
1472we mention above: These are simply shown
for expository purposes.
1474Finally, the last message to take
1475home is that when the data shown above was generated (in 2018), linear
1476systems with 100,000 unknowns are
1477easily solved on a desktop or laptop machine in about a
second, making
1478the solution of relatively simple 2
d problems even to very high
1479accuracy not that big a task as it used to be in the
1480past. At the same time, the situation
for 3
d problems continues to be
1481quite different: A uniform 2
d mesh with 100,000 unknowns corresponds to
1482a grid with about @f$300 \times 300@f$ nodes; the corresponding 3
d mesh has
1483@f$300 \times 300 \times 300@f$ nodes and 30 million unknowns. Because
1484finite element matrices in 3
d have many more
nonzero entries than in 2
d,
1485solving these linear systems will not only take 300 times as much CPU
1486time, but substantially longer. In other words, achieving the
1487same resolution in 3
d *is* quite a large problem, and solving it
1488within a reasonable amount of time *will* require much more work to
1489implement better linear solvers. As mentioned above, multigrid methods
1490and
matrix-
free methods (see,
for example, @ref step_37
"step-37"), along with
1491parallelization (@ref step_40
"step-40") will be necessary, but are then also able
1492to comfortably solve such linear systems.
1496<a name=
"Abettermesh"></a><h4>A better mesh</h4>
1499If you look at the meshes above, you will see even though the domain is the
1500unit disk, and the jump in the coefficient lies along a circle, the cells
1501that make up the mesh
do not track
this geometry well. The reason, already hinted
1502at in @ref step_1
"step-1", is that in the absence of other information,
1504coarse grid cells but has, of course, no real idea what kind of geometry they
1505might represent when looked at together. For this reason, we need to tell
1506the
Triangulation what to do when a cell is refined: where should the new
1507vertices at the edge midpoints and the cell midpoint be located so that the
1508child cells better represent the desired geometry than the parent cell.
1510To visualize what the
triangulation actually knows about the geometry,
1511it is not enough to just output the location of
vertices and draw a
1512straight line for each edge; instead, we have to output both interior
1513and boundary lines as multiple segments so that they look
1514curved. We can do this by making one change to the
gnuplot part of
1515<code>output_results</code>:
1519 std::ofstream output(
"grid-" + std::to_string(cycle) +
".gnuplot");
1527In the code above, we already
do this for faces that sit at the boundary:
this
1530<i>interior</i> also track a circular domain, we need to work a bit harder,
1531though. First, recall that our coarse mesh consists of a central square
1532cell and four cells around it. Now
first consider what would happen
if we
1534but also the four cells at the perimeter as well as all of their faces. We can
1535do this by adding the following snippet (testing that the
center of a cell is
1536larger than a small multiple, say one tenth, of the cell diameter away from
1537center of the mesh only fails
for the central square of the mesh):
1543for (
const auto &cell :
triangulation.active_cell_iterators())
1545 cell->set_all_manifold_ids(0);
1550After a few global refinement steps,
this would lead to a mesh of the following
1554 <div
class=
"onecolumn" style=
"width: 80%">
1556 <img src=
"https://www.dealii.org/images/steps/developer/step_6_bad_grid_4.svg"
1557 alt=
"Grid where some central cells are nearly triangular."
1558 width=
"300" height=
"300">
1562This is not a good mesh: the central cell has been refined in such a way that
1563the children located in the four corners of the original central cell
1564<i>degenerate</i>: they all tend towards triangles as mesh refinement
1565continues. This means that the Jacobian
matrix of the transformation from
1566reference cell to actual cell degenerates
for these cells, and because
1567all error estimates
for finite element solutions contain the
norm of the
1568inverse of the Jacobian
matrix, you will get very large errors on these
1569cells and, in the limit as mesh refinement, a loss of convergence order because
1570the cells in these corners become worse and worse under mesh refinement.
1572So we need something smarter. To
this end, consider the following solution
1573originally developed by Konstantin Ladutenko. We will use the following code:
1578const double core_radius = 1.0/5.0,
1579 inner_radius = 1.0/3.0;
1592for (
const auto &cell :
triangulation.active_cell_iterators())
1593 if (mesh_center.distance(cell->
center()) < 1
e-5)
1596 cell->vertex(v) *= core_radius/mesh_center.distance(cell->vertex(v));
1600for (
const auto &cell :
triangulation.active_cell_iterators())
1601 if (mesh_center.distance(cell->
center()) >= 1
e-5)
1602 cell->set_refine_flag();
1613for (
const auto &cell :
triangulation.active_cell_iterators())
1616 const double dist = mesh_center.
distance(cell->vertex(v));
1617 if (dist > core_radius*1.0001 && dist < 0.9999)
1618 cell->vertex(v) *= inner_radius/dist;
1632for (
const auto &cell :
triangulation.active_cell_iterators())
1634 bool is_in_inner_circle =
false;
1636 if (mesh_center.distance(cell->vertex(v)) < inner_radius)
1638 is_in_inner_circle =
true;
1642 if (is_in_inner_circle ==
false)
1647 cell->set_all_manifold_ids(0);
1651This code then generates the following, much better sequence of meshes:
1653<div
class=
"twocolumn" style=
"width: 80%">
1655 <img src=
"https://www.dealii.org/images/steps/developer/step_6_grid_0_ladutenko.svg"
1656 alt=
"Initial grid: the Ladutenko grid with one global refinement."
1657 width=
"300" height=
"300">
1660 <img src=
"https://www.dealii.org/images/steps/developer/step_6_grid_1_ladutenko.svg"
1661 alt=
"First adaptively refined Ladutenko grid."
1662 width=
"300" height=
"300">
1665 <img src=
"https://www.dealii.org/images/steps/developer/step_6_grid_2_ladutenko.svg"
1666 alt=
"Second adaptively refined Ladutenko grid."
1667 width=
"300" height=
"300">
1670 <img src=
"https://www.dealii.org/images/steps/developer/step_6_grid_3_ladutenko.svg"
1671 alt=
"Third adaptively refined Ladutenko grid."
1672 width=
"300" height=
"300">
1675 <img src=
"https://www.dealii.org/images/steps/developer/step_6_grid_4_ladutenko.svg"
1676 alt=
"Fourth adaptively refined Ladutenko grid. The cells are clustered
1677 along the inner circle."
1678 width=
"300" height=
"300">
1681 <img src=
"https://www.dealii.org/images/steps/developer/step_6_grid_5_ladutenko.svg"
1682 alt=
"Fifth adaptively refined Ladutenko grid: the cells are clustered
1683 along the inner circle."
1684 width=
"300" height=
"300">
1688Creating good meshes, and in particular making them fit the geometry you
1689want, is a complex topic in itself. You can find much more on
this in
1690@ref step_49
"step-49", @ref step_53
"step-53", and @ref step_54
"step-54", among other tutorial programs that cover
1691the issue. @ref step_65
"step-65" shows another, less manual way to achieve a mesh
1692well fit to the problem here.
1693Information on curved domains can also be found in the
1694documentation module on @ref manifold
"Manifold descriptions".
1696Why does it make sense to choose a mesh that tracks the
internal
1697interface? There are a number of reasons, but the most essential one
1698comes down to what we actually integrate in our bilinear
1699form. Conceptually, we want to integrate the term @f$A_{ij}^
K=\int_K
1700a(\mathbf x) \nabla \varphi_i(\mathbf x) \nabla \varphi_j(\mathbf x) ;
dx@f$ as the
1701contribution of cell @f$K@f$ to the
matrix entry @f$A_{ij}@f$. We can not
1702compute it exactly and have to resort to quadrature. We know that
1703quadrature is accurate
if the integrand is smooth. That is because
1704quadrature in essence computes a polynomial approximation to the
1705integrand that coincides with the integrand in the quadrature points,
1706and then computes the
volume under
this polynomial as an approximation
1707to the
volume under the original integrand. This polynomial
1708interpolant is accurate
if the integrand is smooth on a cell, but it
1709is usually rather inaccurate
if the integrand is discontinuous on a
1712Consequently, it is worthwhile to align cells in such a way that the
1713interfaces across which the coefficient is discontinuous are aligned
1714with cell interfaces. This way, the coefficient is
constant on each
1715cell, following which the integrand will be smooth, and its polynomial
1716approximation and the quadrature approximation of the integral will
1717both be accurate. Note that such an alignment is common in many
1718practical cases, so deal.II provides a number of
functions (such as
1719@ref GlossMaterialId
"material_id") to help manage such a scenario.
1720Refer to @ref step_28
"step-28" and @ref step_46
"step-46" for examples of how material ids can be
1723Finally, let us consider the
case of a coefficient that has a smooth
1724and non-uniform distribution in space. We can repeat once again all of
1725the above discussion on the representation of such a function with the
1726quadrature. So, to simulate it accurately there are a few readily
1727available options: you could
reduce the cell size, increase the order
1728of the polynomial used in the quadrature formula, select a more
1729appropriate quadrature formula, or perform a combination of these
1730steps. The key is that providing the best fit of the coefficient
's
1731spatial dependence with the quadrature polynomial will lead to a more
1732accurate finite element solution of the PDE.
1734As a final note: The discussion in the previous paragraphs shows, we here
1735have a very concrete way of stating what we think of a good mesh -- it should
1736be aligned with the jump in the coefficient. But one could also have asked
1737this kind of question in a more general setting: Given some equation with
1738a smooth solution and smooth coefficients, can we say what a good mesh
1739would look like? This is a question for which the answer is easier to state
1740in intuitive terms than mathematically: A good mesh has cells that all,
1741by and large, look like squares (or cubes, in 3d). A bad mesh would contain
1742cells that are very elongated in some directions or, more generally, for which
1743there are cells that have both short and long edges. There are many ways
1744in which one could assign a numerical quality index to each cell that measures
1745whether the cell is "good" or "bad"; some of these are often chosen because
1746they are cheap and easy to compute, whereas others are based on what enters
1747into proofs of convergence. An example of the former would be the ratio of
1748the longest to the shortest edge of a cell: In the ideal case, that ratio
1749would be one; bad cells have values much larger than one. An example of the
1750latter kind would consider the gradient (the "Jacobian") of the mapping
1751from the reference cell @f$\hat K=[0,1]^d@f$ to the real cell @f$K@f$; this
1752gradient is a matrix, and a quantity that enters into error estimates
1753is the maximum over all points on the reference cell of the ratio of the
1754largest to the smallest eigenvalue of this matrix. It is again not difficult
1755to see that this ratio is constant if the cell @f$K@f$ is an affine image of
1756@f$\hat K@f$, and that it is one for squares and cubes.
1758In practice, it might be interesting to visualize such quality measures.
1759The function GridTools::compute_aspect_ratio_of_cells() provides one
1760way to get this kind of information. Even better, visualization tools
1761such as VisIt often allow you to visualize this sort of information
1762for a variety of measures from within the visualization software
1763itself; in the case of VisIt, just add a "pseudo-color" plot and select
1764one of the mesh quality measures instead of the solution field.
1767<a name="Playingwiththeregularityofthesolution"></a><h4>Playing with the regularity of the solution</h4>
1770From a mathematical perspective, solutions of the Laplace equation
1774on smoothly bounded, convex domains are known to be smooth themselves. The exact degree
1775of smoothness, i.e., the function space in which the solution lives, depends
1776on how smooth exactly the boundary of the domain is, and how smooth the right
1777hand side is. Some regularity of the solution may be lost at the boundary, but
1778one generally has that the solution is twice more differentiable in
1779compact subsets of the domain than the right hand side.
1780If, in particular, the right hand side satisfies @f$f\in C^\infty(\Omega)@f$, then
1781@f$u \in C^\infty(\Omega_i)@f$ where @f$\Omega_i@f$ is any compact subset of @f$\Omega@f$
1782(@f$\Omega@f$ is an open domain, so a compact subset needs to keep a positive distance
1783from @f$\partial\Omega@f$).
1785The situation we chose for the current example is different, however: we look
1786at an equation with a non-constant coefficient @f$a(\mathbf x)@f$:
1788 -\nabla \cdot (a \nabla u) = f.
1790Here, if @f$a@f$ is not smooth, then the solution will not be smooth either,
1791regardless of @f$f@f$. In particular, we expect that wherever @f$a@f$ is discontinuous
1792along a line (or along a plane in 3d),
1793the solution will have a kink. This is easy to see: if for example @f$f@f$
1794is continuous, then @f$f=-\nabla \cdot (a \nabla u)@f$ needs to be
1795continuous. This means that @f$a \nabla u@f$ must be continuously differentiable
1796(not have a kink). Consequently, if @f$a@f$ has a discontinuity, then @f$\nabla u@f$
1797must have an opposite discontinuity so that the two exactly cancel and their
1798product yields a function without a discontinuity. But for @f$\nabla u@f$ to have
1799a discontinuity, @f$u@f$ must have a kink. This is of course exactly what is
1800happening in the current example, and easy to observe in the pictures of the
1803In general, if the coefficient @f$a(\mathbf x)@f$ is discontinuous along a line in 2d,
1804or a plane in 3d, then the solution may have a kink, but the gradient of the
1805solution will not go to infinity. That means, that the solution is at least
1806still in the <a href="https://en.wikipedia.org/wiki/Sobolev_space">Sobolev space</a>
1807@f$W^{1,\infty}@f$ (i.e., roughly speaking, in the
1808space of functions whose derivatives are bounded). On the other hand,
1809we know that in the most
1810extreme cases -- i.e., where the domain has reentrant corners, the
1811right hand side only satisfies @f$f\in H^{-1}@f$, or the coefficient @f$a@f$ is only in
1812@f$L^\infty@f$ -- all we can expect is that @f$u\in H^1@f$ (i.e., the
1814href="https://en.wikipedia.org/wiki/Sobolev_space#Sobolev_spaces_with_integer_k">Sobolev
1815space</a> of functions whose derivative is square integrable), a much larger space than
1816@f$W^{1,\infty}@f$. It is not very difficult to create cases where
1817the solution is in a space @f$H^{1+s}@f$ where we can get @f$s@f$ to become as small
1818as we want. Such cases are often used to test adaptive finite element
1819methods because the mesh will have to resolve the singularity that causes
1820the solution to not be in @f$W^{1,\infty}@f$ any more.
1822The typical example one uses for this is called the <i>Kellogg problem</i>
1823(referring to @cite Kel74), which in the commonly used form has a coefficient
1824@f$a(\mathbf x)@f$ that has different values in the four quadrants of the plane
1825(or eight different values in the octants of @f${\mathbb R}^3@f$). The exact degree
1826of regularity (the @f$s@f$ in the index of the Sobolev space above) depends on the
1827values of @f$a(\mathbf x)@f$ coming together at the origin, and by choosing the
1828jumps large enough, the regularity of the solution can be made as close as
1829desired to @f$H^1@f$.
1831To implement something like this, one could replace the coefficient
1832function by the following (shown here only for the 2d case):
1835double coefficient (const Point<dim> &p)
1837 if ((p[0] < 0) && (p[1] < 0)) // lower left quadrant
1839 else if ((p[0] >= 0) && (p[1] < 0)) // lower right quadrant
1841 else if ((p[0] < 0) && (p[1] >= 0)) // upper left quadrant
1843 else if ((p[0] >= 0) && (p[1] >= 0)) // upper right quadrant
1847 Assert(false, ExcInternalError());
1852(Adding the <code>Assert</code> at the end ensures that either an exception
1853is thrown or that the program aborts if we ever get to that point
1854-- which of course we shouldn't,
1855but
this is a good way to insure yourself: we all make mistakes by
1856sometimes not thinking of all cases,
for example by checking
1857for <code>p[0]</code> to be less than and greater than zero,
1858rather than greater-or-
equal to zero, and thereby forgetting
1859some cases that would otherwise lead to bugs that are awkward
1860to find. The <code>
return 0;</code> at the
end is only there to
1861avoid compiler warnings that the function does not
end in a
1862<code>
return</code> statement -- the compiler cannot see that the
1863function would never actually get to that
point because of the
1864preceding <code>
Assert</code> statement.)
1866By playing with such cases where four or more sectors come
1867together and on which the coefficient has different values, one can
1868construct cases where the solution has singularities at the
1869origin. One can also see how the meshes are refined in such cases.
1872<a name="PlainProg"></a>
1873<h1> The plain program</h1>
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 set_flags(const GridOutFlags::DX &flags)
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
@ 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 hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=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)
@ 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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
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 run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation