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