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-20.h
Go to the documentation of this file.
1,
1083 *   const unsigned int /*component*/) const
1084 *   {
1085 *   return 0;
1086 *   }
1087 *  
1088 *  
1089 *  
1090 *   template <int dim>
1091 *   double
1092 *   PressureBoundaryValues<dim>::value(const Point<dim> &p,
1093 *   const unsigned int /*component*/) const
1094 *   {
1095 *   return -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
1096 *   alpha * p[0] * p[0] * p[0] / 6);
1097 *   }
1098 *  
1099 *  
1100 *  
1101 *   template <int dim>
1102 *   void ExactSolution<dim>::vector_value(const Point<dim> &p,
1103 *   Vector<double> & values) const
1104 *   {
1105 *   AssertDimension(values.size(), dim + 1);
1106 *  
1107 *   values(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2;
1108 *   values(1) = alpha * p[0] * p[1];
1109 *   values(2) = -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
1110 *   alpha * p[0] * p[0] * p[0] / 6);
1111 *   }
1112 *  
1113 *  
1114 *  
1115 * @endcode
1116 *
1117 *
1118 * <a name="Theinversepermeabilitytensor"></a>
1119 * <h3>The inverse permeability tensor</h3>
1120 *
1121
1122 *
1123 * In addition to the other equation data, we also want to use a
1124 * permeability tensor, or better -- because this is all that appears in the
1125 * weak form -- the inverse of the permeability tensor,
1126 * <code>KInverse</code>. For the purpose of verifying the exactness of the
1127 * solution and determining convergence orders, this tensor is more in the
1128 * way than helpful. We will therefore simply set it to the identity matrix.
1129 *
1130
1131 *
1132 * However, a spatially varying permeability tensor is indispensable in
1133 * real-life porous media flow simulations, and we would like to use the
1134 * opportunity to demonstrate the technique to use tensor valued functions.
1135 *
1136
1137 *
1138 * Possibly unsurprisingly, deal.II also has a base class not only for
1139 * scalar and generally vector-valued functions (the <code>Function</code>
1140 * base class) but also for functions that return tensors of fixed dimension
1141 * and rank, the <code>TensorFunction</code> template. Here, the function
1142 * under consideration returns a dim-by-dim matrix, i.e. a tensor of rank 2
1143 * and dimension <code>dim</code>. We then choose the template arguments of
1144 * the base class appropriately.
1145 *
1146
1147 *
1148 * The interface that the <code>TensorFunction</code> class provides is
1149 * essentially equivalent to the <code>Function</code> class. In particular,
1150 * there exists a <code>value_list</code> function that takes a list of
1151 * points at which to evaluate the function, and returns the values of the
1152 * function in the second argument, a list of tensors:
1153 *
1154 * @code
1155 *   template <int dim>
1156 *   class KInverse : public TensorFunction<2, dim>
1157 *   {
1158 *   public:
1159 *   KInverse()
1160 *   : TensorFunction<2, dim>()
1161 *   {}
1162 *  
1163 *   virtual void
1164 *   value_list(const std::vector<Point<dim>> &points,
1165 *   std::vector<Tensor<2, dim>> & values) const override;
1166 *   };
1167 *  
1168 *  
1169 * @endcode
1170 *
1171 * The implementation is less interesting. As in previous examples, we add a
1172 * check to the beginning of the class to make sure that the sizes of input
1173 * and output parameters are the same (see @ref step_5 "step-5" for a discussion of this
1174 * technique). Then we loop over all evaluation points, and for each one
1175 * set the output tensor to the identity matrix.
1176 *
1177
1178 *
1179 * There is an oddity at the top of the function (the
1180 * `(void)points;` statement) that is worth discussing. The values
1181 * we put into the output `values` array does not actually depend
1182 * on the `points` arrays of coordinates at which the function is
1183 * evaluated. In other words, the `points` argument is in fact
1184 * unused, and we could have just not given it a name if we had
1185 * wanted. But we want to use the `points` object for checking
1186 * that the `values` object has the correct size. The problem is
1187 * that in release mode, `AssertDimension` is defined as a macro
1188 * that expands to nothing; the compiler will then complain that
1189 * the `points` object is unused. The idiomatic approach to
1190 * silencing this warning is to have a statement that evaluates
1191 * (reads) variable but doesn't actually do anything: That's what
1192 * `(void)points;` does: It reads from `points`, and then casts
1193 * the result of the read to `void`, i.e., nothing. This statement
1194 * is, in other words, completely pointless and implies no actual
1195 * action except to explain to the compiler that yes, this
1196 * variable is in fact used even in release mode. (In debug mode,
1197 * the `AssertDimension` macro expands to something that reads
1198 * from the variable, and so the funny statement would not be
1199 * necessary in debug mode.)
1200 *
1201 * @code
1202 *   template <int dim>
1203 *   void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
1204 *   std::vector<Tensor<2, dim>> & values) const
1205 *   {
1206 *   (void)points;
1207 *   AssertDimension(points.size(), values.size());
1208 *  
1209 *   for (auto &value : values)
1210 *   value = unit_symmetric_tensor<dim>();
1211 *   }
1212 *   } // namespace PrescribedSolution
1213 *  
1214 *  
1215 *  
1216 * @endcode
1217 *
1218 *
1219 * <a name="MixedLaplaceProblemclassimplementation"></a>
1220 * <h3>MixedLaplaceProblem class implementation</h3>
1221 *
1222
1223 *
1224 *
1225 * <a name="MixedLaplaceProblemMixedLaplaceProblem"></a>
1226 * <h4>MixedLaplaceProblem::MixedLaplaceProblem</h4>
1227 *
1228
1229 *
1230 * In the constructor of this class, we first store the value that was
1231 * passed in concerning the degree of the finite elements we shall use (a
1232 * degree of zero, for example, means to use RT(0) and DG(0)), and then
1233 * construct the vector valued element belonging to the space @f$X_h@f$ described
1234 * in the introduction. The rest of the constructor is as in the early
1235 * tutorial programs.
1236 *
1237
1238 *
1239 * The only thing worth describing here is the constructor call of the
1240 * <code>fe</code> variable. The <code>FESystem</code> class to which this
1241 * variable belongs has a number of different constructors that all refer to
1242 * binding simpler elements together into one larger element. In the present
1243 * case, we want to couple a single RT(degree) element with a single
1244 * DQ(degree) element. The constructor to <code>FESystem</code> that does
1245 * this requires us to specify first the first base element (the
1246 * <code>FE_RaviartThomas</code> object of given degree) and then the number
1247 * of copies for this base element, and then similarly the kind and number
1248 * of <code>FE_DGQ</code> elements. Note that the Raviart-Thomas element
1249 * already has <code>dim</code> vector components, so that the coupled
1250 * element will have <code>dim+1</code> vector components, the first
1251 * <code>dim</code> of which correspond to the velocity variable whereas the
1252 * last one corresponds to the pressure.
1253 *
1254
1255 *
1256 * It is also worth comparing the way we constructed this element from its
1257 * base elements, with the way we have done so in @ref step_8 "step-8": there, we have
1258 * built it as <code>fe (FE_Q@<dim@>(1), dim)</code>, i.e. we have simply
1259 * used <code>dim</code> copies of the <code>FE_Q(1)</code> element, one
1260 * copy for the displacement in each coordinate direction.
1261 *
1262 * @code
1263 *   template <int dim>
1264 *   MixedLaplaceProblem<dim>::MixedLaplaceProblem(const unsigned int degree)
1265 *   : degree(degree)
1266 *   , fe(FE_RaviartThomas<dim>(degree), 1, FE_DGQ<dim>(degree), 1)
1267 *   , dof_handler(triangulation)
1268 *   {}
1269 *  
1270 *  
1271 *  
1272 * @endcode
1273 *
1274 *
1275 * <a name="MixedLaplaceProblemmake_grid_and_dofs"></a>
1276 * <h4>MixedLaplaceProblem::make_grid_and_dofs</h4>
1277 *
1278
1279 *
1280 * This next function starts out with well-known functions calls that create
1281 * and refine a mesh, and then associate degrees of freedom with it:
1282 *
1283 * @code
1284 *   template <int dim>
1285 *   void MixedLaplaceProblem<dim>::make_grid_and_dofs()
1286 *   {
1288 *   triangulation.refine_global(5);
1289 *  
1290 *   dof_handler.distribute_dofs(fe);
1291 *  
1292 * @endcode
1293 *
1294 * However, then things become different. As mentioned in the
1295 * introduction, we want to subdivide the matrix into blocks corresponding
1296 * to the two different kinds of variables, velocity and pressure. To this
1297 * end, we first have to make sure that the indices corresponding to
1298 * velocities and pressures are not intermingled: First all velocity
1299 * degrees of freedom, then all pressure DoFs. This way, the global matrix
1300 * separates nicely into a @f$2 \times 2@f$ system. To achieve this, we have to
1301 * renumber degrees of freedom based on their vector component, an
1302 * operation that conveniently is already implemented:
1303 *
1304 * @code
1305 *   DoFRenumbering::component_wise(dof_handler);
1306 *  
1307 * @endcode
1308 *
1309 * The next thing is that we want to figure out the sizes of these blocks
1310 * so that we can allocate an appropriate amount of space. To this end, we
1311 * call the DoFTools::count_dofs_per_fe_component() function that
1312 * counts how many shape functions are non-zero for a particular vector
1313 * component. We have <code>dim+1</code> vector components, and
1314 * DoFTools::count_dofs_per_fe_component() will count how many shape
1315 * functions belong to each of these components.
1316 *
1317
1318 *
1319 * There is one problem here. As described in the documentation of that
1320 * function, it <i>wants</i> to put the number of @f$x@f$-velocity shape
1321 * functions into <code>dofs_per_component[0]</code>, the number of
1322 * @f$y@f$-velocity shape functions into <code>dofs_per_component[1]</code>
1323 * (and similar in 3d), and the number of pressure shape functions into
1324 * <code>dofs_per_component[dim]</code>. But, the Raviart-Thomas element
1325 * is special in that it is non-@ref GlossPrimitive "primitive", i.e.,
1326 * for Raviart-Thomas elements all velocity shape functions
1327 * are nonzero in all components. In other words, the function cannot
1328 * distinguish between @f$x@f$ and @f$y@f$ velocity functions because there
1329 * <i>is</i> no such distinction. It therefore puts the overall number
1330 * of velocity into each of <code>dofs_per_component[c]</code>,
1331 * @f$0\le c\le \text{dim}@f$. On the other hand, the number
1332 * of pressure variables equals the number of shape functions that are
1333 * nonzero in the dim-th component.
1334 *
1335
1336 *
1337 * Using this knowledge, we can get the number of velocity shape
1338 * functions from any of the first <code>dim</code> elements of
1339 * <code>dofs_per_component</code>, and then use this below to initialize
1340 * the vector and matrix block sizes, as well as create output.
1341 *
1342
1343 *
1344 * @note If you find this concept difficult to understand, you may
1345 * want to consider using the function DoFTools::count_dofs_per_fe_block()
1346 * instead, as we do in the corresponding piece of code in @ref step_22 "step-22".
1347 * You might also want to read up on the difference between
1348 * @ref GlossBlock "blocks" and @ref GlossComponent "components"
1349 * in the glossary.
1350 *
1351 * @code
1352 *   const std::vector<types::global_dof_index> dofs_per_component =
1353 *   DoFTools::count_dofs_per_fe_component(dof_handler);
1354 *   const unsigned int n_u = dofs_per_component[0],
1355 *   n_p = dofs_per_component[dim];
1356 *  
1357 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
1358 *   << std::endl
1359 *   << "Total number of cells: " << triangulation.n_cells()
1360 *   << std::endl
1361 *   << "Number of degrees of freedom: " << dof_handler.n_dofs()
1362 *   << " (" << n_u << '+' << n_p << ')' << std::endl;
1363 *  
1364 * @endcode
1365 *
1366 * The next task is to allocate a sparsity pattern for the matrix that we
1367 * will create. We use a compressed sparsity pattern like in the previous
1368 * steps, but as <code>system_matrix</code> is a block matrix we use the
1369 * class <code>BlockDynamicSparsityPattern</code> instead of just
1370 * <code>DynamicSparsityPattern</code>. This block sparsity pattern has
1371 * four blocks in a @f$2 \times 2@f$ pattern. The blocks' sizes depend on
1372 * <code>n_u</code> and <code>n_p</code>, which hold the number of velocity
1373 * and pressure variables.
1374 *
1375 * @code
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);
1379 *  
1380 * @endcode
1381 *
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
1384 * matrix:
1385 *
1386 * @code
1387 *   sparsity_pattern.copy_from(dsp);
1388 *   system_matrix.reinit(sparsity_pattern);
1389 *  
1390 * @endcode
1391 *
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:
1394 *
1395 * @code
1396 *   solution.reinit(block_sizes);
1397 *   system_rhs.reinit(block_sizes);
1398 *   }
1399 *  
1400 *  
1401 * @endcode
1402 *
1403 *
1404 * <a name="MixedLaplaceProblemassemble_system"></a>
1405 * <h4>MixedLaplaceProblem::assemble_system</h4>
1406 *
1407
1408 *
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.
1417 *
1418 * @code
1419 *   template <int dim>
1420 *   void MixedLaplaceProblem<dim>::assemble_system()
1421 *   {
1422 *   QGauss<dim> quadrature_formula(degree + 2);
1423 *   QGauss<dim - 1> face_quadrature_formula(degree + 2);
1424 *  
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);
1434 *  
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();
1438 *  
1439 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1440 *   Vector<double> local_rhs(dofs_per_cell);
1441 *  
1442 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1443 *  
1444 * @endcode
1445 *
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.
1452 *
1453 * @code
1454 *   const PrescribedSolution::RightHandSide<dim> right_hand_side;
1455 *   const PrescribedSolution::PressureBoundaryValues<dim>
1456 *   pressure_boundary_values;
1457 *   const PrescribedSolution::KInverse<dim> k_inverse;
1458 *  
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);
1462 *  
1463 * @endcode
1464 *
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>):
1474 *
1475 * @code
1476 *   const FEValuesExtractors::Vector velocities(0);
1477 *   const FEValuesExtractors::Scalar pressure(dim);
1478 *  
1479 * @endcode
1480 *
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:
1484 *
1485 * @code
1486 *   for (const auto &cell : dof_handler.active_cell_iterators())
1487 *   {
1488 *   fe_values.reinit(cell);
1489 *   local_matrix = 0;
1490 *   local_rhs = 0;
1491 *  
1492 *   right_hand_side.value_list(fe_values.get_quadrature_points(),
1493 *   rhs_values);
1494 *   k_inverse.value_list(fe_values.get_quadrature_points(),
1495 *   k_inverse_values);
1496 *  
1497 *   for (unsigned int q = 0; q < n_q_points; ++q)
1498 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1499 *   {
1500 *   const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
1501 *   const double div_phi_i_u = fe_values[velocities].divergence(i, q);
1502 *   const double phi_i_p = fe_values[pressure].value(i, q);
1503 *  
1504 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1505 *   {
1506 *   const Tensor<1, dim> phi_j_u =
1507 *   fe_values[velocities].value(j, q);
1508 *   const double div_phi_j_u =
1509 *   fe_values[velocities].divergence(j, q);
1510 *   const double phi_j_p = fe_values[pressure].value(j, q);
1511 *  
1512 *   local_matrix(i, j) +=
1513 *   (phi_i_u * k_inverse_values[q] * phi_j_u
1514 *   - phi_i_p * div_phi_j_u
1515 *   - div_phi_i_u * phi_j_p)
1516 *   * fe_values.JxW(q);
1517 *   }
1518 *  
1519 *   local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
1520 *   }
1521 *  
1522 *   for (const auto &face : cell->face_iterators())
1523 *   if (face->at_boundary())
1524 *   {
1525 *   fe_face_values.reinit(cell, face);
1526 *  
1527 *   pressure_boundary_values.value_list(
1528 *   fe_face_values.get_quadrature_points(), boundary_values);
1529 *  
1530 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1531 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1532 *   local_rhs(i) += -(fe_face_values[velocities].value(i, q) *
1533 *   fe_face_values.normal_vector(q) *
1534 *   boundary_values[q] *
1535 *   fe_face_values.JxW(q));
1536 *   }
1537 *  
1538 * @endcode
1539 *
1540 * The final step in the loop over all cells is to transfer local
1541 * contributions into the global matrix and right hand side
1542 * vector. Note that we use exactly the same interface as in previous
1543 * examples, although we now use block matrices and vectors instead of
1544 * the regular ones. In other words, to the outside world, block
1545 * objects have the same interface as matrices and vectors, but they
1546 * additionally allow to access individual blocks.
1547 *
1548 * @code
1549 *   cell->get_dof_indices(local_dof_indices);
1550 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1551 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1552 *   system_matrix.add(local_dof_indices[i],
1553 *   local_dof_indices[j],
1554 *   local_matrix(i, j));
1555 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1556 *   system_rhs(local_dof_indices[i]) += local_rhs(i);
1557 *   }
1558 *   }
1559 *  
1560 *  
1561 * @endcode
1562 *
1563 *
1564 * <a name="Implementationoflinearsolversandpreconditioners"></a>
1565 * <h3>Implementation of linear solvers and preconditioners</h3>
1566 *
1567
1568 *
1569 * The linear solvers and preconditioners we use in this example have
1570 * been discussed in significant detail already in the introduction. We
1571 * will therefore not discuss the rationale for our approach here any
1572 * more, but rather only comment on some remaining implementational
1573 * aspects.
1574 *
1575
1576 *
1577 *
1578 * <a name="MixedLaplacesolve"></a>
1579 * <h4>MixedLaplace::solve</h4>
1580 *
1581
1582 *
1583 * As already outlined in the introduction, the solve function consists
1584 * essentially of two steps. First, we have to form the first equation
1585 * involving the Schur complement and solve for the pressure (component 1
1586 * of the solution). Then, we can reconstruct the velocities from the
1587 * second equation (component 0 of the solution).
1588 *
1589 * @code
1590 *   template <int dim>
1591 *   void MixedLaplaceProblem<dim>::solve()
1592 *   {
1593 * @endcode
1594 *
1595 * As a first step we declare references to all block components of the
1596 * matrix, the right hand side and the solution vector that we will
1597 * need.
1598 *
1599 * @code
1600 *   const auto &M = system_matrix.block(0, 0);
1601 *   const auto &B = system_matrix.block(0, 1);
1602 *  
1603 *   const auto &F = system_rhs.block(0);
1604 *   const auto &G = system_rhs.block(1);
1605 *  
1606 *   auto &U = solution.block(0);
1607 *   auto &P = solution.block(1);
1608 *  
1609 * @endcode
1610 *
1611 * Then, we will create corresponding LinearOperator objects and create
1612 * the <code>op_M_inv</code> operator:
1613 *
1614 * @code
1615 *   const auto op_M = linear_operator(M);
1616 *   const auto op_B = linear_operator(B);
1617 *  
1618 *   ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10);
1619 *   SolverCG<Vector<double>> solver_M(reduction_control_M);
1620 *   PreconditionJacobi<SparseMatrix<double>> preconditioner_M;
1621 *  
1622 *   preconditioner_M.initialize(M);
1623 *  
1624 *   const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
1625 *  
1626 * @endcode
1627 *
1628 * This allows us to declare the Schur complement <code>op_S</code> and
1629 * the approximate Schur complement <code>op_aS</code>:
1630 *
1631 * @code
1632 *   const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
1633 *   const auto op_aS =
1634 *   transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;
1635 *  
1636 * @endcode
1637 *
1638 * We now create a preconditioner out of <code>op_aS</code> that
1639 * applies a fixed number of 30 (inexpensive) CG iterations:
1640 *
1641 * @code
1642 *   IterationNumberControl iteration_number_control_aS(30, 1.e-18);
1643 *   SolverCG<Vector<double>> solver_aS(iteration_number_control_aS);
1644 *  
1645 *   const auto preconditioner_S =
1646 *   inverse_operator(op_aS, solver_aS, PreconditionIdentity());
1647 *  
1648 * @endcode
1649 *
1650 * Now on to the first equation. The right hand side of it is
1651 * @f$B^TM^{-1}F-G@f$, which is what we compute in the first few lines. We
1652 * then solve the first equation with a CG solver and the
1653 * preconditioner we just declared.
1654 *
1655 * @code
1656 *   const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
1657 *  
1658 *   SolverControl solver_control_S(2000, 1.e-12);
1659 *   SolverCG<Vector<double>> solver_S(solver_control_S);
1660 *  
1661 *   const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);
1662 *  
1663 *   P = op_S_inv * schur_rhs;
1664 *  
1665 *   std::cout << solver_control_S.last_step()
1666 *   << " CG Schur complement iterations to obtain convergence."
1667 *   << std::endl;
1668 *  
1669 * @endcode
1670 *
1671 * After we have the pressure, we can compute the velocity. The equation
1672 * reads @f$MU=-BP+F@f$, and we solve it by first computing the right hand
1673 * side, and then multiplying it with the object that represents the
1674 * inverse of the @ref GlossMassMatrix "mass matrix":
1675 *
1676 * @code
1677 *   U = op_M_inv * (F - op_B * P);
1678 *   }
1679 *  
1680 *  
1681 * @endcode
1682 *
1683 *
1684 * <a name="MixedLaplaceProblemclassimplementationcontinued"></a>
1685 * <h3>MixedLaplaceProblem class implementation (continued)</h3>
1686 *
1687
1688 *
1689 *
1690 * <a name="MixedLaplacecompute_errors"></a>
1691 * <h4>MixedLaplace::compute_errors</h4>
1692 *
1693
1694 *
1695 * After we have dealt with the linear solver and preconditioners, we
1696 * continue with the implementation of our main class. In particular, the
1697 * next task is to compute the errors in our numerical solution, in both the
1698 * pressures as well as velocities.
1699 *
1700
1701 *
1702 * To compute errors in the solution, we have already introduced the
1703 * <code>VectorTools::integrate_difference</code> function in @ref step_7 "step-7" and
1704 * @ref step_11 "step-11". However, there we only dealt with scalar solutions, whereas here
1705 * we have a vector-valued solution with components that even denote
1706 * different quantities and may have different orders of convergence (this
1707 * isn't the case here, by choice of the used finite elements, but is
1708 * frequently the case in mixed finite element applications). What we
1709 * therefore have to do is to `mask' the components that we are interested
1710 * in. This is easily done: the
1711 * <code>VectorTools::integrate_difference</code> function takes as one of its
1712 * arguments a pointer to a weight function (the parameter defaults to the
1713 * null pointer, meaning unit weights). What we have to do is to pass
1714 * a function object that equals one in the components we are interested in,
1715 * and zero in the other ones. For example, to compute the pressure error,
1716 * we should pass a function that represents the constant vector with a unit
1717 * value in component <code>dim</code>, whereas for the velocity the
1718 * constant vector should be one in the first <code>dim</code> components,
1719 * and zero in the location of the pressure.
1720 *
1721
1722 *
1723 * In deal.II, the <code>ComponentSelectFunction</code> does exactly this:
1724 * it wants to know how many vector components the function it is to
1725 * represent should have (in our case this would be <code>dim+1</code>, for
1726 * the joint velocity-pressure space) and which individual or range of
1727 * components should be equal to one. We therefore define two such masks at
1728 * the beginning of the function, following by an object representing the
1729 * exact solution and a vector in which we will store the cellwise errors as
1730 * computed by <code>integrate_difference</code>:
1731 *
1732 * @code
1733 *   template <int dim>
1734 *   void MixedLaplaceProblem<dim>::compute_errors() const
1735 *   {
1736 *   const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
1737 *   const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
1738 *   dim + 1);
1739 *  
1740 *   PrescribedSolution::ExactSolution<dim> exact_solution;
1741 *   Vector<double> cellwise_errors(triangulation.n_active_cells());
1742 *  
1743 * @endcode
1744 *
1745 * As already discussed in @ref step_7 "step-7", we have to realize that it is
1746 * impossible to integrate the errors exactly. All we can do is
1747 * approximate this integral using quadrature. This actually presents a
1748 * slight twist here: if we naively chose an object of type
1749 * <code>QGauss@<dim@>(degree+1)</code> as one may be inclined to do (this
1750 * is what we used for integrating the linear system), one realizes that
1751 * the error is very small and does not follow the expected convergence
1752 * curves at all. What is happening is that for the mixed finite elements
1753 * used here, the Gauss points happen to be superconvergence points in
1754 * which the pointwise error is much smaller (and converges with higher
1755 * order) than anywhere else. These are therefore not particularly good
1756 * points for integration. To avoid this problem, we simply use a
1757 * trapezoidal rule and iterate it <code>degree+2</code> times in each
1758 * coordinate direction (again as explained in @ref step_7 "step-7"):
1759 *
1760 * @code
1761 *   QTrapezoid<1> q_trapez;
1762 *   QIterated<dim> quadrature(q_trapez, degree + 2);
1763 *  
1764 * @endcode
1765 *
1766 * With this, we can then let the library compute the errors and output
1767 * them to the screen:
1768 *
1769 * @code
1770 *   VectorTools::integrate_difference(dof_handler,
1771 *   solution,
1772 *   exact_solution,
1773 *   cellwise_errors,
1774 *   quadrature,
1775 *   VectorTools::L2_norm,
1776 *   &pressure_mask);
1777 *   const double p_l2_error =
1778 *   VectorTools::compute_global_error(triangulation,
1779 *   cellwise_errors,
1780 *   VectorTools::L2_norm);
1781 *  
1782 *   VectorTools::integrate_difference(dof_handler,
1783 *   solution,
1784 *   exact_solution,
1785 *   cellwise_errors,
1786 *   quadrature,
1787 *   VectorTools::L2_norm,
1788 *   &velocity_mask);
1789 *   const double u_l2_error =
1790 *   VectorTools::compute_global_error(triangulation,
1791 *   cellwise_errors,
1792 *   VectorTools::L2_norm);
1793 *  
1794 *   std::cout << "Errors: ||e_p||_L2 = " << p_l2_error
1795 *   << ", ||e_u||_L2 = " << u_l2_error << std::endl;
1796 *   }
1797 *  
1798 *  
1799 * @endcode
1800 *
1801 *
1802 * <a name="MixedLaplaceoutput_results"></a>
1803 * <h4>MixedLaplace::output_results</h4>
1804 *
1805
1806 *
1807 * The last interesting function is the one in which we generate graphical
1808 * output. Note that all velocity components get the same solution name
1809 * "u". Together with using
1810 * DataComponentInterpretation::component_is_part_of_vector this will
1811 * cause DataOut<dim>::write_vtu() to generate a vector representation of
1812 * the individual velocity components, see @ref step_22 "step-22" or the
1813 * @ref VVOutput "Generating graphical output"
1814 * section of the
1815 * @ref vector_valued
1816 * module for more information. Finally, it seems inappropriate for higher
1817 * order elements to only show a single bilinear quadrilateral per cell in
1818 * the graphical output. We therefore generate patches of size
1819 * (degree+1)x(degree+1) to capture the full information content of the
1820 * solution. See the @ref step_7 "step-7" tutorial program for more information on this.
1821 *
1822 * @code
1823 *   template <int dim>
1824 *   void MixedLaplaceProblem<dim>::output_results() const
1825 *   {
1826 *   std::vector<std::string> solution_names(dim, "u");
1827 *   solution_names.emplace_back("p");
1828 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1829 *   interpretation(dim,
1830 *   DataComponentInterpretation::component_is_part_of_vector);
1831 *   interpretation.push_back(DataComponentInterpretation::component_is_scalar);
1832 *  
1833 *   DataOut<dim> data_out;
1834 *   data_out.add_data_vector(dof_handler,
1835 *   solution,
1836 *   solution_names,
1837 *   interpretation);
1838 *  
1839 *   data_out.build_patches(degree + 1);
1840 *  
1841 *   std::ofstream output("solution.vtu");
1842 *   data_out.write_vtu(output);
1843 *   }
1844 *  
1845 *  
1846 *  
1847 * @endcode
1848 *
1849 *
1850 * <a name="MixedLaplacerun"></a>
1851 * <h4>MixedLaplace::run</h4>
1852 *
1853
1854 *
1855 * This is the final function of our main class. It's only job is to call
1856 * the other functions in their natural order:
1857 *
1858 * @code
1859 *   template <int dim>
1860 *   void MixedLaplaceProblem<dim>::run()
1861 *   {
1862 *   make_grid_and_dofs();
1863 *   assemble_system();
1864 *   solve();
1865 *   compute_errors();
1866 *   output_results();
1867 *   }
1868 *   } // namespace Step20
1869 *  
1870 *  
1871 * @endcode
1872 *
1873 *
1874 * <a name="Thecodemaincodefunction"></a>
1875 * <h3>The <code>main</code> function</h3>
1876 *
1877
1878 *
1879 * The main function we stole from @ref step_6 "step-6" instead of @ref step_4 "step-4". It is almost
1880 * equal to the one in @ref step_6 "step-6" (apart from the changed class names, of course),
1881 * the only exception is that we pass the degree of the finite element space
1882 * to the constructor of the mixed Laplace problem (here, we use zero-th order
1883 * elements).
1884 *
1885 * @code
1886 *   int main()
1887 *   {
1888 *   try
1889 *   {
1890 *   using namespace Step20;
1891 *  
1892 *   const unsigned int fe_degree = 0;
1893 *   MixedLaplaceProblem<2> mixed_laplace_problem(fe_degree);
1894 *   mixed_laplace_problem.run();
1895 *   }
1896 *   catch (std::exception &exc)
1897 *   {
1898 *   std::cerr << std::endl
1899 *   << std::endl
1900 *   << "----------------------------------------------------"
1901 *   << std::endl;
1902 *   std::cerr << "Exception on processing: " << std::endl
1903 *   << exc.what() << std::endl
1904 *   << "Aborting!" << std::endl
1905 *   << "----------------------------------------------------"
1906 *   << std::endl;
1907 *  
1908 *   return 1;
1909 *   }
1910 *   catch (...)
1911 *   {
1912 *   std::cerr << std::endl
1913 *   << std::endl
1914 *   << "----------------------------------------------------"
1915 *   << std::endl;
1916 *   std::cerr << "Unknown exception!" << std::endl
1917 *   << "Aborting!" << std::endl
1918 *   << "----------------------------------------------------"
1919 *   << std::endl;
1920 *   return 1;
1921 *   }
1922 *  
1923 *   return 0;
1924 *   }
1925 * @endcode
1926<a name="Results"></a><h1>Results</h1>
1927
1928
1929<a name="Outputoftheprogramandgraphicalvisualization"></a><h3>Output of the program and graphical visualization</h3>
1930
1931
1932
1933If we run the program as is, we get this output for the @f$32\times 32@f$
1934mesh we use (for a total of 1024 cells with 1024 pressure degrees of
1935freedom since we use piecewise constants, and 2112 velocities because
1936the Raviart-Thomas element defines one degree per freedom per face and
1937there are @f$1024 + 32 = 1056@f$ faces parallel to the @f$x@f$-axis and the same
1938number parallel to the @f$y@f$-axis):
1939@verbatim
1940\$ make run
1941[ 66%] Built target step-20
1942Scanning dependencies of target run
1943[100%] Run step-20 with Release configuration
1944Number of active cells: 1024
1945Total number of cells: 1365
1946Number of degrees of freedom: 3136 (2112+1024)
194724 CG Schur complement iterations to obtain convergence.
1948Errors: ||e_p||_L2 = 0.0445032, ||e_u||_L2 = 0.010826
1949[100%] Built target run
1950@endverbatim
1951
1952The fact that the number of iterations is so small, of course, is due to
1953the good (but expensive!) preconditioner we have developed. To get
1954confidence in the solution, let us take a look at it. The following three
1955images show (from left to right) the x-velocity, the y-velocity, and the
1956pressure:
1957
1958<table style="width:60%" align="center">
1959 <tr>
1960 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u_new.jpg" width="400" alt=""></td>
1961 <td><img src="https://www.dealii.org/images/steps/developer/step-20.v_new.jpg" width="400" alt=""></td>
1962 <td><img src="https://www.dealii.org/images/steps/developer/step-20.p_new.jpg" width="400" alt=""></td>
1963 </tr>
1964</table>
1965
1966
1967
1968Let us start with the pressure: it is highest at the left and lowest at the
1969right, so flow will be from left to right. In addition, though hardly visible
1970in the graph, we have chosen the pressure field such that the flow left-right
1971flow first channels towards the center and then outward again. Consequently,
1972the x-velocity has to increase to get the flow through the narrow part,
1973something that can easily be seen in the left image. The middle image
1974represents inward flow in y-direction at the left end of the domain, and
1975outward flow in y-direction at the right end of the domain.
1976
1977
1978
1979As an additional remark, note how the x-velocity in the left image is only
1980continuous in x-direction, whereas the y-velocity is continuous in
1981y-direction. The flow fields are discontinuous in the other directions. This
1982very obviously reflects the continuity properties of the Raviart-Thomas
1983elements, which are, in fact, only in the space H(div) and not in the space
1984@f$H^1@f$. Finally, the pressure field is completely discontinuous, but
1985that should not surprise given that we have chosen <code>FE_DGQ(0)</code> as
1986the finite element for that solution component.
1987
1988
1989
1990<a name="Convergence"></a><h3>Convergence</h3>
1991
1992
1993
1994The program offers two obvious places where playing and observing convergence
1995is in order: the degree of the finite elements used (passed to the constructor
1996of the <code>MixedLaplaceProblem</code> class from <code>main()</code>), and
1997the refinement level (determined in
1998<code>MixedLaplaceProblem::make_grid_and_dofs</code>). What one can do is to
1999change these values and observe the errors computed later on in the course of
2000the program run.
2001
2002
2003
2004If one does this, one finds the following pattern for the @f$L_2@f$ error
2005in the pressure variable:
2006<table align="center" class="doxtable">
2007 <tr>
2008 <th></th>
2009 <th colspan="3" align="center">Finite element order</th>
2010 </tr>
2011 <tr>
2012 <th>Refinement level</th>
2013 <th>0</th>
2014 <th>1</th>
2015 <th>2</th>
2016 </tr>
2017 <tr>
2018 <th>0</th> <td>1.45344</td> <td>0.0831743</td> <td>0.0235186</td>
2019 </tr>
2020 <tr>
2021 <th>1</th> <td>0.715099</td> <td>0.0245341</td> <td>0.00293983</td>
2022 </tr>
2023 <tr>
2024 <th>2</th> <td>0.356383</td> <td>0.0063458</td> <td>0.000367478</td>
2025 </tr>
2026 <tr>
2027 <th>3</th> <td>0.178055</td> <td>0.00159944</td> <td>4.59349e-05</td>
2028 </tr>
2029 <tr>
2030 <th>4</th> <td>0.0890105</td> <td>0.000400669</td> <td>5.74184e-06</td>
2031 </tr>
2032 <tr>
2033 <th>5</th> <td>0.0445032</td> <td>0.000100218</td> <td>7.17799e-07</td>
2034 </tr>
2035 <tr>
2036 <th>6</th> <td>0.0222513</td> <td>2.50576e-05</td> <td>9.0164e-08</td>
2037 </tr>
2038 <tr>
2039 <th></th> <th>@f$O(h)@f$</th> <th>@f$O(h^2)@f$</th> <th>@f$O(h^3)@f$</th>
2040 </tr>
2041</table>
2042
2043The theoretically expected convergence orders are very nicely reflected by the
2044experimentally observed ones indicated in the last row of the table.
2045
2046
2047
2048One can make the same experiment with the @f$L_2@f$ error
2049in the velocity variables:
2050<table align="center" class="doxtable">
2051 <tr>
2052 <th></th>
2053 <th colspan="3" align="center">Finite element order</th>
2054 </tr>
2055 <tr>
2056 <th>Refinement level</th>
2057 <th>0</th>
2058 <th>1</th>
2059 <th>2</th>
2060 </tr>
2061 <tr>
2062 <th>0</th> <td>0.367423</td> <td>0.127657</td> <td>5.10388e-14</td>
2063 </tr>
2064 <tr>
2065 <th>1</th> <td>0.175891</td> <td>0.0319142</td> <td>9.04414e-15</td>
2066 </tr>
2067 <tr>
2068 <th>2</th> <td>0.0869402</td> <td>0.00797856</td> <td>1.23723e-14</td>
2069 </tr>
2070 <tr>
2071 <th>3</th> <td>0.0433435</td> <td>0.00199464</td> <td>1.86345e-07</td>
2072 </tr>
2073 <tr>
2074 <th>4</th> <td>0.0216559</td> <td>0.00049866</td> <td>2.72566e-07</td>
2075 </tr>
2076 <tr>
2077 <th>5</th> <td>0.010826</td> <td>0.000124664</td> <td>3.57141e-07</td>
2078 </tr>
2079 <tr>
2080 <th>6</th> <td>0.00541274</td> <td>3.1166e-05</td> <td>4.46124e-07</td>
2081 </tr>
2082 <tr>
2083 <th></th> <td>@f$O(h)@f$</td> <td>@f$O(h^2)@f$</td> <td>@f$O(h^3)@f$</td>
2084 </tr>
2085</table>
2086The result concerning the convergence order is the same here.
2087
2088
2089
2090<a name="extensions"></a>
2091<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2092
2093
2094<a name="Morerealisticpermeabilityfields"></a><h4>More realistic permeability fields</h4>
2095
2096
2097Realistic flow computations for ground water or oil reservoir simulations will
2098not use a constant permeability. Here's a first, rather simple way to change
2099this situation: we use a permeability that decays very rapidly away from a
2100central flowline until it hits a background value of 0.001. This is to mimic
2101the behavior of fluids in sandstone: in most of the domain, the sandstone is
2102homogeneous and, while permeable to fluids, not overly so; on the other stone,
2103the stone has cracked, or faulted, along one line, and the fluids flow much
2104easier along this large crack. Here is how we could implement something like
2105this:
2106@code
2107template <int dim>
2108void
2109KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
2110 std::vector<Tensor<2,dim> > &values) const
2111{
2112 AssertDimension (points.size(), values.size());
2113
2114 for (unsigned int p=0; p<points.size(); ++p)
2115 {
2116 values[p].clear ();
2117
2118 const double distance_to_flowline
2119 = std::fabs(points[p][1]-0.2*std::sin(10*points[p][0]));
2120
2121 const double permeability = std::max(std::exp(-(distance_to_flowline*
2122 distance_to_flowline)
2123 / (0.1 * 0.1)),
2124 0.001);
2125
2126 for (unsigned int d=0; d<dim; ++d)
2127 values[p][d][d] = 1./permeability;
2128 }
2129}
2130@endcode
2131Remember that the function returns the inverse of the permeability tensor.
2132
2133
2134
2135With a significantly higher mesh resolution, we can visualize this, here with
2136x- and y-velocity:
2137
2138<table style="width:60%" align="center">
2139 <tr>
2140 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-wiggle.png" alt=""></td>
2141 <td><img src="https://www.dealii.org/images/steps/developer/step-20.v-wiggle.png" alt=""></td>
2142 </tr>
2143</table>
2144
2145It is obvious how fluids flow essentially only along the middle line, and not
2146anywhere else.
2147
2148
2149
2150Another possibility would be to use a random permeability field. A simple way
2151to achieve this would be to scatter a number of centers around the domain and
2152then use a permeability field that is the sum of (negative) exponentials for
2153each of these centers. Flow would then try to hop from one center of high
2154permeability to the next one. This is an entirely unscientific attempt at
2155describing a random medium, but one possibility to implement this behavior
2156would look like this:
2157@code
2158template <int dim>
2159class KInverse : public TensorFunction<2,dim>
2160{
2161 public:
2162 KInverse ();
2163
2164 virtual void value_list (const std::vector<Point<dim> > &points,
2165 std::vector<Tensor<2,dim> > &values) const;
2166
2167 private:
2168 std::vector<Point<dim> > centers;
2169};
2170
2171
2172template <int dim>
2173KInverse<dim>::KInverse ()
2174{
2175 const unsigned int N = 40;
2176 centers.resize (N);
2177 for (unsigned int i=0; i<N; ++i)
2178 for (unsigned int d=0; d<dim; ++d)
2179 centers[i][d] = 2.*rand()/RAND_MAX-1;
2180}
2181
2182
2183template <int dim>
2184void
2185KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
2186 std::vector<Tensor<2,dim> > &values) const
2187{
2188 AssertDimension (points.size(), values.size());
2189
2190 for (unsigned int p=0; p<points.size(); ++p)
2191 {
2192 values[p].clear ();
2193
2194 double permeability = 0;
2195 for (unsigned int i=0; i<centers.size(); ++i)
2196 permeability += std::exp(-(points[p] - centers[i]).norm_square() / (0.1 * 0.1));
2197
2198 const double normalized_permeability
2199 = std::max(permeability, 0.005);
2200
2201 for (unsigned int d=0; d<dim; ++d)
2202 values[p][d][d] = 1./normalized_permeability;
2203 }
2204}
2205@endcode
2206
2207A piecewise constant interpolation of the diagonal elements of the
2208inverse of this tensor (i.e., of <code>normalized_permeability</code>)
2209looks as follows:
2210
2211<img src="https://www.dealii.org/images/steps/developer/step-20.k-random.png" alt="">
2212
2213
2214With a permeability field like this, we would get x-velocities and pressures as
2215follows:
2216
2217<table style="width:60%" align="center">
2218 <tr>
2219 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-random.png" alt=""></td>
2220 <td><img src="https://www.dealii.org/images/steps/developer/step-20.p-random.png" alt=""></td>
2221 </tr>
2222</table>
2223
2224We will use these permeability fields again in @ref step_21 "step-21" and @ref step_43 "step-43".
2225
2226
2227<a name="Betterlinearsolvers"></a><h4>Better linear solvers</h4>
2228
2229
2230As mentioned in the introduction, the Schur complement solver used here is not
2231the best one conceivable (nor is it intended to be a particularly good
2232one). Better ones can be found in the literature and can be built using the
2233same block matrix techniques that were introduced here. We pick up on this
2234theme again in @ref step_22 "step-22", where we first build a Schur complement solver for the
2235Stokes equation as we did here, and then in the <a
2236href="step_22.html#improved-solver">Improved Solvers</a> section discuss better
2237ways based on solving the system as a whole but preconditioning based on
2238individual blocks. We will also come back to this in @ref step_43 "step-43".
2239 *
2240 *
2241<a name="PlainProg"></a>
2242<h1> The plain program</h1>
2243@include "step-20.cc"
2244*/
Definition fe_q.h:551
Definition point.h:112
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
Point< 3 > center
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 AssertDimension(dim1, dim2)
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
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
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.
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
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()