]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-61. 6455/head
authorZhuoran Wang <zhrwang@rams.colostate.edu>
Fri, 27 Apr 2018 16:31:49 +0000 (16:31 +0000)
committersophy1029 <sophy1029@hotmail.com>
Wed, 30 Jan 2019 17:07:10 +0000 (17:07 +0000)
This tutorial program adds a code that solves the Poisson equation
using the weak Galerkin formulation.

15 files changed:
examples/step-61/CMakeLists.txt [new file with mode: 0644]
examples/step-61/doc/builds-on [new file with mode: 0644]
examples/step-61/doc/intro.dox [new file with mode: 0644]
examples/step-61/doc/kind [new file with mode: 0644]
examples/step-61/doc/results.dox [new file with mode: 0644]
examples/step-61/doc/step-61.wg000_2d_2.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg000_2d_4.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg000_3d_2.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg000_3d_4.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg111_2d_4.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg111_3d_4.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg222_2d_5.png [new file with mode: 0644]
examples/step-61/doc/step-61.wg222_3d_5.png [new file with mode: 0644]
examples/step-61/doc/tooltip [new file with mode: 0644]
examples/step-61/step-61.cc [new file with mode: 0644]

diff --git a/examples/step-61/CMakeLists.txt b/examples/step-61/CMakeLists.txt
new file mode 100644 (file)
index 0000000..83a33d0
--- /dev/null
@@ -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 (file)
index 0000000..2bb1b6e
--- /dev/null
@@ -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 (file)
index 0000000..f5ca096
--- /dev/null
@@ -0,0 +1,279 @@
+<br>
+
+<i>
+This program was contributed by Zhuoran Wang.
+</i>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+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.
+
+<h3>  Weak Galerkin finite element methods</h3>
+
+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: <i>A weak Galerkin finite element method for second order elliptic problems</i>, 
+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.
+
+
+<h3> WGFEM applied to the Poisson equation </h3>
+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.
+
+
+<h3> Weak Galerkin scheme </h3>
+
+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$. 
+
+<h3> Assembling the linear system </h3>
+
+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 <code>cell_matrix_C</code>, 
+$\mathbf{M}_E$ as <code>cell_matrix_rt</code>, 
+$\mathbf{F}_E$ as <code>cell_matrix_F</code>.
+
+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.
+
+<h3> Post-processing and $L_2$-errors </h3>
+
+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 <code>cell_matrix_D</code>, 
+$
+\left[
+    \begin{matrix}
+     \left( \mathbf{Kw}_j,\mathbf{w}_k \right)
+    \end{matrix}
+    \right]
+$
+is named <code>cell_matrix_E</code>.
+
+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 
+<code>beta</code> 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 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-61/doc/results.dox b/examples/step-61/doc/results.dox
new file mode 100644 (file)
index 0000000..b71c57a
--- /dev/null
@@ -0,0 +1,88 @@
+<h1>Results</h1>
+
+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]})$.
+
+<h3>Test results on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$</h3>
+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. 
+
+<table align="center">
+  <tr>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_2d_2.png" alt=""></td>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_3d_2.png" alt=""></td>
+  </tr>
+  <tr>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_2d_4.png" alt=""></td>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_3d_4.png" alt=""></td>
+  </tr>
+</table>
+
+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.
+
+<h4>Convergence table</h4>
+
+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.
+
+
+<h3>Test results on $\mbox{WG}(Q_1,Q_1;RT_{[1]})$</h3>
+
+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.
+
+<table align="center">
+  <tr>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg111_2d_4.png" alt=""></td>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg111_3d_4.png" alt=""></td>
+  </tr>
+</table>
+
+<h4>Convergence table</h4>
+
+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.
+
+<h3>Test results on $WG(Q_2,Q_2;RT_{[2]})$</h3>
+
+These are interior pressures and face pressures implemented on $WG(Q_2,Q_2;RT_{[2]})$, with mesh size $h = 1/32$.
+
+<table align="center">
+  <tr>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg222_2d_5.png" alt=""></td>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg222_3d_5.png" alt=""></td>
+  </tr>
+</table>
+
+<h4>Convergence table</h4>
+
+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 (file)
index 0000000..bf3be23
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 (file)
index 0000000..fe5e5bc
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 (file)
index 0000000..73b6c7b
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 (file)
index 0000000..a76e480
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 (file)
index 0000000..23258e0
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 (file)
index 0000000..c3f0f2f
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 (file)
index 0000000..1be5805
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 (file)
index 0000000..c79c7fb
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 (file)
index 0000000..1f69a06
--- /dev/null
@@ -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 (file)
index 0000000..90ed9c0
--- /dev/null
@@ -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 <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/point.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/component_mask.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/data_out_faces.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+
+#include <fstream>
+#include <iostream>
+
+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 <int dim>
+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<dim> triangulation;
+
+  AffineConstraints<double> constraints;
+
+  FE_RaviartThomas<dim> fe_rt;
+  DoFHandler<dim>       dof_handler_rt;
+
+  // The finite element system is used for interior and face solutions.
+  FESystem<dim>   fe;
+  DoFHandler<dim> dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  Vector<double> solution;
+  Vector<double> 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 <int dim>
+class Coefficient : public TensorFunction<2, dim>
+{
+public:
+  Coefficient()
+    : TensorFunction<2, dim>()
+  {}
+  virtual void value_list(const std::vector<Point<dim>> &points,
+                          std::vector<Tensor<2, dim>> &  values) const override;
+};
+
+template <int dim>
+void Coefficient<dim>::value_list(const std::vector<Point<dim>> &points,
+                                  std::vector<Tensor<2, dim>> &  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 <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues()
+    : Function<dim>(2)
+  {}
+  virtual double value(const Point<dim> & p,
+                       const unsigned int component = 0) const override;
+};
+
+template <int dim>
+double BoundaryValues<dim>::value(const Point<dim> & /*p*/,
+                                  const unsigned int /*component*/) const
+{
+  return 0;
+}
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide()
+    : Function<dim>()
+  {}
+  virtual double value(const Point<dim> & p,
+                       const unsigned int component = 0) const;
+};
+
+template <int dim>
+double RightHandSide<dim>::value(const Point<dim> &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 <int dim>
+class Solution : public Function<dim>
+{
+public:
+  Solution()
+    : Function<dim>(1)
+  {}
+  virtual double value(const Point<dim> &p, const unsigned int) const;
+};
+
+template <int dim>
+double Solution<dim>::value(const Point<dim> &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 <int dim>
+class Velocity : public TensorFunction<1, dim>
+{
+public:
+  Velocity()
+    : TensorFunction<1, dim>()
+  {}
+  virtual Tensor<1, dim> value(const Point<dim> &p) const override;
+};
+
+template <int dim>
+Tensor<1, dim> Velocity<dim>::value(const Point<dim> &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, <code>FE_RaviartThomas</code>. We will need shape functions in
+// this space to approximate discrete weak gradients.
+
+// <code>FESystem</code> 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 <int dim>
+WGDarcyEquation<dim>::WGDarcyEquation()
+  : fe_rt(0)
+  , dof_handler_rt(triangulation)
+  ,
+
+  fe(FE_DGQ<dim>(0), 1, FE_FaceQ<dim>(0), 1)
+  , dof_handler(triangulation)
+
+{}
+
+// @sect4{WGDarcyEquation::make_grid}
+
+// We generate a mesh on the unit square domain and refine it.
+
+template <int dim>
+void WGDarcyEquation<dim>::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
+// <code>DoFHandler</code> objects.
+
+template <int dim>
+void WGDarcyEquation<dim>::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<dim>(), 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
+  // <code>DoFTools::make_sparsity_pattern</code> 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 <code>FEValue</code> 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 <int dim>
+void WGDarcyEquation<dim>::assemble_system()
+{
+  QGauss<dim>              quadrature_formula(fe_rt.degree + 1);
+  QGauss<dim - 1>          face_quadrature_formula(fe_rt.degree + 1);
+  const RightHandSide<dim> 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
+  // <code>FEFaceValues</code>.
+  FEValues<dim> fe_values_rt(fe_rt,
+                             quadrature_formula,
+                             update_values | update_gradients |
+                               update_quadrature_points | update_JxW_values);
+
+  FEValues<dim> fe_values(fe,
+                          quadrature_formula,
+                          update_values | update_quadrature_points |
+                            update_JxW_values);
+
+  FEFaceValues<dim> fe_face_values(fe,
+                                   face_quadrature_formula,
+                                   update_values | update_normal_vectors |
+                                     update_quadrature_points |
+                                     update_JxW_values);
+
+  FEFaceValues<dim> 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<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  // We will construct these cell matrices to solve for the pressure.
+  FullMatrix<double> cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt);
+  FullMatrix<double> cell_matrix_F(dofs_per_cell, dofs_per_cell_rt);
+  FullMatrix<double> cell_matrix_C(dofs_per_cell, dofs_per_cell_rt);
+  FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+  Vector<double>     cell_solution(dofs_per_cell);
+
+  const Coefficient<dim>      coefficient;
+  std::vector<Tensor<2, dim>> coefficient_values(n_q_points_rt);
+
+  // We need <code>FEValuesExtractors</code> 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<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  typename DoFHandler<dim>::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
+      // <code>gauss_jordan()</code>. 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<dim>::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<dim>::solve}
+
+// Solving the system of the Darcy equation. Now, we have pressures in the
+// interior and on the faces of all the cells.
+template <int dim>
+void WGDarcyEquation<dim>::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<dim>::process_solution}
+
+// This part is to calculate the $L_2$ error of the pressure.
+template <int dim>
+void WGDarcyEquation<dim>::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 <code>FESystem</code>, we have two
+  // components, the first one is for interior, the second one is for skeletons.
+  // <code>fe.base_element(0)</code> shows we only need degrees of freedom
+  // defined in cells.
+  DoFHandler<dim> 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<double> interior_solution(interior_dof_handler.n_dofs());
+  {
+    // <code>types::global_dof_index</code> 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<types::global_dof_index> local_dof_indices(fe.dofs_per_cell);
+    std::vector<types::global_dof_index> interior_local_dof_indices(
+      fe.base_element(0).dofs_per_cell);
+    typename DoFHandler<dim>::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 <code>VectorTool::integrate_difference</code>
+  // to compute the error in the $L_2$ norm on each cell.
+  // Finally, we get the global $L_2$ norm.
+  Vector<float> difference_per_cell(triangulation.n_active_cells());
+  VectorTools::integrate_difference(interior_dof_handler,
+                                    interior_solution,
+                                    Solution<dim>(),
+                                    difference_per_cell,
+                                    QGauss<dim>(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<dim>::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 <int dim>
+void WGDarcyEquation<dim>::postprocess()
+{
+  QGauss<dim>     quadrature_formula(fe_rt.degree + 1);
+  QGauss<dim - 1> face_quadrature_formula(fe_rt.degree + 1);
+
+  FEValues<dim> fe_values_rt(fe_rt,
+                             quadrature_formula,
+                             update_values | update_gradients |
+                               update_quadrature_points | update_JxW_values);
+
+  FEValues<dim> fe_values(fe,
+                          quadrature_formula,
+                          update_values | update_quadrature_points |
+                            update_JxW_values);
+
+  FEFaceValues<dim> fe_face_values(fe,
+                                   face_quadrature_formula,
+                                   update_values | update_normal_vectors |
+                                     update_quadrature_points |
+                                     update_JxW_values);
+
+  FEFaceValues<dim> 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<types::global_dof_index> local_dof_indices(dofs_per_cell);
+  FullMatrix<double> cell_matrix_rt(dofs_per_cell_rt, dofs_per_cell_rt);
+  FullMatrix<double> cell_matrix_F(dofs_per_cell, dofs_per_cell_rt);
+  FullMatrix<double> cell_matrix_C(dofs_per_cell, dofs_per_cell_rt);
+  FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> cell_matrix_D(dofs_per_cell_rt, dofs_per_cell_rt);
+  FullMatrix<double> cell_matrix_E(dofs_per_cell_rt, dofs_per_cell_rt);
+  Vector<double>     cell_rhs(dofs_per_cell);
+  Vector<double>     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<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+
+  typename DoFHandler<dim>::active_cell_iterator cell_rt =
+    dof_handler_rt.begin_active();
+
+  const Coefficient<dim>           coefficient;
+  std::vector<Tensor<2, dim>>      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<dim> 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 <code>cell_matrix_E</code> is the integral of
+      // $(\mathbf{K} \mathbf{w}, \mathbf{w})$. <code>cell_matrix_rt</code> 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<dim>::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<double> 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<dim>::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 <code>DataOut</code> to visualize interior
+// results. The graphical output for the skeleton results is done by using the
+// <code>DataOutFaces</code> class.
+template <int dim>
+void WGDarcyEquation<dim>::output_results() const
+{
+  DataOut<dim> 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<dim> data_out_face(false);
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    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 <int dim>
+void WGDarcyEquation<dim>::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 <code>main</code> 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;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.