From: Zhuoran Wang Date: Fri, 27 Apr 2018 16:31:49 +0000 (+0000) Subject: Add step-61. X-Git-Tag: v9.1.0-rc1~357^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6455%2Fhead;p=dealii.git Add step-61. This tutorial program adds a code that solves the Poisson equation using the weak Galerkin formulation. --- diff --git a/examples/step-61/CMakeLists.txt b/examples/step-61/CMakeLists.txt new file mode 100644 index 0000000000..83a33d091a --- /dev/null +++ b/examples/step-61/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-61 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-61") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +FIND_PACKAGE(deal.II 9.0.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-61/doc/builds-on b/examples/step-61/doc/builds-on new file mode 100644 index 0000000000..2bb1b6e1da --- /dev/null +++ b/examples/step-61/doc/builds-on @@ -0,0 +1 @@ +step-51 step-7 step-20 diff --git a/examples/step-61/doc/intro.dox b/examples/step-61/doc/intro.dox new file mode 100644 index 0000000000..f5ca096683 --- /dev/null +++ b/examples/step-61/doc/intro.dox @@ -0,0 +1,279 @@ +
+ + +This program was contributed by Zhuoran Wang. + + + +

Introduction

+ +This tutorial program presents an implementation of the "weak Galerkin" +finite element method for the Poisson equation. In some sense, the motivation for +considering this method starts from the same point as in step-51: We would like to +consider discontinuous shape functions, but then need to address the fact that +the resulting problem has a large number of degrees of freedom (because, for +example, each vertex carries as many degrees of freedom as there are adjacent cells). +We also have to address the fact that, in general, every degree of freedom +on one cell couples with all of the degrees of freedom on each of its face neighbor +cells. Consequently, the matrix one gets from the "traditional" discontinuous +Galerkin methods are both large and relatively dense. + +Both the hybridized discontinuous Galerkin method (HDG) in step-51 and the weak +Galerkin (WG) method in this tutorial address the issue of coupling by introducing +additional degrees of freedom whose shape functions only live on a face between +cells (i.e., on the "skeleton" of the mesh), and which therefore "insulate" the +degrees of freedom on the adjacent cells from each other: cell degrees of freedom +only couple with other cell degrees of freedom on the same cell, as well as face +degrees of freedom, but not with cell degrees of freedom on neighboring cells. +Consequently, the shape functions for these cell degrees of freedom are +discontinuous and only "live" on exactly one cell. + +For a given equation, say the second order Poisson equation, +the difference between the HDG and the WG method is how exactly one formulates +the problem that connects all of these different shape functions. The HDG does +things by reformulating second order problems in terms of a system of first +order equations and then conceptually considers the face degrees of freedom +to be "fluxes" of this first order system. In contrast, the WG method keeps things +in second order form and considers the face degrees of freedom as of the same +type as the primary solution variable, just restricted to the lower-dimensional +faces. For the purposes of the equation, one then needs to somehow "extend" +these shape functions into the interior of the cell when defining what it means +to apply a differential operator to them. Compared to the HDG, the method +has the advantage that it does not lead to a proliferation of unknowns due +to rewriting the equation as a first-order system, but it is also not quite +as easy to implement. However, as we will see in the following, this +additional effort is not prohibitive. + + +

Weak Galerkin finite element methods

+ +Weak Galerkin Finite Element Methods (WGFEMs) use discrete weak functions +to approximate scalar unknowns and discrete weak gradients to +approximate classical gradients. +It was introduced by Junping Wang and Xiu Ye +in the paper: A weak Galerkin finite element method for second order elliptic problems, +J. Comput. Appl. Math., 2013, 103-115. +Compared to the continuous Galerkin method, +the weak Galerkin method satisfies important physical properties, namely +local mass conservation and bulk normal flux continuity. +It results in a SPD linear system, and expected convergence rates can +be obtained with mesh refinement. + + +

WGFEM applied to the Poisson equation

+This program solves the Poisson equation +using the weak Galerkin finite element method: +@f{eqnarray*} + \nabla \cdot \left( -\mathbf{K} \nabla p \right) + = f, + \quad \mathbf{x} \in \Omega, \\ + p = p_D,\quad \mathbf{x} \in \Gamma^D, \\ + \mathbf{u} \cdot \mathbf{n} = u_N, + \quad \mathbf{x} \in \Gamma^N, +@f} +where $\Omega \subset \mathbb{R}^n (n=2,3)$ is a bounded domain. +In the context of the flow of a fluid through a porous medium, +$p$ is the pressure, $\mathbf{K}$ is a permeability tensor, +$ f $ is the source term, and +$ p_D, u_N $ represent Dirichlet and Neumann boundary conditions. +We can introduce a flux, $\mathbf{u} = -\mathbf{K} \nabla p$, that corresponds +to the Darcy velocity (in the way we did in step-20) and this variable will +be important in the considerations below. + +In this program, we will consider a test case where the exact pressure +is $p = \sin \left( \pi x)\sin(\pi y \right)$ on the unit square domain, +with homogenous Dirichelet boundary conditions and identity matrix $\mathbf{K}$. +Then we will calculate $L_2$ errors of pressure, velocity and flux. + + +

Weak Galerkin scheme

+ +Via integration by parts, the weak Galerkin scheme for the Poisson equation is +@f{equation*} +\mathcal{A}_h\left(p_h,q \right) = \mathcal{F} \left(q \right), +@f} +where +@f{equation*} +\mathcal{A}_h\left(p_h,q\right) + := \sum_{T \in \mathcal{T}_h} + \int_T \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q \mathrm{d}x, +@f} +and +@f{equation*} +\mathcal{F}\left(q\right) + := \sum_{T \in \mathcal{T}_h} \int_T f \, q^\circ \mathrm{d}x + - \sum_{\gamma \in \Gamma_h^N} \int_\gamma u_N q^\partial \mathrm{d}x, +@f} +$q^\circ$ is the shape function of the polynomial space in the interior, +$q^\partial$ is the shape function of the polynomial space on faces, $ \nabla_{w,d} $ means the discrete weak gradient used to approximate the classical gradient. +We use FE_DGQ as the interior polynomial space, +FE_FaceQ as the face polynomial space, and Raviart-Thomas elements for the velocity +$\mathbf{u} = -{\mathbf{K}} \nabla_{w,d} p$. + +

Assembling the linear system

+ +First, we solve for the pressure. +We collect two local spaces together in one FESystem, +the first component in this finite element system denotes +the space for interior pressure, and the second denotes +the space for face pressure. +For the interior component, we use the polynomial space FE_DGQ. +For the face component, we use FE_FaceQ. + +We use shape functions defined on spaces FE_DGQ and FE_FaceQ to +approximate pressures, i.e., $p_h = \sum a_i \phi_i,$ +where $\phi_i$ are shape functions of FESystem. +We construct the local system by using discrete weak gradients of +shape functions of FE_DGQ and FE_FaceQ. +The discrete weak gradients of shape functions $\nabla_{w,d} \phi$ are defined as +$\nabla_{w,d} \phi = \sum_{i=1}^m c_i \mathbf{w}_i,$ +where $\mathbf{w}_i$ is the basis function of $RT(k)$. + +Using integration by parts, we have a small linear system +on each element $T$, +@f{equation*} +\int_{T} \left(\nabla_{w,d} \phi \right) \cdot \mathbf{w} \mathrm{d}x= +\int_{T^\partial} \phi^{\partial} \left(\mathbf{w} \cdot \mathbf{n}\right) \mathrm{d}x- +\int_{T^\circ} \phi^{\circ} \left(\nabla \cdot \mathbf{w}\right) \mathrm{d}x, +\quad \forall \mathbf{w} \in RT_{[k]}(E), +@f} + +@f{equation*} +\sum_{i=1}^m c_i \int_T \mathbf{w}_i \cdot \mathbf{w}_j \mathrm{d}x = +\int_{T^{\partial}} \phi_i^{\partial} +\left(\mathbf{w}_j \cdot \mathbf{n} \right) \mathrm{d}x - +\int_{T^{\circ}} \phi_i^{\circ} \left (\nabla \cdot \mathbf{w}_j \right)\mathrm{d}x, +@f} +which can be simplified to be +@f{equation*} +\mathbf{C}_{E}\mathbf{M}_{E} = \mathbf{F}_{E}, +@f} +where $\mathbf{C}_E$ is the matrix with unknown coefficients $c$, +$\mathbf{M}_E$ is the Gram matrix +$\left[ \int_T \mathbf{w}_i \cdot \mathbf{w}_j \right] \mathrm{d}x$, +$\mathbf{F}_E$ is the matrix of right hand side, +$\mathbf{w}$ and $\phi_i^{\circ}$ are in FEValues, +$\phi_i^{\partial}$ is in FEFaceValues. +Then we solve for $\mathbf{C}_E = \mathbf{F}_E \mathbf{M}_E^{-1}$. +Now, discrete weak gradients of shape functions are written as +linear combinations of basis functions of the $RT$ space. +In our code, we name $\mathbf{C}_E$ as cell_matrix_C, +$\mathbf{M}_E$ as cell_matrix_rt, +$\mathbf{F}_E$ as cell_matrix_F. + +The components of the local cell matrices $\mathbf{A}$ are +@f{equation*} +\mathbf{A}_{ij} = +\int_{T} \mathbf{K} \nabla_{w,d} \phi_i \cdot \nabla_{w,d} \phi_j \mathrm{d}x. +@f} +From previous steps, we know $\nabla_{w,d} \phi_i = \sum_{k=1}^m c_{ik} \mathbf{w}_k,$ +and $\nabla_{w,d} \phi_j = \sum_{l=1}^m c_{jl} \mathbf{w}_l.$ +Then combining the coefficients we have calculated, components of $\mathbf{A}$ are calculated as +@f{equation*} +\int_T \sum_{k,l = 1}^{m}c_{ik} c_{jl} \left(\mathbf{K} \mathbf{w}_i \cdot \mathbf{w}_j\right) \mathrm{d}x += \sum_{k,l = 1}^{m}c_{ik} c_{jl} \int_{T} \mathbf{K} \mathbf{w}_i \cdot \mathbf{w}_j \mathrm{d}x. +@f} + +Next, we use ConstraintMatrix::distribute_local_to_global to +distribute contributions from local matrices $\mathbf{A}$ to the system matrix. + +In the scheme +$\mathcal{A}_h\left(p_h,q \right) = \mathcal{F} \left( q \right),$ +we have system matrix and system right hand side, +we can solve for the coefficients of the system matrix. +The solution vector of the scheme represents the pressure values in interiors and on faces. + +

Post-processing and $L_2$-errors

+ +After we have calculated the numerical pressure $p$, +we use discrete weak gradients of $p$ to calculate the velocity on each element. + +On each element the gradient of the numerical pressure $\nabla p$ can be +approximated by discrete weak gradients $ \nabla_{w,d}\phi_i$, so +@f{equation*} +\nabla_{w,d} p_h = \sum_{i} a_i \nabla_{w,d}\phi_i. +@f} + +The numerical velocity $ \mathbf{u}_h = -\mathbf{K} \nabla_{w,d}p_h$ can be written as +@f{equation*} +\mathbf{u}_h = -\mathbf{K} \nabla_{w,d} p = +-\sum_{i} \sum_{j} a_ic_{ij}\mathbf{K}\mathbf{w}_j, +@f} +where $c_{ij}$ is the coefficient of Gram matrix, +$\mathbf{w}_j$ is the basis function of the $RT$ space. +$\mathbf{K} \mathbf{w}_j$ may not be in the $RT$ space. +So we need $L_2$-projection to project it back to the $RT$ space. + +We define the projection as +$ \mathbf{Q}_h \left( \mathbf{K}\mathbf{w}_j \right) = +\sum_{k} d_{jk}\mathbf{w}_k$. +For any $j$, +$\left( \mathbf{Q}_h \left( \mathbf{Kw}_j \right),\mathbf{w}_k \right)_E = +\left( \mathbf{Kw}_j,\mathbf{w}_k \right)_E.$ +So the numerical velocity becomes +@f{equation*} +\mathbf{u}_h = \mathbf{Q}_h \left( -\mathbf{K}\nabla_{w,d}p_h \right) = +-\sum_{i=0}^{4} \sum_{j=1}^{4}a_ib_{ij}\mathbf{Q}_h \left( \mathbf{K}\mathbf{w}_j \right), +@f} +and we have the following system to solve for the coefficients $d_{jk}$, +@f{equation*} + \left[ + \begin{matrix} + \left(\mathbf{w}_i,\mathbf{w}_j \right) + \end{matrix} + \right] + \left[ + \begin{matrix} + d_{jk} + \end{matrix} + \right] + = + \left[ + \begin{matrix} + \left( \mathbf{Kw}_j,\mathbf{w}_k \right) + \end{matrix} + \right]. +@f} +$ + \left[ + \begin{matrix} + d_{jk} + \end{matrix} + \right] +$ +is named cell_matrix_D, +$ +\left[ + \begin{matrix} + \left( \mathbf{Kw}_j,\mathbf{w}_k \right) + \end{matrix} + \right] +$ +is named cell_matrix_E. + +Then the elementwise velocity is +@f{equation*} +\mathbf{u}_h = -\sum_{i} \sum_{j}a_ic_{ij}\sum_{k}d_{jk}\mathbf{w}_k = +\sum_{k}- \left(\sum_{j} \sum_{i} a_ic_{ij}d_{jk} \right)\mathbf{w}_k, +@f} +where $-\sum_{j} \sum_{i} a_ic_{ij}d_{jk}$ is named +beta in the code. + +We calculate the $L_2$-errors of pressure, velocity and flux +by the following formulas, +@f{eqnarray*} +\|p-p_h^\circ\|^2 + = \sum_{T \in \mathcal{T}_h} \|p-p_h^\circ\|_{L^2(E)}^2, \\ + \|\mathbf{u}-\mathbf{u}_h\|^2 + = \sum_{T \in \mathcal{T}_h} \|\mathbf{u}-\mathbf{u}_h\|_{L^2(E)^2}^2,\\ +\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|^2 + = \sum_{T \in \mathcal{T}_h} \sum_{\gamma \subset T^\partial} + \frac{|T|}{|\gamma|} \|\mathbf{u} \cdot \mathbf{n} - \mathbf{u}_h \cdot \mathbf{n}\|_{L^2(\gamma)}^2, +@f} +where $| T |$ is the area of the element, +$\gamma$ are faces of the element, +$\mathbf{n}$ are unit normal vectors of each face. + +We will extract interior pressure solutions of each cell +from the global solution and calculate the $L_2$ error +by using function VectorTools::integrate_difference. diff --git a/examples/step-61/doc/kind b/examples/step-61/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-61/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-61/doc/results.dox b/examples/step-61/doc/results.dox new file mode 100644 index 0000000000..b71c57abb1 --- /dev/null +++ b/examples/step-61/doc/results.dox @@ -0,0 +1,88 @@ +

Results

+ +We run the test example $p = \sin(\pi x) \sin(\pi y)$ with homogenous Dirichelet boundary conditions in the domain $\Omega = (0,1)^2$. And $\mathbf{K}$ is the identity matrix. We test it on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$, $\mbox{WG}(Q_1,Q_1;RT_{[1]})$ and $\mbox{WG}(Q_2,Q_2;RT_{[2]})$. We will visualize pressure values in interiors and on faces. We want to see the pressure maximum is around 1 and the minimum is around 0. With the mesh refinement, the convergence rates of pressure, velocity and flux should be around 1 on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$ , 2 on $\mbox{WG}(Q_1,Q_1;RT_{[1]})$, and 3 on $\mbox{WG}(Q_2,Q_2;RT_{[2]})$. + +

Test results on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$

+The following figures are interior pressures and face pressures implemented on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$. The mesh is refined 2 times and 4 times separately. + + + + + + + + + + +
+ +From the figures, we can see that with the mesh refinement, the maximum and minimum are approaching to what we expect. +Since the mesh is a rectangular mesh and numbers of refinement are even, we have symmetric solutions. From the 3d figures, we can see that on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$, pressure is a constant in the interior of the cell. + +

Convergence table

+ +We run the code with finer meshes and get the following convergence rates of pressure, +velocity and flux. + +@code +number of refinements $\|p-p_h^\circ\|$ $\|\mathbf{u}-\mathbf{u}_h\|$ $\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|$ + 2 1.587e-01 5.113e-01 7.062e-01 + 3 8.000e-02 2.529e-01 3.554e-01 + 4 4.006e-02 1.260e-01 1.780e-01 + 5 2.004e-02 6.297e-02 8.902e-02 +Conv.rate 1.00 1.00 1.00 +@endcode +We can see that the convergence rates of $\mbox{WG}(Q_0,Q_0;RT_{[0]})$ are around 1. + + +

Test results on $\mbox{WG}(Q_1,Q_1;RT_{[1]})$

+ +The following figures are interior pressures and face pressures implemented on $\mbox{WG}(Q_1,Q_1;RT_{[1]})$. The mesh is refined 4 times. Compared to the previous figures on +$\mbox{WG}(Q_0,Q_0;RT_{[0]})$, on each cell, the result is not a constant. Because we use higher order polynomials to do approximation. So there are 4 pressure values in one interior, 2 pressure values on each face. We use data_out_face.build_patches (fe.degree) +to divide each cell interior into 4 subcells. + + + + + + +
+ +

Convergence table

+ +These are the convergence rates of pressure, velocity and flux on $\mbox{WG}(Q_1,Q_1;RT_{[1]})$ + +@code +number of refinements $\|p-p_h^\circ\|$ $\|\mathbf{u}-\mathbf{u}_h\|$ $\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|$ + 2 1.613e-02 5.093e-02 7.167e-02 + 3 4.056e-03 1.276e-02 1.802e-02 + 4 1.015e-03 3.191e-03 4.512e-03 + 5 2.540e-04 7.979e-04 1.128e-03 +Conv.rate 2.00 2.00 2.00 +@endcode +The convergence rates of $WG(Q_1,Q_1;RT_{[1]})$ are around 2. + +

Test results on $WG(Q_2,Q_2;RT_{[2]})$

+ +These are interior pressures and face pressures implemented on $WG(Q_2,Q_2;RT_{[2]})$, with mesh size $h = 1/32$. + + + + + + +
+ +

Convergence table

+ +This is the convergence table of $L_2$ errors of pressure, velocity and flux on $\mbox{WG}(Q_2,Q_2;RT_{[2]})$ + +@code +number of refinements $\|p-p_h^\circ\|$ $\|\mathbf{u}-\mathbf{u}_h\|$ $\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|$ + 2 1.072e-03 3.375e-03 4.762e-03 + 3 1.347e-04 4.233e-04 5.982e-04 + 4 1.685e-05 5.295e-05 7.487e-05 + 5 2.107e-06 6.620e-06 9.362e-06 +Conv.rate 3.00 3.00 3.00 +@endcode +The convergence rates of $\mbox{WG}(Q_2,Q_2;RT_{[2]})$ are around 3. diff --git a/examples/step-61/doc/step-61.wg000_2d_2.png b/examples/step-61/doc/step-61.wg000_2d_2.png new file mode 100644 index 0000000000..bf3be238e6 Binary files /dev/null and b/examples/step-61/doc/step-61.wg000_2d_2.png differ diff --git a/examples/step-61/doc/step-61.wg000_2d_4.png b/examples/step-61/doc/step-61.wg000_2d_4.png new file mode 100644 index 0000000000..fe5e5bc9ff Binary files /dev/null and b/examples/step-61/doc/step-61.wg000_2d_4.png differ diff --git a/examples/step-61/doc/step-61.wg000_3d_2.png b/examples/step-61/doc/step-61.wg000_3d_2.png new file mode 100644 index 0000000000..73b6c7b67e Binary files /dev/null and b/examples/step-61/doc/step-61.wg000_3d_2.png differ diff --git a/examples/step-61/doc/step-61.wg000_3d_4.png b/examples/step-61/doc/step-61.wg000_3d_4.png new file mode 100644 index 0000000000..a76e4801d2 Binary files /dev/null and b/examples/step-61/doc/step-61.wg000_3d_4.png differ diff --git a/examples/step-61/doc/step-61.wg111_2d_4.png b/examples/step-61/doc/step-61.wg111_2d_4.png new file mode 100644 index 0000000000..23258e0b02 Binary files /dev/null and b/examples/step-61/doc/step-61.wg111_2d_4.png differ diff --git a/examples/step-61/doc/step-61.wg111_3d_4.png b/examples/step-61/doc/step-61.wg111_3d_4.png new file mode 100644 index 0000000000..c3f0f2f59c Binary files /dev/null and b/examples/step-61/doc/step-61.wg111_3d_4.png differ diff --git a/examples/step-61/doc/step-61.wg222_2d_5.png b/examples/step-61/doc/step-61.wg222_2d_5.png new file mode 100644 index 0000000000..1be5805591 Binary files /dev/null and b/examples/step-61/doc/step-61.wg222_2d_5.png differ diff --git a/examples/step-61/doc/step-61.wg222_3d_5.png b/examples/step-61/doc/step-61.wg222_3d_5.png new file mode 100644 index 0000000000..c79c7fbffd Binary files /dev/null and b/examples/step-61/doc/step-61.wg222_3d_5.png differ diff --git a/examples/step-61/doc/tooltip b/examples/step-61/doc/tooltip new file mode 100644 index 0000000000..1f69a06061 --- /dev/null +++ b/examples/step-61/doc/tooltip @@ -0,0 +1 @@ +Weak Galerkin method applied to the Poisson equation diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc new file mode 100644 index 0000000000..90ed9c019b --- /dev/null +++ b/examples/step-61/step-61.cc @@ -0,0 +1,915 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2018 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * Author: Zhuoran Wang + */ + +// @sect3{Include files} +// This program is based on step-7, step-20 and step-51, +// we add these include files. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + +// @sect3{The WGDarcyEquation class template} + +// We will solve for the numerical pressure in the interior and on faces and +// calculate its $L_2$ error of pressure. In the post-processing step, we will +// calculate $L_2$-errors of velocity and flux. +template +class WGDarcyEquation +{ +public: + WGDarcyEquation(); + void run(); + +private: + void make_grid(); + void setup_system(); + void assemble_system(); + void solve(); + void postprocess(); + void process_solution(); + void output_results() const; + + Triangulation triangulation; + + AffineConstraints constraints; + + FE_RaviartThomas fe_rt; + DoFHandler dof_handler_rt; + + // The finite element system is used for interior and face solutions. + FESystem fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + +// @sect3{Right hand side, boundary values, and exact solution} + +// Next, we define the coefficient matrix $\mathbf{K}$, +// Dirichlet boundary conditions, the right-hand side $f = 2\pi^2 \sin(\pi x) +// \sin(\pi y)$, and the reference solutions $p = \sin(\pi x) \sin(\pi y) $. +// +// The coefficient matrix $\mathbf{K}$ is the identity matrix as a test example. +template +class Coefficient : public TensorFunction<2, dim> +{ +public: + Coefficient() + : TensorFunction<2, dim>() + {} + virtual void value_list(const std::vector> &points, + std::vector> & values) const override; +}; + +template +void Coefficient::value_list(const std::vector> &points, + std::vector> & values) const +{ + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); + for (unsigned int p = 0; p < points.size(); ++p) + { + values[p].clear(); + for (unsigned int d = 0; d < dim; ++d) + values[p][d][d] = 1; + } +} + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues() + : Function(2) + {} + virtual double value(const Point & p, + const unsigned int component = 0) const override; +}; + +template +double BoundaryValues::value(const Point & /*p*/, + const unsigned int /*component*/) const +{ + return 0; +} + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + virtual double value(const Point & p, + const unsigned int component = 0) const; +}; + +template +double RightHandSide::value(const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0.0; + return_value = 2 * M_PI * M_PI * sin(M_PI * p[0]) * sin(M_PI * p[1]); + return return_value; +} + +template +class Solution : public Function +{ +public: + Solution() + : Function(1) + {} + virtual double value(const Point &p, const unsigned int) const; +}; + +template +double Solution::value(const Point &p, const unsigned int) const +{ + double return_value = 0; + return_value = sin(M_PI * p[0]) * sin(M_PI * p[1]); + return return_value; +} + +template +class Velocity : public TensorFunction<1, dim> +{ +public: + Velocity() + : TensorFunction<1, dim>() + {} + virtual Tensor<1, dim> value(const Point &p) const override; +}; + +template +Tensor<1, dim> Velocity::value(const Point &p) const +{ + Tensor<1, dim> return_value; + return_value[0] = -M_PI * cos(M_PI * p[0]) * sin(M_PI * p[1]); + return_value[1] = -M_PI * sin(M_PI * p[0]) * cos(M_PI * p[1]); + return return_value; +} + +// @sect3{WGDarcyEquation class implementation} + +// @sect4{WGDarcyEquation::WGDarcyEquation} + +// In this constructor, we create a finite element space for vector valued +// functions, FE_RaviartThomas. We will need shape functions in +// this space to approximate discrete weak gradients. + +// FESystem defines finite element spaces in the interior and on +// edges of elements. Each of them gets an individual component. Others are the +// same as previous tutorial programs. +template +WGDarcyEquation::WGDarcyEquation() + : fe_rt(0) + , dof_handler_rt(triangulation) + , + + fe(FE_DGQ(0), 1, FE_FaceQ(0), 1) + , dof_handler(triangulation) + +{} + +// @sect4{WGDarcyEquation::make_grid} + +// We generate a mesh on the unit square domain and refine it. + +template +void WGDarcyEquation::make_grid() +{ + GridGenerator::hyper_cube(triangulation, 0, 1); + triangulation.refine_global(1); + + std::cout << " Number of active cells: " << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " << triangulation.n_cells() + << std::endl; +} + +// @sect4{WGDarcyEquation::setup_system} + +// After we create the mesh, we distribute degrees of freedom for the two +// DoFHandler objects. + +template +void WGDarcyEquation::setup_system() +{ + dof_handler_rt.distribute_dofs(fe_rt); + dof_handler.distribute_dofs(fe); + + std::cout << " Number of flux degrees of freedom: " + << dof_handler_rt.n_dofs() << std::endl; + + std::cout << " Number of pressure degrees of freedom: " + << dof_handler.n_dofs() << std::endl; + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + + { + constraints.clear(); + FEValuesExtractors::Scalar face(1); + ComponentMask face_pressure_mask = fe.component_mask(face); + VectorTools::interpolate_boundary_values( + dof_handler, 0, BoundaryValues(), constraints, face_pressure_mask); + constraints.close(); + } + + + // In the bilinear form, there is no integration term over faces + // between two neighboring cells, so we can just use + // DoFTools::make_sparsity_pattern to calculate the sparse + // matrix. + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + // solution.reinit(dof_handler.n_dofs()); + // system_rhs.reinit(dof_handler.n_dofs()); +} + +// @sect4{WGDarcyEquation::assemble_system} + +// First, we create quadrature points and FEValue objects for cells +// and faces. Then we allocate space for all cell matrices and the right-hand +// side vector. The following definitions have been explained in previous +// tutorials. +template +void WGDarcyEquation::assemble_system() +{ + QGauss quadrature_formula(fe_rt.degree + 1); + QGauss face_quadrature_formula(fe_rt.degree + 1); + const RightHandSide right_hand_side; + + // We define objects to evaluate values and + // gradients of shape functions at the quadrature points. + // Since we need shape functions and normal vectors on faces, we need + // FEFaceValues. + FEValues fe_values_rt(fe_rt, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values_rt(fe_rt, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + + const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + const unsigned int n_q_points = fe_values.get_quadrature().size(); + const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); + const unsigned int n_face_q_points = fe_face_values.get_quadrature().size(); + + std::vector local_dof_indices(dofs_per_cell); + + // We will construct these cell matrices to solve for the pressure. + FullMatrix cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt); + FullMatrix cell_matrix_F(dofs_per_cell, dofs_per_cell_rt); + FullMatrix cell_matrix_C(dofs_per_cell, dofs_per_cell_rt); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + Vector cell_solution(dofs_per_cell); + + const Coefficient coefficient; + std::vector> coefficient_values(n_q_points_rt); + + // We need FEValuesExtractors to access the @p interior and + // @p face component of the FESystem shape functions. + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar interior(0); + const FEValuesExtractors::Scalar face(1); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator cell_rt = + dof_handler_rt.begin_active(); + + // Here, we will calculate cell matrices used to construct the local matrix on + // each cell. We need shape functions for the Raviart-Thomas space as well, so + // we also loop over the corresponding velocity cell iterators. + for (; cell != endc; ++cell, ++cell_rt) + { + // On each cell, cell matrices are different, so in every loop, they need + // to be re-computed. + fe_values_rt.reinit(cell_rt); + fe_values.reinit(cell); + coefficient.value_list(fe_values_rt.get_quadrature_points(), + coefficient_values); + + // This cell matrix is the mass matrix for the Raviart-Thomas space. + // Hence, we need to loop over all the quadrature points + // for the velocity FEValues object. + cell_matrix_rt = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell_rt; ++i) + { + const Tensor<1, dim> phi_i_u = + fe_values_rt[velocities].value(i, q); + for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) + { + const Tensor<1, dim> phi_j_u = + fe_values_rt[velocities].value(j, q); + cell_matrix_rt(i, j) += + (phi_i_u * phi_j_u * fe_values_rt.JxW(q)); + } + } + } + // Next we take the inverse of this matrix by using + // gauss_jordan(). It will be used to calculate the + // coefficient matrix later. + cell_matrix_rt.gauss_jordan(); + + // From the introduction, we know that the right hand side + // is the difference between a face integral and a cell integral. + // Here, we approximate the negative of the contribution in the interior. + // Each component of this matrix is the integral of a product between a + // basis function of the polynomial space and the divergence of a basis + // function of the Raviart-Thomas space. These basis functions are defined + // in the interior. + cell_matrix_F = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const double phi_k_u_div = + fe_values_rt[velocities].divergence(k, q); + cell_matrix_F(i, k) -= (fe_values[interior].value(i, q) * + phi_k_u_div * fe_values.JxW(q)); + } + } + } + + // Now, we approximate the integral on faces. + // Each component is the integral of a product between a basis function of + // the polynomial space and the dot product of a basis function of the + // Raviart-Thomas space and the normal vector. So we loop over all the + // faces of the element and obtain the normal vector. + for (unsigned int face_n = 0; face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + fe_face_values.reinit(cell, face_n); + fe_face_values_rt.reinit(cell_rt, face_n); + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + const Tensor<1, dim> normal = fe_face_values.normal_vector(q); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_face_values_rt[velocities].value(k, q); + cell_matrix_F(i, k) += + (fe_face_values[face].value(i, q) * (phi_k_u * normal) * + fe_face_values.JxW(q)); + } + } + } + } + + // @p cell_matrix_C is matrix product between the inverse of mass matrix @p cell_matrix_rt and @p cell_matrix_F. + cell_matrix_C = 0; + cell_matrix_F.mmult(cell_matrix_C, cell_matrix_rt); + + // Element $a_{ij}$ of the local cell matrix $A$ is given by + // $\int_{E} \sum_{k,l} c_{ik} c_{jl} (\mathbf{K} \mathbf{w}_k) \cdot + // \mathbf{w}_l \mathrm{d}x.$ We have calculated coefficients $c$ in the + // previous step. + local_matrix = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_values_rt[velocities].value(k, q); + for (unsigned int l = 0; l < dofs_per_cell_rt; ++l) + { + const Tensor<1, dim> phi_l_u = + fe_values_rt[velocities].value(l, q); + local_matrix(i, j) += coefficient_values[q] * + cell_matrix_C[i][k] * phi_k_u * + cell_matrix_C[j][l] * phi_l_u * + fe_values_rt.JxW(q); + } + } + } + } + } + + // Next, we calculate the right hand side, $\int_{E} f q \mathrm{d}x$. + cell_rhs = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + cell_rhs(i) += + (fe_values[interior].value(i, q) * + right_hand_side.value(fe_values.quadrature_point(q)) * + fe_values.JxW(q)); + } + } + + // In this part, we distribute components of this local matrix into the + // system matrix and transfer components of the cell right-hand side into + // the system right hand side. + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global( + local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } +} + +// @sect4{WGDarcyEquation::solve} + +// Solving the system of the Darcy equation. Now, we have pressures in the +// interior and on the faces of all the cells. +template +void WGDarcyEquation::solve() +{ + SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm()); + SolverCG<> solver(solver_control); + solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); + constraints.distribute(solution); +} + +// @sect4{WGDarcyEquation::process_solution} + +// This part is to calculate the $L_2$ error of the pressure. +template +void WGDarcyEquation::process_solution() +{ + // Since we have two different spaces for finite elements in interior and on + // faces, if we want to calculate $L_2$ errors in interior, we need degrees of + // freedom only defined in cells. In FESystem, we have two + // components, the first one is for interior, the second one is for skeletons. + // fe.base_element(0) shows we only need degrees of freedom + // defined in cells. + DoFHandler interior_dof_handler(triangulation); + interior_dof_handler.distribute_dofs(fe.base_element(0)); + // We define a vector to extract pressures in cells. + // The size of the vector is the collective number of all degrees of freedom + // in the interior of all the elements. + Vector interior_solution(interior_dof_handler.n_dofs()); + { + // types::global_dof_index is used to know the global indices + // of degrees of freedom. So here, we get the global indices of local + // degrees of freedom and the global indices of interior degrees of freedom. + std::vector local_dof_indices(fe.dofs_per_cell); + std::vector interior_local_dof_indices( + fe.base_element(0).dofs_per_cell); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + interior_cell = interior_dof_handler.begin_active(); + + // In the loop of all cells and interior of the cell, + // we extract interior solutions from the global solution. + for (; cell != endc; ++cell, ++interior_cell) + { + cell->get_dof_indices(local_dof_indices); + interior_cell->get_dof_indices(interior_local_dof_indices); + + for (unsigned int i = 0; i < fe.base_element(0).dofs_per_cell; ++i) + interior_solution(interior_local_dof_indices[i]) = + solution(local_dof_indices[fe.component_to_system_index(0, i)]); + } + } + + // We define a vector that holds the norm of the error on each cell. + // Next, we use VectorTool::integrate_difference + // to compute the error in the $L_2$ norm on each cell. + // Finally, we get the global $L_2$ norm. + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(interior_dof_handler, + interior_solution, + Solution(), + difference_per_cell, + QGauss(fe.degree + 2), + VectorTools::L2_norm); + + const double L2_error = difference_per_cell.l2_norm(); + std::cout << "L2_error_pressure " << L2_error << std::endl; +} + +// @sect4{WGDarcyEquation::postprocess} + +// After we calculated the numerical pressure, we evaluate $L_2$ errors for the +// velocity on each cell and $L_2$ errors for the flux on faces. + +// We are going to evaluate velocities on each cell and calculate the difference +// between numerical and exact velocities. To calculate velocities, we need +// interior and face pressure values of each element, and some other cell +// matrices. + +template +void WGDarcyEquation::postprocess() +{ + QGauss quadrature_formula(fe_rt.degree + 1); + QGauss face_quadrature_formula(fe_rt.degree + 1); + + FEValues fe_values_rt(fe_rt, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_face_values_rt(fe_rt, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); + const unsigned int n_q_points = fe_values.get_quadrature().size(); + const unsigned int n_face_q_points = fe_face_values.get_quadrature().size(); + const unsigned int n_face_q_points_rt = + fe_face_values_rt.get_quadrature().size(); + + + std::vector local_dof_indices(dofs_per_cell); + FullMatrix cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt); + FullMatrix cell_matrix_F(dofs_per_cell, dofs_per_cell_rt); + FullMatrix cell_matrix_C(dofs_per_cell, dofs_per_cell_rt); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix_D(dofs_per_cell_rt, dofs_per_cell_rt); + FullMatrix cell_matrix_E(dofs_per_cell_rt, dofs_per_cell_rt); + Vector cell_rhs(dofs_per_cell); + Vector cell_solution(dofs_per_cell); + Tensor<1, dim> velocity_cell; + Tensor<1, dim> velocity_face; + Tensor<1, dim> exact_velocity_face; + double L2_err_velocity_cell_sqr_global; + L2_err_velocity_cell_sqr_global = 0; + double L2_err_flux_sqr; + L2_err_flux_sqr = 0; + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + + typename DoFHandler::active_cell_iterator cell_rt = + dof_handler_rt.begin_active(); + + const Coefficient coefficient; + std::vector> coefficient_values(n_q_points_rt); + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + const FEValuesExtractors::Scalar interior(0); + const FEValuesExtractors::Scalar face(1); + + Velocity exact_velocity; + + // In the loop over all cells, we will calculate $L_2$ errors of velocity and + // flux. + + // First, we calculate the $L_2$ velocity error. + // In the introduction, we explained how to calculate the numerical velocity + // on the cell. We need the pressure solution values on each cell, + // coefficients of the Gram matrix and coefficients of the $L_2$ projection. + // We have already calculated the global solution, so we will extract the cell + // solution from the global solution. The coefficients of the Gram matrix have + // been calculated when we assembled the system matrix for the pressures. We + // will do the same way here. For the coefficients of the projection, we do + // matrix multiplication, i.e., the inverse of the Gram matrix times the + // matrix with $(\mathbf{K} \mathbf{w}, \mathbf{w})$ as components. Then, we + // multiply all these coefficients and call them beta. The numerical velocity + // is the product of beta and the basis functions of the Raviart-Thomas space. + for (; cell != endc; ++cell, ++cell_rt) + { + fe_values_rt.reinit(cell_rt); + fe_values.reinit(cell); + coefficient.value_list(fe_values_rt.get_quadrature_points(), + coefficient_values); + + // The component of this cell_matrix_E is the integral of + // $(\mathbf{K} \mathbf{w}, \mathbf{w})$. cell_matrix_rt is + // the Gram matrix. + cell_matrix_E = 0; + cell_matrix_rt = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell_rt; ++i) + { + const Tensor<1, dim> phi_i_u = + fe_values_rt[velocities].value(i, q); + + for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) + { + const Tensor<1, dim> phi_j_u = + fe_values_rt[velocities].value(j, q); + + cell_matrix_E(i, j) += (coefficient_values[q] * phi_j_u * + phi_i_u * fe_values_rt.JxW(q)); + cell_matrix_rt(i, j) += + (phi_i_u * phi_j_u * fe_values_rt.JxW(q)); + } + } + } + + // We take the inverse of the Gram matrix, take matrix multiplication and + // get the matrix with coefficients of projection. + cell_matrix_D = 0; + cell_matrix_rt.gauss_jordan(); + cell_matrix_rt.mmult(cell_matrix_D, cell_matrix_E); + + // This cell matrix will be used to calculate the coefficients of the Gram + // matrix. This part is the same as the part in evaluating pressure. + cell_matrix_F = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const double phi_k_u_div = + fe_values_rt[velocities].divergence(k, q); + cell_matrix_F(i, k) -= (fe_values[interior].value(i, q) * + phi_k_u_div * fe_values.JxW(q)); + } + } + } + + for (unsigned int face_n = 0; face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + fe_face_values.reinit(cell, face_n); + fe_face_values_rt.reinit(cell_rt, face_n); + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + const Tensor<1, dim> normal = fe_face_values.normal_vector(q); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_face_values_rt[velocities].value(k, q); + cell_matrix_F(i, k) += + (fe_face_values[face].value(i, q) * (phi_k_u * normal) * + fe_face_values.JxW(q)); + } + } + } + } + cell_matrix_C = 0; + cell_matrix_F.mmult(cell_matrix_C, cell_matrix_rt); + + // This is to extract pressure values of the element. + cell->get_dof_indices(local_dof_indices); + cell_solution = 0; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + cell_solution(i) = solution(local_dof_indices[i]); + } + + // From previous calculations we obtained all the coefficients needed to + // calculate beta. + Vector beta(dofs_per_cell_rt); + beta = 0; + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + for (unsigned int j = 0; j < dofs_per_cell_rt; ++j) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + beta(k) += -(cell_solution(i) * cell_matrix_C(i, j) * + cell_matrix_D(k, j)); + } + } + } + + // Now, we can calculate the numerical velocity at each quadrature point + // and compute the $L_2$ error on each cell. + double L2_err_velocity_cell_sqr_local; + double difference_velocity_cell_sqr; + L2_err_velocity_cell_sqr_local = 0; + velocity_cell = 0; + for (unsigned int q = 0; q < n_q_points_rt; ++q) + { + difference_velocity_cell_sqr = 0; + velocity_cell = 0; + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_values_rt[velocities].value(k, q); + velocity_cell += beta(k) * phi_k_u; + } + difference_velocity_cell_sqr = + (velocity_cell - + exact_velocity.value(fe_values_rt.quadrature_point(q))) * + (velocity_cell - + exact_velocity.value(fe_values_rt.quadrature_point(q))); + L2_err_velocity_cell_sqr_local += + difference_velocity_cell_sqr * fe_values_rt.JxW(q); + } + + L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local; + + // For reconstructing the flux we need the size of cells and faces. Since + // fluxes are calculated on faces, we have the loop over all four faces of + // each cell. To calculate face velocity, we use the coefficient beta we + // have calculated previously. Then, we calculate the squared velocity + // error in normal direction. Finally, we calculate $L_2$ flux error on + // the cell and add it to the global error. + double difference_velocity_face_sqr; + double L2_err_flux_face_sqr_local; + double err_flux_each_face; + double err_flux_face; + L2_err_flux_face_sqr_local = 0; + err_flux_face = 0; + const double cell_area = cell->measure(); + for (unsigned int face_n = 0; face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + const double face_length = cell->face(face_n)->measure(); + fe_face_values.reinit(cell, face_n); + fe_face_values_rt.reinit(cell_rt, face_n); + L2_err_flux_face_sqr_local = 0; + err_flux_each_face = 0; + for (unsigned int q = 0; q < n_face_q_points_rt; ++q) + { + difference_velocity_face_sqr = 0; + velocity_face = 0; + const Tensor<1, dim> normal = fe_face_values.normal_vector(q); + for (unsigned int k = 0; k < dofs_per_cell_rt; ++k) + { + const Tensor<1, dim> phi_k_u = + fe_face_values_rt[velocities].value(k, q); + velocity_face += beta(k) * phi_k_u; + } + exact_velocity_face = + exact_velocity.value(fe_face_values_rt.quadrature_point(q)); + difference_velocity_face_sqr = + (velocity_face * normal - exact_velocity_face * normal) * + (velocity_face * normal - exact_velocity_face * normal); + L2_err_flux_face_sqr_local += + difference_velocity_face_sqr * fe_face_values_rt.JxW(q); + } + err_flux_each_face = + L2_err_flux_face_sqr_local / (face_length) * (cell_area); + err_flux_face += err_flux_each_face; + } + L2_err_flux_sqr += err_flux_face; + } + + // After adding up errors over all cells, we take square root and get the + // $L_2$ errors of velocity and flux. + const double L2_err_velocity_cell = + std::sqrt(L2_err_velocity_cell_sqr_global); + std::cout << "L2_error_vel " << L2_err_velocity_cell << std::endl; + const double L2_err_flux_face = std::sqrt(L2_err_flux_sqr); + std::cout << "L2_error_flux " << L2_err_flux_face << std::endl; +} + + +// @sect4{WGDarcyEquation::output_results} + +// We have 2 sets of results to output: the interior solution +// and the skeleton solution. We use DataOut to visualize interior +// results. The graphical output for the skeleton results is done by using the +// DataOutFaces class. +template +void WGDarcyEquation::output_results() const +{ + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "Pressure_Interior"); + data_out.build_patches(fe.degree); + std::ofstream output("Pressure_Interior.vtk"); + data_out.write_vtk(output); + + DataOutFaces data_out_face(false); + std::vector + face_component_type(2, DataComponentInterpretation::component_is_scalar); + data_out_face.add_data_vector(dof_handler, + solution, + "Pressure_Edge", + face_component_type); + data_out_face.build_patches(fe.degree); + std::ofstream face_output("Pressure_Edge.vtk"); + data_out_face.write_vtk(face_output); +} + + +// @sect4{WGDarcyEquation::run} + +// This is the final function of the main class. It calls the other functions of +// our class. +template +void WGDarcyEquation::run() +{ + std::cout << "Solving problem in " << dim << " space dimensions." + << std::endl; + make_grid(); + setup_system(); + assemble_system(); + solve(); + process_solution(); + postprocess(); + output_results(); +} + +// @sect3{The main function} + +// This is the main function. We can change the dimension here to run in 3d. +int main() +{ + deallog.depth_console(2); + WGDarcyEquation<2> WGDarcyEquationTest; + WGDarcyEquationTest.run(); + + return 0; +}