732 *
const unsigned int )
const
740 *
class RightHandSide :
public Function<dim>
744 *
const unsigned int component = 0)
const override;
750 *
double RightHandSide<dim>::value(
const Point<dim> &p,
751 *
const unsigned int )
const
761 * The
class that
implements the exact pressure solution has an
762 * oddity in that we implement it as a vector-valued one with two
763 * components. (We say that it has two components in the constructor
764 * where we call the constructor of the base Function class.) In the
765 * `
value()` function, we do not test for the
value of the
766 * `component` argument, which implies that we return the same
value
767 * for both components of the vector-valued function. We do this
768 * because we describe the finite element in use in this program as
769 * a vector-valued system that contains the interior and the
770 * interface pressures, and when we compute errors, we will want to
771 * use the same pressure solution to test both of these components.
775 * class ExactPressure :
public Function<dim>
783 *
const unsigned int component)
const override;
789 *
double ExactPressure<dim>::value(
const Point<dim> &p,
790 *
const unsigned int )
const
818 *
return return_value;
826 * <a name=
"step_61-WGDarcyEquationclassimplementation"></a>
827 * <h3>WGDarcyEquation
class implementation</h3>
832 * <a name=
"step_61-WGDarcyEquationWGDarcyEquation"></a>
833 * <h4>WGDarcyEquation::WGDarcyEquation</h4>
837 * In
this constructor, we create a finite element space
for vector valued
838 *
functions, which will here include the ones used
for the interior and
839 *
interface pressures, @f$p^\circ@f$ and @f$p^\partial@f$.
843 * WGDarcyEquation<dim>::WGDarcyEquation(
const unsigned int degree)
855 * <a name=
"step_61-WGDarcyEquationmake_grid"></a>
856 * <h4>WGDarcyEquation::make_grid</h4>
860 * We generate a mesh on the unit square domain and
refine it.
864 *
void WGDarcyEquation<dim>::make_grid()
869 * std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
880 * <a name=
"step_61-WGDarcyEquationsetup_system"></a>
881 * <h4>WGDarcyEquation::setup_system</h4>
885 * After we have created the mesh above, we distribute degrees of
886 * freedom and resize matrices and vectors. The only piece of
887 * interest in
this function is how we
interpolate the boundary
888 *
values for the pressure. Since the pressure consists of interior
889 * and
interface components, we need to make sure that we only
890 *
interpolate onto that component of the vector-valued solution
891 * space that corresponds to the
interface pressures (as these are
892 * the only ones that are defined on the boundary of the domain). We
893 *
do this via a component mask
object for only the interface
898 *
void WGDarcyEquation<dim>::setup_system()
900 * dof_handler.distribute_dofs(fe);
901 * dof_handler_dgrt.distribute_dofs(fe_dgrt);
903 * std::cout <<
" Number of pressure degrees of freedom: "
904 * << dof_handler.n_dofs() << std::endl;
906 * solution.reinit(dof_handler.n_dofs());
907 * system_rhs.reinit(dof_handler.n_dofs());
911 * constraints.clear();
917 * BoundaryValues<dim>(),
919 * interface_pressure_mask);
920 * constraints.close();
926 * In the bilinear form, there is no integration term over faces
927 * between two neighboring cells, so we can just use
934 * sparsity_pattern.copy_from(dsp);
936 * system_matrix.reinit(sparsity_pattern);
944 * <a name=
"step_61-WGDarcyEquationassemble_system"></a>
945 * <h4>WGDarcyEquation::assemble_system</h4>
949 * This function is more interesting. As detailed in the
950 * introduction, the assembly of the linear system
requires us to
952 * element in the Raviart-Thomas space. As a consequence, we need to
953 * define a Raviart-Thomas finite element object, and have
FEValues
954 * objects that evaluate it at quadrature points. We then need to
955 * compute the
matrix @f$C^
K@f$ on every cell @f$K@f$,
for which we need the
956 * matrices @f$M^
K@f$ and @f$G^
K@f$ mentioned in the introduction.
960 * A
point that may not be obvious is that in all previous tutorial
962 * iterator from a
DoFHandler. This is so that one can call
963 * functions such as
FEValuesBase::get_function_values() that
964 * extract the values of a finite element function (represented by a
965 * vector of DoF values) on the quadrature points of a cell. For
966 * this operation to work, one needs to know which vector elements
967 * correspond to the degrees of freedom on a given cell -- i.e.,
968 * exactly the kind of information and operation provided by the
973 * We could create a
DoFHandler object for the "broken" Raviart-Thomas space
975 * least in the current function, we have no need for any globally defined
976 * degrees of freedom associated with this broken space, but really only
977 * need to reference the shape functions of such a space on the current
978 * cell. As a consequence, we use the fact that one can call
981 * can of course only provide us with information that only
982 * references information about cells, rather than degrees of freedom
983 * enumerated on these cells. So we can't use
985 *
FEValues::shape_value() to obtain the values of shape functions
986 * at quadrature points on the current cell. It is this kind of
987 * functionality we will make use of below. The variable that will
988 * give us this information about the Raviart-Thomas functions below
989 * is then the `fe_values_rt` (and corresponding `fe_face_values_rt`)
994 * Given this introduction, the following declarations should be
999 *
void WGDarcyEquation<dim>::assemble_system()
1001 *
const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1002 *
const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1005 * quadrature_formula,
1009 * face_quadrature_formula,
1015 * quadrature_formula,
1020 * face_quadrature_formula,
1027 *
const unsigned int dofs_per_cell_dgrt = fe_dgrt.n_dofs_per_cell();
1029 *
const unsigned int n_q_points = fe_values.get_quadrature().size();
1030 *
const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1032 *
const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
1034 * RightHandSide<dim> right_hand_side;
1035 * std::vector<double> right_hand_side_values(n_q_points);
1037 *
const Coefficient<dim> coefficient;
1038 * std::vector<Tensor<2, dim>> coefficient_values(n_q_points);
1040 * std::vector<types::global_dof_index> local_dof_indices(
dofs_per_cell);
1045 * Next, let us declare the various cell matrices discussed in the
1059 * @p face component of the shape
functions.
1068 * This
finally gets us in position to
loop over all cells. On
1069 * each cell, we will
first calculate the various cell matrices
1070 * used to construct the local
matrix -- as they depend on the
1071 * cell in question, they need to be re-computed on each cell. We
1072 * need shape
functions for the Raviart-Thomas space as well,
for
1073 * which we need to create
first an iterator to the cell of the
1074 *
triangulation, which we can obtain by assignment from the cell
1078 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1080 * fe_values.
reinit(cell);
1084 * fe_values_dgrt.reinit(cell_dgrt);
1086 * right_hand_side.value_list(fe_values.get_quadrature_points(),
1087 * right_hand_side_values);
1088 * coefficient.value_list(fe_values.get_quadrature_points(),
1089 * coefficient_values);
1093 * The
first cell
matrix we will compute is the @ref GlossMassMatrix
"mass matrix"
1094 *
for the Raviart-Thomas space. Hence, we need to
loop over
1098 * cell_matrix_M = 0;
1099 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1100 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1102 *
const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
1103 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1106 * fe_values_dgrt[velocities].value(k, q);
1107 * cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
1112 * Next we take the inverse of
this matrix by
using
1114 * the coefficient
matrix @f$C^
K@f$ later. It is worth recalling
1115 * later that `cell_matrix_M` actually contains the *inverse*
1116 * of @f$M^
K@f$ after
this call.
1119 * cell_matrix_M.gauss_jordan();
1123 * From the introduction, we know that the right hand side
1124 * @f$G^
K@f$ of the equation that defines @f$C^
K@f$ is the difference
1125 * between a face integral and a cell integral. Here, we
1127 * interior. Each component of
this matrix is the integral of
1128 * a product between a basis function of the polynomial space
1129 * and the divergence of a basis function of the
1130 * Raviart-Thomas space. These basis
functions are defined in
1134 * cell_matrix_G = 0;
1135 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1136 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1138 *
const double div_v_i =
1139 * fe_values_dgrt[velocities].divergence(i, q);
1142 *
const double phi_j_interior =
1143 * fe_values[pressure_interior].value(j, q);
1145 * cell_matrix_G(i, j) -=
1146 * (div_v_i * phi_j_interior * fe_values.JxW(q));
1154 * Each component is the integral of a product between a basis function
1155 * of the polynomial space and the dot product of a basis function of
1156 * the Raviart-Thomas space and the normal vector. So we
loop over all
1157 * the faces of the element and obtain the normal vector.
1160 *
for (
const auto &face : cell->face_iterators())
1162 * fe_face_values.
reinit(cell, face);
1163 * fe_face_values_dgrt.reinit(cell_dgrt, face);
1165 *
for (
unsigned int q = 0; q < n_face_q_points; ++q)
1167 *
const Tensor<1, dim> &normal = fe_face_values.normal_vector(q);
1169 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1172 * fe_face_values_dgrt[velocities].value(i, q);
1175 *
const double phi_j_face =
1176 * fe_face_values[pressure_face].value(j, q);
1178 * cell_matrix_G(i, j) +=
1179 * ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
1187 * @p cell_matrix_C is then the
matrix product between the
1189 * (where
this inverse is stored in @p cell_matrix_M):
1192 * cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
1196 * Finally we can compute the local
matrix @f$A^
K@f$. Element
1197 * @f$A^K_{ij}@f$ is given by @f$\int_{
E} \sum_{k,
l} C_{ik} C_{jl}
1198 * (\mathbf{
K} \mathbf{v}_k) \cdot \mathbf{v}_l
1199 * \mathrm{
d}x@f$. We have calculated the coefficients @f$C@f$ in
1200 * the previous step, and so obtain the following after
1201 * suitably re-arranging the loops:
1205 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1207 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1210 * fe_values_dgrt[velocities].value(k, q);
1211 *
for (
unsigned int l = 0;
l < dofs_per_cell_dgrt; ++
l)
1214 * fe_values_dgrt[velocities].value(l, q);
1216 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1217 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1218 * local_matrix(i, j) +=
1219 * (coefficient_values[q] * cell_matrix_C[i][k] * v_k) *
1220 * cell_matrix_C[j][
l] * v_l * fe_values_dgrt.JxW(q);
1227 * Next, we calculate the right hand side, @f$\int_{
K} f q \mathrm{
d}x@f$:
1231 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1232 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1234 * cell_rhs(i) += (fe_values[pressure_interior].value(i, q) *
1235 * right_hand_side_values[q] * fe_values.JxW(q));
1240 * The last step is to distribute components of the local
1241 *
matrix into the system
matrix and transfer components of
1242 * the cell right hand side into the system right hand side:
1245 * cell->get_dof_indices(local_dof_indices);
1246 * constraints.distribute_local_to_global(
1247 * local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1256 * <a name=
"step_61-WGDarcyEquationdimsolve"></a>
1257 * <h4>WGDarcyEquation<dim>::solve</h4>
1261 * This step is rather trivial and the same as in many previous
1262 * tutorial programs:
1265 *
template <
int dim>
1266 *
void WGDarcyEquation<dim>::solve()
1268 *
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
1271 * constraints.distribute(solution);
1278 * <a name=
"step_61-WGDarcyEquationdimcompute_postprocessed_velocity"></a>
1279 * <h4>WGDarcyEquation<dim>::compute_postprocessed_velocity</h4>
1283 * In
this function, compute the velocity field from the pressure
1284 * solution previously computed. The
1285 * velocity is defined as @f$\mathbf{u}_h = \mathbf{Q}_h \left(
1286 * -\mathbf{
K}\nabla_{
w,
d}p_h \right)@f$, which
requires us to compute
1287 * many of the same terms as in the assembly of the system
matrix.
1288 * There are also the matrices @f$E^K,D^K@f$ we need to
assemble (see
1289 * the introduction) but they really just follow the same kind of
1294 * Computing the same matrices here as we have already done in the
1295 * `assemble_system()` function is of course wasteful in terms of
1296 * CPU time. Likewise, we
copy some of the code from there to
this
1297 * function, and
this is also generally a poor idea. A better
1298 * implementation might provide
for a function that encapsulates
1299 *
this duplicated code. One could also think of
using the classic
1300 * trade-off between computing efficiency and memory efficiency to
1301 * only compute the @f$C^
K@f$ matrices once per cell during the
1302 * assembly, storing them somewhere on the side, and re-
using them
1303 * here. (This is what @ref step_51
"step-51" does,
for example, where the
1304 * `assemble_system()` function takes an argument that determines
1305 * whether the local matrices are recomputed, and a similar approach
1306 * -- maybe with storing local matrices elsewhere -- could be
1307 * adapted
for the current program.)
1310 *
template <
int dim>
1311 *
void WGDarcyEquation<dim>::compute_postprocessed_velocity()
1313 * darcy_velocity.reinit(dof_handler_dgrt.n_dofs());
1315 *
const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1316 *
const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1319 * quadrature_formula,
1324 * face_quadrature_formula,
1330 * quadrature_formula,
1336 * face_quadrature_formula,
1342 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1343 *
const unsigned int dofs_per_cell_dgrt = fe_dgrt.n_dofs_per_cell();
1345 *
const unsigned int n_q_points = fe_values.get_quadrature().size();
1346 *
const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1348 *
const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
1351 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1352 * std::vector<types::global_dof_index> local_dof_indices_dgrt(
1353 * dofs_per_cell_dgrt);
1364 *
const Coefficient<dim> coefficient;
1365 * std::vector<Tensor<2, dim>> coefficient_values(n_q_points_dgrt);
1373 * In the introduction, we explained how to calculate the numerical velocity
1374 * on the cell. We need the pressure solution
values on each cell,
1375 * coefficients of the Gram
matrix and coefficients of the @f$L_2@f$ projection.
1376 * We have already calculated the global solution, so we will
extract the
1377 * cell solution from the global solution. The coefficients of the Gram
1378 *
matrix have been calculated when we assembled the system
matrix for the
1379 * pressures. We will
do the same way here. For the coefficients of the
1380 * projection, we
do matrix multiplication, i.e., the inverse of the Gram
1381 *
matrix times the
matrix with @f$(\mathbf{
K} \mathbf{
w}, \mathbf{
w})@f$ as
1382 * components. Then, we multiply all these coefficients and call them beta.
1383 * The numerical velocity is the product of beta and the basis functions of
1384 * the Raviart-Thomas space.
1388 * cell = dof_handler.begin_active(),
1389 * endc = dof_handler.end(), cell_dgrt = dof_handler_dgrt.begin_active();
1390 *
for (; cell != endc; ++cell, ++cell_dgrt)
1392 * fe_values.reinit(cell);
1393 * fe_values_dgrt.reinit(cell_dgrt);
1395 * coefficient.value_list(fe_values_dgrt.get_quadrature_points(),
1396 * coefficient_values);
1400 * The component of this <code>cell_matrix_E</code> is the integral of
1401 * @f$(\mathbf{
K} \mathbf{
w}, \mathbf{
w})@f$. <code>cell_matrix_M</code> is
1405 * cell_matrix_M = 0;
1406 * cell_matrix_E = 0;
1407 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1408 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1410 *
const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
1411 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1414 * fe_values_dgrt[velocities].value(k, q);
1416 * cell_matrix_E(i, k) +=
1417 * (coefficient_values[q] * v_i * v_k * fe_values_dgrt.JxW(q));
1419 * cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
1425 * To compute the
matrix @f$D@f$ mentioned in the introduction, we
1426 * then need to evaluate @f$D=M^{-1}
E@f$ as explained in the
1430 * cell_matrix_M.gauss_jordan();
1431 * cell_matrix_M.mmult(cell_matrix_D, cell_matrix_E);
1435 * Then we also need, again, to compute the
matrix @f$C@f$ that is
1436 * used to evaluate the weak discrete
gradient. This is the
1437 * exact same code as used in the assembly of the system
1441 * cell_matrix_G = 0;
1442 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1443 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1445 *
const double div_v_i =
1446 * fe_values_dgrt[velocities].divergence(i, q);
1447 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1449 *
const double phi_j_interior =
1450 * fe_values[pressure_interior].value(j, q);
1452 * cell_matrix_G(i, j) -=
1453 * (div_v_i * phi_j_interior * fe_values.JxW(q));
1457 *
for (
const auto &face : cell->face_iterators())
1459 * fe_face_values.
reinit(cell, face);
1460 * fe_face_values_dgrt.reinit(cell_dgrt, face);
1462 *
for (
unsigned int q = 0; q < n_face_q_points; ++q)
1464 *
const Tensor<1, dim> &normal = fe_face_values.normal_vector(q);
1466 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1469 * fe_face_values_dgrt[velocities].value(i, q);
1470 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1472 *
const double phi_j_face =
1473 * fe_face_values[pressure_face].value(j, q);
1475 * cell_matrix_G(i, j) +=
1476 * ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
1481 * cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
1485 * Finally, we need to
extract the pressure unknowns that
1486 * correspond to the current cell:
1489 * cell->get_dof_values(solution, cell_solution);
1493 * We are now in a position to compute the local velocity
1494 * unknowns (with respect to the Raviart-Thomas space we are
1495 * projecting the term @f$-\mathbf K \nabla_{
w,
d} p_h@f$ into):
1498 * cell_velocity = 0;
1499 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1500 *
for (
unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
1501 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1502 * cell_velocity(k) +=
1503 * -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
1507 * We compute Darcy velocity.
1508 * This is same as cell_velocity but used to graph Darcy velocity.
1511 * cell_dgrt->get_dof_indices(local_dof_indices_dgrt);
1512 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1513 *
for (
unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
1514 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1515 * darcy_velocity(local_dof_indices_dgrt[k]) +=
1516 * -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
1525 * <a name=
"step_61-WGDarcyEquationdimcompute_pressure_error"></a>
1526 * <h4>WGDarcyEquation<dim>::compute_pressure_error</h4>
1530 * This part is to calculate the @f$L_2@f$ error of the pressure. We
1531 * define a vector that holds the
norm of the error on each cell.
1532 * Next, we use VectorTool::integrate_difference() to compute the
1533 * error in the @f$L_2@f$ norm on each cell. However, we really only
1534 * care about the error in the interior component of the solution
1535 * vector (we can't even evaluate the interface pressures at the
1536 * quadrature points because these are all located in the interior
1537 * of cells) and consequently have to use a weight function that
1538 * ensures that the interface component of the solution variable is
1540 * arguments indicate which component we want to select (component
1541 * zero, i.e., the interior pressures) and how many components there
1542 * are in total (two).
1545 * template <
int dim>
1546 *
void WGDarcyEquation<dim>::compute_pressure_error()
1552 * ExactPressure<dim>(),
1553 * difference_per_cell,
1556 * &select_interior_pressure);
1558 *
const double L2_error = difference_per_cell.l2_norm();
1559 * std::cout <<
"L2_error_pressure " << L2_error << std::endl;
1567 * <a name=
"step_61-WGDarcyEquationdimcompute_velocity_error"></a>
1568 * <h4>WGDarcyEquation<dim>::compute_velocity_error</h4>
1572 * In
this function, we evaluate @f$L_2@f$ errors
for the velocity on
1573 * each cell, and @f$L_2@f$ errors
for the flux on faces. The function
1574 * relies on the `compute_postprocessed_velocity()` function having
1575 * previous computed, which computes the velocity field based on the
1576 * pressure solution that has previously been computed.
1580 * We are going to evaluate velocities on each cell and calculate
1581 * the difference between numerical and exact velocities.
1584 *
template <
int dim>
1585 *
void WGDarcyEquation<dim>::compute_velocity_errors()
1587 *
const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1588 *
const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1591 * quadrature_formula,
1597 * face_quadrature_formula,
1603 *
const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1604 *
const unsigned int n_face_q_points_dgrt =
1605 * fe_face_values_dgrt.get_quadrature().size();
1607 * std::vector<Tensor<1, dim>> velocity_values(n_q_points_dgrt);
1608 * std::vector<Tensor<1, dim>> velocity_face_values(n_face_q_points_dgrt);
1612 *
const ExactVelocity<dim> exact_velocity;
1614 *
double L2_err_velocity_cell_sqr_global = 0;
1615 *
double L2_err_flux_sqr = 0;
1619 * Having previously computed the postprocessed velocity, we here
1620 * only have to
extract the corresponding
values on each cell and
1621 * face and compare it to the exact
values.
1624 *
for (
const auto &cell_dgrt : dof_handler_dgrt.active_cell_iterators())
1626 * fe_values_dgrt.
reinit(cell_dgrt);
1630 * First compute the @f$L_2@f$ error between the postprocessed velocity
1631 * field and the exact one:
1634 * fe_values_dgrt[velocities].get_function_values(darcy_velocity,
1636 *
double L2_err_velocity_cell_sqr_local = 0;
1637 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1641 * exact_velocity.value(fe_values_dgrt.quadrature_point(q));
1643 * L2_err_velocity_cell_sqr_local +=
1644 * ((velocity - true_velocity) * (velocity - true_velocity) *
1645 * fe_values_dgrt.JxW(q));
1647 * L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local;
1651 * For reconstructing the flux we need the size of cells and
1652 * faces. Since fluxes are calculated on faces, we have the
1653 *
loop over all four faces of each cell. To calculate the
1654 * face velocity, we
extract values at the quadrature points from the
1655 * `darcy_velocity` which we have computed previously. Then, we
1656 * calculate the squared velocity error in normal direction. Finally, we
1657 * calculate the @f$L_2@f$ flux error on the cell by appropriately scaling
1658 * with face and cell areas and add it to the global error.
1661 *
const double cell_area = cell_dgrt->measure();
1662 *
for (
const auto &face_dgrt : cell_dgrt->face_iterators())
1664 * const double face_length = face_dgrt->measure();
1665 * fe_face_values_dgrt.reinit(cell_dgrt, face_dgrt);
1666 * fe_face_values_dgrt[velocities].get_function_values(
1667 * darcy_velocity, velocity_face_values);
1669 *
double L2_err_flux_face_sqr_local = 0;
1670 *
for (
unsigned int q = 0; q < n_face_q_points_dgrt; ++q)
1674 * exact_velocity.value(fe_face_values_dgrt.quadrature_point(q));
1677 * fe_face_values_dgrt.normal_vector(q);
1679 * L2_err_flux_face_sqr_local +=
1680 * ((velocity * normal - true_velocity * normal) *
1681 * (velocity * normal - true_velocity * normal) *
1682 * fe_face_values_dgrt.JxW(q));
1684 *
const double err_flux_each_face =
1685 * L2_err_flux_face_sqr_local / face_length * cell_area;
1686 * L2_err_flux_sqr += err_flux_each_face;
1692 * After adding up errors over all cells and faces, we take the
1693 * square root and get the @f$L_2@f$ errors of velocity and
1694 * flux. These we output to screen.
1697 *
const double L2_err_velocity_cell =
1698 *
std::sqrt(L2_err_velocity_cell_sqr_global);
1699 *
const double L2_err_flux_face =
std::sqrt(L2_err_flux_sqr);
1701 * std::cout <<
"L2_error_vel: " << L2_err_velocity_cell << std::endl
1702 * <<
"L2_error_flux: " << L2_err_flux_face << std::endl;
1709 * <a name=
"step_61-WGDarcyEquationoutput_results"></a>
1710 * <h4>WGDarcyEquation::output_results</h4>
1714 * We have two sets of results to output: the interior solution and
1715 * the skeleton solution. We use <code>
DataOut</code> to visualize
1716 * interior results. The graphical output
for the skeleton results
1721 * In both of the output files, both the interior and the face
1722 * variables are stored. For the
interface output, the output file
1723 * simply contains the interpolation of the interior pressures onto
1724 * the faces, but because it is undefined which of the two interior
1725 * pressure variables you get from the two adjacent cells, it is
1726 * best to ignore the interior pressure in the
interface output
1727 * file. Conversely,
for the cell interior output file, it is of
1728 * course impossible to show any
interface pressures @f$p^\partial@f$,
1729 * because these are only available on interfaces and not cell
1730 * interiors. Consequently, you will see them shown as an
invalid
1731 *
value (such as an infinity).
1735 * For the cell interior output, we also want to output the velocity
1736 * variables. This is a bit tricky since it lives on the same mesh
1737 * but uses a different
DoFHandler object (the pressure variables live
1738 * on the `dof_handler`
object, the Darcy velocity on the `dof_handler_dgrt`
1739 *
object). Fortunately, there are variations of the
1741 *
DoFHandler a vector corresponds to, and consequently we can visualize
1742 * the data from both
DoFHandler objects within the same file.
1745 * template <
int dim>
1746 *
void WGDarcyEquation<dim>::output_results() const
1753 * First attach the pressure solution to the
DataOut object:
1756 *
const std::vector<std::string> solution_names = {
"interior_pressure",
1757 *
"interface_pressure"};
1762 * Then
do the same with the Darcy velocity field, and
continue
1763 * with writing everything out into a file.
1766 *
const std::vector<std::string> velocity_names(dim,
"velocity");
1767 *
const std::vector<
1769 * velocity_component_interpretation(
1771 * data_out.add_data_vector(dof_handler_dgrt,
1774 * velocity_component_interpretation);
1776 * data_out.build_patches(fe.degree);
1777 * std::ofstream output(
"solution_interior.vtu");
1778 * data_out.write_vtu(output);
1783 * data_out_faces.attach_dof_handler(dof_handler);
1784 * data_out_faces.add_data_vector(solution,
"Pressure_Face");
1785 * data_out_faces.build_patches(fe.degree);
1786 * std::ofstream face_output(
"solution_interface.vtu");
1787 * data_out_faces.write_vtu(face_output);
1795 * <a name=
"step_61-WGDarcyEquationrun"></a>
1796 * <h4>WGDarcyEquation::run</h4>
1800 * This is the
final function of the main
class. It calls the other
functions
1804 *
template <
int dim>
1805 *
void WGDarcyEquation<dim>::run()
1807 * std::cout <<
"Solving problem in " << dim <<
" space dimensions."
1811 * assemble_system();
1813 * compute_postprocessed_velocity();
1814 * compute_pressure_error();
1815 * compute_velocity_errors();
1825 * <a name=
"step_61-Thecodemaincodefunction"></a>
1826 * <h3>The <code>main</code> function</h3>
1830 * This is the main function. We can change the dimension here to
run in 3
d.
1837 * Step61::WGDarcyEquation<2> wg_darcy(0);
1840 *
catch (std::exception &exc)
1842 * std::cerr << std::endl
1844 * <<
"----------------------------------------------------"
1846 * std::cerr <<
"Exception on processing: " << std::endl
1847 * << exc.what() << std::endl
1848 * <<
"Aborting!" << std::endl
1849 * <<
"----------------------------------------------------"
1855 * std::cerr << std::endl
1857 * <<
"----------------------------------------------------"
1859 * std::cerr <<
"Unknown exception!" << std::endl
1860 * <<
"Aborting!" << std::endl
1861 * <<
"----------------------------------------------------"
1869<a name=
"step_61-Results"></a><h1>Results</h1>
1872We
run the program with a right hand side that will produce the
1873solution @f$p = \sin(\pi x) \sin(\pi y)@f$ and with homogeneous Dirichlet
1874boundary conditions in the domain @f$\Omega = (0,1)^2@f$. In addition, we
1875choose the coefficient matrix in the differential
operator
1877@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ and
1878@f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ element combinations, which one can
1879select by
using the appropriate constructor argument
for the
1880`WGDarcyEquation`
object in `main()`. We will then visualize pressure
1881values in interiors of cells and on faces. We want to see that the
1882pressure maximum is around 1 and the minimum is around 0. With mesh
1883refinement, the convergence rates of pressure, velocity and flux
1884should then be around 1
for @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ , 2
for
1885@f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$, and 3
for
1886@f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$.
1889<a name=
"step_61-TestresultsoniWGQsub0subQsub0subRTsub0subi"></a><h3>Test results on <i>WG(Q<sub>0</sub>,Q<sub>0</sub>;RT<sub>[0]</sub>)</i></h3>
1892The following figures show interior pressures and face pressures
using the
1893@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ element. The mesh is refined 2 times (top)
1894and 4 times (bottom), respectively. (This number can be adjusted in the
1895`make_grid()` function.) When the mesh is coarse, one can see
1896the face pressures @f$p^\partial@f$ neatly between the
values of the interior
1897pressures @f$p^\circ@f$ on the two adjacent cells.
1899<table align=
"center">
1901 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_2d_2.png" alt=
""></td>
1902 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_3d_2.png" alt=
""></td>
1905 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_2d_4.png" alt=
""></td>
1906 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_3d_4.png" alt=
""></td>
1910From the figures, we can see that with the mesh refinement, the maximum and
1911minimum pressure
values are approaching the
values we expect.
1912Since the mesh is a rectangular mesh and
numbers of cells in each direction is even, we
1913have
symmetric solutions. From the 3
d figures on the right,
1914we can see that on @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, the pressure is a constant
1915in the interior of the cell, as expected.
1917<a name=
"step_61-Convergencetableforik0i"></a><h4>Convergence table for <i>k=0</i></h4>
1920We run the code with differently refined meshes (chosen in the `make_grid()` function)
1921and get the following convergence rates of pressure,
1922velocity, and flux (as defined in the introduction).
1924<table align=
"center" class=
"doxtable">
1926 <th>number of refinements </th><th> @f$\|p-p_h^\circ\|@f$ </th><th> @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
1929 <td> 2 </td><td> 1.587e-01 </td><td> 5.113e-01 </td><td> 7.062e-01 </td>
1932 <td> 3 </td><td> 8.000e-02 </td><td> 2.529e-01 </td><td> 3.554e-01 </td>
1935 <td> 4 </td><td> 4.006e-02 </td><td> 1.260e-01 </td><td> 1.780e-01 </td>
1938 <td> 5 </td><td> 2.004e-02 </td><td> 6.297e-02 </td><td> 8.902e-02 </td>
1941 <th>Conv.rate </th><th> 1.00 </th><th> 1.00 </th><th> 1.00 </th>
1945We can see that the convergence rates of @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ are around 1.
1946This, of course, matches our theoretical expectations.
1949<a name=
"step_61-TestresultsoniWGQsub1subQsub1subRTsub1subi"></a><h3>Test results on <i>WG(Q<sub>1</sub>,Q<sub>1</sub>;RT<sub>[1]</sub>)</i></h3>
1952We can repeat the experiment from above
using the next higher polynomial
1954The following figures are interior pressures and face pressures implemented using
1955@f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$. The mesh is refined 4 times. Compared to the
1956previous figures
using
1957@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, on each cell, the solution is no longer constant
1958on each cell, as we now use bilinear polynomials to
do the approximation.
1959Consequently, there are 4 pressure values in one interior, 2 pressure values on
1962<table align=
"center">
1964 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg111_2d_4.png" alt=
""></td>
1965 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg111_3d_4.png" alt=
""></td>
1969Compared to the corresponding image
for the @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$
1970combination, the solution is now substantially more accurate and, in
1971particular so close to being continuous at the interfaces that we can
1972no longer distinguish the interface pressures @f$p^\partial@f$ from the
1973interior pressures @f$p^\circ@f$ on the adjacent cells.
1975<a name=
"step_61-Convergencetableforik1i"></a><h4>Convergence table for <i>k=1</i></h4>
1978The following are the convergence rates of pressure, velocity, and flux
1979we obtain from
using the @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ element combination:
1981<table align=
"center" class=
"doxtable">
1983 <th>number of refinements </th><th> @f$\|p-p_h^\circ\|@f$ </th><th> @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
1986 <td> 2 </td><td> 1.613e-02 </td><td> 5.093e-02 </td><td> 7.167e-02 </td>
1989 <td> 3 </td><td> 4.056e-03 </td><td> 1.276e-02 </td><td> 1.802e-02 </td>
1992 <td> 4 </td><td> 1.015e-03 </td><td> 3.191e-03 </td><td> 4.512e-03 </td>
1995 <td> 5 </td><td> 2.540e-04 </td><td> 7.979e-04 </td><td> 1.128e-03 </td>
1998 <th>Conv.rate </th><th> 2.00 </th><th> 2.00 </th><th> 2.00 </th>
2002The convergence rates of @f$WG(Q_1,Q_1;RT_{[1]})@f$ are around 2, as expected.
2006<a name=
"step_61-TestresultsoniWGQsub2subQsub2subRTsub2subi"></a><h3>Test results on <i>WG(Q<sub>2</sub>,Q<sub>2</sub>;RT<sub>[2]</sub>)</i></h3>
2009Let us go one polynomial degree higher.
2010The following are interior pressures and face pressures implemented
using
2011@f$WG(Q_2,Q_2;RT_{[2]})@f$, with mesh size @f$h = 1/32@f$ (i.e., 5 global mesh
2012refinement steps). In the program, we use
2013`data_out_face.build_patches(fe.degree)` when generating graphical output
2015we divide each 2
d cell interior into 4 subcells in order to provide a better
2016visualization of the quadratic polynomials.
2017<table align=
"center">
2019 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg222_2d_5.png" alt=
""></td>
2020 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg222_3d_5.png" alt=
""></td>
2025<a name=
"step_61-Convergencetableforik2i"></a><h4>Convergence table for <i>k=2</i></h4>
2028As before, we can generate convergence data
for the
2029@f$L_2@f$ errors of pressure, velocity, and flux
2030using the @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ combination:
2032<table align=
"center" class=
"doxtable">
2034 <th>number of refinements </th><th> @f$\|p-p_h^\circ\|@f$ </th><th> @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
2037 <td> 2 </td><td> 1.072e-03 </td><td> 3.375e-03 </td><td> 4.762e-03 </td>
2040 <td> 3 </td><td> 1.347e-04 </td><td> 4.233e-04 </td><td> 5.982e-04 </td>
2043 <td> 4 </td><td> 1.685e-05 </td><td> 5.295e-05 </td><td> 7.487e-05 </td>
2046 <td> 5 </td><td> 2.107e-06 </td><td> 6.620e-06 </td><td> 9.362e-06 </td>
2049 <th>Conv.rate </th><th> 3.00 </th><th> 3.00 </th><th> 3.00 </th>
2053Once more, the convergence rates of @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ is
2054as expected, with values around 3.
2057<a name=
"step_61-PlainProg"></a>
2058<h1> The plain program</h1>
2059@include
"step-61.cc"
std::vector< bool > component_mask
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
const unsigned int dofs_per_cell
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
const Quadrature< dim > quadrature
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual value_type value(const Point< dim > &p) const
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
DataComponentInterpretation
@ component_is_part_of_vector
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static constexpr double PI
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation