729 *
const unsigned int )
const
737 *
class RightHandSide :
public Function<dim>
741 *
const unsigned int component = 0)
const override;
748 *
const unsigned int )
const
758 * The
class that
implements the exact pressure solution has an
759 * oddity in that we implement it as a vector-valued
one with two
760 * components. (We say that it has two components in the constructor
761 * where we call the constructor of the base Function class.) In the
762 * `
value()` function, we do not test for the
value of the
763 * `component` argument, which implies that we return the same
value
764 * for both components of the vector-valued function. We do this
765 * because we describe the finite element in use in this program as
766 * a vector-valued system that contains the interior and the
767 * interface pressures, and when we compute errors, we will want to
768 * use the same pressure solution to test both of these components.
772 * class ExactPressure :
public Function<dim>
780 *
const unsigned int component)
const override;
787 *
const unsigned int )
const
815 *
return return_value;
823 * <a name=
"WGDarcyEquationclassimplementation"></a>
824 * <h3>WGDarcyEquation
class implementation</h3>
829 * <a name=
"WGDarcyEquationWGDarcyEquation"></a>
830 * <h4>WGDarcyEquation::WGDarcyEquation</h4>
834 * In
this constructor, we create a finite element space
for vector valued
835 *
functions, which will here include the ones used
for the interior and
836 *
interface pressures, @f$p^\circ@f$ and @f$p^\partial@f$.
840 * WGDarcyEquation<dim>::WGDarcyEquation(
const unsigned int degree)
852 * <a name=
"WGDarcyEquationmake_grid"></a>
853 * <h4>WGDarcyEquation::make_grid</h4>
857 * We generate a mesh on the unit square domain and
refine it.
861 *
void WGDarcyEquation<dim>::make_grid()
866 * std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
877 * <a name=
"WGDarcyEquationsetup_system"></a>
878 * <h4>WGDarcyEquation::setup_system</h4>
882 * After we have created the mesh above, we distribute degrees of
883 * freedom and resize matrices and vectors. The only piece of
884 * interest in
this function is how we
interpolate the boundary
885 * values
for the pressure. Since the pressure consists of interior
886 * and
interface components, we need to make sure that we only
887 *
interpolate onto that component of the vector-valued solution
888 * space that corresponds to the
interface pressures (as these are
889 * the only ones that are defined on the boundary of the domain). We
890 *
do this via a component mask
object for only the interface
895 *
void WGDarcyEquation<dim>::setup_system()
897 * dof_handler.distribute_dofs(fe);
898 * dof_handler_dgrt.distribute_dofs(fe_dgrt);
900 * std::cout <<
" Number of pressure degrees of freedom: "
901 * << dof_handler.n_dofs() << std::endl;
903 * solution.reinit(dof_handler.n_dofs());
904 * system_rhs.reinit(dof_handler.n_dofs());
908 * constraints.clear();
914 * BoundaryValues<dim>(),
916 * interface_pressure_mask);
917 * constraints.close();
923 * In the bilinear form, there is no integration term over faces
924 * between two neighboring cells, so we can just use
931 * sparsity_pattern.copy_from(dsp);
933 * system_matrix.reinit(sparsity_pattern);
941 * <a name=
"WGDarcyEquationassemble_system"></a>
942 * <h4>WGDarcyEquation::assemble_system</h4>
946 * This
function is more interesting. As detailed in the
947 * introduction, the assembly of the linear system requires us to
948 * evaluate the weak gradient of the shape
functions, which is an
949 * element in the Raviart-Thomas space. As a consequence, we need to
950 * define a Raviart-Thomas finite element object, and have
FEValues
951 * objects that evaluate it at quadrature points. We then need to
952 * compute the
matrix @f$C^K@f$ on every cell @f$K@f$,
for which we need the
953 * matrices @f$M^K@f$ and @f$G^K@f$ mentioned in the introduction.
957 *
A point that may not be obvious is that in all previous tutorial
961 *
extract the values of a finite element function (represented by a
962 * vector of DoF values) on the quadrature points of a cell. For
963 * this operation to work,
one needs to know which vector elements
964 * correspond to the degrees of freedom on a given cell -- i.
e.,
965 * exactly the kind of information and operation provided by the
970 * We could create a
DoFHandler object for the "broken" Raviart-Thomas space
971 * (using the FE_DGRT class), but we really don't want to here: At
972 * least in the current function, we have no need for any globally defined
973 * degrees of freedom associated with this broken space, but really only
974 * need to reference the shape
functions of such a space on the current
975 * cell. As a consequence, we use the fact that
one can
call
978 * can of course only provide us with information that only
979 * references information about cells, rather than degrees of freedom
980 * enumerated on these cells. So we can't use
983 * at quadrature points on the current cell. It is this kind of
984 * functionality we will make use of below. The variable that will
985 * give us this information about the Raviart-Thomas
functions below
986 * is then the `fe_values_rt` (and corresponding `fe_face_values_rt`)
991 * Given this introduction, the following declarations should be
996 *
void WGDarcyEquation<dim>::assemble_system()
998 *
const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
999 *
const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1002 * quadrature_formula,
1006 * face_quadrature_formula,
1012 * quadrature_formula,
1017 * face_quadrature_formula,
1023 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1024 *
const unsigned int dofs_per_cell_dgrt = fe_dgrt.dofs_per_cell;
1026 *
const unsigned int n_q_points = fe_values.get_quadrature().size();
1027 *
const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1029 *
const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
1031 * RightHandSide<dim> right_hand_side;
1032 * std::vector<double> right_hand_side_values(n_q_points);
1034 *
const Coefficient<dim> coefficient;
1035 * std::vector<Tensor<2, dim>> coefficient_values(n_q_points);
1037 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1042 * Next, let us declare the various cell matrices discussed in the
1056 * @p face component of the shape
functions.
1065 * This
finally gets us in position to
loop over all cells. On
1066 * each cell, we will
first calculate the various cell matrices
1067 * used to construct the local
matrix -- as they depend on the
1068 * cell in question, they need to be re-computed on each cell. We
1069 * need shape
functions for the Raviart-Thomas space as well,
for
1070 * which we need to create
first an iterator to the cell of the
1071 *
triangulation, which we can obtain by assignment from the cell
1075 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1077 * fe_values.reinit(cell);
1081 * fe_values_dgrt.reinit(cell_dgrt);
1083 * right_hand_side.value_list(fe_values.get_quadrature_points(),
1084 * right_hand_side_values);
1085 * coefficient.value_list(fe_values.get_quadrature_points(),
1086 * coefficient_values);
1091 *
for the Raviart-Thomas space. Hence, we need to
loop over
1095 * cell_matrix_M = 0;
1096 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1097 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1099 *
const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
1100 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1103 * fe_values_dgrt[velocities].value(k, q);
1104 * cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
1109 * Next we take the inverse of
this matrix by
using
1111 * the coefficient
matrix @f$C^K@f$ later. It is worth recalling
1112 * later that `cell_matrix_M` actually contains the *inverse*
1113 * of @f$M^K@f$ after
this call.
1116 * cell_matrix_M.gauss_jordan();
1120 * From the introduction, we know that the right hand side
1121 * @f$G^K@f$ of the equation that defines @f$C^K@f$ is the difference
1122 * between a face integral and a cell integral. Here, we
1123 *
approximate the negative of the contribution in the
1124 * interior. Each component of
this matrix is the integral of
1125 * a product between a basis
function of the polynomial space
1126 * and the divergence of a basis
function of the
1127 * Raviart-Thomas space. These basis
functions are defined in
1131 * cell_matrix_G = 0;
1132 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1133 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1135 *
const double div_v_i =
1136 * fe_values_dgrt[velocities].divergence(i, q);
1137 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1139 *
const double phi_j_interior =
1140 * fe_values[pressure_interior].value(j, q);
1142 * cell_matrix_G(i, j) -=
1143 * (div_v_i * phi_j_interior * fe_values.JxW(q));
1151 * Each component is the integral of a product between a basis
function
1152 * of the polynomial space and the dot product of a basis
function of
1153 * the Raviart-Thomas space and the normal vector. So we
loop over all
1154 * the faces of the element and obtain the normal vector.
1157 *
for (
const auto &face : cell->face_iterators())
1159 * fe_face_values.reinit(cell, face);
1160 * fe_face_values_dgrt.reinit(cell_dgrt, face);
1162 *
for (
unsigned int q = 0; q < n_face_q_points; ++q)
1166 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1169 * fe_face_values_dgrt[velocities].value(i, q);
1170 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1172 *
const double phi_j_face =
1173 * fe_face_values[pressure_face].value(j, q);
1175 * cell_matrix_G(i, j) +=
1176 * ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
1184 * @p cell_matrix_C is then the
matrix product between the
1186 * (where
this inverse is stored in @p cell_matrix_M):
1189 * cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
1193 * Finally we can compute the local
matrix @f$A^K@f$. Element
1194 * @f$A^K_{ij}@f$ is given by @f$\int_{
E} \sum_{k,
l} C_{ik} C_{jl}
1195 * (\mathbf{K} \mathbf{v}_k) \cdot \mathbf{v}_l
1196 * \mathrm{
d}x@f$. We have calculated the coefficients @f$C@f$ in
1197 * the previous step, and so obtain the following after
1198 * suitably re-arranging the loops:
1202 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1204 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1207 * fe_values_dgrt[velocities].value(k, q);
1208 *
for (
unsigned int l = 0;
l < dofs_per_cell_dgrt; ++
l)
1211 * fe_values_dgrt[velocities].value(
l, q);
1213 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1214 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1215 * local_matrix(i, j) +=
1216 * (coefficient_values[q] * cell_matrix_C[i][k] * v_k) *
1217 * cell_matrix_C[j][
l] * v_l * fe_values_dgrt.JxW(q);
1224 * Next, we calculate the right hand side, @f$\int_{K} f q \mathrm{
d}x@f$:
1228 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1229 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1231 * cell_rhs(i) += (fe_values[pressure_interior].value(i, q) *
1232 * right_hand_side_values[q] * fe_values.JxW(q));
1237 * The last step is to distribute components of the local
1238 *
matrix into the system
matrix and transfer components of
1239 * the cell right hand side into the system right hand side:
1242 * cell->get_dof_indices(local_dof_indices);
1243 * constraints.distribute_local_to_global(
1244 * local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1253 * <a name=
"WGDarcyEquationdimsolve"></a>
1254 * <h4>WGDarcyEquation<dim>::solve</h4>
1258 * This step is rather trivial and the same as in many previous
1259 * tutorial programs:
1262 *
template <
int dim>
1263 *
void WGDarcyEquation<dim>::solve()
1265 *
SolverControl solver_control(1000, 1
e-8 * system_rhs.l2_norm());
1268 * constraints.distribute(solution);
1275 * <a name=
"WGDarcyEquationdimcompute_postprocessed_velocity"></a>
1276 * <h4>WGDarcyEquation<dim>::compute_postprocessed_velocity</h4>
1280 * In
this function, compute the velocity field from the pressure
1281 * solution previously computed. The
1282 * velocity is defined as @f$\mathbf{u}_h = \mathbf{Q}_h \left(
1283 * -\mathbf{K}\nabla_{
w,
d}p_h \right)@f$, which requires us to compute
1284 * many of the same terms as in the assembly of the system
matrix.
1285 * There are also the matrices @f$E^K,D^K@f$ we need to
assemble (see
1286 * the introduction) but they really just follow the same kind of
1291 * Computing the same matrices here as we have already done in the
1292 * `assemble_system()`
function is of course wasteful in terms of
1293 * CPU time. Likewise, we
copy some of the code from there to
this
1294 *
function, and
this is also generally a poor idea.
A better
1295 * implementation might provide
for a
function that encapsulates
1296 *
this duplicated code. One could also think of
using the classic
1297 * trade-off between computing efficiency and memory efficiency to
1298 * only compute the @f$C^K@f$ matrices once per cell during the
1299 * assembly, storing them somewhere on the side, and re-
using them
1300 * here. (This is what @ref step_51
"step-51" does,
for example, where the
1301 * `assemble_system()`
function takes an argument that determines
1302 * whether the local matrices are recomputed, and a similar approach
1303 * -- maybe with storing local matrices elsewhere -- could be
1304 * adapted
for the current program.)
1307 *
template <
int dim>
1308 *
void WGDarcyEquation<dim>::compute_postprocessed_velocity()
1310 * darcy_velocity.reinit(dof_handler_dgrt.n_dofs());
1312 *
const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1313 *
const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1316 * quadrature_formula,
1321 * face_quadrature_formula,
1327 * quadrature_formula,
1333 * face_quadrature_formula,
1339 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1340 *
const unsigned int dofs_per_cell_dgrt = fe_dgrt.dofs_per_cell;
1342 *
const unsigned int n_q_points = fe_values.get_quadrature().size();
1343 *
const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1345 *
const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
1348 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1349 * std::vector<types::global_dof_index> local_dof_indices_dgrt(
1350 * dofs_per_cell_dgrt);
1361 *
const Coefficient<dim> coefficient;
1362 * std::vector<Tensor<2, dim>> coefficient_values(n_q_points_dgrt);
1370 * In the introduction, we explained how to calculate the numerical velocity
1371 * on the cell. We need the pressure solution values on each cell,
1372 * coefficients of the Gram
matrix and coefficients of the @f$L_2@f$ projection.
1373 * We have already calculated the global solution, so we will
extract the
1374 * cell solution from the global solution. The coefficients of the Gram
1375 *
matrix have been calculated when we assembled the system
matrix for the
1376 * pressures. We will
do the same way here. For the coefficients of the
1377 * projection, we
do matrix multiplication, i.e., the inverse of the Gram
1378 *
matrix times the
matrix with @f$(\mathbf{K} \mathbf{
w}, \mathbf{
w})@f$ as
1379 * components. Then, we multiply all these coefficients and
call them beta.
1380 * The numerical velocity is the product of beta and the basis
functions of
1381 * the Raviart-Thomas space.
1386 * endc = dof_handler.end(), cell_dgrt = dof_handler_dgrt.
begin_active();
1387 *
for (; cell != endc; ++cell, ++cell_dgrt)
1389 * fe_values.reinit(cell);
1390 * fe_values_dgrt.reinit(cell_dgrt);
1392 * coefficient.value_list(fe_values_dgrt.get_quadrature_points(),
1393 * coefficient_values);
1397 * The component of this <code>cell_matrix_E</code> is the integral of
1398 * @f$(\mathbf{K} \mathbf{
w}, \mathbf{
w})@f$. <code>cell_matrix_M</code> is
1402 * cell_matrix_M = 0;
1403 * cell_matrix_E = 0;
1404 * for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1405 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1407 *
const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
1408 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1411 * fe_values_dgrt[velocities].value(k, q);
1413 * cell_matrix_E(i, k) +=
1414 * (coefficient_values[q] * v_i * v_k * fe_values_dgrt.JxW(q));
1416 * cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
1422 * To compute the
matrix @f$D@f$ mentioned in the introduction, we
1423 * then need to evaluate @f$D=M^{-1}
E@f$ as explained in the
1427 * cell_matrix_M.gauss_jordan();
1428 * cell_matrix_M.mmult(cell_matrix_D, cell_matrix_E);
1432 * Then we also need, again, to compute the
matrix @f$C@f$ that is
1433 * used to evaluate the weak discrete gradient. This is the
1434 * exact same code as used in the assembly of the system
1438 * cell_matrix_G = 0;
1439 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1440 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1442 *
const double div_v_i =
1443 * fe_values_dgrt[velocities].divergence(i, q);
1444 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1446 *
const double phi_j_interior =
1447 * fe_values[pressure_interior].value(j, q);
1449 * cell_matrix_G(i, j) -=
1450 * (div_v_i * phi_j_interior * fe_values.JxW(q));
1454 *
for (
const auto &face : cell->face_iterators())
1456 * fe_face_values.reinit(cell, face);
1457 * fe_face_values_dgrt.reinit(cell_dgrt, face);
1459 *
for (
unsigned int q = 0; q < n_face_q_points; ++q)
1463 *
for (
unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1466 * fe_face_values_dgrt[velocities].value(i, q);
1467 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1469 *
const double phi_j_face =
1470 * fe_face_values[pressure_face].value(j, q);
1472 * cell_matrix_G(i, j) +=
1473 * ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
1478 * cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
1482 * Finally, we need to
extract the pressure unknowns that
1483 * correspond to the current cell:
1486 * cell->get_dof_values(solution, cell_solution);
1490 * We are now in a position to compute the local velocity
1491 * unknowns (with respect to the Raviart-Thomas space we are
1492 * projecting the term @f$-\mathbf K \nabla_{
w,
d} p_h@f$ into):
1495 * cell_velocity = 0;
1496 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1497 *
for (
unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
1498 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1499 * cell_velocity(k) +=
1500 * -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
1504 * We compute Darcy velocity.
1505 * This is same as cell_velocity but used to graph Darcy velocity.
1508 * cell_dgrt->get_dof_indices(local_dof_indices_dgrt);
1509 *
for (
unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1510 *
for (
unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
1511 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1512 * darcy_velocity(local_dof_indices_dgrt[k]) +=
1513 * -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
1522 * <a name=
"WGDarcyEquationdimcompute_pressure_error"></a>
1523 * <h4>WGDarcyEquation<dim>::compute_pressure_error</h4>
1527 * This part is to calculate the @f$L_2@f$ error of the pressure. We
1528 * define a vector that holds the
norm of the error on each cell.
1530 * error in the @f$L_2@f$
norm on each cell. However, we really only
1531 * care about the error in the interior component of the solution
1532 * vector (we can't even evaluate the interface pressures at the
1533 * quadrature points because these are all located in the interior
1534 * of cells) and consequently have to use a weight function that
1535 * ensures that the interface component of the solution variable is
1537 * arguments indicate which component we want to select (component
1538 *
zero, i.
e., the interior pressures) and how many components there
1539 * are in total (two).
1542 * template <
int dim>
1543 *
void WGDarcyEquation<dim>::compute_pressure_error()
1549 * ExactPressure<dim>(),
1550 * difference_per_cell,
1553 * &select_interior_pressure);
1555 *
const double L2_error = difference_per_cell.l2_norm();
1556 * std::cout <<
"L2_error_pressure " << L2_error << std::endl;
1564 * <a name=
"WGDarcyEquationdimcompute_velocity_error"></a>
1565 * <h4>WGDarcyEquation<dim>::compute_velocity_error</h4>
1569 * In
this function, we evaluate @f$L_2@f$ errors
for the velocity on
1570 * each cell, and @f$L_2@f$ errors
for the flux on faces. The
function
1571 * relies on the `compute_postprocessed_velocity()`
function having
1572 * previous computed, which computes the velocity field based on the
1573 * pressure solution that has previously been computed.
1577 * We are going to evaluate velocities on each cell and calculate
1578 * the difference between numerical and exact velocities.
1581 *
template <
int dim>
1582 *
void WGDarcyEquation<dim>::compute_velocity_errors()
1584 *
const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1585 *
const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1588 * quadrature_formula,
1594 * face_quadrature_formula,
1600 *
const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1601 *
const unsigned int n_face_q_points_dgrt =
1602 * fe_face_values_dgrt.get_quadrature().size();
1604 * std::vector<Tensor<1, dim>> velocity_values(n_q_points_dgrt);
1605 * std::vector<Tensor<1, dim>> velocity_face_values(n_face_q_points_dgrt);
1609 *
const ExactVelocity<dim> exact_velocity;
1611 *
double L2_err_velocity_cell_sqr_global = 0;
1612 *
double L2_err_flux_sqr = 0;
1616 * Having previously computed the postprocessed velocity, we here
1617 * only have to
extract the corresponding values on each cell and
1618 * face and compare it to the exact values.
1621 *
for (
const auto &cell_dgrt : dof_handler_dgrt.active_cell_iterators())
1623 * fe_values_dgrt.reinit(cell_dgrt);
1627 * First compute the @f$L_2@f$ error between the postprocessed velocity
1628 * field and the exact
one:
1631 * fe_values_dgrt[velocities].get_function_values(darcy_velocity,
1633 *
double L2_err_velocity_cell_sqr_local = 0;
1634 *
for (
unsigned int q = 0; q < n_q_points_dgrt; ++q)
1638 * exact_velocity.value(fe_values_dgrt.quadrature_point(q));
1640 * L2_err_velocity_cell_sqr_local +=
1641 * ((velocity - true_velocity) * (velocity - true_velocity) *
1642 * fe_values_dgrt.JxW(q));
1644 * L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local;
1648 * For reconstructing the flux we need the size of cells and
1649 * faces. Since fluxes are calculated on faces, we have the
1650 *
loop over all four faces of each cell. To calculate the
1651 * face velocity, we
extract values at the quadrature points from the
1652 * `darcy_velocity` which we have computed previously. Then, we
1653 * calculate the squared velocity error in normal direction. Finally, we
1654 * calculate the @f$L_2@f$ flux error on the cell by appropriately scaling
1655 * with face and cell areas and add it to the global error.
1658 *
const double cell_area = cell_dgrt->measure();
1659 *
for (
const auto &face_dgrt : cell_dgrt->face_iterators())
1661 *
const double face_length = face_dgrt->measure();
1662 * fe_face_values_dgrt.reinit(cell_dgrt, face_dgrt);
1663 * fe_face_values_dgrt[velocities].get_function_values(
1664 * darcy_velocity, velocity_face_values);
1666 *
double L2_err_flux_face_sqr_local = 0;
1667 *
for (
unsigned int q = 0; q < n_face_q_points_dgrt; ++q)
1671 * exact_velocity.value(fe_face_values_dgrt.quadrature_point(q));
1674 * fe_face_values_dgrt.normal_vector(q);
1676 * L2_err_flux_face_sqr_local +=
1677 * ((velocity * normal - true_velocity * normal) *
1678 * (velocity * normal - true_velocity * normal) *
1679 * fe_face_values_dgrt.JxW(q));
1681 *
const double err_flux_each_face =
1682 * L2_err_flux_face_sqr_local / face_length * cell_area;
1683 * L2_err_flux_sqr += err_flux_each_face;
1689 * After adding up errors over all cells and faces, we take the
1690 * square root and get the @f$L_2@f$ errors of velocity and
1691 * flux. These we output to screen.
1694 *
const double L2_err_velocity_cell =
1695 *
std::sqrt(L2_err_velocity_cell_sqr_global);
1696 *
const double L2_err_flux_face =
std::sqrt(L2_err_flux_sqr);
1698 * std::cout <<
"L2_error_vel: " << L2_err_velocity_cell << std::endl
1699 * <<
"L2_error_flux: " << L2_err_flux_face << std::endl;
1706 * <a name=
"WGDarcyEquationoutput_results"></a>
1707 * <h4>WGDarcyEquation::output_results</h4>
1711 * We have two sets of results to output: the interior solution and
1712 * the skeleton solution. We use <code>
DataOut</code> to visualize
1713 * interior results. The graphical output
for the skeleton results
1718 * In both of the output files, both the interior and the face
1719 * variables are stored. For the
interface output, the output file
1720 * simply contains the interpolation of the interior pressures onto
1721 * the faces, but because it is undefined which of the two interior
1722 * pressure variables you get from the two adjacent cells, it is
1723 * best to ignore the interior pressure in the
interface output
1724 * file. Conversely,
for the cell interior output file, it is of
1725 * course impossible to show any
interface pressures @f$p^\partial@f$,
1726 * because these are only available on interfaces and not cell
1727 * interiors. Consequently, you will see them shown as an
invalid
1728 *
value (such as an infinity).
1732 * For the cell interior output, we also want to output the velocity
1733 * variables. This is a bit tricky since it lives on the same mesh
1734 * but uses a different
DoFHandler object (the pressure variables live
1735 * on the `dof_handler`
object, the Darcy velocity on the `dof_handler_dgrt`
1736 *
object). Fortunately, there are variations of the
1738 *
DoFHandler a vector corresponds to, and consequently we can visualize
1739 * the data from both
DoFHandler objects within the same file.
1742 * template <
int dim>
1743 *
void WGDarcyEquation<dim>::output_results() const
1750 * First attach the pressure solution to the
DataOut object:
1753 *
const std::vector<std::string> solution_names = {
"interior_pressure",
1754 *
"interface_pressure"};
1759 * Then
do the same with the Darcy velocity field, and
continue
1760 * with writing everything out into a file.
1763 *
const std::vector<std::string> velocity_names(dim,
"velocity");
1764 *
const std::vector<
1766 * velocity_component_interpretation(
1771 * velocity_component_interpretation);
1774 * std::ofstream output(
"solution_interior.vtu");
1780 * data_out_faces.attach_dof_handler(dof_handler);
1781 * data_out_faces.add_data_vector(solution,
"Pressure_Face");
1782 * data_out_faces.build_patches(fe.degree);
1783 * std::ofstream face_output(
"solution_interface.vtu");
1784 * data_out_faces.write_vtu(face_output);
1792 * <a name=
"WGDarcyEquationrun"></a>
1797 * This is the
final function of the main
class. It calls the other
functions
1801 *
template <
int dim>
1804 * std::cout <<
"Solving problem in " << dim <<
" space dimensions."
1808 * assemble_system();
1810 * compute_postprocessed_velocity();
1811 * compute_pressure_error();
1812 * compute_velocity_errors();
1822 * <a name=
"Thecodemaincodefunction"></a>
1823 * <h3>The <code>main</code>
function</h3>
1827 * This is the main
function. We can change the dimension here to
run in 3
d.
1834 * Step61::WGDarcyEquation<2> wg_darcy(0);
1837 *
catch (std::exception &exc)
1839 * std::cerr << std::endl
1841 * <<
"----------------------------------------------------"
1843 * std::cerr <<
"Exception on processing: " << std::endl
1844 * << exc.what() << std::endl
1845 * <<
"Aborting!" << std::endl
1846 * <<
"----------------------------------------------------"
1852 * std::cerr << std::endl
1854 * <<
"----------------------------------------------------"
1856 * std::cerr <<
"Unknown exception!" << std::endl
1857 * <<
"Aborting!" << std::endl
1858 * <<
"----------------------------------------------------"
1866 <a name=
"Results"></a><h1>Results</h1>
1869 We
run the program with a right hand side that will produce the
1870 solution @f$p =
\sin(\pi x)
\sin(\pi y)@f$ and with homogeneous Dirichlet
1871 boundary conditions in the domain @f$\Omega = (0,1)^2@f$. In addition, we
1872 choose the coefficient
matrix in the differential
operator
1873 @f$\mathbf{K}@f$ as the
identity matrix. We test
this setup
using
1874 @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ and
1875 @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ element combinations, which
one can
1876 select by
using the appropriate constructor argument
for the
1877 `WGDarcyEquation`
object in `main()`. We will then visualize pressure
1878 values in interiors of cells and on faces. We want to see that the
1879 pressure maximum is around 1 and the minimum is around 0. With mesh
1880 refinement, the convergence rates of pressure, velocity and flux
1881 should then be around 1
for @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ , 2
for
1882 @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$, and 3
for
1883 @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$.
1886 <a name=
"TestresultsoniWGQsub0subQsub0subRTsub0subi"></a><h3>Test results on <i>WG(Q<sub>0</sub>,Q<sub>0</sub>;RT<sub>[0]</sub>)</i></h3>
1889 The following figures show interior pressures and face pressures
using the
1890 @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ element. The mesh is refined 2 times (top)
1891 and 4 times (bottom), respectively. (This number can be adjusted in the
1892 `make_grid()`
function.) When the mesh is coarse,
one can see
1893 the face pressures @f$p^\partial@f$ neatly between the values of the interior
1894 pressures @f$p^\circ@f$ on the two adjacent cells.
1896 <table align=
"center">
1898 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_2d_2.png" alt=
""></td>
1899 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_3d_2.png" alt=
""></td>
1902 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_2d_4.png" alt=
""></td>
1903 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg000_3d_4.png" alt=
""></td>
1907 From the figures, we can see that with the mesh refinement, the maximum and
1908 minimum pressure values are approaching the values we expect.
1909 Since the mesh is a rectangular mesh and
numbers of cells in each direction is even, we
1910 have
symmetric solutions. From the 3
d figures on the right,
1911 we can see that on @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, the pressure is a constant
1912 in the interior of the cell, as expected.
1914 <a name=
"Convergencetableforik0i"></a><h4>Convergence table for <i>k=0</i></h4>
1917 We
run the code with differently refined meshes (chosen in the `make_grid()`
function)
1918 and get the following convergence rates of pressure,
1919 velocity, and flux (as defined in the introduction).
1921 <table align=
"center" class=
"doxtable">
1923 <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>
1926 <td> 2 </td><td> 1.587e-01 </td><td> 5.113e-01 </td><td> 7.062e-01 </td>
1929 <td> 3 </td><td> 8.000e-02 </td><td> 2.529e-01 </td><td> 3.554e-01 </td>
1932 <td> 4 </td><td> 4.006e-02 </td><td> 1.260e-01 </td><td> 1.780e-01 </td>
1935 <td> 5 </td><td> 2.004e-02 </td><td> 6.297e-02 </td><td> 8.902e-02 </td>
1938 <th>Conv.rate </th><th> 1.00 </th><th> 1.00 </th><th> 1.00 </th>
1942 We can see that the convergence rates of @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ are around 1.
1943 This, of course, matches our theoretical expectations.
1946 <a name=
"TestresultsoniWGQsub1subQsub1subRTsub1subi"></a><h3>Test results on <i>WG(Q<sub>1</sub>,Q<sub>1</sub>;RT<sub>[1]</sub>)</i></h3>
1949 We can repeat the experiment from above
using the next higher polynomial
1951 The following figures are interior pressures and face pressures implemented
using
1952 @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$. The mesh is refined 4 times. Compared to the
1953 previous figures
using
1954 @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, on each cell, the solution is no longer constant
1955 on each cell, as we now use bilinear polynomials to
do the approximation.
1956 Consequently, there are 4 pressure values in
one interior, 2 pressure values on
1959 <table align=
"center">
1961 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg111_2d_4.png" alt=
""></td>
1962 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg111_3d_4.png" alt=
""></td>
1966 Compared to the corresponding image
for the @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$
1967 combination, the solution is now substantially more accurate and, in
1968 particular so close to being continuous at the interfaces that we can
1969 no longer distinguish the interface pressures @f$p^\partial@f$ from the
1970 interior pressures @f$p^\circ@f$ on the adjacent cells.
1972 <a name=
"Convergencetableforik1i"></a><h4>Convergence table for <i>k=1</i></h4>
1975 The following are the convergence rates of pressure, velocity, and flux
1976 we obtain from
using the @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ element combination:
1978 <table align=
"center" class=
"doxtable">
1980 <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>
1983 <td> 2 </td><td> 1.613e-02 </td><td> 5.093e-02 </td><td> 7.167e-02 </td>
1986 <td> 3 </td><td> 4.056e-03 </td><td> 1.276e-02 </td><td> 1.802e-02 </td>
1989 <td> 4 </td><td> 1.015e-03 </td><td> 3.191e-03 </td><td> 4.512e-03 </td>
1992 <td> 5 </td><td> 2.540e-04 </td><td> 7.979e-04 </td><td> 1.128e-03 </td>
1995 <th>Conv.rate </th><th> 2.00 </th><th> 2.00 </th><th> 2.00 </th>
1999 The convergence rates of @f$WG(Q_1,Q_1;RT_{[1]})@f$ are around 2, as expected.
2003 <a name=
"TestresultsoniWGQsub2subQsub2subRTsub2subi"></a><h3>Test results on <i>WG(Q<sub>2</sub>,Q<sub>2</sub>;RT<sub>[2]</sub>)</i></h3>
2006 Let us go
one polynomial degree higher.
2007 The following are interior pressures and face pressures implemented
using
2008 @f$WG(Q_2,Q_2;RT_{[2]})@f$, with mesh size @f$h = 1/32@f$ (i.e., 5 global mesh
2009 refinement steps). In the program, we use
2010 `data_out_face.build_patches(fe.degree)` when generating graphical output
2012 we divide each 2
d cell interior into 4 subcells in order to provide a better
2013 visualization of the quadratic polynomials.
2014 <table align=
"center">
2016 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg222_2d_5.png" alt=
""></td>
2017 <td><img src=
"https://www.dealii.org/images/steps/developer/step-61.wg222_3d_5.png" alt=
""></td>
2022 <a name=
"Convergencetableforik2i"></a><h4>Convergence table for <i>k=2</i></h4>
2025 As before, we can generate convergence data
for the
2026 @f$L_2@f$ errors of pressure, velocity, and flux
2027 using the @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ combination:
2029 <table align=
"center" class=
"doxtable">
2031 <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>
2034 <td> 2 </td><td> 1.072e-03 </td><td> 3.375e-03 </td><td> 4.762e-03 </td>
2037 <td> 3 </td><td> 1.347e-04 </td><td> 4.233e-04 </td><td> 5.982e-04 </td>
2040 <td> 4 </td><td> 1.685e-05 </td><td> 5.295e-05 </td><td> 7.487e-05 </td>
2043 <td> 5 </td><td> 2.107e-06 </td><td> 6.620e-06 </td><td> 9.362e-06 </td>
2046 <th>Conv.rate </th><th> 3.00 </th><th> 3.00 </th><th> 3.00 </th>
2050 Once more, the convergence rates of @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ is
2051 as expected, with values around 3.
2054 <a name=
"PlainProg"></a>
2055 <h1> The plain program</h1>
2056 @include
"step-61.cc"