Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
step-9.h
Go to the documentation of this file.
1false);
1085 *   sparsity_pattern.copy_from(dsp);
1086 *  
1087 *   system_matrix.reinit(sparsity_pattern);
1088 *  
1089 *   solution.reinit(dof_handler.n_dofs());
1090 *   system_rhs.reinit(dof_handler.n_dofs());
1091 *   }
1092 *  
1093 *  
1094 *  
1095 * @endcode
1096 *
1097 * In the following function, the matrix and right hand side are
1098 * assembled. As stated in the documentation of the main class above, it
1099 * does not do this itself, but rather delegates to the function following
1100 * next, utilizing the WorkStream concept discussed in @ref threads .
1101 *
1102
1103 *
1104 * If you have looked through the @ref threads module, you will have
1105 * seen that assembling in parallel does not take an incredible
1106 * amount of extra code as long as you diligently describe what the
1107 * scratch and copy data objects are, and if you define suitable
1108 * functions for the local assembly and the copy operation from local
1109 * contributions to global objects. This done, the following will do
1110 * all the heavy lifting to get these operations done on multiple
1111 * threads on as many cores as you have in your system:
1112 *
1113 * @code
1114 *   template <int dim>
1115 *   void AdvectionProblem<dim>::assemble_system()
1116 *   {
1117 *   WorkStream::run(dof_handler.begin_active(),
1118 *   dof_handler.end(),
1119 *   *this,
1120 *   &AdvectionProblem::local_assemble_system,
1121 *   &AdvectionProblem::copy_local_to_global,
1122 *   AssemblyScratchData(fe),
1123 *   AssemblyCopyData());
1124 *   }
1125 *  
1126 *  
1127 *  
1128 * @endcode
1129 *
1130 * As already mentioned above, we need to have scratch objects for
1131 * the parallel computation of local contributions. These objects
1132 * contain FEValues and FEFaceValues objects (as well as some arrays), and so
1133 * we will need to have constructors and copy constructors that allow us to
1134 * create them. For the cell terms we need the values
1135 * and gradients of the shape functions, the quadrature points in
1136 * order to determine the source density and the advection field at
1137 * a given point, and the weights of the quadrature points times the
1138 * determinant of the Jacobian at these points. In contrast, for the
1139 * boundary integrals, we don't need the gradients, but rather the
1140 * normal vectors to the cells. This determines which update flags
1141 * we will have to pass to the constructors of the members of the
1142 * class:
1143 *
1144 * @code
1145 *   template <int dim>
1146 *   AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
1147 *   const FiniteElement<dim> &fe)
1148 *   : fe_values(fe,
1149 *   QGauss<dim>(fe.degree + 1),
1151 *   update_JxW_values)
1152 *   , fe_face_values(fe,
1153 *   QGauss<dim - 1>(fe.degree + 1),
1156 *   , rhs_values(fe_values.get_quadrature().size())
1157 *   , advection_directions(fe_values.get_quadrature().size())
1158 *   , face_boundary_values(fe_face_values.get_quadrature().size())
1159 *   , face_advection_directions(fe_face_values.get_quadrature().size())
1160 *   {}
1161 *  
1162 *  
1163 *  
1164 *   template <int dim>
1165 *   AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
1166 *   const AssemblyScratchData &scratch_data)
1167 *   : fe_values(scratch_data.fe_values.get_fe(),
1168 *   scratch_data.fe_values.get_quadrature(),
1170 *   update_JxW_values)
1171 *   , fe_face_values(scratch_data.fe_face_values.get_fe(),
1172 *   scratch_data.fe_face_values.get_quadrature(),
1175 *   , rhs_values(scratch_data.rhs_values.size())
1176 *   , advection_directions(scratch_data.advection_directions.size())
1177 *   , face_boundary_values(scratch_data.face_boundary_values.size())
1178 *   , face_advection_directions(scratch_data.face_advection_directions.size())
1179 *   {}
1180 *  
1181 *  
1182 *  
1183 * @endcode
1184 *
1185 * Now, this is the function that does the actual work. It is not very
1186 * different from the <code>assemble_system</code> functions of previous
1187 * example programs, so we will again only comment on the differences. The
1188 * mathematical stuff closely follows what we have said in the introduction.
1189 *
1190
1191 *
1192 * There are a number of points worth mentioning here, though. The
1193 * first one is that we have moved the FEValues and FEFaceValues
1194 * objects into the ScratchData object. We have done so because the
1195 * alternative would have been to simply create one every time we
1196 * get into this function -- i.e., on every cell. It now turns out
1197 * that the FEValues classes were written with the explicit goal of
1198 * moving everything that remains the same from cell to cell into
1199 * the construction of the object, and only do as little work as
1200 * possible in FEValues::reinit() whenever we move to a new
1201 * cell. What this means is that it would be very expensive to
1202 * create a new object of this kind in this function as we would
1203 * have to do it for every cell -- exactly the thing we wanted to
1204 * avoid with the FEValues class. Instead, what we do is create it
1205 * only once (or a small number of times) in the scratch objects and
1206 * then re-use it as often as we can.
1207 *
1208
1209 *
1210 * This begs the question of whether there are other objects we
1211 * create in this function whose creation is expensive compared to
1212 * its use. Indeed, at the top of the function, we declare all sorts
1213 * of objects. The <code>AdvectionField</code>,
1214 * <code>RightHandSide</code> and <code>BoundaryValues</code> do not
1215 * cost much to create, so there is no harm here. However,
1216 * allocating memory in creating the <code>rhs_values</code> and
1217 * similar variables below typically costs a significant amount of
1218 * time, compared to just accessing the (temporary) values we store
1219 * in them. Consequently, these would be candidates for moving into
1220 * the <code>AssemblyScratchData</code> class. We will leave this as
1221 * an exercise.
1222 *
1223 * @code
1224 *   template <int dim>
1225 *   void AdvectionProblem<dim>::local_assemble_system(
1226 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1227 *   AssemblyScratchData & scratch_data,
1228 *   AssemblyCopyData & copy_data)
1229 *   {
1230 * @endcode
1231 *
1232 * We define some abbreviations to avoid unnecessarily long lines:
1233 *
1234 * @code
1235 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1236 *   const unsigned int n_q_points =
1237 *   scratch_data.fe_values.get_quadrature().size();
1238 *   const unsigned int n_face_q_points =
1239 *   scratch_data.fe_face_values.get_quadrature().size();
1240 *  
1241 * @endcode
1242 *
1243 * We declare cell matrix and cell right hand side...
1244 *
1245 * @code
1246 *   copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1247 *   copy_data.cell_rhs.reinit(dofs_per_cell);
1248 *  
1249 * @endcode
1250 *
1251 * ... an array to hold the global indices of the degrees of freedom of
1252 * the cell on which we are presently working...
1253 *
1254 * @code
1255 *   copy_data.local_dof_indices.resize(dofs_per_cell);
1256 *  
1257 * @endcode
1258 *
1259 * ... then initialize the <code>FEValues</code> object...
1260 *
1261 * @code
1262 *   scratch_data.fe_values.reinit(cell);
1263 *  
1264 * @endcode
1265 *
1266 * ... obtain the values of right hand side and advection directions
1267 * at the quadrature points...
1268 *
1269 * @code
1270 *   scratch_data.advection_field.value_list(
1271 *   scratch_data.fe_values.get_quadrature_points(),
1272 *   scratch_data.advection_directions);
1273 *   scratch_data.right_hand_side.value_list(
1274 *   scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
1275 *  
1276 * @endcode
1277 *
1278 * ... set the value of the streamline diffusion parameter as
1279 * described in the introduction...
1280 *
1281 * @code
1282 *   const double delta = 0.1 * cell->diameter();
1283 *  
1284 * @endcode
1285 *
1286 * ... and assemble the local contributions to the system matrix and
1287 * right hand side as also discussed above:
1288 *
1289 * @code
1290 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1291 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1292 *   {
1293 * @endcode
1294 *
1295 * Alias the AssemblyScratchData object to keep the lines from
1296 * getting too long:
1297 *
1298 * @code
1299 *   const auto &sd = scratch_data;
1300 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1301 *   copy_data.cell_matrix(i, j) +=
1302 *   ((sd.fe_values.shape_value(i, q_point) + // (phi_i +
1303 *   delta * (sd.advection_directions[q_point] * // delta beta
1304 *   sd.fe_values.shape_grad(i, q_point))) * // grad phi_i)
1305 *   sd.advection_directions[q_point] * // beta
1306 *   sd.fe_values.shape_grad(j, q_point)) * // grad phi_j
1307 *   sd.fe_values.JxW(q_point); // dx
1308 *  
1309 *   copy_data.cell_rhs(i) +=
1310 *   (sd.fe_values.shape_value(i, q_point) + // (phi_i +
1311 *   delta * (sd.advection_directions[q_point] * // delta beta
1312 *   sd.fe_values.shape_grad(i, q_point))) * // grad phi_i)
1313 *   sd.rhs_values[q_point] * // f
1314 *   sd.fe_values.JxW(q_point); // dx
1315 *   }
1316 *  
1317 * @endcode
1318 *
1319 * Besides the cell terms which we have built up now, the bilinear
1320 * form of the present problem also contains terms on the boundary of
1321 * the domain. Therefore, we have to check whether any of the faces of
1322 * this cell are on the boundary of the domain, and if so assemble the
1323 * contributions of this face as well. Of course, the bilinear form
1324 * only contains contributions from the <code>inflow</code> part of
1325 * the boundary, but to find out whether a certain part of a face of
1326 * the present cell is part of the inflow boundary, we have to have
1327 * information on the exact location of the quadrature points and on
1328 * the direction of flow at this point; we obtain this information
1329 * using the FEFaceValues object and only decide within the main loop
1330 * whether a quadrature point is on the inflow boundary.
1331 *
1332 * @code
1333 *   for (const auto &face : cell->face_iterators())
1334 *   if (face->at_boundary())
1335 *   {
1336 * @endcode
1337 *
1338 * Ok, this face of the present cell is on the boundary of the
1339 * domain. Just as for the usual FEValues object which we have
1340 * used in previous examples and also above, we have to
1341 * reinitialize the FEFaceValues object for the present face:
1342 *
1343 * @code
1344 *   scratch_data.fe_face_values.reinit(cell, face);
1345 *  
1346 * @endcode
1347 *
1348 * For the quadrature points at hand, we ask for the values of
1349 * the inflow function and for the direction of flow:
1350 *
1351 * @code
1352 *   scratch_data.boundary_values.value_list(
1353 *   scratch_data.fe_face_values.get_quadrature_points(),
1354 *   scratch_data.face_boundary_values);
1355 *   scratch_data.advection_field.value_list(
1356 *   scratch_data.fe_face_values.get_quadrature_points(),
1357 *   scratch_data.face_advection_directions);
1358 *  
1359 * @endcode
1360 *
1361 * Now loop over all quadrature points and see whether this face is on
1362 * the inflow or outflow part of the boundary. The normal
1363 * vector points out of the cell: since the face is at
1364 * the boundary, the normal vector points out of the domain,
1365 * so if the advection direction points into the domain, its
1366 * scalar product with the normal vector must be negative (to see why
1367 * this is true, consider the scalar product definition that uses a
1368 * cosine):
1369 *
1370 * @code
1371 *   for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
1372 *   if (scratch_data.fe_face_values.normal_vector(q_point) *
1373 *   scratch_data.face_advection_directions[q_point] <
1374 *   0.)
1375 * @endcode
1376 *
1377 * If the face is part of the inflow boundary, then compute the
1378 * contributions of this face to the global matrix and right
1379 * hand side, using the values obtained from the
1380 * FEFaceValues object and the formulae discussed in the
1381 * introduction:
1382 *
1383 * @code
1384 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1385 *   {
1386 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1387 *   copy_data.cell_matrix(i, j) -=
1388 *   (scratch_data.face_advection_directions[q_point] *
1389 *   scratch_data.fe_face_values.normal_vector(q_point) *
1390 *   scratch_data.fe_face_values.shape_value(i, q_point) *
1391 *   scratch_data.fe_face_values.shape_value(j, q_point) *
1392 *   scratch_data.fe_face_values.JxW(q_point));
1393 *  
1394 *   copy_data.cell_rhs(i) -=
1395 *   (scratch_data.face_advection_directions[q_point] *
1396 *   scratch_data.fe_face_values.normal_vector(q_point) *
1397 *   scratch_data.face_boundary_values[q_point] *
1398 *   scratch_data.fe_face_values.shape_value(i, q_point) *
1399 *   scratch_data.fe_face_values.JxW(q_point));
1400 *   }
1401 *   }
1402 *  
1403 * @endcode
1404 *
1405 * The final piece of information the copy routine needs is the global
1406 * indices of the degrees of freedom on this cell, so we end by writing
1407 * them to the local array:
1408 *
1409 * @code
1410 *   cell->get_dof_indices(copy_data.local_dof_indices);
1411 *   }
1412 *  
1413 *  
1414 *  
1415 * @endcode
1416 *
1417 * The second function we needed to write was the one that copies
1418 * the local contributions the previous function computed (and
1419 * put into the AssemblyCopyData object) into the global matrix and right
1420 * hand side vector objects. This is essentially what we always had
1421 * as the last block of code when assembling something on every
1422 * cell. The following should therefore be pretty obvious:
1423 *
1424 * @code
1425 *   template <int dim>
1426 *   void
1427 *   AdvectionProblem<dim>::copy_local_to_global(const AssemblyCopyData &copy_data)
1428 *   {
1429 *   hanging_node_constraints.distribute_local_to_global(
1430 *   copy_data.cell_matrix,
1431 *   copy_data.cell_rhs,
1432 *   copy_data.local_dof_indices,
1433 *   system_matrix,
1434 *   system_rhs);
1435 *   }
1436 *  
1437 * @endcode
1438 *
1439 * Here comes the linear solver routine. As the system is no longer
1440 * symmetric positive definite as in all the previous examples, we cannot
1441 * use the Conjugate Gradient method anymore. Rather, we use a solver that
1442 * is more general and does not rely on any special properties of the
1443 * matrix: the GMRES method. GMRES, like the conjugate gradient method,
1444 * requires a decent preconditioner: we use a Jacobi preconditioner here,
1445 * which works well enough for this problem.
1446 *
1447 * @code
1448 *   template <int dim>
1449 *   void AdvectionProblem<dim>::solve()
1450 *   {
1451 *   SolverControl solver_control(std::max<std::size_t>(1000,
1452 *   system_rhs.size() / 10),
1453 *   1e-10 * system_rhs.l2_norm());
1454 *   SolverGMRES<Vector<double>> solver(solver_control);
1455 *   PreconditionJacobi<SparseMatrix<double>> preconditioner;
1456 *   preconditioner.initialize(system_matrix, 1.0);
1457 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
1458 *  
1459 *   Vector<double> residual(dof_handler.n_dofs());
1460 *  
1461 *   system_matrix.vmult(residual, solution);
1462 *   residual -= system_rhs;
1463 *   std::cout << " Iterations required for convergence: "
1464 *   << solver_control.last_step() << '\n'
1465 *   << " Max norm of residual: "
1466 *   << residual.linfty_norm() << '\n';
1467 *  
1468 *   hanging_node_constraints.distribute(solution);
1469 *   }
1470 *  
1471 * @endcode
1472 *
1473 * The following function refines the grid according to the quantity
1474 * described in the introduction. The respective computations are made in
1475 * the class <code>GradientEstimation</code>.
1476 *
1477 * @code
1478 *   template <int dim>
1479 *   void AdvectionProblem<dim>::refine_grid()
1480 *   {
1481 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1482 *  
1483 *   GradientEstimation::estimate(dof_handler,
1484 *   solution,
1485 *   estimated_error_per_cell);
1486 *  
1488 *   estimated_error_per_cell,
1489 *   0.3,
1490 *   0.03);
1491 *  
1492 *   triangulation.execute_coarsening_and_refinement();
1493 *   }
1494 *  
1495 * @endcode
1496 *
1497 * This function is similar to the one in step 6, but since we use a higher
1498 * degree finite element we save the solution in a different
1499 * way. Visualization programs like VisIt and Paraview typically only
1500 * understand data that is associated with nodes: they cannot plot
1501 * fifth-degree basis functions, which results in a very inaccurate picture
1502 * of the solution we computed. To get around this we save multiple
1503 * <em>patches</em> per cell: in 2d we save 64 bilinear `cells' to the VTU
1504 * file for each cell, and in 3d we save 512. The end result is that the
1505 * visualization program will use a piecewise linear interpolation of the
1506 * cubic basis functions: this captures the solution detail and, with most
1507 * screen resolutions, looks smooth. We save the grid in a separate step
1508 * with no extra patches so that we have a visual representation of the cell
1509 * faces.
1510 *
1511
1512 *
1513 * Version 9.1 of deal.II gained the ability to write higher degree
1514 * polynomials (i.e., write piecewise bicubic visualization data for our
1515 * piecewise bicubic solution) VTK and VTU output: however, not all recent
1516 * versions of ParaView and VisIt (as of 2018) can read this format, so we
1517 * use the older, more general (but less efficient) approach here.
1518 *
1519 * @code
1520 *   template <int dim>
1521 *   void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1522 *   {
1523 *   {
1524 *   GridOut grid_out;
1525 *   std::ofstream output("grid-" + std::to_string(cycle) + ".vtu");
1526 *   grid_out.write_vtu(triangulation, output);
1527 *   }
1528 *  
1529 *   {
1530 *   DataOut<dim> data_out;
1531 *   data_out.attach_dof_handler(dof_handler);
1532 *   data_out.add_data_vector(solution, "solution");
1533 *   data_out.build_patches(8);
1534 *  
1535 * @endcode
1536 *
1537 * VTU output can be expensive, both to compute and to write to
1538 * disk. Here we ask ZLib, a compression library, to compress the data
1539 * in a way that maximizes throughput.
1540 *
1541 * @code
1542 *   DataOutBase::VtkFlags vtk_flags;
1543 *   vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
1544 *   data_out.set_flags(vtk_flags);
1545 *  
1546 *   std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
1547 *   data_out.write_vtu(output);
1548 *   }
1549 *   }
1550 *  
1551 *  
1552 * @endcode
1553 *
1554 * ... as is the main loop (setup -- solve -- refine), aside from the number
1555 * of cycles and the initial grid:
1556 *
1557 * @code
1558 *   template <int dim>
1559 *   void AdvectionProblem<dim>::run()
1560 *   {
1561 *   for (unsigned int cycle = 0; cycle < 10; ++cycle)
1562 *   {
1563 *   std::cout << "Cycle " << cycle << ':' << std::endl;
1564 *  
1565 *   if (cycle == 0)
1566 *   {
1567 *   GridGenerator::hyper_cube(triangulation, -1, 1);
1568 *   triangulation.refine_global(3);
1569 *   }
1570 *   else
1571 *   {
1572 *   refine_grid();
1573 *   }
1574 *  
1575 *  
1576 *   std::cout << " Number of active cells: "
1577 *   << triangulation.n_active_cells() << std::endl;
1578 *  
1579 *   setup_system();
1580 *  
1581 *   std::cout << " Number of degrees of freedom: "
1582 *   << dof_handler.n_dofs() << std::endl;
1583 *  
1584 *   assemble_system();
1585 *   solve();
1586 *   output_results(cycle);
1587 *   }
1588 *   }
1589 *  
1590 *  
1591 *  
1592 * @endcode
1593 *
1594 *
1595 * <a name="GradientEstimationclassimplementation"></a>
1596 * <h3>GradientEstimation class implementation</h3>
1597 *
1598
1599 *
1600 * Now for the implementation of the <code>GradientEstimation</code> class.
1601 * Let us start by defining constructors for the
1602 * <code>EstimateScratchData</code> class used by the
1603 * <code>estimate_cell()</code> function:
1604 *
1605 * @code
1606 *   template <int dim>
1607 *   GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1608 *   const FiniteElement<dim> &fe,
1609 *   const Vector<double> & solution,
1610 *   Vector<float> & error_per_cell)
1611 *   : fe_midpoint_value(fe,
1612 *   QMidpoint<dim>(),
1613 *   update_values | update_quadrature_points)
1614 *   , solution(solution)
1615 *   , error_per_cell(error_per_cell)
1616 *   , cell_midpoint_value(1)
1617 *   , neighbor_midpoint_value(1)
1618 *   {
1619 * @endcode
1620 *
1621 * We allocate a vector to hold iterators to all active neighbors of
1622 * a cell. We reserve the maximal number of active neighbors in order to
1623 * avoid later reallocations. Note how this maximal number of active
1624 * neighbors is computed here.
1625 *
1626 * @code
1627 *   active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
1628 *   GeometryInfo<dim>::max_children_per_face);
1629 *   }
1630 *  
1631 *  
1632 *   template <int dim>
1633 *   GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1634 *   const EstimateScratchData &scratch_data)
1635 *   : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(),
1636 *   scratch_data.fe_midpoint_value.get_quadrature(),
1637 *   update_values | update_quadrature_points)
1638 *   , solution(scratch_data.solution)
1639 *   , error_per_cell(scratch_data.error_per_cell)
1640 *   , cell_midpoint_value(1)
1641 *   , neighbor_midpoint_value(1)
1642 *   {}
1643 *  
1644 *  
1645 * @endcode
1646 *
1647 * Next comes the implementation of the <code>GradientEstimation</code>
1648 * class. The first function does not much except for delegating work to the
1649 * other function, but there is a bit of setup at the top.
1650 *
1651
1652 *
1653 * Before starting with the work, we check that the vector into
1654 * which the results are written has the right size. Programming
1655 * mistakes in which one forgets to size arguments correctly at the
1656 * calling site are quite common. Because the resulting damage from
1657 * not catching such errors is often subtle (e.g., corruption of
1658 * data somewhere in memory, or non-reproducible results), it is
1659 * well worth the effort to check for such things.
1660 *
1661 * @code
1662 *   template <int dim>
1663 *   void GradientEstimation::estimate(const DoFHandler<dim> &dof_handler,
1664 *   const Vector<double> & solution,
1665 *   Vector<float> & error_per_cell)
1666 *   {
1667 *   Assert(
1668 *   error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(),
1669 *   ExcInvalidVectorLength(error_per_cell.size(),
1670 *   dof_handler.get_triangulation().n_active_cells()));
1671 *  
1672 *   WorkStream::run(dof_handler.begin_active(),
1673 *   dof_handler.end(),
1674 *   &GradientEstimation::template estimate_cell<dim>,
1675 *   std::function<void(const EstimateCopyData &)>(),
1676 *   EstimateScratchData<dim>(dof_handler.get_fe(),
1677 *   solution,
1678 *   error_per_cell),
1679 *   EstimateCopyData());
1680 *   }
1681 *  
1682 *  
1683 * @endcode
1684 *
1685 * Here comes the function that estimates the local error by computing the
1686 * finite difference approximation of the gradient. The function first
1687 * computes the list of active neighbors of the present cell and then
1688 * computes the quantities described in the introduction for each of
1689 * the neighbors. The reason for this order is that it is not a one-liner
1690 * to find a given neighbor with locally refined meshes. In principle, an
1691 * optimized implementation would find neighbors and the quantities
1692 * depending on them in one step, rather than first building a list of
1693 * neighbors and in a second step their contributions but we will gladly
1694 * leave this as an exercise. As discussed before, the worker function
1695 * passed to WorkStream::run works on "scratch" objects that keep all
1696 * temporary objects. This way, we do not need to create and initialize
1697 * objects that are expensive to initialize within the function that does
1698 * the work every time it is called for a given cell. Such an argument is
1699 * passed as the second argument. The third argument would be a "copy-data"
1700 * object (see @ref threads for more information) but we do not actually use
1701 * any of these here. Since WorkStream::run() insists on passing three
1702 * arguments, we declare this function with three arguments, but simply
1703 * ignore the last one.
1704 *
1705
1706 *
1707 * (This is unsatisfactory from an aesthetic perspective. It can be avoided
1708 * by using an anonymous (lambda) function. If you allow, let us here show
1709 * how. First, assume that we had declared this function to only take two
1710 * arguments by omitting the unused last one. Now, WorkStream::run still
1711 * wants to call this function with three arguments, so we need to find a
1712 * way to "forget" the third argument in the call. Simply passing
1713 * WorkStream::run the pointer to the function as we do above will not do
1714 * this -- the compiler will complain that a function declared to have two
1715 * arguments is called with three arguments. However, we can do this by
1716 * passing the following as the third argument to WorkStream::run():
1717 * <div class=CodeFragmentInTutorialComment>
1718 * @code
1719 * [](const typename DoFHandler<dim>::active_cell_iterator &cell,
1720 * EstimateScratchData<dim> & scratch_data,
1721 * EstimateCopyData &)
1722 * {
1723 * GradientEstimation::estimate_cell<dim>(cell, scratch_data);
1724 * }
1725 * @endcode
1726 * </div>
1727 * This is not much better than the solution implemented below: either the
1728 * routine itself must take three arguments or it must be wrapped by
1729 * something that takes three arguments. We don't use this since adding the
1730 * unused argument at the beginning is simpler.
1731 *
1732
1733 *
1734 * Now for the details:
1735 *
1736 * @code
1737 *   template <int dim>
1738 *   void GradientEstimation::estimate_cell(
1739 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1740 *   EstimateScratchData<dim> & scratch_data,
1741 *   const EstimateCopyData &)
1742 *   {
1743 * @endcode
1744 *
1745 * We need space for the tensor <code>Y</code>, which is the sum of
1746 * outer products of the y-vectors.
1747 *
1748 * @code
1749 *   Tensor<2, dim> Y;
1750 *  
1751 * @endcode
1752 *
1753 * First initialize the <code>FEValues</code> object, as well as the
1754 * <code>Y</code> tensor:
1755 *
1756 * @code
1757 *   scratch_data.fe_midpoint_value.reinit(cell);
1758 *  
1759 * @endcode
1760 *
1761 * Now, before we go on, we first compute a list of all active neighbors
1762 * of the present cell. We do so by first looping over all faces and see
1763 * whether the neighbor there is active, which would be the case if it
1764 * is on the same level as the present cell or one level coarser (note
1765 * that a neighbor can only be once coarser than the present cell, as
1766 * we only allow a maximal difference of one refinement over a face in
1767 * deal.II). Alternatively, the neighbor could be on the same level
1768 * and be further refined; then we have to find which of its children
1769 * are next to the present cell and select these (note that if a child
1770 * of a neighbor of an active cell that is next to this active cell,
1771 * needs necessarily be active itself, due to the one-refinement rule
1772 * cited above).
1773 *
1774
1775 *
1776 * Things are slightly different in one space dimension, as there the
1777 * one-refinement rule does not exist: neighboring active cells may
1778 * differ in as many refinement levels as they like. In this case, the
1779 * computation becomes a little more difficult, but we will explain
1780 * this below.
1781 *
1782
1783 *
1784 * Before starting the loop over all neighbors of the present cell, we
1785 * have to clear the array storing the iterators to the active
1786 * neighbors, of course.
1787 *
1788 * @code
1789 *   scratch_data.active_neighbors.clear();
1790 *   for (const auto face_n : cell->face_indices())
1791 *   if (!cell->at_boundary(face_n))
1792 *   {
1793 * @endcode
1794 *
1795 * First define an abbreviation for the iterator to the face and
1796 * the neighbor
1797 *
1798 * @code
1799 *   const auto face = cell->face(face_n);
1800 *   const auto neighbor = cell->neighbor(face_n);
1801 *  
1802 * @endcode
1803 *
1804 * Then check whether the neighbor is active. If it is, then it
1805 * is on the same level or one level coarser (if we are not in
1806 * 1d), and we are interested in it in any case.
1807 *
1808 * @code
1809 *   if (neighbor->is_active())
1810 *   scratch_data.active_neighbors.push_back(neighbor);
1811 *   else
1812 *   {
1813 * @endcode
1814 *
1815 * If the neighbor is not active, then check its children.
1816 *
1817 * @code
1818 *   if (dim == 1)
1819 *   {
1820 * @endcode
1821 *
1822 * To find the child of the neighbor which bounds to the
1823 * present cell, successively go to its right child if
1824 * we are left of the present cell (n==0), or go to the
1825 * left child if we are on the right (n==1), until we
1826 * find an active cell.
1827 *
1828 * @code
1829 *   auto neighbor_child = neighbor;
1830 *   while (neighbor_child->has_children())
1831 *   neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0);
1832 *  
1833 * @endcode
1834 *
1835 * As this used some non-trivial geometrical intuition,
1836 * we might want to check whether we did it right,
1837 * i.e., check whether the neighbor of the cell we found
1838 * is indeed the cell we are presently working
1839 * on. Checks like this are often useful and have
1840 * frequently uncovered errors both in algorithms like
1841 * the line above (where it is simple to involuntarily
1842 * exchange <code>n==1</code> for <code>n==0</code> or
1843 * the like) and in the library (the assumptions
1844 * underlying the algorithm above could either be wrong,
1845 * wrongly documented, or are violated due to an error
1846 * in the library). One could in principle remove such
1847 * checks after the program works for some time, but it
1848 * might be a good things to leave it in anyway to check
1849 * for changes in the library or in the algorithm above.
1850 *
1851
1852 *
1853 * Note that if this check fails, then this is certainly
1854 * an error that is irrecoverable and probably qualifies
1855 * as an internal error. We therefore use a predefined
1856 * exception class to throw here.
1857 *
1858 * @code
1859 *   Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell,
1860 *   ExcInternalError());
1861 *  
1862 * @endcode
1863 *
1864 * If the check succeeded, we push the active neighbor
1865 * we just found to the stack we keep:
1866 *
1867 * @code
1868 *   scratch_data.active_neighbors.push_back(neighbor_child);
1869 *   }
1870 *   else
1871 * @endcode
1872 *
1873 * If we are not in 1d, we collect all neighbor children
1874 * `behind' the subfaces of the current face and move on:
1875 *
1876 * @code
1877 *   for (unsigned int subface_n = 0; subface_n < face->n_children();
1878 *   ++subface_n)
1879 *   scratch_data.active_neighbors.push_back(
1880 *   cell->neighbor_child_on_subface(face_n, subface_n));
1881 *   }
1882 *   }
1883 *  
1884 * @endcode
1885 *
1886 * OK, now that we have all the neighbors, lets start the computation
1887 * on each of them. First we do some preliminaries: find out about the
1888 * center of the present cell and the solution at this point. The
1889 * latter is obtained as a vector of function values at the quadrature
1890 * points, of which there are only one, of course. Likewise, the
1891 * position of the center is the position of the first (and only)
1892 * quadrature point in real space.
1893 *
1894 * @code
1895 *   const Point<dim> this_center =
1896 *   scratch_data.fe_midpoint_value.quadrature_point(0);
1897 *  
1898 *   scratch_data.fe_midpoint_value.get_function_values(
1899 *   scratch_data.solution, scratch_data.cell_midpoint_value);
1900 *  
1901 * @endcode
1902 *
1903 * Now loop over all active neighbors and collect the data we
1904 * need.
1905 *
1906 * @code
1907 *   Tensor<1, dim> projected_gradient;
1908 *   for (const auto &neighbor : scratch_data.active_neighbors)
1909 *   {
1910 * @endcode
1911 *
1912 * Then get the center of the neighbor cell and the value of the
1913 * finite element function at that point. Note that for this
1914 * information we have to reinitialize the <code>FEValues</code>
1915 * object for the neighbor cell.
1916 *
1917 * @code
1918 *   scratch_data.fe_midpoint_value.reinit(neighbor);
1919 *   const Point<dim> neighbor_center =
1920 *   scratch_data.fe_midpoint_value.quadrature_point(0);
1921 *  
1922 *   scratch_data.fe_midpoint_value.get_function_values(
1923 *   scratch_data.solution, scratch_data.neighbor_midpoint_value);
1924 *  
1925 * @endcode
1926 *
1927 * Compute the vector <code>y</code> connecting the centers of the
1928 * two cells. Note that as opposed to the introduction, we denote
1929 * by <code>y</code> the normalized difference vector, as this is
1930 * the quantity used everywhere in the computations.
1931 *
1932 * @code
1933 *   Tensor<1, dim> y = neighbor_center - this_center;
1934 *   const double distance = y.norm();
1935 *   y /= distance;
1936 *  
1937 * @endcode
1938 *
1939 * Then add up the contribution of this cell to the Y matrix...
1940 *
1941 * @code
1942 *   for (unsigned int i = 0; i < dim; ++i)
1943 *   for (unsigned int j = 0; j < dim; ++j)
1944 *   Y[i][j] += y[i] * y[j];
1945 *  
1946 * @endcode
1947 *
1948 * ... and update the sum of difference quotients:
1949 *
1950 * @code
1951 *   projected_gradient += (scratch_data.neighbor_midpoint_value[0] -
1952 *   scratch_data.cell_midpoint_value[0]) /
1953 *   distance * y;
1954 *   }
1955 *  
1956 * @endcode
1957 *
1958 * If now, after collecting all the information from the neighbors, we
1959 * can determine an approximation of the gradient for the present
1960 * cell, then we need to have passed over vectors <code>y</code> which
1961 * span the whole space, otherwise we would not have all components of
1962 * the gradient. This is indicated by the invertibility of the matrix.
1963 *
1964
1965 *
1966 * If the matrix is not invertible, then the present
1967 * cell had an insufficient number of active neighbors. In contrast to
1968 * all previous cases (where we raised exceptions) this is, however,
1969 * not a programming error: it is a runtime error that can happen in
1970 * optimized mode even if it ran well in debug mode, so it is
1971 * reasonable to try to catch this error also in optimized mode. For
1972 * this case, there is the <code>AssertThrow</code> macro: it checks
1973 * the condition like the <code>Assert</code> macro, but not only in
1974 * debug mode; it then outputs an error message, but instead of
1975 * aborting the program as in the case of the <code>Assert</code>
1976 * macro, the exception is thrown using the <code>throw</code> command
1977 * of C++. This way, one has the possibility to catch this error and
1978 * take reasonable counter actions. One such measure would be to
1979 * refine the grid globally, as the case of insufficient directions
1980 * can not occur if every cell of the initial grid has been refined at
1981 * least once.
1982 *
1983 * @code
1984 *   AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
1985 *  
1986 * @endcode
1987 *
1988 * If, on the other hand, the matrix is invertible, then invert it,
1989 * multiply the other quantity with it, and compute the estimated error
1990 * using this quantity and the correct powers of the mesh width:
1991 *
1992 * @code
1993 *   const Tensor<2, dim> Y_inverse = invert(Y);
1994 *  
1995 *   const Tensor<1, dim> gradient = Y_inverse * projected_gradient;
1996 *  
1997 * @endcode
1998 *
1999 * The last part of this function is the one where we write into
2000 * the element of the output vector what we have just
2001 * computed. The address of this vector has been stored in the
2002 * scratch data object, and all we have to do is know how to get
2003 * at the correct element inside this vector -- but we can ask the
2004 * cell we're on the how-manyth active cell it is for this:
2005 *
2006 * @code
2007 *   scratch_data.error_per_cell(cell->active_cell_index()) =
2008 *   (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) * gradient.norm());
2009 *   }
2010 *   } // namespace Step9
2011 *  
2012 *  
2013 * @endcode
2014 *
2015 *
2016 * <a name="Mainfunction"></a>
2017 * <h3>Main function</h3>
2018 *
2019
2020 *
2021 * The <code>main</code> function is similar to the previous examples. The
2022 * primary difference is that we use MultithreadInfo to set the maximum
2023 * number of threads (see the documentation module @ref threads
2024 * "Parallel computing with multiple processors accessing shared memory"
2025 * for more information). The number of threads used is the minimum of the
2026 * environment variable DEAL_II_NUM_THREADS and the parameter of
2027 * <code>set_thread_limit</code>. If no value is given to
2028 * <code>set_thread_limit</code>, the default value from the Intel Threading
2029 * Building Blocks (TBB) library is used. If the call to
2030 * <code>set_thread_limit</code> is omitted, the number of threads will be
2031 * chosen by TBB independently of DEAL_II_NUM_THREADS.
2032 *
2033 * @code
2034 *   int main()
2035 *   {
2036 *   using namespace dealii;
2037 *   try
2038 *   {
2040 *  
2041 *   Step9::AdvectionProblem<2> advection_problem_2d;
2042 *   advection_problem_2d.run();
2043 *   }
2044 *   catch (std::exception &exc)
2045 *   {
2046 *   std::cerr << std::endl
2047 *   << std::endl
2048 *   << "----------------------------------------------------"
2049 *   << std::endl;
2050 *   std::cerr << "Exception on processing: " << std::endl
2051 *   << exc.what() << std::endl
2052 *   << "Aborting!" << std::endl
2053 *   << "----------------------------------------------------"
2054 *   << std::endl;
2055 *   return 1;
2056 *   }
2057 *   catch (...)
2058 *   {
2059 *   std::cerr << std::endl
2060 *   << std::endl
2061 *   << "----------------------------------------------------"
2062 *   << std::endl;
2063 *   std::cerr << "Unknown exception!" << std::endl
2064 *   << "Aborting!" << std::endl
2065 *   << "----------------------------------------------------"
2066 *   << std::endl;
2067 *   return 1;
2068 *   }
2069 *  
2070 *   return 0;
2071 *   }
2072 * @endcode
2073<a name="Results"></a><h1>Results</h1>
2074
2075
2076
2077The results of this program are not particularly spectacular. They
2078consist of the console output, some grid files, and the solution on
2079each of these grids. First for the console output:
2080@code
2081Cycle 0:
2082 Number of active cells: 64
2083 Number of degrees of freedom: 1681
2084 Iterations required for convergence: 298
2085 Max norm of residual: 3.60316e-12
2086Cycle 1:
2087 Number of active cells: 124
2088 Number of degrees of freedom: 3537
2089 Iterations required for convergence: 415
2090 Max norm of residual: 3.70682e-12
2091Cycle 2:
2092 Number of active cells: 247
2093 Number of degrees of freedom: 6734
2094 Iterations required for convergence: 543
2095 Max norm of residual: 7.19716e-13
2096Cycle 3:
2097 Number of active cells: 502
2098 Number of degrees of freedom: 14105
2099 Iterations required for convergence: 666
2100 Max norm of residual: 3.45628e-13
2101Cycle 4:
2102 Number of active cells: 1003
2103 Number of degrees of freedom: 27462
2104 Iterations required for convergence: 1064
2105 Max norm of residual: 1.86495e-13
2106Cycle 5:
2107 Number of active cells: 1993
2108 Number of degrees of freedom: 55044
2109 Iterations required for convergence: 1251
2110 Max norm of residual: 1.28765e-13
2111Cycle 6:
2112 Number of active cells: 3985
2113 Number of degrees of freedom: 108492
2114 Iterations required for convergence: 2035
2115 Max norm of residual: 6.78085e-14
2116Cycle 7:
2117 Number of active cells: 7747
2118 Number of degrees of freedom: 210612
2119 Iterations required for convergence: 2187
2120 Max norm of residual: 2.61457e-14
2121Cycle 8:
2122 Number of active cells: 15067
2123 Number of degrees of freedom: 406907
2124 Iterations required for convergence: 3079
2125 Max norm of residual: 2.9932e-14
2126Cycle 9:
2127 Number of active cells: 29341
2128 Number of degrees of freedom: 780591
2129 Iterations required for convergence: 3913
2130 Max norm of residual: 8.15689e-15
2131@endcode
2132
2133Quite a number of cells are used on the finest level to resolve the features of
2134the solution. Here are the fourth and tenth grids:
2135<div class="twocolumn" style="width: 80%">
2136 <div>
2137 <img src="https://www.dealii.org/images/steps/developer/step-9-grid-3.png"
2138 alt="Fourth grid in the refinement cycle, showing some adaptivity to features."
2139 width="400" height="400">
2140 </div>
2141 <div>
2142 <img src="https://www.dealii.org/images/steps/developer/step-9-grid-9.png"
2143 alt="Tenth grid in the refinement cycle, showing that the waves are fully captured."
2144 width="400" height="400">
2145 </div>
2146</div>
2147and the fourth and tenth solutions:
2148<div class="twocolumn" style="width: 80%">
2149 <div>
2150 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-3.png"
2151 alt="Fourth solution, showing that we resolve most features but some
2152 are sill unresolved and appear blury."
2153 width="400" height="400">
2154 </div>
2155 <div>
2156 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-9.png"
2157 alt="Tenth solution, showing a fully resolved flow."
2158 width="400" height="400">
2159 </div>
2160</div>
2161and both the grid and solution zoomed in:
2162<div class="twocolumn" style="width: 80%">
2163 <div>
2164 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-3-zoom.png"
2165 alt="Detail of the fourth solution, showing that we resolve most
2166 features but some are sill unresolved and appear blury. In particular,
2167 the larger cells need to be refined."
2168 width="400" height="400">
2169 </div>
2170 <div>
2171 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-9-zoom.png"
2172 alt="Detail of the tenth solution, showing that we needed a lot more
2173 cells than were present in the fourth solution."
2174 width="400" height="400">
2175 </div>
2176</div>
2177
2178The solution is created by that part that is transported along the wiggly
2179advection field from the left and lower boundaries to the top right, and the
2180part that is created by the source in the lower left corner, and the results of
2181which are also transported along. The grid shown above is well-adapted to
2182resolve these features. The comparison between plots shows that, even though we
2183are using a high-order approximation, we still need adaptive mesh refinement to
2184fully resolve the wiggles.
2185 *
2186 *
2187<a name="PlainProg"></a>
2188<h1> The plain program</h1>
2189@include "step-9.cc"
2190*/
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:439
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
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.
@ symmetric
Matrix is symmetric.
@ general
No special properties.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
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:71
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)