1070 *
sparsity_pattern.copy_from(
dsp);
1072 *
system_matrix.reinit(sparsity_pattern);
1074 *
solution.reinit(dof_handler.n_dofs());
1103 *
dof_handler.
end(),
1134 *
QGauss<dim>(fe.degree + 1),
1137 *
, fe_face_values(fe,
1138 *
QGauss<dim - 1>(fe.degree + 1),
1152 *
: fe_values(scratch_data.fe_values.get_fe(),
1153 *
scratch_data.fe_values.get_quadrature(),
1156 *
, fe_face_values(scratch_data.fe_face_values.get_fe(),
1157 *
scratch_data.fe_face_values.get_quadrature(),
1220 *
const unsigned
int dofs_per_cell = fe.n_dofs_per_cell();
1221 *
const unsigned int n_q_points =
1222 *
scratch_data.fe_values.get_quadrature().size();
1224 *
scratch_data.fe_face_values.get_quadrature().size();
1231 *
copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1232 *
copy_data.cell_rhs.reinit(dofs_per_cell);
1240 *
copy_data.local_dof_indices.resize(dofs_per_cell);
1247 *
scratch_data.fe_values.reinit(cell);
1252 * at
the quadrature points...
1255 *
scratch_data.advection_field.value_list(
1256 *
scratch_data.fe_values.get_quadrature_points(),
1257 *
scratch_data.advection_directions);
1258 *
scratch_data.right_hand_side.value_list(
1259 *
scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
1267 *
const double delta = 0.1 * cell->diameter();
1276 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1284 *
const auto &
sd = scratch_data;
1285 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1286 *
copy_data.cell_matrix(i,
j) +=
1287 *
((
sd.fe_values.shape_value(i,
q_point) +
1288 *
delta * (
sd.advection_directions[
q_point] *
1289 *
sd.fe_values.shape_grad(i,
q_point))) *
1294 *
copy_data.cell_rhs(i) +=
1295 *
(
sd.fe_values.shape_value(i,
q_point) +
1296 *
delta * (
sd.advection_directions[
q_point] *
1297 *
sd.fe_values.shape_grad(i,
q_point))) *
1318 *
for (
const auto &face : cell->face_iterators())
1319 *
if (face->at_boundary())
1329 *
scratch_data.fe_face_values.
reinit(cell, face);
1337 *
scratch_data.boundary_values.value_list(
1338 *
scratch_data.fe_face_values.get_quadrature_points(),
1339 *
scratch_data.face_boundary_values);
1340 *
scratch_data.advection_field.value_list(
1341 *
scratch_data.fe_face_values.get_quadrature_points(),
1342 *
scratch_data.face_advection_directions);
1357 *
if (scratch_data.fe_face_values.normal_vector(
q_point) *
1358 *
scratch_data.face_advection_directions[
q_point] <
1369 *
for (unsigned
int i = 0; i < dofs_per_cell; ++i)
1371 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1372 *
copy_data.cell_matrix(i,
j) -=
1373 *
(scratch_data.face_advection_directions[
q_point] *
1374 *
scratch_data.fe_face_values.normal_vector(
q_point) *
1375 *
scratch_data.fe_face_values.shape_value(i,
q_point) *
1376 *
scratch_data.fe_face_values.shape_value(
j,
q_point) *
1377 *
scratch_data.fe_face_values.JxW(
q_point));
1379 *
copy_data.cell_rhs(i) -=
1380 *
(scratch_data.face_advection_directions[
q_point] *
1381 *
scratch_data.fe_face_values.normal_vector(
q_point) *
1382 *
scratch_data.face_boundary_values[
q_point] *
1383 *
scratch_data.fe_face_values.shape_value(i,
q_point) *
1384 *
scratch_data.fe_face_values.JxW(
q_point));
1395 *
cell->get_dof_indices(copy_data.local_dof_indices);
1410 *
template <
int dim>
1415 *
copy_data.cell_matrix,
1416 *
copy_data.cell_rhs,
1417 *
copy_data.local_dof_indices,
1429 *
requires a
decent preconditioner:
we use a Jacobi preconditioner
here,
1433 *
template <
int dim>
1441 *
preconditioner.initialize(system_matrix, 1.0);
1442 *
solver.solve(system_matrix, solution,
system_rhs, preconditioner);
1446 *
std::cout <<
" Iterations required for convergence: "
1447 *
<< solver_control.last_step() <<
'\n'
1448 *
<<
" Norm of residual at convergence: "
1449 *
<< solver_control.last_value() <<
'\n';
1459 *
template <
int dim>
1464 *
GradientEstimation::estimate(dof_handler,
1485 * `
sub-cells
' to the VTU file for each cell, and in 3d we save
1486 * @f$8\times 8\times 8 = 512@f$. The end result is that the
1487 * visualization program will use a piecewise linear interpolation of the
1488 * cubic basis functions on a 3 times refined mesh:
1489 * This captures the solution detail and, with most
1490 * screen resolutions, looks smooth. We save the grid in a separate step
1491 * with no extra patches so that we have a visual representation of the cell
1496 * @note Version 9.1 of deal.II gained the ability to write higher degree
1497 * polynomials (i.e., write piecewise bicubic visualization data for our
1498 * piecewise bicubic solution) VTK and VTU output: however, not all recent
1499 * versions of ParaView and VisIt (as of 2018) can read this format, so we
1500 * use the older, more general (but less efficient) approach here.
1503 * template <int dim>
1504 * void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1508 * const std::string filename = "grid-" + std::to_string(cycle) + ".vtu";
1509 * std::ofstream output(filename);
1510 * grid_out.write_vtu(triangulation, output);
1511 * std::cout << " Grid written to " << filename << std::endl;
1515 * DataOut<dim> data_out;
1516 * data_out.attach_dof_handler(dof_handler);
1517 * data_out.add_data_vector(solution, "solution");
1518 * data_out.build_patches(8);
1522 * VTU output can be expensive, both to compute and to write to
1523 * disk. Here we ask ZLib, a compression library, to compress the data
1524 * in a way that maximizes throughput.
1527 * DataOutBase::VtkFlags vtk_flags;
1528 * vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
1529 * data_out.set_flags(vtk_flags);
1531 * const std::string filename = "solution-" + std::to_string(cycle) + ".vtu";
1532 * std::ofstream output(filename);
1533 * data_out.write_vtu(output);
1534 * std::cout << " Solution written to " << filename << std::endl;
1541 * ... as is the main loop (setup -- solve -- refine), aside from the number
1542 * of cycles and the initial grid:
1545 * template <int dim>
1546 * void AdvectionProblem<dim>::run()
1548 * for (unsigned int cycle = 0; cycle < 10; ++cycle)
1550 * std::cout << "Cycle " << cycle << ':
' << std::endl;
1554 * GridGenerator::hyper_cube(triangulation, -1, 1);
1555 * triangulation.refine_global(3);
1563 * std::cout << " Number of active cells: "
1564 * << triangulation.n_active_cells() << std::endl;
1568 * std::cout << " Number of degrees of freedom: "
1569 * << dof_handler.n_dofs() << std::endl;
1571 * assemble_system();
1573 * output_results(cycle);
1582 * <a name="step_9-GradientEstimationclassimplementation"></a>
1583 * <h3>GradientEstimation class implementation</h3>
1587 * Now for the implementation of the <code>GradientEstimation</code> class.
1588 * Let us start by defining constructors for the
1589 * <code>EstimateScratchData</code> class used by the
1590 * <code>estimate_cell()</code> function:
1593 * template <int dim>
1594 * GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1595 * const FiniteElement<dim> &fe,
1596 * const Vector<double> &solution,
1597 * Vector<float> &error_per_cell)
1598 * : fe_midpoint_value(fe,
1600 * update_values | update_quadrature_points)
1601 * , solution(solution)
1602 * , error_per_cell(error_per_cell)
1603 * , cell_midpoint_value(1)
1604 * , neighbor_midpoint_value(1)
1608 * We allocate a vector to hold iterators to all active neighbors of
1609 * a cell. We reserve the maximal number of active neighbors in order to
1610 * avoid later reallocations. Note how this maximal number of active
1611 * neighbors is computed here.
1614 * active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
1615 * GeometryInfo<dim>::max_children_per_face);
1619 * template <int dim>
1620 * GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1621 * const EstimateScratchData &scratch_data)
1622 * : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(),
1623 * scratch_data.fe_midpoint_value.get_quadrature(),
1624 * update_values | update_quadrature_points)
1625 * , solution(scratch_data.solution)
1626 * , error_per_cell(scratch_data.error_per_cell)
1627 * , cell_midpoint_value(1)
1628 * , neighbor_midpoint_value(1)
1634 * Next comes the implementation of the <code>GradientEstimation</code>
1635 * class. The first function does not much except for delegating work to the
1636 * other function, but there is a bit of setup at the top.
1640 * Before starting with the work, we check that the vector into
1641 * which the results are written has the right size, using the `Assert`
1642 * macro and the exception class we declared above. Programming
1643 * mistakes in which one forgets to size arguments correctly at the
1644 * calling site are quite common. Because the resulting damage from
1645 * not catching such errors is often subtle (e.g., corruption of
1646 * data somewhere in memory, or non-reproducible results), it is
1647 * well worth the effort to check for such things.
1650 * template <int dim>
1651 * void GradientEstimation::estimate(const DoFHandler<dim> &dof_handler,
1652 * const Vector<double> &solution,
1653 * Vector<float> &error_per_cell)
1656 * error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(),
1657 * ExcInvalidVectorLength(error_per_cell.size(),
1658 * dof_handler.get_triangulation().n_active_cells()));
1660 * WorkStream::run(dof_handler.begin_active(),
1661 * dof_handler.end(),
1662 * &GradientEstimation::template estimate_cell<dim>,
1663 * std::function<void(const EstimateCopyData &)>(),
1664 * EstimateScratchData<dim>(dof_handler.get_fe(),
1667 * EstimateCopyData());
1673 * Here comes the function that estimates the local error by computing the
1674 * finite difference approximation of the gradient. The function first
1675 * computes the list of active neighbors of the present cell and then
1676 * computes the quantities described in the introduction for each of
1677 * the neighbors. The reason for this order is that it is not a one-liner
1678 * to find a given neighbor with locally refined meshes. In principle, an
1679 * optimized implementation would find neighbors and the quantities
1680 * depending on them in one step, rather than first building a list of
1681 * neighbors and in a second step their contributions but we will gladly
1682 * leave this as an exercise. As discussed before, the worker function
1683 * passed to WorkStream::run works on "scratch" objects that keep all
1684 * temporary objects. This way, we do not need to create and initialize
1685 * objects that are expensive to initialize within the function that does
1686 * the work every time it is called for a given cell. Such an argument is
1687 * passed as the second argument. The third argument would be a "copy-data"
1688 * object (see @ref threads for more information) but we do not actually use
1689 * any of these here. Since WorkStream::run() insists on passing three
1690 * arguments, we declare this function with three arguments, but simply
1691 * ignore the last one.
1695 * (This is unsatisfactory from an aesthetic perspective. It can be avoided
1696 * by using an anonymous (lambda) function. If you allow, let us here show
1697 * how. First, assume that we had declared this function to only take two
1698 * arguments by omitting the unused last one. Now, WorkStream::run still
1699 * wants to call this function with three arguments, so we need to find a
1700 * way to "forget" the third argument in the call. Simply passing
1701 * WorkStream::run the pointer to the function as we do above will not do
1702 * this -- the compiler will complain that a function declared to have two
1703 * arguments is called with three arguments. However, we can do this by
1704 * passing the following as the third argument to WorkStream::run():
1705 * <div class=CodeFragmentInTutorialComment>
1707 * [](const typename DoFHandler<dim>::active_cell_iterator &cell,
1708 * EstimateScratchData<dim> & scratch_data,
1709 * EstimateCopyData &)
1711 * GradientEstimation::estimate_cell<dim>(cell, scratch_data);
1715 * This is not much better than the solution implemented below: either the
1716 * routine itself must take three arguments or it must be wrapped by
1725 *
template <
int dim>
1745 *
scratch_data.fe_midpoint_value.reinit(cell);
1777 *
scratch_data.active_neighbors.
clear();
1778 *
for (
const auto face_n : cell->face_indices())
1788 *
const auto neighbor = cell->neighbor(
face_n);
1797 *
if (neighbor->is_active())
1798 *
scratch_data.active_neighbors.push_back(neighbor);
1862 * `
behind' the subfaces of the current face and move on:
1865 * for (unsigned int subface_n = 0; subface_n < face->n_children();
1867 * scratch_data.active_neighbors.push_back(
1868 * cell->neighbor_child_on_subface(face_n, subface_n));
1874 * OK, now that we have all the neighbors, lets start the computation
1875 * on each of them. First we do some preliminaries: find out about the
1876 * center of the present cell and the solution at this point. The
1877 * latter is obtained as a vector of function values at the quadrature
1878 * points, of which there are only one, of course. Likewise, the
1879 * position of the center is the position of the first (and only)
1880 * quadrature point in real space.
1883 * const Point<dim> this_center =
1884 * scratch_data.fe_midpoint_value.quadrature_point(0);
1886 * scratch_data.fe_midpoint_value.get_function_values(
1887 * scratch_data.solution, scratch_data.cell_midpoint_value);
1891 * Now loop over all active neighbors and collect the data we
1895 * Tensor<1, dim> projected_gradient;
1896 * for (const auto &neighbor : scratch_data.active_neighbors)
1900 * Then get the center of the neighbor cell and the value of the
1901 * finite element function at that point. Note that for this
1902 * information we have to reinitialize the <code>FEValues</code>
1903 * object for the neighbor cell.
1906 * scratch_data.fe_midpoint_value.reinit(neighbor);
1907 * const Point<dim> neighbor_center =
1908 * scratch_data.fe_midpoint_value.quadrature_point(0);
1910 * scratch_data.fe_midpoint_value.get_function_values(
1911 * scratch_data.solution, scratch_data.neighbor_midpoint_value);
1915 * Compute the vector <code>y</code> connecting the centers of the
1916 * two cells. Note that as opposed to the introduction, we denote
1917 * by <code>y</code> the normalized difference vector, as this is
1918 * the quantity used everywhere in the computations.
1921 * Tensor<1, dim> y = neighbor_center - this_center;
1922 * const double distance = y.norm();
1927 * Then add up the contribution of this cell to the Y matrix...
1930 * for (unsigned int i = 0; i < dim; ++i)
1931 * for (unsigned int j = 0; j < dim; ++j)
1932 * Y[i][j] += y[i] * y[j];
1936 * ... and update the sum of difference quotients:
1939 * projected_gradient += (scratch_data.neighbor_midpoint_value[0] -
1940 * scratch_data.cell_midpoint_value[0]) /
1946 * If now, after collecting all the information from the neighbors, we
1947 * can determine an approximation of the gradient for the present
1948 * cell, then we need to have passed over vectors <code>y</code> which
1949 * span the whole space, otherwise we would not have all components of
1950 * the gradient. This is indicated by the invertibility of the matrix.
1954 * If the matrix is not invertible, then the present
1955 * cell had an insufficient number of active neighbors. In contrast to
1956 * all previous cases (where we raised exceptions) this is, however,
1957 * not a programming error: it is a runtime error that can happen in
1958 * optimized mode even if it ran well in debug mode, so it is
1959 * reasonable to try to catch this error also in optimized mode. For
1960 * this case, there is the <code>AssertThrow</code> macro: it checks
1961 * the condition like the <code>Assert</code> macro, but not only in
1962 * debug mode; it then outputs an error message, but instead of
1963 * aborting the program as in the case of the <code>Assert</code>
1964 * macro, the exception is thrown using the <code>throw</code> command
1965 * of C++. This way, one has the possibility to catch this error and
1966 * take reasonable counter actions. One such measure would be to
1967 * refine the grid globally, as the case of insufficient directions
1968 * can not occur if every cell of the initial grid has been refined at
1972 * AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
1976 * If, on the other hand, the matrix is invertible, then invert it,
1977 * multiply the other quantity with it, and compute the estimated error
1978 * using this quantity and the correct powers of the mesh width:
1981 * const Tensor<2, dim> Y_inverse = invert(Y);
1983 * const Tensor<1, dim> gradient = Y_inverse * projected_gradient;
1987 * The last part of this function is the one where we write into
1988 * the element of the output vector what we have just
1989 * computed. The address of this vector has been stored in the
1990 * scratch data object, and all we have to do is know how to get
1991 * at the correct element inside this vector -- but we can ask the
1995 *
scratch_data.error_per_cell(cell->active_cell_index()) =
2004 * <a name=
"step_9-Mainfunction"></a>
2012 *
"Parallel computing with multiple processors accessing shared memory"
2024 *
using namespace dealii;
2032 *
catch (std::exception &exc)
2034 *
std::cerr << std::endl
2036 *
<<
"----------------------------------------------------"
2038 *
std::cerr <<
"Exception on processing: " << std::endl
2039 *
<< exc.what() << std::endl
2040 *
<<
"Aborting!" << std::endl
2041 *
<<
"----------------------------------------------------"
2047 *
std::cerr << std::endl
2049 *
<<
"----------------------------------------------------"
2051 *
std::cerr <<
"Unknown exception!" << std::endl
2052 *
<<
"Aborting!" << std::endl
2053 *
<<
"----------------------------------------------------"
2143<
div class=
"twocolumn" style=
"width: 80%">
2145 <
img src=
"https://www.dealii.org/images/steps/developer/step-9-grid-3.png"
2146 alt=
"Fourth grid in the refinement cycle, showing some adaptivity to features."
2147 width=
"400" height=
"400">
2150 <
img src=
"https://www.dealii.org/images/steps/developer/step-9-grid-9.png"
2151 alt=
"Tenth grid in the refinement cycle, showing that the waves are fully captured."
2152 width=
"400" height=
"400">
2156<
div class=
"twocolumn" style=
"width: 80%">
2158 <
img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-3.png"
2159 alt=
"Fourth solution, showing that we resolve most features but some
2160 are sill unresolved and appear blury."
2161 width=
"400" height=
"400">
2164 <
img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-9.png"
2165 alt=
"Tenth solution, showing a fully resolved flow."
2166 width=
"400" height=
"400">
2170<
div class=
"twocolumn" style=
"width: 80%">
2172 <
img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-3-zoom.png"
2173 alt=
"Detail of the fourth solution, showing that we resolve most
2174 features but some are sill unresolved and appear blury. In particular,
2175 the larger cells need to be refined."
2176 width=
"400" height=
"400">
2179 <
img src=
"https://www.dealii.org/images/steps/developer/step-9-solution-9-zoom.png"
2180 alt=
"Detail of the tenth solution, showing that we needed a lot more
2181 cells than were present in the fourth solution."
2182 width=
"400" height=
"400">
2195<a name=
"step_9-PlainProg"></a>
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)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
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.
std::vector< index_type > data
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.
constexpr types::blas_int one
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)
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 > &)