Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-6.h
Go to the documentation of this file.
1  false);
714  *
715  * @endcode
716  *
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
720  * sparsity pattern:
721  *
722  * @code
723  * sparsity_pattern.copy_from(dsp);
724  *
725  * @endcode
726  *
727  * We may now, finally, initialize the sparse matrix:
728  *
729  * @code
730  * system_matrix.reinit(sparsity_pattern);
731  * }
732  *
733  *
734  * @endcode
735  *
736  *
737  * <a name="Step6assemble_system"></a>
738  * <h4>Step6::assemble_system</h4>
739  *
740 
741  *
742  * Next, we have to assemble the matrix. However, to copy the local matrix and
743  * vector on each cell into the global system, we are no longer using a
744  * hand-written loop. Instead, we use
745  * AffineConstraints::distribute_local_to_global() that internally executes
746  * this loop while performing Gaussian elimination on rows and columns
747  * corresponding to constrained degrees on freedom.
748  *
749 
750  *
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.
761  *
762  * @code
763  * template <int dim>
764  * void Step6<dim>::assemble_system()
765  * {
766  * const QGauss<dim> quadrature_formula(fe.degree + 1);
767  *
768  * FEValues<dim> fe_values(fe,
769  * quadrature_formula,
772  *
773  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
774  *
775  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
776  * Vector<double> cell_rhs(dofs_per_cell);
777  *
778  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
779  *
780  * for (const auto &cell : dof_handler.active_cell_iterators())
781  * {
782  * cell_matrix = 0;
783  * cell_rhs = 0;
784  *
785  * fe_values.reinit(cell);
786  *
787  * for (const unsigned int q_index : fe_values.quadrature_point_indices())
788  * {
789  * const double current_coefficient =
790  * coefficient<dim>(fe_values.quadrature_point(q_index));
791  * for (const unsigned int i : fe_values.dof_indices())
792  * {
793  * for (const unsigned int j : fe_values.dof_indices())
794  * cell_matrix(i, j) +=
795  * (current_coefficient * // a(x_q)
796  * fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
797  * fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
798  * fe_values.JxW(q_index)); // dx
799  *
800  * cell_rhs(i) += (1.0 * // f(x)
801  * fe_values.shape_value(i, q_index) * // phi_i(x_q)
802  * fe_values.JxW(q_index)); // dx
803  * }
804  * }
805  *
806  * @endcode
807  *
808  * Finally, transfer the contributions from @p cell_matrix and
809  * @p cell_rhs into the global objects.
810  *
811  * @code
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);
815  * }
816  * @endcode
817  *
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>
826  * function.
827  *
828  * @code
829  * }
830  *
831  *
832  * @endcode
833  *
834  *
835  * <a name="Step6solve"></a>
836  * <h4>Step6::solve</h4>
837  *
838 
839  *
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,
843  * the degrees of freedom from the AffineConstraints object corresponding to
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
851  * unconstrained ones, and requires only a single additional function call
852  * that you find at the end of this function:
853  *
854 
855  *
856  *
857  * @code
858  * template <int dim>
859  * void Step6<dim>::solve()
860  * {
861  * SolverControl solver_control(1000, 1e-12);
862  * SolverCG<Vector<double>> solver(solver_control);
863  *
864  * PreconditionSSOR<SparseMatrix<double>> preconditioner;
865  * preconditioner.initialize(system_matrix, 1.2);
866  *
867  * solver.solve(system_matrix, solution, system_rhs, preconditioner);
868  *
869  * constraints.distribute(solution);
870  * }
871  *
872  *
873  * @endcode
874  *
875  *
876  * <a name="Step6refine_grid"></a>
877  * <h4>Step6::refine_grid</h4>
878  *
879 
880  *
881  * We use a sophisticated error estimation scheme to refine the mesh instead
882  * of global refinement. We will use the KellyErrorEstimator class which
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
887  * grids.
888  *
889 
890  *
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.
903  *
904 
905  *
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.)
922  *
923 
924  *
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.
934  *
935 
936  *
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.
943  *
944  * @code
945  * template <int dim>
946  * void Step6<dim>::refine_grid()
947  * {
948  * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
949  *
950  * KellyErrorEstimator<dim>::estimate(dof_handler,
951  * QGauss<dim - 1>(fe.degree + 1),
952  * {},
953  * solution,
954  * estimated_error_per_cell);
955  *
956  * @endcode
957  *
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.
962  *
963 
964  *
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.
972  *
973 
974  *
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.
983  *
984 
985  *
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
990  * indicators.
991  *
992  * @code
993  * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
994  * estimated_error_per_cell,
995  * 0.3,
996  * 0.03);
997  *
998  * @endcode
999  *
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:
1006  *
1007  * @code
1008  * triangulation.execute_coarsening_and_refinement();
1009  * }
1010  *
1011  *
1012  * @endcode
1013  *
1014  *
1015  * <a name="Step6output_results"></a>
1016  * <h4>Step6::output_results</h4>
1017  *
1018 
1019  *
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
1022  * cycle.
1023  *
1024 
1025  *
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:
1028  * <ol>
1029  * <li>We use two different formats: gnuplot and VTU.</li>
1030  * <li>We embed the cycle number in the output file name.</li>
1031  * <li>For gnuplot output, we set up a GridOutFlags::Gnuplot object to
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>
1034  * </ol>
1035  *
1036  * @code
1037  * template <int dim>
1038  * void Step6<dim>::output_results(const unsigned int cycle) const
1039  * {
1040  * {
1041  * GridOut grid_out;
1042  * std::ofstream output("grid-" + std::to_string(cycle) + ".gnuplot");
1043  * GridOutFlags::Gnuplot gnuplot_flags(false, 5);
1044  * grid_out.set_flags(gnuplot_flags);
1045  * MappingQGeneric<dim> mapping(3);
1046  * grid_out.write_gnuplot(triangulation, output, &mapping);
1047  * }
1048  *
1049  * {
1050  * DataOut<dim> data_out;
1051  * data_out.attach_dof_handler(dof_handler);
1052  * data_out.add_data_vector(solution, "solution");
1053  * data_out.build_patches();
1054  *
1055  * std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
1056  * data_out.write_vtu(output);
1057  * }
1058  * }
1059  *
1060  *
1061  * @endcode
1062  *
1063  *
1064  * <a name="Step6run"></a>
1065  * <h4>Step6::run</h4>
1066  *
1067 
1068  *
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.
1074  *
1075 
1076  *
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).
1083  *
1084 
1085  *
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.
1090  *
1091 
1092  *
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.
1096  *
1097 
1098  *
1099  * The rest of the loop looks as before:
1100  *
1101  * @code
1102  * template <int dim>
1103  * void Step6<dim>::run()
1104  * {
1105  * for (unsigned int cycle = 0; cycle < 8; ++cycle)
1106  * {
1107  * std::cout << "Cycle " << cycle << ':' << std::endl;
1108  *
1109  * if (cycle == 0)
1110  * {
1112  * triangulation.refine_global(1);
1113  * }
1114  * else
1115  * refine_grid();
1116  *
1117  *
1118  * std::cout << " Number of active cells: "
1119  * << triangulation.n_active_cells() << std::endl;
1120  *
1121  * setup_system();
1122  *
1123  * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1124  * << std::endl;
1125  *
1126  * assemble_system();
1127  * solve();
1128  * output_results(cycle);
1129  * }
1130  * }
1131  *
1132  *
1133  * @endcode
1134  *
1135  *
1136  * <a name="Thecodemaincodefunction"></a>
1137  * <h3>The <code>main</code> function</h3>
1138  *
1139 
1140  *
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.
1159  *
1160  * @code
1161  * int main()
1162  * {
1163  * @endcode
1164  *
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...
1167  *
1168  * @code
1169  * try
1170  * {
1171  * Step6<2> laplace_problem_2d;
1172  * laplace_problem_2d.run();
1173  * }
1174  * @endcode
1175  *
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
1181  * thrown.
1182  *
1183 
1184  *
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
1192  * print.
1193  *
1194 
1195  *
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>
1198  * does):
1199  *
1200  * @code
1201  * catch (std::exception &exc)
1202  * {
1203  * std::cerr << std::endl
1204  * << std::endl
1205  * << "----------------------------------------------------"
1206  * << std::endl;
1207  * std::cerr << "Exception on processing: " << std::endl
1208  * << exc.what() << std::endl
1209  * << "Aborting!" << std::endl
1210  * << "----------------------------------------------------"
1211  * << std::endl;
1212  *
1213  * return 1;
1214  * }
1215  * @endcode
1216  *
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.
1220  *
1221  * @code
1222  * catch (...)
1223  * {
1224  * std::cerr << std::endl
1225  * << std::endl
1226  * << "----------------------------------------------------"
1227  * << std::endl;
1228  * std::cerr << "Unknown exception!" << std::endl
1229  * << "Aborting!" << std::endl
1230  * << "----------------------------------------------------"
1231  * << std::endl;
1232  * return 1;
1233  * }
1234  *
1235  * @endcode
1236  *
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.
1241  *
1242  * @code
1243  * return 0;
1244  * }
1245  * @endcode
1246 <a name="Results"></a><h1>Results</h1>
1247 
1248 
1249 
1250 The output of the program looks as follows:
1251 @code
1252 Cycle 0:
1253  Number of active cells: 20
1254  Number of degrees of freedom: 89
1255 Cycle 1:
1256  Number of active cells: 44
1257  Number of degrees of freedom: 209
1258 Cycle 2:
1259  Number of active cells: 92
1260  Number of degrees of freedom: 449
1261 Cycle 3:
1262  Number of active cells: 200
1263  Number of degrees of freedom: 921
1264 Cycle 4:
1265  Number of active cells: 440
1266  Number of degrees of freedom: 2017
1267 Cycle 5:
1268  Number of active cells: 956
1269  Number of degrees of freedom: 4425
1270 Cycle 6:
1271  Number of active cells: 1916
1272  Number of degrees of freedom: 8993
1273 Cycle 7:
1274  Number of active cells: 3860
1275  Number of degrees of freedom: 18353
1276 @endcode
1277 
1278 
1279 
1280 As intended, the number of cells roughly doubles in each cycle. The
1281 number of degrees is slightly more than four times the number of
1282 cells; one would expect a factor of exactly four in two spatial
1283 dimensions on an infinite grid (since the spacing between the degrees
1284 of freedom is half the cell width: one additional degree of freedom on
1285 each edge and one in the middle of each cell), but it is larger than
1286 that factor due to the finite size of the mesh and due to additional
1287 degrees of freedom which are introduced by hanging nodes and local
1288 refinement.
1289 
1290 
1291 
1292 The program outputs the solution and mesh in each cycle of the
1293 refinement loop. The solution looks as follows:
1294 
1295 <img src="https://www.dealii.org/images/steps/developer/step-6.solution.9.2.png" alt="">
1296 
1297 It is interesting to follow how the program arrives at the final mesh:
1298 
1299 <div class="twocolumn" style="width: 80%">
1300  <div>
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">
1304  </div>
1305  <div>
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">
1309  </div>
1310  <div>
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">
1314  </div>
1315  <div>
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">
1320  </div>
1321  <div>
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">
1326  </div>
1327  <div>
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">
1332  </div>
1333  <div>
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">
1338  </div>
1339  <div>
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">
1344  </div>
1345 </div>
1346 
1347 
1348 It is clearly visible that the region where the solution has a kink,
1349 i.e. the circle at radial distance 0.5 from the center, is
1350 refined most. Furthermore, the central region where the solution is
1351 very smooth and almost flat, is almost not refined at all, but this
1352 results from the fact that we did not take into account that the
1353 coefficient is large there. The region outside is refined rather
1354 arbitrarily, since the second derivative is constant there and refinement
1355 is therefore mostly based on the size of the cells and their deviation
1356 from the optimal square.
1357 
1358 
1359 
1360 <a name="extensions"></a>
1361 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1362 
1363 
1364 <a name="Solversandpreconditioners"></a><h4>Solvers and preconditioners</h4>
1365 
1366 
1367 
1368 One thing that is always worth playing around with if one solves
1369 problems of appreciable size (much bigger than the one we have here)
1370 is to try different solvers or preconditioners. In the current case,
1371 the linear system is symmetric and positive definite, which makes the
1372 CG algorithm pretty much the canonical choice for solving. However,
1373 the SSOR preconditioner we use in the <code>solve()</code> function is
1374 up for grabs.
1375 
1376 In deal.II, it is relatively simple to change the preconditioner. For
1377 example, by changing the existing lines of code
1378 @code
1379  PreconditionSSOR<SparseMatrix<double>> preconditioner;
1380  preconditioner.initialize(system_matrix, 1.2);
1381 @endcode
1382 into
1383 @code
1384  PreconditionSSOR<SparseMatrix<double>> preconditioner;
1385  preconditioner.initialize(system_matrix, 1.0);
1386 @endcode
1387 we can try out different relaxation parameters for SSOR. By using
1388 @code
1389  PreconditionJacobi<SparseMatrix<double>> preconditioner;
1390  preconditioner.initialize(system_matrix);
1391 @endcode
1392 we can use Jacobi as a preconditioner. And by using
1393 @code
1394  SparseILU<double> preconditioner;
1395  preconditioner.initialize(system_matrix);
1396 @endcode
1397 we can use a simple incomplete LU decomposition without any thresholding or
1398 strengthening of the diagonal (to use this preconditioner, you have to also
1399 add the header file <code>deal.II/lac/sparse_ilu.h</code> to the include list
1400 at the top of the file).
1401 
1402 Using these various different preconditioners, we can compare the
1403 number of CG iterations needed (available through the
1404 <code>solver_control.last_step()</code> call, see
1405 @ref step_4 "step-4") as well as CPU time needed (using the Timer class,
1406 discussed, for example, in @ref step_28 "step-28") and get the
1407 following results (left: iterations; right: CPU time):
1408 
1409 <table width="60%" align="center">
1410  <tr>
1411  <td align="center">
1412  <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_iterations.png" alt="">
1413  </td>
1414  <td align="center">
1415  <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_time.png" alt="">
1416  </td>
1417  </tr>
1418 </table>
1419 
1420 As we can see, all preconditioners behave pretty much the same on this
1421 simple problem, with the number of iterations growing like @f${\cal
1422 O}(N^{1/2})@f$ and because each iteration requires around @f${\cal
1423 O}(N)@f$ operations the total CPU time grows like @f${\cal
1424 O}(N^{3/2})@f$ (for the few smallest meshes, the CPU time is so small
1425 that it doesn't record). Note that even though it is the simplest
1426 method, Jacobi is the fastest for this problem.
1427 
1428 The situation changes slightly when the finite element is not a
1429 bi-quadratic one as set in the constructor of this program, but a
1430 bi-linear one. If one makes this change, the results are as follows:
1431 
1432 <table width="60%" align="center">
1433  <tr>
1434  <td align="center">
1435  <img src="https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_iterations.png" alt="">
1436  </td>
1437  <td align="center">
1438  <img src="https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_time.png" alt="">
1439  </td>
1440  </tr>
1441 </table>
1442 
1443 In other words, while the increase in iterations and CPU time is as
1444 before, Jacobi is now the method that requires the most iterations; it
1445 is still the fastest one, however, owing to the simplicity of the
1446 operations it has to perform. This is not to say that Jacobi
1447 is actually a good preconditioner -- for problems of appreciable size, it is
1448 definitely not, and other methods will be substantially better -- but really
1449 only that it is fast because its implementation is so simple that it can
1450 compensate for a larger number of iterations.
1451 
1452 The message to take away from this is not that simplicity in
1453 preconditioners is always best. While this may be true for the current
1454 problem, it definitely is not once we move to more complicated
1455 problems (elasticity or Stokes, for examples @ref step_8 "step-8" or
1456 @ref step_22 "step-22"). Secondly, all of these preconditioners still
1457 lead to an increase in the number of iterations as the number @f$N@f$ of
1458 degrees of freedom grows, for example @f${\cal O}(N^\alpha)@f$; this, in
1459 turn, leads to a total growth in effort as @f${\cal O}(N^{1+\alpha})@f$
1460 since each iteration takes @f${\cal O}(N)@f$ work. This behavior is
1461 undesirable: we would really like to solve linear systems with @f$N@f$
1462 unknowns in a total of @f${\cal O}(N)@f$ work; there is a class
1463 of preconditioners that can achieve this, namely geometric (@ref step_16 "step-16",
1464 @ref step_37 "step-37", @ref step_39 "step-39")
1465 or algebraic multigrid (@ref step_31 "step-31", @ref step_40 "step-40", and several others)
1466 preconditioners. They are, however, significantly more complex than
1467 the preconditioners outlined above.
1468 
1469 Finally, the last message to take
1470 home is that when the data shown above was generated (in 2018), linear
1471 systems with 100,000 unknowns are
1472 easily solved on a desktop machine in about a second, making
1473 the solution of relatively simple 2d problems even to very high
1474 accuracy not that big a task as it used to be even in the
1475 past. At the time, the situation for 3d problems was entirely different,
1476 but even that has changed substantially in the intervening time -- though
1477 solving problems in 3d to high accuracy remains a challenge.
1478 
1479 
1480 <a name="Abettermesh"></a><h4>A better mesh</h4>
1481 
1482 
1483 If you look at the meshes above, you will see even though the domain is the
1484 unit disk, and the jump in the coefficient lies along a circle, the cells
1485 that make up the mesh do not track this geometry well. The reason, already hinted
1486 at in @ref step_1 "step-1", is that in the absence of other information,
1487 the Triangulation class only sees a bunch of
1488 coarse grid cells but has, of course, no real idea what kind of geometry they
1489 might represent when looked at together. For this reason, we need to tell
1490 the Triangulation what to do when a cell is refined: where should the new
1491 vertices at the edge midpoints and the cell midpoint be located so that the
1492 child cells better represent the desired geometry than the parent cell.
1493 
1494 To visualize what the triangulation actually knows about the geometry,
1495 it is not enough to just output the location of vertices and draw a
1496 straight line for each edge; instead, we have to output both interior
1497 and boundary lines as multiple segments so that they look
1498 curved. We can do this by making one change to the gnuplot part of
1499 <code>output_results</code>:
1500 @code
1501 {
1502  GridOut grid_out;
1503  std::ofstream output("grid-" + std::to_string(cycle) + ".gnuplot");
1504  GridOutFlags::Gnuplot gnuplot_flags(false, 5, /*curved_interior_cells*/true);
1505  grid_out.set_flags(gnuplot_flags);
1506  MappingQGeneric<dim> mapping(3);
1507  grid_out.write_gnuplot(triangulation, output, &mapping);
1508 }
1509 @endcode
1510 
1511 In the code above, we already do this for faces that sit at the boundary: this
1512 happens automatically since we use GridGenerator::hyper_ball, which attaches a
1513 SphericalManifold to the boundary of the domain. To make the mesh
1514 <i>interior</i> also track a circular domain, we need to work a bit harder,
1515 though. First, recall that our coarse mesh consists of a central square
1516 cell and four cells around it. Now first consider what would happen if we
1517 also attached the SphericalManifold object not only to the four exterior faces
1518 but also the four cells at the perimeter as well as all of their faces. We can
1519 do this by adding the following snippet (testing that the center of a cell is
1520 larger than a small multiple, say one tenth, of the cell diameter away from
1521 center of the mesh only fails for the central square of the mesh):
1522 @code
1524 // after GridGenerator::hyper_ball is called the Triangulation has
1525 // a SphericalManifold with id 0. We can use it again on the interior.
1526 const Point<dim> mesh_center;
1527 for (const auto &cell : triangulation.active_cell_iterators())
1528  if (mesh_center.distance (cell->center()) > cell->diameter()/10)
1529  cell->set_all_manifold_ids(0);
1530 
1531 triangulation.refine_global(1);
1532 @endcode
1533 
1534 After a few global refinement steps, this would lead to a mesh of the following
1535 kind:
1536 
1537 
1538  <div class="onecolumn" style="width: 80%">
1539  <div>
1540  <img src="https://www.dealii.org/images/steps/developer/step_6_bad_grid_4.svg"
1541  alt="Grid where some central cells are nearly triangular."
1542  width="300" height="300">
1543  </div>
1544  </div>
1545 
1546 This is not a good mesh: the central cell has been refined in such a way that
1547 the children located in the four corners of the original central cell
1548 <i>degenerate</i>: they all tend towards triangles as mesh refinement
1549 continues. This means that the Jacobian matrix of the transformation from
1550 reference cell to actual cell degenerates for these cells, and because
1551 all error estimates for finite element solutions contain the norm of the
1552 inverse of the Jacobian matrix, you will get very large errors on these
1553 cells and, in the limit as mesh refinement, a loss of convergence order because
1554 the cells in these corners become worse and worse under mesh refinement.
1555 
1556 So we need something smarter. To this end, consider the following solution
1557 originally developed by Konstantin Ladutenko. We will use the following code:
1558 @code
1560 
1561 const Point<dim> mesh_center;
1562 const double core_radius = 1.0/5.0,
1563  inner_radius = 1.0/3.0;
1564 
1565 // Step 1: Shrink the inner cell
1566 //
1567 // We cannot get a circle out of the inner cell because of
1568 // the degeneration problem mentioned above. Rather, shrink
1569 // the inner cell to a core radius of 1/5 that stays
1570 // sufficiently far away from the place where the
1571 // coefficient will have a discontinuity and where we want
1572 // to have cell interfaces that actually lie on a circle.
1573 // We do this shrinking by just scaling the location of each
1574 // of the vertices, given that the center of the circle is
1575 // simply the origin of the coordinate system.
1576 for (const auto &cell : triangulation.active_cell_iterators())
1577  if (mesh_center.distance(cell->center()) < 1e-5)
1578  {
1579  for (unsigned int v=0;
1580  v < GeometryInfo<dim>::vertices_per_cell;
1581  ++v)
1582  cell->vertex(v) *= core_radius/mesh_center.distance(cell->vertex(v));
1583  }
1584 
1585 // Step 2: Refine all cells except the central one
1586 for (const auto &cell : triangulation.active_cell_iterators())
1587  if (mesh_center.distance(cell->center()) >= 1e-5)
1588  cell->set_refine_flag();
1589 triangulation.execute_coarsening_and_refinement();
1590 
1591 // Step 3: Resize the inner children of the outer cells
1592 //
1593 // The previous step replaced each of the four outer cells
1594 // by its four children, but the radial distance at which we
1595 // have intersected is not what we want to later refinement
1596 // steps. Consequently, move the vertices that were just
1597 // created in radial direction to a place where we need
1598 // them.
1599 for (const auto &cell : triangulation.active_cell_iterators())
1600  for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1601  {
1602  const double dist = mesh_center.distance(cell->vertex(v));
1603  if (dist > core_radius*1.0001 && dist < 0.9999)
1604  cell->vertex(v) *= inner_radius/dist;
1605  }
1606 
1607 // Step 4: Apply curved manifold description
1608 //
1609 // As discussed above, we can not expect to subdivide the
1610 // inner four cells (or their faces) onto concentric rings,
1611 // but we can do so for all other cells that are located
1612 // outside the inner radius. To this end, we loop over all
1613 // cells and determine whether it is in this zone. If it
1614 // isn't, then we set the manifold description of the cell
1615 // and all of its bounding faces to the one that describes
1616 // the spherical manifold already introduced above and that
1617 // will be used for all further mesh refinement.
1618 for (const auto &cell : triangulation.active_cell_iterators())
1619  {
1620  bool is_in_inner_circle = false;
1621  for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v)
1622  if (mesh_center.distance(cell->vertex(v)) < inner_radius)
1623  {
1624  is_in_inner_circle = true;
1625  break;
1626  }
1627 
1628  if (is_in_inner_circle == false)
1629  // The Triangulation already has a SphericalManifold with
1630  // manifold id 0 (see the documentation of
1631  // GridGenerator::hyper_ball) so we just attach it to the outer
1632  // ring here:
1633  cell->set_all_manifold_ids(0);
1634  }
1635 @endcode
1636 
1637 This code then generates the following, much better sequence of meshes:
1638 
1639 <div class="twocolumn" style="width: 80%">
1640  <div>
1641  <img src="https://www.dealii.org/images/steps/developer/step_6_grid_0_ladutenko.svg"
1642  alt="Initial grid: the Ladutenko grid with one global refinement."
1643  width="300" height="300">
1644  </div>
1645  <div>
1646  <img src="https://www.dealii.org/images/steps/developer/step_6_grid_1_ladutenko.svg"
1647  alt="First adaptively refined Ladutenko grid."
1648  width="300" height="300">
1649  </div>
1650  <div>
1651  <img src="https://www.dealii.org/images/steps/developer/step_6_grid_2_ladutenko.svg"
1652  alt="Second adaptively refined Ladutenko grid."
1653  width="300" height="300">
1654  </div>
1655  <div>
1656  <img src="https://www.dealii.org/images/steps/developer/step_6_grid_3_ladutenko.svg"
1657  alt="Third adaptively refined Ladutenko grid."
1658  width="300" height="300">
1659  </div>
1660  <div>
1661  <img src="https://www.dealii.org/images/steps/developer/step_6_grid_4_ladutenko.svg"
1662  alt="Fourth adaptively refined Ladutenko grid. The cells are clustered
1663  along the inner circle."
1664  width="300" height="300">
1665  </div>
1666  <div>
1667  <img src="https://www.dealii.org/images/steps/developer/step_6_grid_5_ladutenko.svg"
1668  alt="Fifth adaptively refined Ladutenko grid: the cells are clustered
1669  along the inner circle."
1670  width="300" height="300">
1671  </div>
1672 </div>
1673 
1674 Creating good meshes, and in particular making them fit the geometry you
1675 want, is a complex topic in itself. You can find much more on this in
1676 @ref step_49 "step-49", @ref step_53 "step-53", and @ref step_54 "step-54", among other tutorial programs that cover
1677 the issue. @ref step_65 "step-65" shows another, less manual way to achieve a mesh
1678 well fit to the problem here.
1679 Information on curved domains can also be found in the
1680 documentation module on @ref manifold "Manifold descriptions".
1681 
1682 Why does it make sense to choose a mesh that tracks the internal
1683 interface? There are a number of reasons, but the most essential one
1684 comes down to what we actually integrate in our bilinear
1685 form. Conceptually, we want to integrate the term @f$A_{ij}^K=\int_K
1686 a(\mathbf x) \nabla \varphi_i(\mathbf x) \nabla \varphi_j(\mathbf x) ; dx@f$ as the
1687 contribution of cell @f$K@f$ to the matrix entry @f$A_{ij}@f$. We can not
1688 compute it exactly and have to resort to quadrature. We know that
1689 quadrature is accurate if the integrand is smooth. That is because
1690 quadrature in essence computes a polynomial approximation to the
1691 integrand that coincides with the integrand in the quadrature points,
1692 and then computes the volume under this polynomial as an approximation
1693 to the volume under the original integrand. This polynomial
1694 interpolant is accurate if the integrand is smooth on a cell, but it
1695 is usually rather inaccurate if the integrand is discontinuous on a
1696 cell.
1697 
1698 Consequently, it is worthwhile to align cells in such a way that the
1699 interfaces across which the coefficient is discontinuous are aligned
1700 with cell interfaces. This way, the coefficient is constant on each
1701 cell, following which the integrand will be smooth, and its polynomial
1702 approximation and the quadrature approximation of the integral will
1703 both be accurate. Note that such an alignment is common in many
1704 practical cases, so deal.II provides a number of functions (such as
1705 @ref GlossMaterialId "material_id") to help manage such a scenario.
1706 Refer to @ref step_28 "step-28" and @ref step_46 "step-46" for examples of how material ids can be
1707 applied.
1708 
1709 Finally, let us consider the case of a coefficient that has a smooth
1710 and non-uniform distribution in space. We can repeat once again all of
1711 the above discussion on the representation of such a function with the
1712 quadrature. So, to simulate it accurately there are a few readily
1713 available options: you could reduce the cell size, increase the order
1714 of the polynomial used in the quadrature formula, select a more
1715 appropriate quadrature formula, or perform a combination of these
1716 steps. The key is that providing the best fit of the coefficient's
1717 spatial dependence with the quadrature polynomial will lead to a more
1718 accurate finite element solution of the PDE.
1719 
1720 As a final note: The discussion in the previous paragraphs shows, we here
1721 have a very concrete way of stating what we think of a good mesh -- it should
1722 be aligned with the jump in the coefficient. But one could also have asked
1723 this kind of question in a more general setting: Given some equation with
1724 a smooth solution and smooth coefficients, can we say what a good mesh
1725 would look like? This is a question for which the answer is easier to state
1726 in intuitive terms than mathematically: A good mesh has cells that all,
1727 by and large, look like squares (or cubes, in 3d). A bad mesh would contain
1728 cells that are very elongated in some directions or, more generally, for which
1729 there are cells that have both short and long edges. There are many ways
1730 in which one could assign a numerical quality index to each cell that measures
1731 whether the cell is "good" or "bad"; some of these are often chosen because
1732 they are cheap and easy to compute, whereas others are based on what enters
1733 into proofs of convergence. An example of the former would be the ratio of
1734 the longest to the shortest edge of a cell: In the ideal case, that ratio
1735 would be one; bad cells have values much larger than one. An example of the
1736 latter kind would consider the gradient (the "Jacobian") of the mapping
1737 from the reference cell @f$\hat K=[0,1]^d@f$ to the real cell @f$K@f$; this
1738 gradient is a matrix, and a quantity that enters into error estimates
1739 is the maximum over all points on the reference cell of the ratio of the
1740 largest to the smallest eigenvalue of this matrix. It is again not difficult
1741 to see that this ratio is constant if the cell @f$K@f$ is an affine image of
1742 @f$\hat K@f$, and that it is one for squares and cubes.
1743 
1744 In practice, it might be interesting to visualize such quality measures.
1745 The function GridTools::compute_aspect_ratio_of_cells() provides one
1746 way to get this kind of information. Even better, visualization tools
1747 such as VisIt often allow you to visualize this sort of information
1748 for a variety of measures from within the visualization software
1749 itself; in the case of VisIt, just add a "pseudo-color" plot and select
1750 one of the mesh quality measures instead of the solution field.
1751 
1752 
1753 <a name="Playingwiththeregularityofthesolution"></a><h4>Playing with the regularity of the solution</h4>
1754 
1755 
1756 From a mathematical perspective, solutions of the Laplace equation
1757 @f[
1758  -\Delta u = f
1759 @f]
1760 on smoothly bounded, convex domains are known to be smooth themselves. The exact degree
1761 of smoothness, i.e., the function space in which the solution lives, depends
1762 on how smooth exactly the boundary of the domain is, and how smooth the right
1763 hand side is. Some regularity of the solution may be lost at the boundary, but
1764 one generally has that the solution is twice more differentiable in
1765 compact subsets of the domain than the right hand side.
1766 If, in particular, the right hand side satisfies @f$f\in C^\infty(\Omega)@f$, then
1767 @f$u \in C^\infty(\Omega_i)@f$ where @f$\Omega_i@f$ is any compact subset of @f$\Omega@f$
1768 (@f$\Omega@f$ is an open domain, so a compact subset needs to keep a positive distance
1769 from @f$\partial\Omega@f$).
1770 
1771 The situation we chose for the current example is different, however: we look
1772 at an equation with a non-constant coefficient @f$a(\mathbf x)@f$:
1773 @f[
1774  -\nabla \cdot (a \nabla u) = f.
1775 @f]
1776 Here, if @f$a@f$ is not smooth, then the solution will not be smooth either,
1777 regardless of @f$f@f$. In particular, we expect that wherever @f$a@f$ is discontinuous
1778 along a line (or along a plane in 3d),
1779 the solution will have a kink. This is easy to see: if for example @f$f@f$
1780 is continuous, then @f$f=-\nabla \cdot (a \nabla u)@f$ needs to be
1781 continuous. This means that @f$a \nabla u@f$ must be continuously differentiable
1782 (not have a kink). Consequently, if @f$a@f$ has a discontinuity, then @f$\nabla u@f$
1783 must have an opposite discontinuity so that the two exactly cancel and their
1784 product yields a function without a discontinuity. But for @f$\nabla u@f$ to have
1785 a discontinuity, @f$u@f$ must have a kink. This is of course exactly what is
1786 happening in the current example, and easy to observe in the pictures of the
1787 solution.
1788 
1789 In general, if the coefficient @f$a(\mathbf x)@f$ is discontinuous along a line in 2d,
1790 or a plane in 3d, then the solution may have a kink, but the gradient of the
1791 solution will not go to infinity. That means, that the solution is at least
1792 still in the <a href="https://en.wikipedia.org/wiki/Sobolev_space">Sobolev space</a>
1793 @f$W^{1,\infty}@f$ (i.e., roughly speaking, in the
1794 space of functions whose derivatives are bounded). On the other hand,
1795 we know that in the most
1796 extreme cases -- i.e., where the domain has reentrant corners, the
1797 right hand side only satisfies @f$f\in H^{-1}@f$, or the coefficient @f$a@f$ is only in
1798 @f$L^\infty@f$ -- all we can expect is that @f$u\in H^1@f$ (i.e., the
1799 <a
1800 href="https://en.wikipedia.org/wiki/Sobolev_space#Sobolev_spaces_with_integer_k">Sobolev
1801 space</a> of functions whose derivative is square integrable), a much larger space than
1802 @f$W^{1,\infty}@f$. It is not very difficult to create cases where
1803 the solution is in a space @f$H^{1+s}@f$ where we can get @f$s@f$ to become as small
1804 as we want. Such cases are often used to test adaptive finite element
1805 methods because the mesh will have to resolve the singularity that causes
1806 the solution to not be in @f$W^{1,\infty}@f$ any more.
1807 
1808 The typical example one uses for this is called the <i>Kellogg problem</i>
1809 (referring to @cite Kel74), which in the commonly used form has a coefficient
1810 @f$a(\mathbf x)@f$ that has different values in the four quadrants of the plane
1811 (or eight different values in the octants of @f${\mathbb R}^3@f$). The exact degree
1812 of regularity (the @f$s@f$ in the index of the Sobolev space above) depends on the
1813 values of @f$a(\mathbf x)@f$ coming together at the origin, and by choosing the
1814 jumps large enough, the regularity of the solution can be made as close as
1815 desired to @f$H^1@f$.
1816 
1817 To implement something like this, one could replace the coefficient
1818 function by the following (shown here only for the 2d case):
1819 @code
1820 template <int dim>
1821 double coefficient (const Point<dim> &p)
1822 {
1823  if ((p[0] < 0) && (p[1] < 0)) // lower left quadrant
1824  return 1;
1825  else if ((p[0] >= 0) && (p[1] < 0)) // lower right quadrant
1826  return 10;
1827  else if ((p[0] < 0) && (p[1] >= 0)) // upper left quadrant
1828  return 100;
1829  else if ((p[0] >= 0) && (p[1] >= 0)) // upper right quadrant
1830  return 1000;
1831  else
1832  {
1833  Assert(false, ExcInternalError());
1834  return 0;
1835  }
1836 }
1837 @endcode
1838 (Adding the <code>Assert</code> at the end ensures that either an exception
1839 is thrown or that the program aborts if we ever get to that point
1840 -- which of course we shouldn't,
1841 but this is a good way to insure yourself: we all make mistakes by
1842 sometimes not thinking of all cases, for example by checking
1843 for <code>p[0]</code> to be less than and greater than zero,
1844 rather than greater-or-equal to zero, and thereby forgetting
1845 some cases that would otherwise lead to bugs that are awkward
1846 to find. The <code>return 0;</code> at the end is only there to
1847 avoid compiler warnings that the function does not end in a
1848 <code>return</code> statement -- the compiler cannot see that the
1849 function would never actually get to that point because of the
1850 preceding <code>Assert</code> statement.)
1851 
1852 By playing with such cases where four or more sectors come
1853 together and on which the coefficient has different values, one can
1854 construct cases where the solution has singularities at the
1855 origin. One can also see how the meshes are refined in such cases.
1856  *
1857  *
1858 <a name="PlainProg"></a>
1859 <h1> The plain program</h1>
1860 @include "step-6.cc"
1861 */
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition: lapack_support.h:121
GridOutFlags::Gnuplot
Definition: grid_out.h:234
DoFTools::nonzero
@ nonzero
Definition: dof_tools.h:244
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
DataOutInterface::write_vtu
void write_vtu(std::ostream &out) const
Definition: data_out_base.cc:6864
SolverCG
Definition: solver_cg.h:98
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
Triangulation
Definition: tria.h:1109
MappingQGeneric
Definition: mapping_q_generic.h:135
DataOutBase::dx
@ dx
Definition: data_out_base.h:1562
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
SphericalManifold
Definition: manifold_lib.h:231
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
DoFTools::always
@ always
Definition: dof_tools.h:236
LocalIntegrators::Advection::cell_matrix
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:80
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
KellyErrorEstimator
Definition: error_estimator.h:262
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
GridOut::set_flags
void set_flags(const GridOutFlags::DX &flags)
Definition: grid_out.cc:491
AffineConstraints::distribute_local_to_global
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
Definition: affine_constraints.h:1859
TimeStepping::invalid
@ invalid
Definition: time_stepping.h:71
GridOut::write_gnuplot
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:4337
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
LAPACKSupport::general
@ general
No special properties.
Definition: lapack_support.h:113
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
GridGenerator
Definition: grid_generator.h:57
QGauss
Definition: quadrature_lib.h:40
SIMDComparison::equal
@ equal
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
value
static const bool value
Definition: dof_tools_constraints.cc:433
GridGenerator::hyper_ball
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)
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
AffineConstraints
Definition: affine_constraints.h:180
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
GridTools::diameter
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:76
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
DataOutBase::gnuplot
@ gnuplot
Definition: data_out_base.h:1572
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
Point< dim >
PreconditionSSOR
Definition: precondition.h:665
FullMatrix< double >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
LAPACKSupport::N
static const char N
Definition: lapack_support.h:159
SolverControl
Definition: solver_control.h:67
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type 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())
Definition: loop.h:443
first
Point< 2 > first
Definition: grid_out.cc:4352
GridOut
Definition: grid_out.h:1020
DataOut< dim >
Vector< double >
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41
center
Point< 3 > center
Definition: data_out_base.cc:169
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
DataOut_DoFData::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090