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