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-27.h
Go to the documentation of this file.
1) const
941 *   {
942 *   double product = 1;
943 *   for (unsigned int d = 0; d < dim; ++d)
944 *   product *= (p[d] + 1);
945 *   return product;
946 *   }
947 *  
948 *  
949 *  
950 * @endcode
951 *
952 *
953 * <a name="step_27-Implementationofthemainclass"></a>
954 * <h3>Implementation of the main class</h3>
955 *
956
957 *
958 *
959 * <a name="step_27-LaplaceProblemLaplaceProblemconstructor"></a>
961 *
962
963 *
966 * maximal polynomial degree to 7 (in 1d and 2d) or 5 (in 3d and higher). We
967 * do so because using higher order polynomial degrees becomes prohibitively
968 * expensive, especially in higher space dimensions.
969 *
970
971 *
972 * Following this, we fill the collections of finite element, and cell and
973 * face quadrature objects. We start with quadratic elements, and each
974 * quadrature formula is chosen so that it is appropriate for the matching
975 * finite element in the hp::FECollection object.
976 *
977 * @code
978 *   template <int dim>
980 *   : dof_handler(triangulation)
981 *   , max_degree(dim <= 2 ? 7 : 5)
982 *   {
983 *   for (unsigned int degree = 2; degree <= max_degree; ++degree)
984 *   {
985 *   fe_collection.push_back(FE_Q<dim>(degree));
986 *   quadrature_collection.push_back(QGauss<dim>(degree + 1));
987 *   face_quadrature_collection.push_back(QGauss<dim - 1>(degree + 1));
988 *   }
989 *   }
990 *  
991 *  
992 * @endcode
993 *
994 *
995 * <a name="step_27-LaplaceProblemLaplaceProblemdestructor"></a>
997 *
998
999 *
1000 * The destructor is unchanged from what we already did in @ref step_6 "step-6":
1001 *
1002 * @code
1003 *   template <int dim>
1005 *   {
1006 *   dof_handler.clear();
1007 *   }
1008 *  
1009 *  
1010 * @endcode
1011 *
1012 *
1013 * <a name="step_27-LaplaceProblemsetup_system"></a>
1015 *
1016
1017 *
1018 * This function is again a verbatim copy of what we already did in
1019 * @ref step_6 "step-6". Despite function calls with exactly the same names and arguments,
1021 * dof_handler variable here is in @f$hp@f$-mode.
1022 *
1023 * @code
1024 *   template <int dim>
1026 *   {
1027 *   dof_handler.distribute_dofs(fe_collection);
1028 *  
1029 *   solution.reinit(dof_handler.n_dofs());
1030 *   system_rhs.reinit(dof_handler.n_dofs());
1031 *  
1032 *   constraints.clear();
1033 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1035 *   0,
1037 *   constraints);
1038 *   constraints.close();
1039 *  
1040 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1041 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
1042 *   sparsity_pattern.copy_from(dsp);
1043 *  
1044 *   system_matrix.reinit(sparsity_pattern);
1045 *   }
1046 *  
1047 *  
1048 *  
1049 * @endcode
1050 *
1051 *
1052 * <a name="step_27-LaplaceProblemassemble_system"></a>
1054 *
1055
1056 *
1057 * This is the function that assembles the global matrix and right hand side
1058 * vector from the local contributions of each cell. Its main working is as
1064 * local contributions into global objects. Both of these are explained in
1065 * detail in the introduction of this program.
1066 *
1067
1068 *
1070 * polynomial degrees on different cells, the matrices and vectors holding
1071 * local contributions do not have the same size on all cells. At the
1072 * beginning of the loop over all cells, we therefore each time have to
1075 * of a matrix or vector does not release the currently allocated memory
1076 * (unless the new size is zero), the process of resizing at the beginning of
1078 * iterations. Once we have found in a cell with the maximal finite element
1082 * expensive, and doing so every time we visit a new cell would take
1083 * significant compute time.
1084 *
1085 * @code
1088 *   {
1089 *   hp::FEValues<dim> hp_fe_values(fe_collection,
1094 *  
1096 *  
1099 *  
1100 *   std::vector<types::global_dof_index> local_dof_indices;
1101 *  
1102 *   for (const auto &cell : dof_handler.active_cell_iterators())
1103 *   {
1104 *   const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1105 *  
1106 *   cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1107 *   cell_matrix = 0;
1108 *  
1109 *   cell_rhs.reinit(dofs_per_cell);
1110 *   cell_rhs = 0;
1111 *  
1112 *   hp_fe_values.reinit(cell);
1113 *  
1114 *   const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
1115 *  
1116 *   std::vector<double> rhs_values(fe_values.n_quadrature_points);
1117 *   rhs_function.value_list(fe_values.get_quadrature_points(), rhs_values);
1118 *  
1119 *   for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
1120 *   ++q_point)
1121 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1122 *   {
1123 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1124 *   cell_matrix(i, j) +=
1125 *   (fe_values.shape_grad(i, q_point) * // grad phi_i(x_q)
1126 *   fe_values.shape_grad(j, q_point) * // grad phi_j(x_q)
1127 *   fe_values.JxW(q_point)); // dx
1128 *  
1129 *   cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
1130 *   rhs_values[q_point] * // f(x_q)
1131 *   fe_values.JxW(q_point)); // dx
1132 *   }
1133 *  
1134 *   local_dof_indices.resize(dofs_per_cell);
1135 *   cell->get_dof_indices(local_dof_indices);
1136 *  
1137 *   constraints.distribute_local_to_global(
1138 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1139 *   }
1140 *   }
1141 *  
1142 *  
1143 *  
1144 * @endcode
1145 *
1146 *
1147 * <a name="step_27-LaplaceProblemsolve"></a>
1148 * <h4>LaplaceProblem::solve</h4>
1149 *
1150
1151 *
1152 * The function solving the linear system is entirely unchanged from
1154 * equals the @f$l_2@f$ norm of the right hand side) by a certain factor:
1155 *
1156 * @code
1157 *   template <int dim>
1159 *   {
1160 *   SolverControl solver_control(system_rhs.size(),
1161 *   1e-12 * system_rhs.l2_norm());
1162 *   SolverCG<Vector<double>> cg(solver_control);
1163 *  
1164 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1165 *   preconditioner.initialize(system_matrix, 1.2);
1166 *  
1167 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
1168 *  
1169 *   constraints.distribute(solution);
1170 *   }
1171 *  
1172 *  
1173 *  
1174 * @endcode
1175 *
1176 *
1177 * <a name="step_27-LaplaceProblempostprocess"></a>
1178 * <h4>LaplaceProblem::postprocess</h4>
1179 *
1180
1181 *
1182 * After solving the linear system, we will want to postprocess the
1183 * solution. Here, all we do is to estimate the error, estimate the local
1184 * smoothness of the solution as described in the introduction, then write
1185 * graphical output, and finally refine the mesh in both @f$h@f$ and @f$p@f$
1186 * according to the indicators computed before. We do all this in the same
1188 * not only for refinement, but also include them in the graphical output.
1189 *
1190 * @code
1191 *   template <int dim>
1192 *   void LaplaceProblem<dim>::postprocess(const unsigned int cycle)
1193 *   {
1194 * @endcode
1195 *
1197 * which each are one number for each active cell of our
1199 * class as always.
1200 *
1201 * @code
1204 *   dof_handler,
1205 *   face_quadrature_collection,
1206 *   std::map<types::boundary_id, const Function<dim> *>(),
1207 *   solution,
1209 *  
1210 * @endcode
1211 *
1213 * expansion coefficients as outlined in the introduction. We will first
1214 * need to create an object capable of transforming the finite element
1215 * solution on every single cell into a sequence of Fourier series
1216 * coefficients. The SmoothnessEstimator namespace offers a factory function
1217 * for such a FESeries::Fourier object that is optimized for the process of
1219 * coefficients on every individual cell then happens in the last function.
1220 *
1221 * @code
1226 *   dof_handler,
1227 *   solution,
1229 *  
1230 * @endcode
1231 *
1232 * Next we want to generate graphical output. In addition to the two
1234 * polynomial degree of the finite elements used on each of the elements
1235 * on the mesh.
1236 *
1237
1238 *
1239 * The way to do that requires that we loop over all cells and poll the
1240 * active finite element index of them using
1241 * <code>cell-@>active_fe_index()</code>. We then use the result of this
1242 * operation and query the finite element collection for the finite
1243 * element with that index, and finally determine the polynomial degree of
1244 * that element. The result we put into a vector with one element per
1245 * cell. The DataOut class requires this to be a vector of
1246 * <code>float</code> or <code>double</code>, even though our values are
1247 * all integers, so that is what we use:
1248 *
1249 * @code
1250 *   {
1251 *   Vector<float> fe_degrees(triangulation.n_active_cells());
1252 *   for (const auto &cell : dof_handler.active_cell_iterators())
1253 *   fe_degrees(cell->active_cell_index()) =
1254 *   fe_collection[cell->active_fe_index()].degree;
1255 *  
1256 * @endcode
1257 *
1258 * With now all data vectors available -- solution, estimated errors and
1259 * smoothness indicators, and finite element degrees --, we create a
1260 * DataOut object for graphical output and attach all data:
1261 *
1262 * @code
1263 *   DataOut<dim> data_out;
1264 *  
1265 *   data_out.attach_dof_handler(dof_handler);
1266 *   data_out.add_data_vector(solution, "solution");
1267 *   data_out.add_data_vector(estimated_error_per_cell, "error");
1268 *   data_out.add_data_vector(smoothness_indicators, "smoothness");
1269 *   data_out.add_data_vector(fe_degrees, "fe_degree");
1270 *   data_out.build_patches();
1271 *  
1272 * @endcode
1273 *
1274 * The final step in generating output is to determine a file name, open
1275 * the file, and write the data into it (here, we use VTK format):
1276 *
1277 * @code
1278 *   const std::string filename =
1279 *   "solution-" + Utilities::int_to_string(cycle, 2) + ".vtk";
1280 *   std::ofstream output(filename);
1281 *   data_out.write_vtk(output);
1282 *   }
1283 *  
1284 * @endcode
1285 *
1286 * After this, we would like to actually refine the mesh, in both @f$h@f$ and
1287 * @f$p@f$. The way we are going to do this is as follows: first, we use the
1288 * estimated error to flag those cells for refinement that have the
1290 *
1291 * @code
1292 *   {
1295 *   0.3,
1296 *   0.03);
1297 *  
1298 * @endcode
1299 *
1302 * @f$h@f$ decreased. The strategy we choose here is that we look at the
1305 * relative threshold. In other words, for every cell for which (i) the
1307 * the threshold, and (iii) we still have a finite element with a
1308 * polynomial degree higher than the current one in the finite element
1309 * collection, we will assign a future FE index that corresponds to a
1310 * polynomial with degree one higher than it currently is. The following
1311 * function is capable of doing exactly this. Absent any better
1312 * strategies, we will set the threshold via interpolation between the
1316 * with a low threshold by setting a small interpolation factor of 0.2. In
1318 * decrease their polynomial degree when their smoothness indicator is
1319 * below the corresponding threshold determined on cells to be coarsened.
1320 *
1321 * @code
1323 *   dof_handler, smoothness_indicators, 0.2, 0.2);
1324 *  
1325 * @endcode
1326 *
1327 * The above function only determines whether the polynomial degree will
1328 * change via future FE indices, but does not manipulate the
1329 * @f$h@f$-refinement flags. So for cells that are flagged for both refinement
1332 * not both at once.
1333 *
1334 * @code
1335 *   hp::Refinement::choose_p_over_h(dof_handler);
1336 *  
1337 * @endcode
1338 *
1339 * For grid adaptive refinement, we ensure a 2:1 mesh balance by limiting
1340 * the difference of refinement levels of neighboring cells to one by
1343 * cells: levels of future finite elements are not allowed to differ by
1344 * more than a specified difference. With its default parameters, a call
1346 * difference is limited to one. This will not necessarily decrease the
1347 * number of hanging nodes in the domain, but makes sure that high order
1348 * polynomials are not constrained to much lower polynomials on faces,
1349 * e.g., fifth order to second order polynomials.
1350 *
1351 * @code
1352 *   triangulation.prepare_coarsening_and_refinement();
1353 *   hp::Refinement::limit_p_level_difference(dof_handler);
1354 *  
1355 * @endcode
1356 *
1357 * At the end of this procedure, we then refine the mesh. During this
1358 * process, children of cells undergoing bisection inherit their mother
1359 * cell's finite element index. Further, future finite element indices
1361 * assigned to cells after the next call of DoFHandler::distribute_dofs().
1362 *
1363 * @code
1364 *   triangulation.execute_coarsening_and_refinement();
1365 *   }
1366 *   }
1367 *  
1368 *  
1369 * @endcode
1370 *
1371 *
1372 * <a name="step_27-LaplaceProblemcreate_coarse_grid"></a>
1374 *
1375
1376 *
1377 * The following function is used when creating the initial grid. The grid we
1378 * would like to create is actually similar to the one from @ref step_14 "step-14", i.e., the
1379 * square domain with the square hole in the middle. It can be generated by
1383 *
1384
1385 *
1387 * already holds our desired domain @f$[-1,1]^d@f$, subdivided into @f$4^d@f$ cells.
1388 * We then remove those cells in the center of the domain by testing the
1389 * coordinate values of the vertices on each cell. In the end, we refine the
1390 * so created grid globally as usual.
1391 *
1392 * @code
1393 *   template <int dim>
1394 *   void LaplaceProblem<dim>::create_coarse_grid()
1395 *   {
1398 *  
1399 *   std::set<typename Triangulation<dim>::active_cell_iterator> cells_to_remove;
1400 *   for (const auto &cell : cube.active_cell_iterators())
1401 *   for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1402 *   if (cell->vertex(v).square() < .1)
1403 *   cells_to_remove.insert(cell);
1404 *  
1407 *   triangulation);
1408 *  
1409 *   triangulation.refine_global(3);
1410 *   }
1411 *  
1412 *  
1413 *  
1414 * @endcode
1415 *
1416 *
1417 * <a name="step_27-LaplaceProblemrun"></a>
1418 * <h4>LaplaceProblem::run</h4>
1419 *
1420
1421 *
1423 * function in most of the previous programs already, see for example @ref step_6 "step-6".
1424 *
1425
1426 *
1427 * Basically, it contains the adaptive loop: in the first iteration create a
1428 * coarse grid, and then set up the linear system, assemble it, solve, and
1429 * postprocess the solution including mesh refinement. Then start over
1431 * the screen trying to figure out what the program does:
1432 *
1433 * @code
1434 *   template <int dim>
1436 *   {
1437 *   for (unsigned int cycle = 0; cycle < 6; ++cycle)
1438 *   {
1439 *   std::cout << "Cycle " << cycle << ':' << std::endl;
1440 *  
1441 *   if (cycle == 0)
1443 *  
1444 *   setup_system();
1445 *  
1446 *   std::cout << " Number of active cells : "
1447 *   << triangulation.n_active_cells() << std::endl
1448 *   << " Number of degrees of freedom: " << dof_handler.n_dofs()
1449 *   << std::endl
1450 *   << " Number of constraints : "
1451 *   << constraints.n_constraints() << std::endl;
1452 *  
1453 *   assemble_system();
1454 *   solve();
1455 *   postprocess(cycle);
1456 *   }
1457 *   }
1458 *   } // namespace Step27
1459 *  
1460 *  
1461 * @endcode
1462 *
1463 *
1464 * <a name="step_27-Themainfunction"></a>
1465 * <h3>The main function</h3>
1466 *
1467
1468 *
1469 * The main function is again verbatim what we had before: wrap creating and
1470 * running an object of the main class into a <code>try</code> block and catch
1473 *
1474 * @code
1475 *   int main()
1476 *   {
1477 *   try
1478 *   {
1479 *   using namespace Step27;
1480 *  
1482 *   laplace_problem.run();
1483 *   }
1484 *   catch (std::exception &exc)
1485 *   {
1486 *   std::cerr << std::endl
1487 *   << std::endl
1488 *   << "----------------------------------------------------"
1489 *   << std::endl;
1490 *   std::cerr << "Exception on processing: " << std::endl
1491 *   << exc.what() << std::endl
1492 *   << "Aborting!" << std::endl
1493 *   << "----------------------------------------------------"
1494 *   << std::endl;
1495 *  
1496 *   return 1;
1497 *   }
1498 *   catch (...)
1499 *   {
1500 *   std::cerr << std::endl
1501 *   << std::endl
1502 *   << "----------------------------------------------------"
1503 *   << std::endl;
1504 *   std::cerr << "Unknown exception!" << std::endl
1505 *   << "Aborting!" << std::endl
1506 *   << "----------------------------------------------------"
1507 *   << std::endl;
1508 *   return 1;
1509 *   }
1510 *  
1511 *   return 0;
1512 *   }
1513 * @endcode
1514<a name="step_27-Results"></a><h1>Results</h1>
1515
1516
1517In this section, we discuss a few results produced from running the
1520components of the program take, are given in the @ref hp_paper "hp-paper".
1521
1522When run, this is what the program produces:
1523
1524@code
1525> make run
1526[ 66%] Built target @ref step_27 "step-27"
1527[100%] Run @ref step_27 "step-27" with Release configuration
1528Cycle 0:
1529 Number of active cells : 768
1530 Number of degrees of freedom: 3264
1531 Number of constraints : 384
1532Cycle 1:
1533 Number of active cells : 807
1534 Number of degrees of freedom: 4764
1535 Number of constraints : 756
1536Cycle 2:
1537 Number of active cells : 927
1538 Number of degrees of freedom: 8226
1539 Number of constraints : 1856
1540Cycle 3:
1541 Number of active cells : 978
1542 Number of degrees of freedom: 12146
1543 Number of constraints : 2944
1544Cycle 4:
1545 Number of active cells : 1104
1546 Number of degrees of freedom: 16892
1547 Number of constraints : 3998
1548Cycle 5:
1549 Number of active cells : 1149
1550 Number of degrees of freedom: 22078
1551 Number of constraints : 5230
1552@endcode
1553
1554The first thing we learn from this is that the number of constrained degrees
1555of freedom is on the order of 20-25% of the total number of degrees of
1557high order (in 3d, the fraction of constrained degrees of freedom can be up
1558to 30%). This is, in fact, on the same order of magnitude as for
1559non-@f$hp@f$-discretizations. For example, in the last step of the @ref step_6 "step-6"
1560program, we have 18353 degrees of freedom, 4432 of which are
1566constraints. In cycle 0, all the constraints are actually because of
1567boundary conditions.
1568
1570solution of the problem:
1571
1572<img src="https://www.dealii.org/images/steps/developer/step-27-solution.png"
1573 alt="Elevation plot of the solution, showing the lack of regularity near
1574 the interior (reentrant) corners."
1575 width="200" height="200">
1576
1577Secondly, let us look at the sequence of meshes generated:
1578
1579<div class="threecolumn" style="width: 80%">
1580 <div>
1581 <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-00.svg"
1582 alt="Triangulation containing reentrant corners without adaptive refinement."
1583 width="200" height="200">
1584 </div>
1585 <div>
1586 <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-01.svg"
1587 alt="Triangulation containing reentrant corners with one level of
1588 refinement. New cells are placed near the corners."
1589 width="200" height="200">
1590 </div>
1591 <div>
1592 <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-02.svg"
1593 alt="Triangulation containing reentrant corners with two levels of
1594 refinement. New cells are placed near the corners."
1595 width="200" height="200">
1596 </div>
1597 <div>
1598 <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-03.svg"
1599 alt="Triangulation containing reentrant corners with three levels of
1600 refinement. New cells are placed near the corners."
1601 width="200" height="200">
1602 </div>
1603 <div>
1604 <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-04.svg"
1605 alt="Triangulation containing reentrant corners with four levels of
1606 refinement. New cells are placed near the corners."
1607 width="200" height="200">
1608 </div>
1609 <div>
1610 <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-05.svg"
1611 alt="Triangulation containing reentrant corners with five levels of
1612 refinement. New cells are placed near the corners."
1613 width="200" height="200">
1614 </div>
1615</div>
1616
1619distribution of finite element polynomial degrees to these mesh cells, where
1621to degree seven:
1622
1623<div class="threecolumn" style="width: 80%">
1624 <div>
1625 <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-00.svg"
1626 alt="Initial grid where all cells contain just biquadratic functions."
1627 width="200" height="200">
1628 </div>
1629 <div>
1630 <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-01.svg"
1631 alt="Depiction of local approximation degrees after one refinement."
1632 width="200" height="200">
1633 </div>
1634 <div>
1635 <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-02.svg"
1636 alt="Depiction of local approximation degrees after two refinements."
1637 width="200" height="200">
1638 </div>
1639 <div>
1640 <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-03.svg"
1641 alt="Depiction of local approximation degrees after three refinements."
1642 width="200" height="200">
1643 </div>
1644 <div>
1645 <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-04.svg"
1646 alt="Depiction of local approximation degrees after four refinements."
1647 width="200" height="200">
1648 </div>
1649 <div>
1650 <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-05.svg"
1651 alt="Depiction of local approximation degrees after five refinements."
1652 width="200" height="200">
1653 </div>
1654</div>
1655
1657use low order elements close to boundaries and corners where regularity is
1658low. On the other hand, higher order elements are used where (i) the error was
1660singularities and in the top right corner where the solution is large, and
1661(ii) where the solution is smooth, i.e., far away from the boundary.
1662
1663This arrangement of polynomial degrees of course follows from our smoothness
1666
1667<div class="threecolumn" style="width: 80%">
1668 <div>
1669 <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-00.svg"
1670 alt="Estimated regularity per cell on the initial grid."
1671 width="200" height="200">
1672 </div>
1673 <div>
1674 <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-01.svg"
1675 alt="Depiction of the estimated regularity per cell after one refinement."
1676 width="200" height="200">
1677 </div>
1678 <div>
1679 <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-02.svg"
1680 alt="Depiction of the estimated regularity per cell after two refinements."
1681 width="200" height="200">
1682 </div>
1683 <div>
1684 <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-03.svg"
1685 alt="Depiction of the estimated regularity per cell after three refinements."
1686 width="200" height="200">
1687 </div>
1688 <div>
1689 <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-04.svg"
1690 alt="Depiction of the estimated regularity per cell after four refinements."
1691 width="200" height="200">
1692 </div>
1693 <div>
1694 <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-05.svg"
1695 alt="Depiction of the estimated regularity per cell after five refinements."
1696 width="200" height="200">
1697 </div>
1698</div>
1699
1700The primary conclusion one can draw from this is that the loss of regularity at
1703coloring is no longer visible. Besides the corners, this sequence of plots
1707size of the solution (see the picture of the solution above), as it should be.
1710cells with hanging nodes. This in turn leads to higher polynomial degrees in
1711these areas, skewing the allocation of finite elements onto cells.
1712
1715and therefore not entirely free to explore the function space to get close to
1716the exact solution. This lack of degrees of freedom may manifest itself by
1721
1722The bigger question is, of course, how to avoid this problem. Possibilities
1724patches surrounding each cell. It may also be possible to find simple
1725correction factors for each cell depending on the number of constrained
1726degrees of freedom it has. In either case, there are ample opportunities for
1730sub-optimal refinement criterion.
1731
1732
1733
1734<a name="step-27-extensions"></a>
1735<a name="step_27-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1736
1737
1738<a name="step_27-Differenthpdecisionstrategies"></a><h4>Different hp-decision strategies</h4>
1739
1740
1744<ul>
1745 <li><i>Fourier coefficient decay:</i> This is the strategy currently
1746 implemented in this tutorial. For more information on this strategy, see
1748
1749 <li><i>Legendre coefficient decay:</i> This strategy is quite similar
1751 Fourier one: instead of sinusoids as basis functions, this strategy uses
1752 Legendre polynomials. Of course, since we approximate the solution using a
1753 finite-dimensional polynomial on each cell, the expansion of the solution in
1754 Legendre polynomials is also finite and, consequently, when we talk about the
1756 coefficients of this expansion, rather than think about it in asymptotic terms.
1757 But, if we have enough of these coefficients, we can certainly think of the
1758 decay of these coefficients as characteristic of the decay of the coefficients
1759 of the exact solution (which is, in general, not polynomial and so will have an
1760 infinite Legendre expansion), and considering the coefficients we have should
1761 reveal something about the properties of the exact solution.
1762
1763 The transition from the Fourier strategy to the Legendre one is quite simple:
1766 FESeries::Legendre and SmoothnessEstimator::Legendre. This strategy is used
1767 in @ref step_75 "step-75". For the theoretical background of this strategy, consult the
1770
1771 <li><i>Refinement history:</i> The last strategy is quite different
1775 out in @ref step_7 "step-7". If the solution is sufficiently smooth, though, we
1777 polynomial degree of the finite element. We can compare a proper
1780
1781 The transition to this strategy is a bit more complicated. For this, we need
1788
1790 the next time step in time-dependent problems. Therefore, this strategy
1793 strategies as well: start each time step with a coarse mesh, keep refining
1794 until happy with the result, and only then move on to the next time step.</li>
1795</ul>
1796
1802in @cite fehling2020 .
1803
1804
1805<a name="step_27-Parallelhpadaptivefiniteelements"></a><h4>Parallel hp-adaptive finite elements</h4>
1806
1807
1812it, we recommend reading @ref step_18 "step-18" for the former and @ref step_40 "step-40" for the
1815
1816We go one step further in @ref step_75 "step-75": Here, we combine hp-adaptive and
1818parallel::distributed::Triangulation objects.
1819 *
1820 *
1821<a name="step_27-PlainProg"></a>
1822<h1> The plain program</h1>
1823@include "step-27.cc"
1824*/
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
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
virtual bool prepare_coarsening_and_refinement()
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
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
const Event initial
Definition event.cc:68
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr types::blas_int zero
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 > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
FESeries::Fourier< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
void coefficient_decay(FESeries::Fourier< dim, spacedim > &fe_fourier, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const VectorTools::NormType regression_strategy=VectorTools::Linfty_norm, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
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)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474
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)
bool limit_p_level_difference(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int max_difference=1, const unsigned int contains_fe_index=0)
void predict_error(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &error_indicators, Vector< Number > &predicted_errors, const double gamma_p=std::sqrt(0.4), const double gamma_h=2., const double gamma_n=1.)
void choose_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_from_relative_threshold(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const double p_refine_fraction=0.5, const double p_coarsen_fraction=0.5, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
Definition hp.h:117
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.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation