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