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-27.h
Go to the documentation of this file.
1 ) const
938  * {
939  * double product = 1;
940  * for (unsigned int d = 0; d < dim; ++d)
941  * product *= (p[d] + 1);
942  * return product;
943  * }
944  *
945  *
946  *
947  * @endcode
948  *
949  *
950  * <a name="Implementationofthemainclass"></a>
951  * <h3>Implementation of the main class</h3>
952  *
953 
954  *
955  *
956  * <a name="LaplaceProblemLaplaceProblemconstructor"></a>
957  * <h4>LaplaceProblem::LaplaceProblem constructor</h4>
958  *
959 
960  *
961  * The constructor of this class is fairly straightforward. It associates
962  * the hp::DoFHandler object with the triangulation, and then sets the
963  * maximal polynomial degree to 7 (in 1d and 2d) or 5 (in 3d and higher). We
964  * do so because using higher order polynomial degrees becomes prohibitively
965  * expensive, especially in higher space dimensions.
966  *
967 
968  *
969  * Following this, we fill the collections of finite element, and cell and
970  * face quadrature objects. We start with quadratic elements, and each
971  * quadrature formula is chosen so that it is appropriate for the matching
972  * finite element in the hp::FECollection object.
973  *
974 
975  *
976  * Finally, we initialize FESeries::Fourier object which will be used to
977  * calculate coefficient in Fourier series as described in the introduction.
978  * In addition to the hp::FECollection, we need to provide quadrature rules
979  * hp::QCollection for integration on the reference cell.
980  *
981 
982  *
983  * In order to resize fourier_coefficients Table, we use the following
984  * auxiliary function
985  *
986  * @code
987  * template <int dim, typename T>
988  * void resize(Table<dim, T> &coeff, const unsigned int N)
989  * {
990  * TableIndices<dim> size;
991  * for (unsigned int d = 0; d < dim; d++)
992  * size[d] = N;
993  * coeff.reinit(size);
994  * }
995  *
996  * template <int dim>
997  * LaplaceProblem<dim>::LaplaceProblem()
998  * : dof_handler(triangulation)
999  * , max_degree(dim <= 2 ? 7 : 5)
1000  * {
1001  * for (unsigned int degree = 2; degree <= max_degree; ++degree)
1002  * {
1003  * fe_collection.push_back(FE_Q<dim>(degree));
1004  * quadrature_collection.push_back(QGauss<dim>(degree + 1));
1005  * face_quadrature_collection.push_back(QGauss<dim - 1>(degree + 1));
1006  * }
1007  *
1008  * @endcode
1009  *
1010  * As described in the introduction, we define the Fourier vectors @f${\bf
1011  * k}@f$ for which we want to compute Fourier coefficients of the solution
1012  * on each cell as follows. In 2d, we will need coefficients corresponding
1013  * to vectors @f${\bf k}=(2 \pi i, 2\pi j)^T@f$ for which @f$\sqrt{i^2+j^2}\le N@f$,
1014  * with @f$i,j@f$ integers and @f$N@f$ being the maximal polynomial degree we use
1015  * for the finite elements in this program. The FESeries::Fourier class'
1016  * constructor first parameter @f$N@f$ defines the number of coefficients in 1D
1017  * with the total number of coefficients being @f$N^{dim}@f$. Although we will
1018  * not use coefficients corresponding to
1019  * @f$\sqrt{i^2+j^2}> N@f$ and @f$i+j==0@f$, the overhead of their calculation is
1020  * minimal. The transformation matrices for each FiniteElement will be
1021  * calculated only once the first time they are required in the course of
1022  * hp-adaptive refinement. Because we work on the unit cell, we can do all
1023  * this work without a mapping from reference to real cell and consequently
1024  * can precalculate these matrices. The calculation of expansion
1025  * coefficients for a particular set of local degrees of freedom on a given
1026  * cell then follows as a simple matrix-vector product.
1027  * The 3d case is handled analogously.
1028  *
1029  * @code
1030  * const unsigned int N = max_degree;
1031  *
1032  * @endcode
1033  *
1034  * We will need to assemble the matrices that do the Fourier transforms
1035  * for each of the finite elements we deal with, i.e. the matrices @f${\cal
1036  * F}_{{\bf k},j}@f$ defined in the introduction. We have to do that for
1037  * each of the finite elements in use. To that end we need a quadrature
1038  * rule. In this example we use the same quadrature formula for each
1039  * finite element, namely that is obtained by iterating a
1040  * 2-point Gauss formula as many times as the maximal exponent we use for
1041  * the term @f$e^{i{\bf k}\cdot{\bf x}}@f$:
1042  *
1043  * @code
1044  * QGauss<1> base_quadrature(2);
1045  * QIterated<dim> quadrature(base_quadrature, N);
1046  * for (unsigned int i = 0; i < fe_collection.size(); i++)
1047  * fourier_q_collection.push_back(quadrature);
1048  *
1049  * @endcode
1050  *
1051  * Now we are ready to set-up the FESeries::Fourier object
1052  *
1053  * @code
1054  * const std::vector<unsigned int> n_coefficients_per_direction(
1055  * fe_collection.size(), N);
1056  * fourier = std_cxx14::make_unique<FESeries::Fourier<dim>>(
1057  * n_coefficients_per_direction, fe_collection, fourier_q_collection);
1058  *
1059  * @endcode
1060  *
1061  * We need to resize the matrix of fourier coefficients according to the
1062  * number of modes N.
1063  *
1064  * @code
1065  * resize(fourier_coefficients, N);
1066  * }
1067  *
1068  *
1069  * @endcode
1070  *
1071  *
1072  * <a name="LaplaceProblemLaplaceProblemdestructor"></a>
1073  * <h4>LaplaceProblem::~LaplaceProblem destructor</h4>
1074  *
1075 
1076  *
1077  * The destructor is unchanged from what we already did in @ref step_6 "step-6":
1078  *
1079  * @code
1080  * template <int dim>
1081  * LaplaceProblem<dim>::~LaplaceProblem()
1082  * {
1083  * dof_handler.clear();
1084  * }
1085  *
1086  *
1087  * @endcode
1088  *
1089  *
1090  * <a name="LaplaceProblemsetup_system"></a>
1091  * <h4>LaplaceProblem::setup_system</h4>
1092  *
1093 
1094  *
1095  * This function is again a verbatim copy of what we already did in
1096  * @ref step_6 "step-6". Despite function calls with exactly the same names and arguments,
1097  * the algorithms used internally are different in some aspect since the
1098  * dof_handler variable here is an hp object.
1099  *
1100  * @code
1101  * template <int dim>
1102  * void LaplaceProblem<dim>::setup_system()
1103  * {
1104  * dof_handler.distribute_dofs(fe_collection);
1105  *
1106  * solution.reinit(dof_handler.n_dofs());
1107  * system_rhs.reinit(dof_handler.n_dofs());
1108  *
1109  * constraints.clear();
1110  * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1111  * VectorTools::interpolate_boundary_values(dof_handler,
1112  * 0,
1113  * Functions::ZeroFunction<dim>(),
1114  * constraints);
1115  * constraints.close();
1116  *
1117  * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1118  * DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
1119  * sparsity_pattern.copy_from(dsp);
1120  *
1121  * system_matrix.reinit(sparsity_pattern);
1122  * }
1123  *
1124  *
1125  *
1126  * @endcode
1127  *
1128  *
1129  * <a name="LaplaceProblemassemble_system"></a>
1130  * <h4>LaplaceProblem::assemble_system</h4>
1131  *
1132 
1133  *
1134  * This is the function that assembles the global matrix and right hand side
1135  * vector from the local contributions of each cell. Its main working is as
1136  * has been described in many of the tutorial programs before. The
1137  * significant deviations are the ones necessary for <i>hp</i> finite
1138  * element methods. In particular, that we need to use a collection of
1139  * FEValues object (implemented through the hp::FEValues class), and that we
1140  * have to eliminate constrained degrees of freedom already when copying
1141  * local contributions into global objects. Both of these are explained in
1142  * detail in the introduction of this program.
1143  *
1144 
1145  *
1146  * One other slight complication is the fact that because we use different
1147  * polynomial degrees on different cells, the matrices and vectors holding
1148  * local contributions do not have the same size on all cells. At the
1149  * beginning of the loop over all cells, we therefore each time have to
1150  * resize them to the correct size (given by
1151  * <code>dofs_per_cell</code>). Because these classes are implement in such
1152  * a way that reducing the size of a matrix or vector does not release the
1153  * currently allocated memory (unless the new size is zero), the process of
1154  * resizing at the beginning of the loop will only require re-allocation of
1155  * memory during the first few iterations. Once we have found in a cell with
1156  * the maximal finite element degree, no more re-allocations will happen
1157  * because all subsequent <code>reinit</code> calls will only set the size
1158  * to something that fits the currently allocated memory. This is important
1159  * since allocating memory is expensive, and doing so every time we visit a
1160  * new cell would take significant compute time.
1161  *
1162  * @code
1163  * template <int dim>
1164  * void LaplaceProblem<dim>::assemble_system()
1165  * {
1166  * hp::FEValues<dim> hp_fe_values(fe_collection,
1167  * quadrature_collection,
1168  * update_values | update_gradients |
1169  * update_quadrature_points |
1170  * update_JxW_values);
1171  *
1172  * RightHandSide<dim> rhs_function;
1173  *
1174  * FullMatrix<double> cell_matrix;
1175  * Vector<double> cell_rhs;
1176  *
1177  * std::vector<types::global_dof_index> local_dof_indices;
1178  *
1179  * for (const auto &cell : dof_handler.active_cell_iterators())
1180  * {
1181  * const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
1182  *
1183  * cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1184  * cell_matrix = 0;
1185  *
1186  * cell_rhs.reinit(dofs_per_cell);
1187  * cell_rhs = 0;
1188  *
1189  * hp_fe_values.reinit(cell);
1190  *
1191  * const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
1192  *
1193  * std::vector<double> rhs_values(fe_values.n_quadrature_points);
1194  * rhs_function.value_list(fe_values.get_quadrature_points(), rhs_values);
1195  *
1196  * for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
1197  * ++q_point)
1198  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1199  * {
1200  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1201  * cell_matrix(i, j) +=
1202  * (fe_values.shape_grad(i, q_point) * // grad phi_i(x_q)
1203  * fe_values.shape_grad(j, q_point) * // grad phi_j(x_q)
1204  * fe_values.JxW(q_point)); // dx
1205  *
1206  * cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
1207  * rhs_values[q_point] * // f(x_q)
1208  * fe_values.JxW(q_point)); // dx
1209  * }
1210  *
1211  * local_dof_indices.resize(dofs_per_cell);
1212  * cell->get_dof_indices(local_dof_indices);
1213  *
1214  * constraints.distribute_local_to_global(
1215  * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1216  * }
1217  * }
1218  *
1219  *
1220  *
1221  * @endcode
1222  *
1223  *
1224  * <a name="LaplaceProblemsolve"></a>
1225  * <h4>LaplaceProblem::solve</h4>
1226  *
1227 
1228  *
1229  * The function solving the linear system is entirely unchanged from
1230  * previous examples. We simply try to reduce the initial residual (which
1231  * equals the @f$l_2@f$ norm of the right hand side) by a certain factor:
1232  *
1233  * @code
1234  * template <int dim>
1235  * void LaplaceProblem<dim>::solve()
1236  * {
1237  * SolverControl solver_control(system_rhs.size(),
1238  * 1e-12 * system_rhs.l2_norm());
1239  * SolverCG<Vector<double>> cg(solver_control);
1240  *
1241  * PreconditionSSOR<SparseMatrix<double>> preconditioner;
1242  * preconditioner.initialize(system_matrix, 1.2);
1243  *
1244  * cg.solve(system_matrix, solution, system_rhs, preconditioner);
1245  *
1246  * constraints.distribute(solution);
1247  * }
1248  *
1249  *
1250  *
1251  * @endcode
1252  *
1253  *
1254  * <a name="LaplaceProblempostprocess"></a>
1255  * <h4>LaplaceProblem::postprocess</h4>
1256  *
1257 
1258  *
1259  * After solving the linear system, we will want to postprocess the
1260  * solution. Here, all we do is to estimate the error, estimate the local
1261  * smoothness of the solution as described in the introduction, then write
1262  * graphical output, and finally refine the mesh in both @f$h@f$ and @f$p@f$
1263  * according to the indicators computed before. We do all this in the same
1264  * function because we want the estimated error and smoothness indicators
1265  * not only for refinement, but also include them in the graphical output.
1266  *
1267  * @code
1268  * template <int dim>
1269  * void LaplaceProblem<dim>::postprocess(const unsigned int cycle)
1270  * {
1271  * @endcode
1272  *
1273  * Let us start with computing estimated error and smoothness indicators,
1274  * which each are one number for each active cell of our
1275  * triangulation. For the error indicator, we use the KellyErrorEstimator
1276  * class as always. Estimating the smoothness is done in the respective
1277  * function of this class; that function is discussed further down below:
1278  *
1279  * @code
1280  * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1281  * KellyErrorEstimator<dim>::estimate(
1282  * dof_handler,
1283  * face_quadrature_collection,
1284  * std::map<types::boundary_id, const Function<dim> *>(),
1285  * solution,
1286  * estimated_error_per_cell);
1287  *
1288  *
1289  * Vector<float> smoothness_indicators(triangulation.n_active_cells());
1290  * estimate_smoothness(smoothness_indicators);
1291  *
1292  * @endcode
1293  *
1294  * Next we want to generate graphical output. In addition to the two
1295  * estimated quantities derived above, we would also like to output the
1296  * polynomial degree of the finite elements used on each of the elements
1297  * on the mesh.
1298  *
1299 
1300  *
1301  * The way to do that requires that we loop over all cells and poll the
1302  * active finite element index of them using
1303  * <code>cell-@>active_fe_index()</code>. We then use the result of this
1304  * operation and query the finite element collection for the finite
1305  * element with that index, and finally determine the polynomial degree of
1306  * that element. The result we put into a vector with one element per
1307  * cell. The DataOut class requires this to be a vector of
1308  * <code>float</code> or <code>double</code>, even though our values are
1309  * all integers, so that it what we use:
1310  *
1311  * @code
1312  * {
1313  * Vector<float> fe_degrees(triangulation.n_active_cells());
1314  * for (const auto &cell : dof_handler.active_cell_iterators())
1315  * fe_degrees(cell->active_cell_index()) =
1316  * fe_collection[cell->active_fe_index()].degree;
1317  *
1318  * @endcode
1319  *
1320  * With now all data vectors available -- solution, estimated errors and
1321  * smoothness indicators, and finite element degrees --, we create a
1322  * DataOut object for graphical output and attach all data. Note that
1323  * the DataOut class has a second template argument (which defaults to
1324  * DoFHandler@<dim@>, which is why we have never seen it in previous
1325  * tutorial programs) that indicates the type of DoF handler to be
1326  * used. Here, we have to use the hp::DoFHandler class:
1327  *
1328  * @code
1329  * DataOut<dim, hp::DoFHandler<dim>> data_out;
1330  *
1331  * data_out.attach_dof_handler(dof_handler);
1332  * data_out.add_data_vector(solution, "solution");
1333  * data_out.add_data_vector(estimated_error_per_cell, "error");
1334  * data_out.add_data_vector(smoothness_indicators, "smoothness");
1335  * data_out.add_data_vector(fe_degrees, "fe_degree");
1336  * data_out.build_patches();
1337  *
1338  * @endcode
1339  *
1340  * The final step in generating output is to determine a file name, open
1341  * the file, and write the data into it (here, we use VTK format):
1342  *
1343  * @code
1344  * const std::string filename =
1345  * "solution-" + Utilities::int_to_string(cycle, 2) + ".vtk";
1346  * std::ofstream output(filename);
1347  * data_out.write_vtk(output);
1348  * }
1349  *
1350  * @endcode
1351  *
1352  * After this, we would like to actually refine the mesh, in both @f$h@f$ and
1353  * @f$p@f$. The way we are going to do this is as follows: first, we use the
1354  * estimated error to flag those cells for refinement that have the
1355  * largest error. This is what we have always done:
1356  *
1357  * @code
1358  * {
1359  * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1360  * estimated_error_per_cell,
1361  * 0.3,
1362  * 0.03);
1363  *
1364  * @endcode
1365  *
1366  * Next we would like to figure out which of the cells that have been
1367  * flagged for refinement should actually have @f$p@f$ increased instead of
1368  * @f$h@f$ decreased. The strategy we choose here is that we look at the
1369  * smoothness indicators of those cells that are flagged for refinement,
1370  * and increase @f$p@f$ for those with a smoothness larger than a certain
1371  * threshold. For this, we first have to determine the maximal and
1372  * minimal values of the smoothness indicators of all flagged cells,
1373  * which we do using a loop over all cells and comparing current minimal
1374  * and maximal values. (We start with the minimal and maximal values of
1375  * <i>all</i> cells, a range within which the minimal and maximal values
1376  * on cells flagged for refinement must surely lie.) Absent any better
1377  * strategies, we will then set the threshold above which will increase
1378  * @f$p@f$ instead of reducing @f$h@f$ as the mean value between minimal and
1379  * maximal smoothness indicators on cells flagged for refinement:
1380  *
1381  * @code
1382  * float max_smoothness = *std::min_element(smoothness_indicators.begin(),
1383  * smoothness_indicators.end()),
1384  * min_smoothness = *std::max_element(smoothness_indicators.begin(),
1385  * smoothness_indicators.end());
1386  * for (const auto &cell : dof_handler.active_cell_iterators())
1387  * if (cell->refine_flag_set())
1388  * {
1389  * max_smoothness =
1390  * std::max(max_smoothness,
1391  * smoothness_indicators(cell->active_cell_index()));
1392  * min_smoothness =
1393  * std::min(min_smoothness,
1394  * smoothness_indicators(cell->active_cell_index()));
1395  * }
1396  * const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;
1397  *
1398  * @endcode
1399  *
1400  * With this, we can go back, loop over all cells again, and for those
1401  * cells for which (i) the refinement flag is set, (ii) the smoothness
1402  * indicator is larger than the threshold, and (iii) we still have a
1403  * finite element with a polynomial degree higher than the current one
1404  * in the finite element collection, we then increase the polynomial
1405  * degree and in return remove the flag indicating that the cell should
1406  * undergo bisection. For all other cells, the refinement flags remain
1407  * untouched:
1408  *
1409  * @code
1410  * for (const auto &cell : dof_handler.active_cell_iterators())
1411  * if (cell->refine_flag_set() &&
1412  * (smoothness_indicators(cell->active_cell_index()) >
1413  * threshold_smoothness) &&
1414  * (cell->active_fe_index() + 1 < fe_collection.size()))
1415  * {
1416  * cell->clear_refine_flag();
1417  * cell->set_active_fe_index(cell->active_fe_index() + 1);
1418  * }
1419  *
1420  * @endcode
1421  *
1422  * At the end of this procedure, we then refine the mesh. During this
1423  * process, children of cells undergoing bisection inherit their mother
1424  * cell's finite element index:
1425  *
1426  * @code
1427  * triangulation.execute_coarsening_and_refinement();
1428  * }
1429  * }
1430  *
1431  *
1432  * @endcode
1433  *
1434  *
1435  * <a name="LaplaceProblemcreate_coarse_grid"></a>
1436  * <h4>LaplaceProblem::create_coarse_grid</h4>
1437  *
1438 
1439  *
1440  * The following function is used when creating the initial grid. It is a
1441  * specialization for the 2d case, i.e. a corresponding function needs to be
1442  * implemented if the program is run in anything other then 2d. The function
1443  * is actually stolen from @ref step_14 "step-14" and generates the same mesh used already
1444  * there, i.e. the square domain with the square hole in the middle. The
1445  * meaning of the different parts of this function are explained in the
1446  * documentation of @ref step_14 "step-14":
1447  *
1448  * @code
1449  * template <>
1450  * void LaplaceProblem<2>::create_coarse_grid()
1451  * {
1452  * const unsigned int dim = 2;
1453  *
1454  * const std::vector<Point<2>> vertices = {
1455  * {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0},
1456  * {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5},
1457  * {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0},
1458  * {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5},
1459  * {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}};
1460  *
1461  * const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
1462  * cell_vertices = {{{0, 1, 5, 6}},
1463  * {{1, 2, 6, 7}},
1464  * {{2, 3, 7, 8}},
1465  * {{3, 4, 8, 9}},
1466  * {{5, 6, 10, 11}},
1467  * {{8, 9, 12, 13}},
1468  * {{10, 11, 14, 15}},
1469  * {{12, 13, 17, 18}},
1470  * {{14, 15, 19, 20}},
1471  * {{15, 16, 20, 21}},
1472  * {{16, 17, 21, 22}},
1473  * {{17, 18, 22, 23}}};
1474  *
1475  * const unsigned int n_cells = cell_vertices.size();
1476  *
1477  * std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
1478  * for (unsigned int i = 0; i < n_cells; ++i)
1479  * {
1480  * for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
1481  * cells[i].vertices[j] = cell_vertices[i][j];
1482  * cells[i].material_id = 0;
1483  * }
1484  *
1485  * triangulation.create_triangulation(vertices, cells, SubCellData());
1486  * triangulation.refine_global(3);
1487  * }
1488  *
1489  *
1490  *
1491  * @endcode
1492  *
1493  *
1494  * <a name="LaplaceProblemrun"></a>
1495  * <h4>LaplaceProblem::run</h4>
1496  *
1497 
1498  *
1499  * This function implements the logic of the program, as did the respective
1500  * function in most of the previous programs already, see for example
1501  * @ref step_6 "step-6".
1502  *
1503 
1504  *
1505  * Basically, it contains the adaptive loop: in the first iteration create a
1506  * coarse grid, and then set up the linear system, assemble it, solve, and
1507  * postprocess the solution including mesh refinement. Then start over
1508  * again. In the meantime, also output some information for those staring at
1509  * the screen trying to figure out what the program does:
1510  *
1511  * @code
1512  * template <int dim>
1513  * void LaplaceProblem<dim>::run()
1514  * {
1515  * for (unsigned int cycle = 0; cycle < 6; ++cycle)
1516  * {
1517  * std::cout << "Cycle " << cycle << ':' << std::endl;
1518  *
1519  * if (cycle == 0)
1520  * create_coarse_grid();
1521  *
1522  * setup_system();
1523  *
1524  * std::cout << " Number of active cells : "
1525  * << triangulation.n_active_cells() << std::endl
1526  * << " Number of degrees of freedom: " << dof_handler.n_dofs()
1527  * << std::endl
1528  * << " Number of constraints : "
1529  * << constraints.n_constraints() << std::endl;
1530  *
1531  * assemble_system();
1532  * solve();
1533  * postprocess(cycle);
1534  * }
1535  * }
1536  *
1537  *
1538  * @endcode
1539  *
1540  *
1541  * <a name="LaplaceProblemestimate_smoothness"></a>
1542  * <h4>LaplaceProblem::estimate_smoothness</h4>
1543  *
1544 
1545  *
1546  * As described in the introduction, we will need to take the maximum
1547  * absolute value of fourier coefficients which correspond to @f$k@f$-vector
1548  * @f$|{\bf k}|= const@f$. To filter the coefficients Table we
1549  * will use the FESeries::process_coefficients() which requires a predicate
1550  * to be specified. The predicate should operate on TableIndices and return
1551  * a pair of <code>bool</code> and <code>unsigned int</code>. The latter
1552  * is the value of the map from TableIndicies to unsigned int. It is
1553  * used to define subsets of coefficients from which we search for the one
1554  * with highest absolute value, i.e. @f$l^\infty@f$-norm. The <code>bool</code>
1555  * parameter defines which indices should be used in processing. In the
1556  * current case we are interested in coefficients which correspond to
1557  * @f$0 < i*i+j*j < N*N@f$ and @f$0 < i*i+j*j+k*k < N*N@f$ in 2D and 3D, respectively.
1558  *
1559  * @code
1560  * template <int dim>
1561  * std::pair<bool, unsigned int>
1562  * LaplaceProblem<dim>::predicate(const TableIndices<dim> &ind)
1563  * {
1564  * unsigned int v = 0;
1565  * for (unsigned int i = 0; i < dim; i++)
1566  * v += ind[i] * ind[i];
1567  * if (v > 0 && v < max_degree * max_degree)
1568  * return std::make_pair(true, v);
1569  * else
1570  * return std::make_pair(false, v);
1571  * }
1572  *
1573  * @endcode
1574  *
1575  * This last function of significance implements the algorithm to estimate
1576  * the smoothness exponent using the algorithms explained in detail in the
1577  * introduction. We will therefore only comment on those points that are of
1578  * implementational importance.
1579  *
1580  * @code
1581  * template <int dim>
1582  * void
1583  * LaplaceProblem<dim>::estimate_smoothness(Vector<float> &smoothness_indicators)
1584  * {
1585  * @endcode
1586  *
1587  * Since most of the hard work is done for us in FESeries::Fourier and
1588  * we set up the object of this class in the constructor, what we are left
1589  * to do here is apply this class to calculate coefficients and then
1590  * perform linear regression to fit their decay slope.
1591  *
1592 
1593  *
1594  *
1595 
1596  *
1597  * First thing to do is to loop over all cells and do our work there, i.e.
1598  * to locally do the Fourier transform and estimate the decay coefficient.
1599  * We will use the following array as a scratch array in the loop to store
1600  * local DoF values:
1601  *
1602  * @code
1603  * Vector<double> local_dof_values;
1604  *
1605  * @endcode
1606  *
1607  * Then here is the loop:
1608  *
1609  * @code
1610  * for (const auto &cell : dof_handler.active_cell_iterators())
1611  * {
1612  * @endcode
1613  *
1614  * Inside the loop, we first need to get the values of the local
1615  * degrees of freedom (which we put into the
1616  * <code>local_dof_values</code> array after setting it to the right
1617  * size) and then need to compute the Fourier transform by multiplying
1618  * this vector with the matrix @f${\cal F}@f$ corresponding to this finite
1619  * element. This is done by calling FESeries::Fourier::calculate(),
1620  * that has to be provided with the <code>local_dof_values</code>,
1621  * <code>cell->active_fe_index()</code> and a Table to store
1622  * coefficients.
1623  *
1624  * @code
1625  * local_dof_values.reinit(cell->get_fe().dofs_per_cell);
1626  * cell->get_dof_values(solution, local_dof_values);
1627  *
1628  * fourier->calculate(local_dof_values,
1629  * cell->active_fe_index(),
1630  * fourier_coefficients);
1631  *
1632  * @endcode
1633  *
1634  * The next thing, as explained in the introduction, is that we wanted
1635  * to only fit our exponential decay of Fourier coefficients to the
1636  * largest coefficients for each possible value of @f$|{\bf k}|@f$. To
1637  * this end, we use FESeries::process_coefficients() to rework
1638  * coefficients into the desired format. We'll only take those Fourier
1639  * coefficients with the largest magnitude for a given value of @f$|{\bf
1640  * k}|@f$ and thereby need to use VectorTools::Linfty_norm:
1641  *
1642  * @code
1643  * std::pair<std::vector<unsigned int>, std::vector<double>> res =
1644  * FESeries::process_coefficients<dim>(
1645  * fourier_coefficients,
1646  * [this](const TableIndices<dim> &indices) {
1647  * return this->predicate(indices);
1648  * },
1650  *
1651  * Assert(res.first.size() == res.second.size(), ExcInternalError());
1652  *
1653  * @endcode
1654  *
1655  * The first vector in the <code>std::pair</code> will store values of
1656  * the predicate, that is @f$i*i+j*j= const@f$ or @f$i*i+j*j+k*k = const@f$ in
1657  * 2D or 3D respectively. This vector will be the same for all the cells
1658  * so we can calculate logarithms of the corresponding Fourier vectors
1659  * @f$|{\bf k}|@f$ only once in the whole hp-refinement cycle:
1660  *
1661  * @code
1662  * if (ln_k.size() == 0)
1663  * {
1664  * ln_k.resize(res.first.size(), 0);
1665  * for (unsigned int f = 0; f < ln_k.size(); f++)
1666  * ln_k[f] =
1667  * std::log(2.0 * numbers::PI * std::sqrt(1. * res.first[f]));
1668  * }
1669  *
1670  * @endcode
1671  *
1672  * We have to calculate the logarithms of absolute values of
1673  * coefficients and use it in a linear regression fit to obtain @f$\mu@f$.
1674  *
1675  * @code
1676  * for (double &residual_element : res.second)
1677  * residual_element = std::log(residual_element);
1678  *
1679  * std::pair<double, double> fit =
1680  * FESeries::linear_regression(ln_k, res.second);
1681  *
1682  * @endcode
1683  *
1684  * The final step is to compute the Sobolev index @f$s=\mu-\frac d2@f$ and
1685  * store it in the vector of estimated values for each cell:
1686  *
1687  * @code
1688  * smoothness_indicators(cell->active_cell_index()) =
1689  * -fit.first - 1. * dim / 2;
1690  * }
1691  * }
1692  * } // namespace Step27
1693  *
1694  *
1695  * @endcode
1696  *
1697  *
1698  * <a name="Themainfunction"></a>
1699  * <h3>The main function</h3>
1700  *
1701 
1702  *
1703  * The main function is again verbatim what we had before: wrap creating and
1704  * running an object of the main class into a <code>try</code> block and catch
1705  * whatever exceptions are thrown, thereby producing meaningful output if
1706  * anything should go wrong:
1707  *
1708  * @code
1709  * int main()
1710  * {
1711  * try
1712  * {
1713  * using namespace Step27;
1714  *
1715  * LaplaceProblem<2> laplace_problem;
1716  * laplace_problem.run();
1717  * }
1718  * catch (std::exception &exc)
1719  * {
1720  * std::cerr << std::endl
1721  * << std::endl
1722  * << "----------------------------------------------------"
1723  * << std::endl;
1724  * std::cerr << "Exception on processing: " << std::endl
1725  * << exc.what() << std::endl
1726  * << "Aborting!" << std::endl
1727  * << "----------------------------------------------------"
1728  * << std::endl;
1729  *
1730  * return 1;
1731  * }
1732  * catch (...)
1733  * {
1734  * std::cerr << std::endl
1735  * << std::endl
1736  * << "----------------------------------------------------"
1737  * << std::endl;
1738  * std::cerr << "Unknown exception!" << std::endl
1739  * << "Aborting!" << std::endl
1740  * << "----------------------------------------------------"
1741  * << std::endl;
1742  * return 1;
1743  * }
1744  *
1745  * return 0;
1746  * }
1747  * @endcode
1748 <a name="Results"></a><h1>Results</h1>
1749 
1750 
1751 In this section, we discuss a few results produced from running the
1752 current tutorial program. More results, in particular the extension to
1753 3d calculations and determining how much compute time the individual
1754 components of the program take, are given in the @ref hp_paper .
1755 
1756 When run, this is what the program produces:
1757 
1758 @code
1759 > make run
1760 [ 66%] Built target @ref step_27 "step-27"
1761 [100%] Run @ref step_27 "step-27" with Release configuration
1762 Cycle 0:
1763  Number of active cells : 768
1764  Number of degrees of freedom: 3264
1765  Number of constraints : 384
1766 Cycle 1:
1767  Number of active cells : 966
1768  Number of degrees of freedom: 5245
1769  Number of constraints : 936
1770 Cycle 2:
1771  Number of active cells : 1143
1772  Number of degrees of freedom: 8441
1773  Number of constraints : 1929
1774 Cycle 3:
1775  Number of active cells : 1356
1776  Number of degrees of freedom: 12349
1777  Number of constraints : 3046
1778 Cycle 4:
1779  Number of active cells : 1644
1780  Number of degrees of freedom: 18178
1781  Number of constraints : 4713
1782 Cycle 5:
1783  Number of active cells : 1728
1784  Number of degrees of freedom: 22591
1785  Number of constraints : 6095
1786 @endcode
1787 
1788 The first thing we learn from this is that the number of constrained degrees
1789 of freedom is on the order of 20-25% of the total number of degrees of
1790 freedom, at least on the later grids when we have elements of relatively
1791 high order (in 3d, the fraction of constrained degrees of freedom can be up
1792 to 30%). This is, in fact, on the same order of magnitude as for non-@f$hp@f$
1793 discretizations. For example, in the last step of the @ref step_6 "step-6"
1794 program, we have 18353 degrees of freedom, 4432 of which are
1795 constrained. The difference is that in the latter program, each constrained
1796 hanging node is constrained against only the two adjacent degrees of
1797 freedom, whereas in the @f$hp@f$ case, constrained nodes are constrained against
1798 many more degrees of freedom. Note also that the current program also
1799 includes nodes subject to Dirichlet boundary conditions in the list of
1800 constraints. In cycle 0, all the constraints are actually because of
1801 boundary conditions.
1802 
1803 Of maybe more interest is to look at the graphical output. First, here is the
1804 solution of the problem:
1805 
1806 <img src="https://www.dealii.org/images/steps/developer/step-27-solution.png"
1807  alt="Elevation plot of the solution, showing the lack of regularity near
1808  the interior (reentrant) corners."
1809  width="200" height="200">
1810 
1811 Secondly, let us look at the sequence of meshes generated:
1812 
1813 <div class="threecolumn" style="width: 80%">
1814  <div>
1815  <img src="https://www.dealii.org/images/steps/developer/step-27-mesh-0.svg"
1816  alt="Triangulation containing reentrant corners without adaptive refinement."
1817  width="200" height="200">
1818  </div>
1819  <div>
1820  <img src="https://www.dealii.org/images/steps/developer/step-27-mesh-1.svg"
1821  alt="Triangulation containing reentrant corners with one level of
1822  refinement. New cells are placed near the corners."
1823  width="200" height="200">
1824  </div>
1825  <div>
1826  <img src="https://www.dealii.org/images/steps/developer/step-27-mesh-2.svg"
1827  alt="Triangulation containing reentrant corners with two levels of
1828  refinement. New cells are placed near the corners."
1829  width="200" height="200">
1830  </div>
1831  <div>
1832  <img src="https://www.dealii.org/images/steps/developer/step-27-mesh-3.svg"
1833  alt="Triangulation containing reentrant corners with three levels of
1834  refinement. New cells are placed near the corners."
1835  width="200" height="200">
1836  </div>
1837  <div>
1838  <img src="https://www.dealii.org/images/steps/developer/step-27-mesh-4.svg"
1839  alt="Triangulation containing reentrant corners with four levels of
1840  refinement. New cells are placed near the corners."
1841  width="200" height="200">
1842  </div>
1843  <div>
1844  <img src="https://www.dealii.org/images/steps/developer/step-27-mesh-5.svg"
1845  alt="Triangulation containing reentrant corners with five levels of
1846  refinement. New cells are placed near the corners."
1847  width="200" height="200">
1848  </div>
1849 </div>
1850 
1851 It is clearly visible how the mesh is refined near the corner singularities,
1852 as one would expect it. More interestingly, we should be curious to see the
1853 distribution of finite element polynomial degrees to these mesh cells, where
1854 grey corresponds to degree two and pink corresponds to degree seven:
1855 
1856 <div class="threecolumn" style="width: 80%">
1857  <div>
1858  <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-0.png"
1859  alt="Initial grid where all cells contain just biquadratic functions."
1860  width="200" height="200">
1861  </div>
1862  <div>
1863  <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-1.png"
1864  alt="Depiction of local approximation degrees after one refinement."
1865  width="200" height="200">
1866  </div>
1867  <div>
1868  <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-2.png"
1869  alt="Depiction of local approximation degrees after two refinements."
1870  width="200" height="200">
1871  </div>
1872  <div>
1873  <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-3.png"
1874  alt="Depiction of local approximation degrees after three refinements."
1875  width="200" height="200">
1876  </div>
1877  <div>
1878  <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-4.png"
1879  alt="Depiction of local approximation degrees after four refinements."
1880  width="200" height="200">
1881  </div>
1882  <div>
1883  <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-5.png"
1884  alt="Depiction of local approximation degrees after five refinements."
1885  width="200" height="200">
1886  </div>
1887 </div>
1888 
1889 While this is certainly not a perfect arrangement, it does make some sense: we
1890 use low order elements close to boundaries and corners where regularity is
1891 low. On the other hand, higher order elements are used where (i) the error was
1892 at one point fairly large, i.e. mainly in the general area around the corner
1893 singularities and in the top right corner where the solution is large, and
1894 (ii) where the solution is smooth, i.e. far away from the boundary.
1895 
1896 This arrangement of polynomial degrees of course follows from our smoothness
1897 estimator. Here is the estimated smoothness of the solution, with darker colors
1898 indicating least smoothness and lighter indicating the smoothest areas:
1899 
1900 <div class="threecolumn" style="width: 80%">
1901  <div>
1902  <img src="https://www.dealii.org/images/steps/developer/step-27-smoothness-0.png"
1903  alt="Estimated regularity per cell on the initial grid."
1904  width="200" height="200">
1905  </div>
1906  <div>
1907  <img src="https://www.dealii.org/images/steps/developer/step-27-smoothness-1.png"
1908  alt="Depiction of the estimated regularity per cell after one refinement."
1909  width="200" height="200">
1910  </div>
1911  <div>
1912  <img src="https://www.dealii.org/images/steps/developer/step-27-smoothness-2.png"
1913  alt="Depiction of the estimated regularity per cell after two refinements."
1914  width="200" height="200">
1915  </div>
1916  <div>
1917  <img src="https://www.dealii.org/images/steps/developer/step-27-smoothness-3.png"
1918  alt="Depiction of the estimated regularity per cell after three refinements."
1919  width="200" height="200">
1920  </div>
1921  <div>
1922  <img src="https://www.dealii.org/images/steps/developer/step-27-smoothness-4.png"
1923  alt="Depiction of the estimated regularity per cell after four refinements."
1924  width="200" height="200">
1925  </div>
1926  <div>
1927  <img src="https://www.dealii.org/images/steps/developer/step-27-smoothness-5.png"
1928  alt="Depiction of the estimated regularity per cell after five refinements."
1929  width="200" height="200">
1930  </div>
1931 </div>
1932 
1933 The primary conclusion one can draw from this is that the loss of regularity at
1934 the internal corners is a highly localized phenomenon; it only seems to impact
1935 the cells adjacent to the corner itself, so when we refine the mesh the black
1936 coloring is no longer visible. Besides the corners, this sequence of plots
1937 implies that the smoothness estimates are somewhat independent of the mesh
1938 refinement, particularly when we are far away from boundaries.
1939 It is also obvious that the smoothness estimates are independent of the actual
1940 size of the solution (see the picture of the solution above), as it should be.
1941 A point of larger concern, however, is that one realizes on closer inspection
1942 that the estimator we have overestimates the smoothness of the solution on
1943 cells with hanging nodes. This in turn leads to higher polynomial degrees in
1944 these areas, skewing the allocation of finite elements onto cells.
1945 
1946 We have no good explanation for this effect at the moment. One theory is that
1947 the numerical solution on cells with hanging nodes is, of course, constrained
1948 and therefore not entirely free to explore the function space to get close to
1949 the exact solution. This lack of degrees of freedom may manifest itself by
1950 yielding numerical solutions on these cells with suppressed oscillation,
1951 meaning a higher degree of smoothness. The estimator picks this signal up and
1952 the estimated smoothness overestimates the actual value. However, a definite
1953 answer to what is going on currently eludes the authors of this program.
1954 
1955 The bigger question is, of course, how to avoid this problem. Possibilities
1956 include estimating the smoothness not on single cells, but cell assemblies or
1957 patches surrounding each cell. It may also be possible to find simple
1958 correction factors for each cell depending on the number of constrained
1959 degrees of freedom it has. In either case, there are ample opportunities for
1960 further research on finding good @f$hp@f$ refinement criteria. On the other hand,
1961 the main point of the current program was to demonstrate using the @f$hp@f$
1962 technology in deal.II, which is unaffected by our use of a possible
1963 sub-optimal refinement criterion.
1964  *
1965  *
1966 <a name="PlainProg"></a>
1967 <h1> The plain program</h1>
1968 @include "step-27.cc"
1969 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
TableIndices
Definition: table_indices.h:45
FE_Q
Definition: fe_q.h:554
CellData
Definition: tria_description.h:67
parallel::transform
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:211
hp::QCollection
Definition: q_collection.h:48
TableBase::reinit
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
hp::DoFHandler
Definition: dof_handler.h:203
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
FESeries::Fourier
Definition: fe_series.h:91
Table
Definition: table.h:699
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
Utilities::CUDA::free
void free(T *&pointer)
Definition: cuda.h:96
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
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
LAPACKSupport::T
static const char T
Definition: lapack_support.h:163
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
FESeries::Fourier::calculate
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
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
QGauss
Definition: quadrature_lib.h:40
value
static const bool value
Definition: dof_tools_constraints.cc:433
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
FESeries::process_coefficients
std::pair< std::vector< unsigned int >, std::vector< double > > process_coefficients(const Table< dim, CoefficientType > &coefficients, const std::function< std::pair< bool, unsigned int >(const TableIndices< dim > &)> &predicate, const VectorTools::NormType norm_type, const double smallest_abs_coefficient=1e-10)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
VectorizedArray::sqrt
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
FESeries::linear_regression
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
hp
Definition: hp.h:117
VectorTools::Linfty_norm
@ Linfty_norm
Definition: vector_tools_common.h:148
std::log
inline ::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5390
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
SubCellData
Definition: tria_description.h:199
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12618
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
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector
Definition: mapping_q1_eulerian.h:32
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
hp::FECollection
Definition: fe_collection.h:54