1069 * sparsity_pattern.copy_from(dsp);
1071 * system_matrix.reinit(sparsity_pattern);
1073 * solution.reinit(dof_handler.n_dofs());
1074 * system_rhs.reinit(dof_handler.n_dofs());
1081 * In the following function, the
matrix and right hand side are
1082 * assembled. As stated in the documentation of the main
class above, it
1083 * does not
do this itself, but rather delegates to the function following
1084 * next, utilizing the
WorkStream concept discussed in @ref threads .
1088 * If you have looked through the @ref threads topic, you will have
1089 * seen that assembling in
parallel does not take an incredible
1090 * amount of extra code as long as you diligently describe what the
1091 * scratch and
copy data objects are, and if you define suitable
1092 *
functions for the local assembly and the
copy operation from local
1093 * contributions to global objects. This done, the following will do
1094 * all the heavy lifting to get these operations done on multiple
1095 * threads on as many cores as you have in your system:
1098 * template <
int dim>
1099 * void AdvectionProblem<dim>::assemble_system()
1102 * dof_handler.
end(),
1104 * &AdvectionProblem::local_assemble_system,
1105 * &AdvectionProblem::copy_local_to_global,
1106 * AssemblyScratchData(fe),
1107 * AssemblyCopyData());
1114 * As already mentioned above, we need to have scratch objects for
1115 * the
parallel computation of local contributions. These objects
1117 * we will need to have constructors and
copy constructors that allow us to
1118 * create them. For the cell terms we need the
values
1120 * order to determine the source density and the advection field at
1121 * a given
point, and the weights of the quadrature points times the
1122 *
determinant of the Jacobian at these points. In contrast, for the
1123 * boundary integrals, we don't need the
gradients, but rather the
1124 * normal vectors to the cells. This determines which update flags
1125 * we will have to pass to the constructors of the members of the
1129 * template <
int dim>
1130 * AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
1133 *
QGauss<dim>(fe.degree + 1),
1136 * , fe_face_values(fe,
1137 *
QGauss<dim - 1>(fe.degree + 1),
1140 * , rhs_values(fe_values.get_quadrature().size())
1141 * , advection_directions(fe_values.get_quadrature().size())
1142 * , face_boundary_values(fe_face_values.get_quadrature().size())
1143 * , face_advection_directions(fe_face_values.get_quadrature().size())
1148 * template <
int dim>
1149 * AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
1150 * const AssemblyScratchData &scratch_data)
1151 * : fe_values(scratch_data.fe_values.get_fe(),
1152 * scratch_data.fe_values.get_quadrature(),
1155 * , fe_face_values(scratch_data.fe_face_values.get_fe(),
1156 * scratch_data.fe_face_values.get_quadrature(),
1159 * , rhs_values(scratch_data.rhs_values.size())
1160 * , advection_directions(scratch_data.advection_directions.size())
1161 * , face_boundary_values(scratch_data.face_boundary_values.size())
1162 * , face_advection_directions(scratch_data.face_advection_directions.size())
1169 * Now, this is the function that does the actual work. It is not very
1170 * different from the <code>assemble_system</code>
functions of previous
1171 * example programs, so we will again only comment on the differences. The
1172 * mathematical stuff closely follows what we have said in the introduction.
1176 * There are a number of points worth mentioning here, though. The
1178 * objects into the ScratchData object. We have done so because the
1179 * alternative would have been to simply create one every time we
1180 * get into this function -- i.
e., on every cell. It now turns out
1181 * that the
FEValues classes were written with the explicit goal of
1182 * moving everything that remains the same from cell to cell into
1183 * the construction of the object, and only do as little work as
1185 * cell. What this means is that it would be very expensive to
1186 * create a new object of this kind in this function as we would
1187 * have to do it for every cell -- exactly the thing we wanted to
1188 * avoid with the
FEValues class. Instead, what we do is create it
1189 * only once (or a small number of times) in the scratch objects and
1190 * then re-use it as often as we can.
1194 * This begs the question of whether there are other objects we
1195 * create in this function whose creation is expensive compared to
1196 * its use. Indeed, at the top of the function, we declare all sorts
1197 * of objects. The <code>AdvectionField</code>,
1198 * <code>RightHandSide</code> and <code>BoundaryValues</code> do not
1199 * cost much to create, so there is no harm here. However,
1200 * allocating memory in creating the <code>rhs_values</code> and
1201 * similar variables below typically costs a significant amount of
1202 * time, compared to just accessing the (temporary)
values we store
1203 * in them. Consequently, these would be candidates for moving into
1204 * the <code>AssemblyScratchData</code> class. We will leave this as
1208 * template <
int dim>
1209 * void AdvectionProblem<dim>::local_assemble_system(
1210 * const typename
DoFHandler<dim>::active_cell_iterator &cell,
1211 * AssemblyScratchData &scratch_data,
1212 * AssemblyCopyData ©_data)
1216 * We define some abbreviations to avoid unnecessarily long lines:
1219 * const unsigned
int dofs_per_cell = fe.n_dofs_per_cell();
1220 *
const unsigned int n_q_points =
1221 * scratch_data.fe_values.get_quadrature().size();
1222 *
const unsigned int n_face_q_points =
1223 * scratch_data.fe_face_values.get_quadrature().size();
1227 * We declare cell
matrix and cell right hand side...
1230 * copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1231 * copy_data.cell_rhs.reinit(dofs_per_cell);
1235 * ... an array to hold the global indices of the degrees of freedom of
1236 * the cell on which we are presently working...
1239 * copy_data.local_dof_indices.resize(dofs_per_cell);
1243 * ... then initialize the <code>
FEValues</code>
object...
1246 * scratch_data.fe_values.
reinit(cell);
1250 * ... obtain the
values of right hand side and advection directions
1251 * at the quadrature points...
1254 * scratch_data.advection_field.value_list(
1255 * scratch_data.fe_values.get_quadrature_points(),
1256 * scratch_data.advection_directions);
1257 * scratch_data.right_hand_side.value_list(
1258 * scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
1262 * ... set the
value of the streamline diffusion parameter as
1263 * described in the introduction...
1266 *
const double delta = 0.1 * cell->diameter();
1270 * ... and
assemble the local contributions to the system
matrix and
1271 * right hand side as also discussed above:
1274 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1275 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1279 * Alias the AssemblyScratchData
object to keep the lines from
1283 *
const auto &sd = scratch_data;
1284 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1285 * copy_data.cell_matrix(i, j) +=
1286 * ((sd.fe_values.shape_value(i, q_point) +
1287 * delta * (sd.advection_directions[q_point] *
1288 * sd.fe_values.shape_grad(i, q_point))) *
1289 * sd.advection_directions[q_point] *
1290 * sd.fe_values.shape_grad(j, q_point)) *
1291 * sd.fe_values.JxW(q_point);
1293 * copy_data.cell_rhs(i) +=
1294 * (sd.fe_values.shape_value(i, q_point) +
1295 * delta * (sd.advection_directions[q_point] *
1296 * sd.fe_values.shape_grad(i, q_point))) *
1297 * sd.rhs_values[q_point] *
1298 * sd.fe_values.JxW(q_point);
1303 * Besides the cell terms which we have built up now, the bilinear
1304 * form of the present problem also contains terms on the boundary of
1305 * the domain. Therefore, we have to
check whether any of the faces of
1306 *
this cell are on the boundary of the domain, and
if so
assemble the
1307 * contributions of
this face as well. Of course, the bilinear form
1308 * only contains contributions from the <code>inflow</code> part of
1309 * the boundary, but to find out whether a certain part of a face of
1310 * the present cell is part of the inflow boundary, we have to have
1311 * information on the exact location of the quadrature points and on
1312 * the direction of flow at
this point; we obtain
this information
1314 * whether a quadrature
point is on the inflow boundary.
1317 *
for (
const auto &face : cell->face_iterators())
1318 * if (face->at_boundary())
1322 * Ok, this face of the present cell is on the boundary of the
1323 * domain. Just as for the usual
FEValues object which we have
1324 * used in previous examples and also above, we have to
1325 * reinitialize the
FEFaceValues object for the present face:
1328 * scratch_data.fe_face_values.
reinit(cell, face);
1332 * For the quadrature points at hand, we ask
for the
values of
1333 * the inflow function and
for the direction of flow:
1336 * scratch_data.boundary_values.value_list(
1337 * scratch_data.fe_face_values.get_quadrature_points(),
1338 * scratch_data.face_boundary_values);
1339 * scratch_data.advection_field.value_list(
1340 * scratch_data.fe_face_values.get_quadrature_points(),
1341 * scratch_data.face_advection_directions);
1345 * Now
loop over all quadrature points and see whether
this face is on
1346 * the inflow or outflow part of the boundary. The normal
1347 * vector points out of the cell: since the face is at
1348 * the boundary, the normal vector points out of the domain,
1349 * so
if the advection direction points into the domain, its
1350 *
scalar product with the normal vector must be
negative (to see why
1351 *
this is
true, consider the scalar product definition that uses a
1355 * for (unsigned
int q_point = 0; q_point < n_face_q_points; ++q_point)
1356 *
if (scratch_data.fe_face_values.normal_vector(q_point) *
1357 * scratch_data.face_advection_directions[q_point] <
1361 * If the face is part of the inflow boundary, then compute the
1362 * contributions of
this face to the global matrix and right
1363 * hand side,
using the values obtained from the
1364 *
FEFaceValues object and the formulae discussed in the
1368 * for (unsigned
int i = 0; i < dofs_per_cell; ++i)
1370 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1371 * copy_data.cell_matrix(i, j) -=
1372 * (scratch_data.face_advection_directions[q_point] *
1373 * scratch_data.fe_face_values.normal_vector(q_point) *
1374 * scratch_data.fe_face_values.shape_value(i, q_point) *
1375 * scratch_data.fe_face_values.shape_value(j, q_point) *
1376 * scratch_data.fe_face_values.JxW(q_point));
1378 * copy_data.cell_rhs(i) -=
1379 * (scratch_data.face_advection_directions[q_point] *
1380 * scratch_data.fe_face_values.normal_vector(q_point) *
1381 * scratch_data.face_boundary_values[q_point] *
1382 * scratch_data.fe_face_values.shape_value(i, q_point) *
1383 * scratch_data.fe_face_values.JxW(q_point));
1389 * The
final piece of information the
copy routine needs is the global
1390 * indices of the degrees of freedom on
this cell, so we
end by writing
1391 * them to the local array:
1394 * cell->get_dof_indices(copy_data.local_dof_indices);
1401 * The
second function we needed to write was the one that copies
1402 * the local contributions the previous function computed (and
1403 * put into the AssemblyCopyData
object) into the global
matrix and right
1404 * hand side vector objects. This is essentially what we
always had
1405 * as the last block of code when assembling something on every
1406 * cell. The following should therefore be pretty obvious:
1409 *
template <
int dim>
1411 * AdvectionProblem<dim>::copy_local_to_global(
const AssemblyCopyData ©_data)
1413 * hanging_node_constraints.distribute_local_to_global(
1414 * copy_data.cell_matrix,
1415 * copy_data.cell_rhs,
1416 * copy_data.local_dof_indices,
1423 * The next function is the linear solver routine. As the system is no longer
1425 * use the Conjugate Gradient method any more. Rather, we use a solver that
1426 * is more
general and does not rely on any special properties of the
1427 *
matrix: the GMRES method. GMRES, like the conjugate
gradient method,
1428 *
requires a decent preconditioner: we use a Jacobi preconditioner here,
1429 * which works well enough
for this problem.
1432 *
template <
int dim>
1433 *
void AdvectionProblem<dim>::solve()
1436 * system_rhs.size() / 10),
1437 * 1e-10 * system_rhs.l2_norm());
1440 * preconditioner.
initialize(system_matrix, 1.0);
1441 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
1443 * hanging_node_constraints.distribute(solution);
1445 * std::cout <<
" Iterations required for convergence: "
1446 * << solver_control.last_step() <<
'\n'
1447 * <<
" Norm of residual at convergence: "
1448 * << solver_control.last_value() <<
'\n';
1453 * The following function refines the grid according to the quantity
1454 * described in the introduction. The respective computations are made in
1455 * the class <code>GradientEstimation</code>.
1458 *
template <
int dim>
1459 *
void AdvectionProblem<dim>::refine_grid()
1463 * GradientEstimation::estimate(dof_handler,
1465 * estimated_error_per_cell);
1468 * estimated_error_per_cell,
1477 * This function is similar to the one in @ref step_6
"step-6", but since we use a higher
1478 * degree finite element we save the solution in a different
1479 * way. Visualization programs like VisIt and Paraview typically only
1480 * understand data that is associated with nodes: they cannot plot
1481 * fifth-degree basis
functions, which results in a very inaccurate picture
1482 * of the solution we computed. To get around
this we save multiple
1483 * <em>patches</em> per cell: in 2
d we save @f$8\times 8=64@f$ bilinear
1484 * `sub-cells
' to the VTU file for each cell, and in 3d we save
1485 * @f$8\times 8\times 8 = 512@f$. The end result is that the
1486 * visualization program will use a piecewise linear interpolation of the
1487 * cubic basis functions on a 3 times refined mesh:
1488 * This captures the solution detail and, with most
1489 * screen resolutions, looks smooth. We save the grid in a separate step
1490 * with no extra patches so that we have a visual representation of the cell
1495 * @note Version 9.1 of deal.II gained the ability to write higher degree
1496 * polynomials (i.e., write piecewise bicubic visualization data for our
1497 * piecewise bicubic solution) VTK and VTU output: however, not all recent
1498 * versions of ParaView and VisIt (as of 2018) can read this format, so we
1499 * use the older, more general (but less efficient) approach here.
1502 * template <int dim>
1503 * void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1507 * const std::string filename = "grid-" + std::to_string(cycle) + ".vtu";
1508 * std::ofstream output(filename);
1509 * grid_out.write_vtu(triangulation, output);
1510 * std::cout << " Grid written to " << filename << std::endl;
1514 * DataOut<dim> data_out;
1515 * data_out.attach_dof_handler(dof_handler);
1516 * data_out.add_data_vector(solution, "solution");
1517 * data_out.build_patches(8);
1521 * VTU output can be expensive, both to compute and to write to
1522 * disk. Here we ask ZLib, a compression library, to compress the data
1523 * in a way that maximizes throughput.
1526 * DataOutBase::VtkFlags vtk_flags;
1527 * vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
1528 * data_out.set_flags(vtk_flags);
1530 * const std::string filename = "solution-" + std::to_string(cycle) + ".vtu";
1531 * std::ofstream output(filename);
1532 * data_out.write_vtu(output);
1533 * std::cout << " Solution written to " << filename << std::endl;
1540 * ... as is the main loop (setup -- solve -- refine), aside from the number
1541 * of cycles and the initial grid:
1544 * template <int dim>
1545 * void AdvectionProblem<dim>::run()
1547 * for (unsigned int cycle = 0; cycle < 10; ++cycle)
1549 * std::cout << "Cycle " << cycle << ':
' << std::endl;
1553 * GridGenerator::hyper_cube(triangulation, -1, 1);
1554 * triangulation.refine_global(3);
1562 * std::cout << " Number of active cells: "
1563 * << triangulation.n_active_cells() << std::endl;
1567 * std::cout << " Number of degrees of freedom: "
1568 * << dof_handler.n_dofs() << std::endl;
1570 * assemble_system();
1572 * output_results(cycle);
1581 * <a name="step_9-GradientEstimationclassimplementation"></a>
1582 * <h3>GradientEstimation class implementation</h3>
1586 * Now for the implementation of the <code>GradientEstimation</code> class.
1587 * Let us start by defining constructors for the
1588 * <code>EstimateScratchData</code> class used by the
1589 * <code>estimate_cell()</code> function:
1592 * template <int dim>
1593 * GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1594 * const FiniteElement<dim> &fe,
1595 * const Vector<double> &solution,
1596 * Vector<float> &error_per_cell)
1597 * : fe_midpoint_value(fe,
1599 * update_values | update_quadrature_points)
1600 * , solution(solution)
1601 * , error_per_cell(error_per_cell)
1602 * , cell_midpoint_value(1)
1603 * , neighbor_midpoint_value(1)
1607 * We allocate a vector to hold iterators to all active neighbors of
1608 * a cell. We reserve the maximal number of active neighbors in order to
1609 * avoid later reallocations. Note how this maximal number of active
1610 * neighbors is computed here.
1613 * active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
1614 * GeometryInfo<dim>::max_children_per_face);
1618 * template <int dim>
1619 * GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1620 * const EstimateScratchData &scratch_data)
1621 * : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(),
1622 * scratch_data.fe_midpoint_value.get_quadrature(),
1623 * update_values | update_quadrature_points)
1624 * , solution(scratch_data.solution)
1625 * , error_per_cell(scratch_data.error_per_cell)
1626 * , cell_midpoint_value(1)
1627 * , neighbor_midpoint_value(1)
1633 * Next comes the implementation of the <code>GradientEstimation</code>
1634 * class. The first function does not much except for delegating work to the
1635 * other function, but there is a bit of setup at the top.
1639 * Before starting with the work, we check that the vector into
1640 * which the results are written has the right size, using the `Assert`
1641 * macro and the exception class we declared above. Programming
1642 * mistakes in which one forgets to size arguments correctly at the
1643 * calling site are quite common. Because the resulting damage from
1644 * not catching such errors is often subtle (e.g., corruption of
1645 * data somewhere in memory, or non-reproducible results), it is
1646 * well worth the effort to check for such things.
1649 * template <int dim>
1650 * void GradientEstimation::estimate(const DoFHandler<dim> &dof_handler,
1651 * const Vector<double> &solution,
1652 * Vector<float> &error_per_cell)
1655 * error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(),
1656 * ExcInvalidVectorLength(error_per_cell.size(),
1657 * dof_handler.get_triangulation().n_active_cells()));
1659 * WorkStream::run(dof_handler.begin_active(),
1660 * dof_handler.end(),
1661 * &GradientEstimation::template estimate_cell<dim>,
1662 * std::function<void(const EstimateCopyData &)>(),
1663 * EstimateScratchData<dim>(dof_handler.get_fe(),
1666 * EstimateCopyData());
1672 * Here comes the function that estimates the local error by computing the
1673 * finite difference approximation of the gradient. The function first
1674 * computes the list of active neighbors of the present cell and then
1675 * computes the quantities described in the introduction for each of
1676 * the neighbors. The reason for this order is that it is not a one-liner
1677 * to find a given neighbor with locally refined meshes. In principle, an
1678 * optimized implementation would find neighbors and the quantities
1679 * depending on them in one step, rather than first building a list of
1680 * neighbors and in a second step their contributions but we will gladly
1681 * leave this as an exercise. As discussed before, the worker function
1682 * passed to WorkStream::run works on "scratch" objects that keep all
1683 * temporary objects. This way, we do not need to create and initialize
1684 * objects that are expensive to initialize within the function that does
1685 * the work every time it is called for a given cell. Such an argument is
1686 * passed as the second argument. The third argument would be a "copy-data"
1687 * object (see @ref threads for more information) but we do not actually use
1688 * any of these here. Since WorkStream::run() insists on passing three
1689 * arguments, we declare this function with three arguments, but simply
1690 * ignore the last one.
1694 * (This is unsatisfactory from an aesthetic perspective. It can be avoided
1695 * by using an anonymous (lambda) function. If you allow, let us here show
1696 * how. First, assume that we had declared this function to only take two
1697 * arguments by omitting the unused last one. Now, WorkStream::run still
1698 * wants to call this function with three arguments, so we need to find a
1699 * way to "forget" the third argument in the call. Simply passing
1700 * WorkStream::run the pointer to the function as we do above will not do
1701 * this -- the compiler will complain that a function declared to have two
1702 * arguments is called with three arguments. However, we can do this by
1703 * passing the following as the third argument to WorkStream::run():
1704 * <div class=CodeFragmentInTutorialComment>
1706 * [](const typename DoFHandler<dim>::active_cell_iterator &cell,
1707 * EstimateScratchData<dim> & scratch_data,
1708 * EstimateCopyData &)
1710 * GradientEstimation::estimate_cell<dim>(cell, scratch_data);
1714 * This is not much better than the solution implemented below: either the
1715 * routine itself must take three arguments or it must be wrapped by
1716 * something that takes three arguments. We don't use
this since adding the
1717 * unused argument at the beginning is simpler.
1721 * Now
for the details:
1724 *
template <
int dim>
1725 *
void GradientEstimation::estimate_cell(
1727 * EstimateScratchData<dim> &scratch_data,
1728 *
const EstimateCopyData &)
1732 * We need space
for the tensor <code>Y</code>, which is the
sum of
1733 * outer products of the y-vectors.
1740 * First initialize the <code>
FEValues</code> object, as well as the
1741 * <code>Y</code> tensor:
1744 * scratch_data.fe_midpoint_value.
reinit(cell);
1748 * Now, before we go on, we
first compute a list of all active neighbors
1749 * of the present cell. We
do so by
first looping over all faces and see
1750 * whether the neighbor there is active, which would be the
case if it
1751 * is on the same
level as the present cell or one
level coarser (note
1752 * that a neighbor can only be once coarser than the present cell, as
1753 * we only allow a maximal difference of one refinement over a face in
1754 * deal.II). Alternatively, the neighbor could be on the same
level
1755 * and be further refined; then we have to find which of its children
1756 * are next to the present cell and select these (note that
if a child
1757 * of a neighbor of an active cell that is next to
this active cell,
1758 * needs necessarily be active itself, due to the one-refinement rule
1763 * Things are slightly different in one space dimension, as there the
1764 * one-refinement rule does not exist: neighboring active cells may
1765 * differ in as many refinement levels as they like. In
this case, the
1766 * computation becomes a little more difficult, but we will explain
1771 * Before starting the
loop over all neighbors of the present cell, we
1772 * have to clear the array storing the iterators to the active
1773 * neighbors, of course.
1776 * scratch_data.active_neighbors.clear();
1777 *
for (
const auto face_n : cell->face_indices())
1778 * if (!cell->at_boundary(face_n))
1782 * First define an abbreviation for the iterator to the face and
1786 * const auto face = cell->face(face_n);
1787 *
const auto neighbor = cell->neighbor(face_n);
1791 * Then
check whether the neighbor is active. If it is, then it
1792 * is on the same
level or one
level coarser (
if we are not in
1793 * 1d), and we are interested in it in any
case.
1796 *
if (neighbor->is_active())
1797 * scratch_data.active_neighbors.push_back(neighbor);
1802 * If the neighbor is not active, then
check its children.
1809 * To find the child of the neighbor which bounds to the
1810 * present cell, successively go to its right child
if
1811 * we are left of the present cell (n==0), or go to the
1812 * left child
if we are on the right (n==1), until we
1813 * find an active cell.
1816 *
auto neighbor_child = neighbor;
1817 *
while (neighbor_child->has_children())
1818 * neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0);
1822 * As
this used some non-trivial geometrical intuition,
1823 * we might want to
check whether we did it right,
1824 * i.e.,
check whether the neighbor of the cell we found
1825 * is indeed the cell we are presently working
1826 * on. Checks like
this are often useful and have
1827 * frequently uncovered errors both in algorithms like
1828 * the line above (where it is simple to involuntarily
1829 * exchange <code>n==1</code> for <code>n==0</code> or
1830 * the like) and in the library (the assumptions
1831 * underlying the algorithm above could either be wrong,
1832 * wrongly documented, or are violated due to an error
1833 * in the library). One could in principle remove such
1834 * checks after the program works
for some time, but it
1835 * might be a good things to leave it in anyway to
check
1836 *
for changes in the library or in the algorithm above.
1840 * Note that
if this check fails, then
this is certainly
1841 * an error that is irrecoverable and probably qualifies
1842 * as an
internal error. We therefore use a predefined
1843 * exception
class to throw here.
1846 *
Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell,
1847 * ExcInternalError());
1851 * If the
check succeeded, we push the active neighbor
1852 * we just found to the stack we keep:
1855 * scratch_data.active_neighbors.push_back(neighbor_child);
1860 * If we are not in 1
d, we collect all neighbor children
1861 * `behind
' the subfaces of the current face and move on:
1864 * for (unsigned int subface_n = 0; subface_n < face->n_children();
1866 * scratch_data.active_neighbors.push_back(
1867 * cell->neighbor_child_on_subface(face_n, subface_n));
1873 * OK, now that we have all the neighbors, lets start the computation
1874 * on each of them. First we do some preliminaries: find out about the
1875 * center of the present cell and the solution at this point. The
1876 * latter is obtained as a vector of function values at the quadrature
1877 * points, of which there are only one, of course. Likewise, the
1878 * position of the center is the position of the first (and only)
1879 * quadrature point in real space.
1882 * const Point<dim> this_center =
1883 * scratch_data.fe_midpoint_value.quadrature_point(0);
1885 * scratch_data.fe_midpoint_value.get_function_values(
1886 * scratch_data.solution, scratch_data.cell_midpoint_value);
1890 * Now loop over all active neighbors and collect the data we
1894 * Tensor<1, dim> projected_gradient;
1895 * for (const auto &neighbor : scratch_data.active_neighbors)
1899 * Then get the center of the neighbor cell and the value of the
1900 * finite element function at that point. Note that for this
1901 * information we have to reinitialize the <code>FEValues</code>
1902 * object for the neighbor cell.
1905 * scratch_data.fe_midpoint_value.reinit(neighbor);
1906 * const Point<dim> neighbor_center =
1907 * scratch_data.fe_midpoint_value.quadrature_point(0);
1909 * scratch_data.fe_midpoint_value.get_function_values(
1910 * scratch_data.solution, scratch_data.neighbor_midpoint_value);
1914 * Compute the vector <code>y</code> connecting the centers of the
1915 * two cells. Note that as opposed to the introduction, we denote
1916 * by <code>y</code> the normalized difference vector, as this is
1917 * the quantity used everywhere in the computations.
1920 * Tensor<1, dim> y = neighbor_center - this_center;
1921 * const double distance = y.norm();
1926 * Then add up the contribution of this cell to the Y matrix...
1929 * for (unsigned int i = 0; i < dim; ++i)
1930 * for (unsigned int j = 0; j < dim; ++j)
1931 * Y[i][j] += y[i] * y[j];
1935 * ... and update the sum of difference quotients:
1938 * projected_gradient += (scratch_data.neighbor_midpoint_value[0] -
1939 * scratch_data.cell_midpoint_value[0]) /
1945 * If now, after collecting all the information from the neighbors, we
1946 * can determine an approximation of the gradient for the present
1947 * cell, then we need to have passed over vectors <code>y</code> which
1948 * span the whole space, otherwise we would not have all components of
1949 * the gradient. This is indicated by the invertibility of the matrix.
1953 * If the matrix is not invertible, then the present
1954 * cell had an insufficient number of active neighbors. In contrast to
1955 * all previous cases (where we raised exceptions) this is, however,
1956 * not a programming error: it is a runtime error that can happen in
1957 * optimized mode even if it ran well in debug mode, so it is
1958 * reasonable to try to catch this error also in optimized mode. For
1959 * this case, there is the <code>AssertThrow</code> macro: it checks
1960 * the condition like the <code>Assert</code> macro, but not only in
1961 * debug mode; it then outputs an error message, but instead of
1962 * aborting the program as in the case of the <code>Assert</code>
1963 * macro, the exception is thrown using the <code>throw</code> command
1964 * of C++. This way, one has the possibility to catch this error and
1965 * take reasonable counter actions. One such measure would be to
1966 * refine the grid globally, as the case of insufficient directions
1967 * can not occur if every cell of the initial grid has been refined at
1971 * AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
1975 * If, on the other hand, the matrix is invertible, then invert it,
1976 * multiply the other quantity with it, and compute the estimated error
1977 * using this quantity and the correct powers of the mesh width:
1980 * const Tensor<2, dim> Y_inverse = invert(Y);
1982 * const Tensor<1, dim> gradient = Y_inverse * projected_gradient;
1986 * The last part of this function is the one where we write into
1987 * the element of the output vector what we have just
1988 * computed. The address of this vector has been stored in the
1989 * scratch data object, and all we have to do is know how to get
1990 * at the correct element inside this vector -- but we can ask the
1991 * cell we're on the how-manyth active cell it is
for this:
1994 * scratch_data.error_per_cell(cell->active_cell_index()) =
2003 * <a name=
"step_9-Mainfunction"></a>
2004 * <h3>Main function</h3>
2008 * The <code>main</code> function is similar to the previous examples. The
2009 * primary difference is that we use
MultithreadInfo to set the maximum
2010 * number of threads (see the documentation topic @ref threads
2011 *
"Parallel computing with multiple processors accessing shared memory"
2012 *
for more information). The number of threads used is the minimum of the
2013 * environment variable DEAL_II_NUM_THREADS and the parameter of
2014 * <code>set_thread_limit</code>. If no
value is given to
2015 * <code>set_thread_limit</code>, the
default value from the Intel Threading
2016 * Building Blocks (TBB) library is used. If the call to
2017 * <code>set_thread_limit</code> is omitted, the number of threads will be
2018 * chosen by TBB independently of DEAL_II_NUM_THREADS.
2023 *
using namespace dealii;
2028 * Step9::AdvectionProblem<2> advection_problem_2d;
2029 * advection_problem_2d.run();
2031 *
catch (std::exception &exc)
2033 * std::cerr << std::endl
2035 * <<
"----------------------------------------------------"
2037 * std::cerr <<
"Exception on processing: " << std::endl
2038 * << exc.what() << std::endl
2039 * <<
"Aborting!" << std::endl
2040 * <<
"----------------------------------------------------"
2046 * std::cerr << std::endl
2048 * <<
"----------------------------------------------------"
2050 * std::cerr <<
"Unknown exception!" << std::endl
2051 * <<
"Aborting!" << std::endl
2052 * <<
"----------------------------------------------------"
2060<a name=
"step_9-Results"></a><h1>Results</h1>
2064The results of
this program are not particularly spectacular. They
2065consist of the console output, some grid files, and the solution on
2066each of these grids. First
for the console output:
2069 Number of active cells: 64
2070 Number of degrees of freedom: 1681
2071 Iterations required
for convergence: 338
2072 Norm of residual at convergence: 2.9177e-11
2073 Grid written to grid-0.
vtu
2074 Solution written to solution-0.
vtu
2076 Number of active cells: 121
2077 Number of degrees of freedom: 3436
2078 Iterations required
for convergence: 448
2079 Norm of residual at convergence: 2.40051e-11
2080 Grid written to grid-1.
vtu
2081 Solution written to solution-1.
vtu
2083 Number of active cells: 238
2084 Number of degrees of freedom: 6487
2085 Iterations required
for convergence: 535
2086 Norm of residual at convergence: 2.06897e-11
2087 Grid written to grid-2.
vtu
2088 Solution written to solution-2.
vtu
2090 Number of active cells: 481
2091 Number of degrees of freedom: 13510
2092 Iterations required
for convergence: 669
2093 Norm of residual at convergence: 1.41502e-11
2094 Grid written to grid-3.
vtu
2095 Solution written to solution-3.
vtu
2097 Number of active cells: 958
2098 Number of degrees of freedom: 26137
2099 Iterations required
for convergence: 1039
2100 Norm of residual at convergence: 1.27042e-11
2101 Grid written to grid-4.
vtu
2102 Solution written to solution-4.
vtu
2104 Number of active cells: 1906
2105 Number of degrees of freedom: 52832
2106 Iterations required
for convergence: 1345
2107 Norm of residual at convergence: 1.00694e-11
2108 Grid written to grid-5.
vtu
2109 Solution written to solution-5.
vtu
2111 Number of active cells: 3829
2112 Number of degrees of freedom: 104339
2113 Iterations required
for convergence: 1976
2114 Norm of residual at convergence: 7.43452e-12
2115 Grid written to grid-6.
vtu
2116 Solution written to solution-6.
vtu
2118 Number of active cells: 7414
2119 Number of degrees of freedom: 201946
2120 Iterations required
for convergence: 2256
2121 Norm of residual at convergence: 7.88403e-12
2122 Grid written to grid-7.
vtu
2123 Solution written to solution-7.
vtu
2125 Number of active cells: 14413
2126 Number of degrees of freedom: 389558
2127 Iterations required
for convergence: 2968
2128 Norm of residual at convergence: 6.72725e-12
2129 Grid written to grid-8.
vtu
2130 Solution written to solution-8.
vtu
2132 Number of active cells: 28141
2133 Number of degrees of freedom: 750187
2134 Iterations required
for convergence: 3885
2135 Norm of residual at convergence: 5.86246e-12
2136 Grid written to grid-9.
vtu
2137 Solution written to solution-9.
vtu
2140Quite a number of cells are used on the finest
level to resolve the features of
2141the solution. Here are the fourth and tenth grids:
2142<div
class=
"twocolumn" style=
"width: 80%">
2144 <img src=
"https://www.dealii.org/images/steps/developer/step-9-grid-3.png"
2145 alt=
"Fourth grid in the refinement cycle, showing some adaptivity to features."
2146 width=
"400" height=
"400">
2149 <img src=
"https://www.dealii.org/images/steps/developer/step-9-grid-9.png"
2150 alt=
"Tenth grid in the refinement cycle, showing that the waves are fully captured."
2151 width=
"400" height=
"400">
2154and the fourth and tenth solutions:
2155<div
class=
"twocolumn" style=
"width: 80%">
2157 <img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-3.png"
2158 alt=
"Fourth solution, showing that we resolve most features but some
2159 are sill unresolved and appear blury."
2160 width=
"400" height=
"400">
2163 <img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-9.png"
2164 alt=
"Tenth solution, showing a fully resolved flow."
2165 width=
"400" height=
"400">
2168and both the grid and solution zoomed in:
2169<div
class=
"twocolumn" style=
"width: 80%">
2171 <img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-3-zoom.png"
2172 alt=
"Detail of the fourth solution, showing that we resolve most
2173 features but some are sill unresolved and appear blury. In particular,
2174 the larger cells need to be refined."
2175 width=
"400" height=
"400">
2178 <img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-9-zoom.png"
2179 alt=
"Detail of the tenth solution, showing that we needed a lot more
2180 cells than were present in the fourth solution."
2181 width=
"400" height=
"400">
2185The solution is created by that part that is transported along the wiggly
2186advection field from the left and lower boundaries to the top right, and the
2187part that is created by the source in the lower left corner, and the results of
2188which are also transported along. The grid shown above is well-adapted to
2189resolve these features. The comparison between plots shows that, even though we
2190are
using a high-order approximation, we still need adaptive mesh refinement to
2191fully resolve the wiggles.
2194<a name=
"step_9-PlainProg"></a>
2195<h1> The plain program</h1>
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 ¶meters=AdditionalData())
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
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())
@ 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.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
T 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)
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 > &)