1083 *
const unsigned int )
const
1090 *
template <
int dim>
1093 *
const unsigned int )
const
1095 *
return -(
alpha * p[0] * p[1] * p[1] / 2 +
beta * p[0] -
1096 *
alpha * p[0] * p[0] * p[0] / 6);
1101 *
template <
int dim>
1110 *
alpha * p[0] * p[0] * p[0] / 6);
1118 * <a name=
"step_20-Theinversepermeabilitytensor"></a>
1155 *
template <
int dim>
1164 *
value_list(
const std::vector<
Point<dim>> &points,
1191 * (
reads) variable
but doesn't actually do anything: That's what
1202 *
template <
int dim>
1209 *
for (
auto &value :
values)
1219 * <a name=
"step_20-MixedLaplaceProblemclassimplementation"></a>
1225 * <a name=
"step_20-MixedLaplaceProblemMixedLaplaceProblem"></a>
1263 *
template <
int dim>
1275 * <a name=
"step_20-MixedLaplaceProblemmake_grid_and_dofs"></a>
1284 *
template <
int dim>
1290 *
dof_handler.distribute_dofs(fe);
1354 *
const unsigned int n_u = dofs_per_component[0],
1355 *
n_p = dofs_per_component[dim];
1357 *
std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
1361 *
<<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1362 *
<<
" (" <<
n_u <<
'+' <<
n_p <<
')' << std::endl;
1372 * <code>n_u</code> and <code>n_p</code>, which hold the number of velocity
1373 * and pressure variables.
1376 * const std::vector<types::global_dof_index> block_sizes = {n_u, n_p};
1377 * BlockDynamicSparsityPattern dsp(block_sizes, block_sizes);
1378 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
1382 * We use the compressed block sparsity pattern in the same way as the
1383 * non-block version to create the sparsity pattern and then the system
1387 * sparsity_pattern.copy_from(dsp);
1388 * system_matrix.reinit(sparsity_pattern);
1392 * Then we have to resize the solution and right hand side vectors in
1393 * exactly the same way as the block compressed sparsity pattern:
1396 * solution.reinit(block_sizes);
1397 * system_rhs.reinit(block_sizes);
1404 * <a name="step_20-MixedLaplaceProblemassemble_system"></a>
1405 * <h4>MixedLaplaceProblem::assemble_system</h4>
1409 * Similarly, the function that assembles the linear system has mostly been
1410 * discussed already in the introduction to this example. At its top, what
1411 * happens are all the usual steps, with the addition that we do not only
1412 * allocate quadrature and <code>FEValues</code> objects for the cell terms,
1413 * but also for face terms. After that, we define the usual abbreviations
1414 * for variables, and the allocate space for the local matrix and right hand
1415 * side contributions, and the array that holds the global numbers of the
1416 * degrees of freedom local to the present cell.
1419 * template <int dim>
1420 * void MixedLaplaceProblem<dim>::assemble_system()
1422 * const QGauss<dim> quadrature_formula(degree + 2);
1423 * const QGauss<dim - 1> face_quadrature_formula(degree + 2);
1425 * FEValues<dim> fe_values(fe,
1426 * quadrature_formula,
1427 * update_values | update_gradients |
1428 * update_quadrature_points | update_JxW_values);
1429 * FEFaceValues<dim> fe_face_values(fe,
1430 * face_quadrature_formula,
1431 * update_values | update_normal_vectors |
1432 * update_quadrature_points |
1433 * update_JxW_values);
1435 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1436 * const unsigned int n_q_points = quadrature_formula.size();
1437 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1439 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1440 * Vector<double> local_rhs(dofs_per_cell);
1442 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1446 * The next step is to declare objects that represent the source term,
1447 * pressure boundary value, and coefficient in the equation. In addition
1448 * to these objects that represent continuous functions, we also need
1449 * arrays to hold their values at the quadrature points of individual
1450 * cells (or faces, for the boundary values). Note that in the case of the
1451 * coefficient, the array has to be one of matrices.
1454 * const PrescribedSolution::RightHandSide<dim> right_hand_side;
1455 * const PrescribedSolution::PressureBoundaryValues<dim>
1456 * pressure_boundary_values;
1457 * const PrescribedSolution::KInverse<dim> k_inverse;
1459 * std::vector<double> rhs_values(n_q_points);
1460 * std::vector<double> boundary_values(n_face_q_points);
1461 * std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1465 * Finally, we need a couple of extractors that we will use to get at the
1466 * velocity and pressure components of vector-valued shape
1467 * functions. Their function and use is described in detail in the @ref
1468 * vector_valued report. Essentially, we will use them as subscripts on
1469 * the FEValues objects below: the FEValues object describes all vector
1470 * components of shape functions, while after subscription, it will only
1471 * refer to the velocities (a set of <code>dim</code> components starting
1472 * at component zero) or the pressure (a scalar component located at
1473 * position <code>dim</code>):
1476 * const FEValuesExtractors::Vector velocities(0);
1477 * const FEValuesExtractors::Scalar pressure(dim);
1481 * With all this in place, we can go on with the loop over all cells. The
1482 * body of this loop has been discussed in the introduction, and will not
1483 * be commented any further here:
1486 * for (const auto &cell : dof_handler.active_cell_iterators())
1488 * fe_values.reinit(cell);
1493 * right_hand_side.value_list(fe_values.get_quadrature_points(),
1495 * k_inverse.value_list(fe_values.get_quadrature_points(),
1496 * k_inverse_values);
1498 * for (unsigned int q = 0; q < n_q_points; ++q)
1499 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1501 * const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
1502 * const double div_phi_i_u = fe_values[velocities].divergence(i, q);
1503 * const double phi_i_p = fe_values[pressure].value(i, q);
1505 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1507 * const Tensor<1, dim> phi_j_u =
1508 * fe_values[velocities].value(j, q);
1509 * const double div_phi_j_u =
1510 * fe_values[velocities].divergence(j, q);
1511 * const double phi_j_p = fe_values[pressure].value(j, q);
1513 * local_matrix(i, j) +=
1514 * (phi_i_u * k_inverse_values[q] * phi_j_u
1515 * - phi_i_p * div_phi_j_u
1516 * - div_phi_i_u * phi_j_p)
1517 * * fe_values.JxW(q);
1520 * local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
1523 * for (const auto &face : cell->face_iterators())
1524 * if (face->at_boundary())
1526 * fe_face_values.reinit(cell, face);
1528 * pressure_boundary_values.value_list(
1529 * fe_face_values.get_quadrature_points(), boundary_values);
1531 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1532 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1533 * local_rhs(i) += -(fe_face_values[velocities].value(i, q) *
1534 * fe_face_values.normal_vector(q) *
1535 * boundary_values[q] *
1536 * fe_face_values.JxW(q));
1541 * The final step in the loop over all cells is to transfer local
1542 * contributions into the global matrix and right hand side
1543 * vector. Note that we use exactly the same interface as in previous
1544 * examples, although we now use block matrices and vectors instead of
1545 * the regular ones. In other words, to the outside world, block
1546 * objects have the same interface as matrices and vectors, but they
1547 * additionally allow to access individual blocks.
1550 * cell->get_dof_indices(local_dof_indices);
1551 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1552 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1553 * system_matrix.add(local_dof_indices[i],
1554 * local_dof_indices[j],
1555 * local_matrix(i, j));
1556 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1557 * system_rhs(local_dof_indices[i]) += local_rhs(i);
1565 * <a name="step_20-Implementationoflinearsolversandpreconditioners"></a>
1566 * <h3>Implementation of linear solvers and preconditioners</h3>
1570 * The linear solvers and preconditioners we use in this example have
1571 * been discussed in significant detail already in the introduction. We
1572 * will therefore not discuss the rationale for our approach here any
1573 * more, but rather only comment on some remaining implementational
1579 * <a name="step_20-MixedLaplacesolve"></a>
1580 * <h4>MixedLaplace::solve</h4>
1584 * As already outlined in the introduction, the solve function consists
1585 * essentially of two steps. First, we have to form the first equation
1586 * involving the Schur complement and solve for the pressure (component 1
1587 * of the solution). Then, we can reconstruct the velocities from the
1588 * second equation (component 0 of the solution).
1591 * template <int dim>
1592 * void MixedLaplaceProblem<dim>::solve()
1596 * As a first step we declare references to all block components of the
1597 * matrix, the right hand side and the solution vector that we will
1601 * const auto &M = system_matrix.block(0, 0);
1602 * const auto &B = system_matrix.block(0, 1);
1604 * const auto &F = system_rhs.block(0);
1605 * const auto &G = system_rhs.block(1);
1607 * auto &U = solution.block(0);
1608 * auto &P = solution.block(1);
1612 * Then, we will create corresponding LinearOperator objects and create
1613 * the <code>op_M_inv</code> operator:
1616 * const auto op_M = linear_operator(M);
1617 * const auto op_B = linear_operator(B);
1619 * ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10);
1620 * SolverCG<Vector<double>> solver_M(reduction_control_M);
1621 * PreconditionJacobi<SparseMatrix<double>> preconditioner_M;
1623 * preconditioner_M.initialize(M);
1625 * const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
1629 * This allows us to declare the Schur complement <code>op_S</code> and
1630 * the approximate Schur complement <code>op_aS</code>:
1633 * const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
1634 * const auto op_aS =
1635 * transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;
1639 * We now create a preconditioner out of <code>op_aS</code> that
1640 * applies a fixed number of 30 (inexpensive) CG iterations:
1643 * IterationNumberControl iteration_number_control_aS(30, 1.e-18);
1644 * SolverCG<Vector<double>> solver_aS(iteration_number_control_aS);
1646 * const auto preconditioner_S =
1647 * inverse_operator(op_aS, solver_aS, PreconditionIdentity());
1651 * Now on to the first equation. The right hand side of it is
1652 * @f$B^TM^{-1}F-G@f$, which is what we compute in the first few lines. We
1653 * then solve the first equation with a CG solver and the
1654 * preconditioner we just declared.
1657 * const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
1659 * SolverControl solver_control_S(2000, 1.e-12);
1660 * SolverCG<Vector<double>> solver_S(solver_control_S);
1662 * const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);
1664 * P = op_S_inv * schur_rhs;
1666 * std::cout << solver_control_S.last_step()
1667 * << " CG Schur complement iterations to obtain convergence."
1672 * After we have the pressure, we can compute the velocity. The equation
1673 * reads @f$MU=-BP+F@f$, and we solve it by first computing the right hand
1674 * side, and then multiplying it with the object that represents the
1675 * inverse of the @ref GlossMassMatrix "mass matrix":
1678 * U = op_M_inv * (F - op_B * P);
1685 * <a name="step_20-MixedLaplaceProblemclassimplementationcontinued"></a>
1686 * <h3>MixedLaplaceProblem class implementation (continued)</h3>
1691 * <a name="step_20-MixedLaplacecompute_errors"></a>
1692 * <h4>MixedLaplace::compute_errors</h4>
1696 * After we have dealt with the linear solver and preconditioners, we
1697 * continue with the implementation of our main class. In particular, the
1698 * next task is to compute the errors in our numerical solution, in both the
1699 * pressures as well as velocities.
1703 * To compute errors in the solution, we have already introduced the
1704 * <code>VectorTools::integrate_difference</code> function in @ref step_7 "step-7" and
1705 * @ref step_11 "step-11". However, there we only dealt with scalar solutions, whereas here
1706 * we have a vector-valued solution with components that even denote
1707 * different quantities and may have different orders of convergence (this
1711 * in. This is easily done: the
1712 * <code>VectorTools::integrate_difference</code> function takes as one of its
1713 * arguments a pointer to a weight function (the parameter defaults to the
1714 * null pointer, meaning unit weights). What we have to do is to pass
1715 * a function object that equals one in the components we are interested in,
1716 * and zero in the other ones. For example, to compute the pressure error,
1717 * we should pass a function that represents the constant vector with a unit
1718 * value in component <code>dim</code>, whereas for the velocity the
1719 * constant vector should be one in the first <code>dim</code> components,
1720 * and zero in the location of the pressure.
1724 * In deal.II, the <code>ComponentSelectFunction</code> does exactly this:
1725 * it wants to know how many vector components the function it is to
1726 * represent should have (in our case this would be <code>dim+1</code>, for
1727 * the joint velocity-pressure space) and which individual or range of
1728 * components should be equal to one. We therefore define two such masks at
1729 * the beginning of the function, following by an object representing the
1730 * exact solution and a vector in which we will store the cellwise errors as
1731 * computed by <code>integrate_difference</code>:
1734 * template <int dim>
1735 * void MixedLaplaceProblem<dim>::compute_errors() const
1737 * const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
1738 * const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
1741 * PrescribedSolution::ExactSolution<dim> exact_solution;
1742 * Vector<double> cellwise_errors(triangulation.n_active_cells());
1746 * As already discussed in @ref step_7 "step-7", we have to realize that it is
1747 * impossible to integrate the errors exactly. All we can do is
1748 * approximate this integral using quadrature. This actually presents a
1749 * slight twist here: if we naively chose an object of type
1750 * <code>QGauss@<dim@>(degree+1)</code> as one may be inclined to do (this
1751 * is what we used for integrating the linear system), one realizes that
1752 * the error is very small and does not follow the expected convergence
1753 * curves at all. What is happening is that for the mixed finite elements
1754 * used here, the Gauss points happen to be superconvergence points in
1755 * which the pointwise error is much smaller (and converges with higher
1756 * order) than anywhere else. These are therefore not particularly good
1757 * points for integration. To avoid this problem, we simply use a
1758 * trapezoidal rule and iterate it <code>degree+2</code> times in each
1759 * coordinate direction (again as explained in @ref step_7 "step-7"):
1762 * const QTrapezoid<1> q_trapez;
1763 * const QIterated<dim> quadrature(q_trapez, degree + 2);
1767 * With this, we can then let the library compute the errors and output
1768 * them to the screen:
1771 * VectorTools::integrate_difference(dof_handler,
1776 * VectorTools::L2_norm,
1778 * const double p_l2_error =
1779 * VectorTools::compute_global_error(triangulation,
1781 * VectorTools::L2_norm);
1783 * VectorTools::integrate_difference(dof_handler,
1788 * VectorTools::L2_norm,
1790 * const double u_l2_error =
1791 * VectorTools::compute_global_error(triangulation,
1793 * VectorTools::L2_norm);
1795 * std::cout << "Errors: ||e_p||_L2 = " << p_l2_error
1796 * << ", ||e_u||_L2 = " << u_l2_error << std::endl;
1803 * <a name="step_20-MixedLaplaceoutput_results"></a>
1804 * <h4>MixedLaplace::output_results</h4>
1808 * The last interesting function is the one in which we generate graphical
1809 * output. Note that all velocity components get the same solution name
1810 * "u". Together with using
1811 * DataComponentInterpretation::component_is_part_of_vector this will
1812 * cause DataOut<dim>::write_vtu() to generate a vector representation of
1813 * the individual velocity components, see @ref step_22 "step-22" or the
1814 * @ref VVOutput "Generating graphical output"
1816 * @ref vector_valued
1817 * topic for more information. Finally, it seems inappropriate for higher
1818 * order elements to only show a single bilinear quadrilateral per cell in
1819 * the graphical output. We therefore generate patches of size
1820 * (degree+1)x(degree+1) to capture the full information content of the
1821 * solution. See the @ref step_7 "step-7" tutorial program for more information on this.
1824 * template <int dim>
1825 * void MixedLaplaceProblem<dim>::output_results() const
1827 * std::vector<std::string> solution_names(dim, "u");
1828 * solution_names.emplace_back("p");
1829 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1830 * interpretation(dim,
1831 * DataComponentInterpretation::component_is_part_of_vector);
1832 * interpretation.push_back(DataComponentInterpretation::component_is_scalar);
1834 * DataOut<dim> data_out;
1835 * data_out.add_data_vector(dof_handler,
1840 * data_out.build_patches(degree + 1);
1842 * std::ofstream output("solution.vtu");
1843 * data_out.write_vtu(output);
1851 * <a name="step_20-MixedLaplacerun"></a>
1852 * <h4>MixedLaplace::run</h4>
1856 * This is the final function of our main class. It's
only job is to call
1875 * <a name=
"step_20-Thecodemaincodefunction"></a>
1891 *
using namespace Step20;
1893 *
const unsigned int fe_degree = 0;
1897 *
catch (std::exception &exc)
1899 *
std::cerr << std::endl
1901 *
<<
"----------------------------------------------------"
1903 *
std::cerr <<
"Exception on processing: " << std::endl
1904 *
<< exc.what() << std::endl
1905 *
<<
"Aborting!" << std::endl
1906 *
<<
"----------------------------------------------------"
1913 *
std::cerr << std::endl
1915 *
<<
"----------------------------------------------------"
1917 *
std::cerr <<
"Unknown exception!" << std::endl
1918 *
<<
"Aborting!" << std::endl
1919 *
<<
"----------------------------------------------------"
1961 <
td><
img src=
"https://www.dealii.org/images/steps/developer/step-20.u_new.jpg" width=
"400" alt=
""></
td>
1962 <
td><
img src=
"https://www.dealii.org/images/steps/developer/step-20.v_new.jpg" width=
"400" alt=
""></
td>
1963 <
td><
img src=
"https://www.dealii.org/images/steps/developer/step-20.p_new.jpg" width=
"400" alt=
""></
td>
2091<a name=
"step-20-extensions"></a>
2092<a name=
"step_20-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
2100this situation: we use a permeability that decays very rapidly away from a
2101central flowline until it hits a background value of 0.001. This is to mimic
2102the behavior of fluids in sandstone: in most of the domain, the sandstone is
2103homogeneous and, while permeable to fluids, not overly so; on the other stone,
2104the stone has cracked, or faulted, along one line, and the fluids flow much
2105easier along this large crack. Here is how we could implement something like
2110KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
2111 std::vector<Tensor<2,dim> > &values) const
2113 AssertDimension (points.size(), values.size());
2115 for (unsigned int p=0; p<points.size(); ++p)
2119 const double distance_to_flowline
2120 = std::fabs(points[p][1]-0.2*std::sin(10*points[p][0]));
2122 const double permeability = std::max(std::exp(-(distance_to_flowline*
2123 distance_to_flowline)
2127 for (unsigned int d=0; d<dim; ++d)
2128 values[p][d][d] = 1./permeability;
2132Remember that the function returns the inverse of the permeability tensor.
2136With a significantly higher mesh resolution, we can visualize this, here with
2139<table style="width:60%" align="center">
2141 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-wiggle.png" alt=""></td>
2142 <td><img src="https://www.dealii.org/images/steps/developer/step-20.v-wiggle.png" alt=""></td>
2146It is obvious how fluids flow essentially only along the middle line, and not
2151Another possibility would be to use a random permeability field. A simple way
2152to achieve this would be to scatter a number of centers around the domain and
2153then use a permeability field that is the sum of (negative) exponentials for
2154each of these centers. Flow would then try to hop from one center of high
2155permeability to the next one. This is an entirely unscientific attempt at
2156describing a random medium, but one possibility to implement this behavior
2157would look like this:
2160class KInverse : public TensorFunction<2,dim>
2165 virtual void value_list (const std::vector<Point<dim> > &points,
2166 std::vector<Tensor<2,dim> > &values) const;
2169 std::vector<Point<dim> > centers;
2174KInverse<dim>::KInverse ()
2176 const unsigned int N = 40;
2178 for (unsigned int i=0; i<N; ++i)
2179 for (unsigned int d=0; d<dim; ++d)
2180 centers[i][d] = 2.*rand()/RAND_MAX-1;
2186KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
2187 std::vector<Tensor<2,dim> > &values) const
2189 AssertDimension (points.size(), values.size());
2191 for (unsigned int p=0; p<points.size(); ++p)
2195 double permeability = 0;
2196 for (unsigned int i=0; i<centers.size(); ++i)
2197 permeability += std::exp(-(points[p] - centers[i]).norm_square() / (0.1 * 0.1));
2199 const double normalized_permeability
2200 = std::max(permeability, 0.005);
2202 for (unsigned int d=0; d<dim; ++d)
2203 values[p][d][d] = 1./normalized_permeability;
2208A piecewise constant interpolation of the diagonal elements of the
2209inverse of this tensor (i.e., of <code>normalized_permeability</code>)
2212<img src="https://www.dealii.org/images/steps/developer/step-20.k-random.png" alt="">
2215With a permeability field like this, we would get x-velocities and pressures as
2218<table style="width:60%" align="center">
2220 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-random.png" alt=""></td>
2221 <td><img src="https://www.dealii.org/images/steps/developer/step-20.p-random.png" alt=""></td>
2225We will use these permeability fields again in @ref step_21 "step-21" and @ref step_43 "step-43".
2228<a name="step_20-Betterlinearsolvers"></a><h4>Better linear solvers</h4>
2231As mentioned in the introduction, the Schur complement solver used here is not
2232the best one conceivable (nor is it intended to be a particularly good
2233one). Better ones can be found in the literature and can be built using the
2234same block matrix techniques that were introduced here. We pick up on this
2235theme again in @ref step_22 "step-22", where we first build a Schur complement solver for the
2236Stokes equation as we did here, and then in the <a
2237href="step_22.html#improved-solver">Improved Solvers</a> section discuss better
2238ways based on solving the system as a whole but preconditioning based on
2239individual blocks. We will also come back to this in @ref step_43 "step-43".
2242<a name="step_20-PlainProg"></a>
2243<h1> The plain program</h1>
2244@include "step-20.cc"
static constexpr unsigned int rank
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define AssertDimension(dim1, dim2)
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())
std::vector< index_type > data
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool 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)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()