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