]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add step-74
authorJiaqi Zhang <zjiaqi@vt.edu>
Wed, 12 Aug 2020 00:13:40 +0000 (20:13 -0400)
committerJiaqi Zhang <zjiaqi@vt.edu>
Mon, 30 Nov 2020 17:06:40 +0000 (12:06 -0500)
doc/doxygen/references.bib
examples/step-74/CMakeLists.txt [new file with mode: 0644]
examples/step-74/doc/builds-on [new file with mode: 0644]
examples/step-74/doc/intro.dox [new file with mode: 0644]
examples/step-74/doc/kind [new file with mode: 0644]
examples/step-74/doc/results.dox [new file with mode: 0644]
examples/step-74/doc/tooltip [new file with mode: 0644]
examples/step-74/step-74.cc [new file with mode: 0644]

index 5ed6044eb3b9d8b92bb27bbccb6b04b47d2ad51b..8daefc86529c9fb547ce85fc59c4668f9168c66b 100644 (file)
@@ -823,6 +823,31 @@ year = {2008},
   title = {Improved energy estimates for interior penalty, constrained and discontinuous {G}alerkin methods for elliptic problems. {P}art {I}}.
 }
 
+
+% ------------------------------------
+% Step 70
+% ------------------------------------
+
+@book{di2011mathematical,
+  title={Mathematical aspects of discontinuous Galerkin methods},
+  author={Di Pietro, Daniele Antonio and Ern, Alexandre},
+  volume={69},
+  year={2011},
+  publisher={Springer Science \& Business Media}
+}
+
+@article{ainsworth2007posteriori,
+  title={A posteriori error estimation for discontinuous Galerkin finite element approximation},
+  author={Ainsworth, Mark},
+  journal={SIAM Journal on Numerical Analysis},
+  volume={45},
+  number={4},
+  pages={1777--1798},
+  year={2007},
+  publisher={SIAM}
+}
+
+
 % ------------------------------------
 % References used elsewhere
 % ------------------------------------
diff --git a/examples/step-74/CMakeLists.txt b/examples/step-74/CMakeLists.txt
new file mode 100644 (file)
index 0000000..6d16682
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-74 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-74")
+
+# 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.3.0
+  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-74/doc/builds-on b/examples/step-74/doc/builds-on
new file mode 100644 (file)
index 0000000..20c3064
--- /dev/null
@@ -0,0 +1 @@
+step-12
diff --git a/examples/step-74/doc/intro.dox b/examples/step-74/doc/intro.dox
new file mode 100644 (file)
index 0000000..7bd4844
--- /dev/null
@@ -0,0 +1,125 @@
+<br>
+
+<i>
+This program was contributed by Timo Heister and Jiaqi Zhang.
+<br>
+This material is based upon work partly supported by the National
+Science Foundation Award DMS-2028346, OAC-2015848, EAR-1925575, by the Computational
+Infrastructure in Geodynamics initiative (CIG), through the NSF under Award
+EAR-0949446 and EAR-1550901 and The University of California -- Davis.
+</i>
+
+
+<a name="Intro"></a>
+<h1><em>Symmetric interior penalty Galerkin</em> (SIPG) method for Poisson's equation</h1>
+
+<h3>Overview</h3>
+In this tutorial, we display the usage of the FEInterfaceValues class,
+which is designed for assembling face terms arising from discontinuous Galerkin (DG) methods.
+The FEInterfaceValues class provides an easy way to obtain the jump and the average of the solution across cell faces.
+This tutorial includes the following topics.
+<ol>
+  <li> The SIPG method for Poisson's equation, which has been used in step-39 and step-59.
+  <li> Assembling of face terms using FEInterfaceValues and the system matrix using MeshWorker::mesh_loop(), which is similar to step-12.
+  <li> Adaptive mesh refinement using an error estimator.
+  <li> Two test cases: convergence test for a smooth function and adaptive mesh refinement test for a singular solution.
+</ol>
+
+<h3>The equation</h3>
+In this example, we consider Poisson's equation with a Dirichlet boundary condition
+@f[
+- \nu  \Delta u = f  \qquad   \mbox{in } \Omega,
+@f]
+subject to the boundary condition
+@f[
+u = g_D \qquad \mbox{on } \partial \Omega.
+@f]
+For simplicity, we assume that the diffusion coefficient $\nu$ is constant here.
+Note that if $\nu$ is discontinuous, we need to take this into account when computing jump terms
+on cell faces.
+
+We denote the mesh by $\Gamma_h$, and $K\in\Gamma_h$ is a mesh cell.
+The sets of interior and boundary faces are denoted by $F^i_h$ and $F^b_h$
+respectively. Let $K^0$ and $K^1$ be the two cells sharing a face $f\in F_h^i$,
+and $\mathbf n$ be the outer normal vector of $K^0$. Then the jump and average
+operators are given by
+@f[
+\jump{v} = v^0 - v^1
+@f]
+and
+@f[
+\average{v} = \frac{v^0 + v^1}{2}
+@f]
+respectively. Note that when $f\in \partial Omega$, we define $\jump{v} = v$ and
+$\average{v}=v$.
+The discretization using the SIPG is given by the following weak formula
+(more details can be found in @cite di2011mathematical and the references therein)
+@f{multline*}
+  \sum_{K\in \Gamma_h} (\nabla v_h, \nu \nabla u_h)_K
+  \\
+  - \sum_{F \in F_h^i} \biggl\{
+    \bigl< \jump{v_h}, \nu\average{ \nabla u_h} \cdot  \mathbf n \bigr>_F
+   +\bigl<\average{ \nabla v_h }\cdot \mathbf n,\nu\jump{u_h}\bigr>_F
+   -\bigl<\jump{v_h},\nu \sigma \jump{u_h} \bigr>_F
+  \biggr\}
+  \\
+  - \sum_{F \in F_h^b} \biggl\{
+    \bigl<v_h, \nv \nabla u_h\cdot \mathbf n \bigr>_F
+  + \bigl< \nabla v_h \cdot \mathbb n , \nu u_h\bigr>_F
+  - \bigl< v_h,\nu \sigma u_h\bigr>_F
+  \biggr\}
+  \\
+  = (v_h, f)_\Omega
+  - \sum_{F \in F_h^b} \biggl\{
+    \bigl< \nabla v_h \cdot \mathbf n, \nu g_D\bigr>_F - \bigl<v_h,\nu \sigma g_D\bigr>_F
+  \biggr\}.
+@f}
+
+<h3>The penalty parameter</h3>
+The parameter is defined as $\sigma = \gamma/h_f$, where $h_f$ a local length scale associated
+with the cell face, here we choose the length of the cell in the direction normal to the face,
+and $\gamma$ is the penalization constant. The lower bound of $\gamma$ @cite ainsworth2007posteriori
+is given by
+@f[
+\gamma > 4\max_{K\in \Gamma_h}\rho\left( S_K \right).
+@f]
+Here $[S_K]_{i,j} = (\nu \nabla \phi_i, \nabla\phi_j)_K$.
+
+To ensure the discrete coercivity, the penalization constant has to be large enough.
+There is no consensus in the literature on how to determine $\gamma$ in practice. One can just pick a large constant,
+while other options could be the multiples of $(p+1)^2$ or $p(p+1)$. In this code,
+we follow step-39 and use $\gamma = p(p+1)$.
+
+<h3>Posteriori error estimator</h3>
+In this example, we use the error estimator by Karakashian and Pascal @cite karakashian2003posteriori
+@f[
+\eta^2 = \sum_{K \in \Gamma_h} \eta^2_{K} +  \sum_{f_i \in F^i_h}  \eta^2_{f_i} + \sum_{f_b \in F^i_b}\eta^2_{f_b}
+@f]
+where
+@f[
+\eta^2_{K} = h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2
+@f]
+@f[
+\eta^2_{f_i} = \sigma \left\| \jump{u_h}  \right\|_f^2   +  h_f \left\|  \jump{\nu \nabla u_h} \cdot \mathbf n   \right\|_f^2
+@f]
+@f[
+\eta_{f_b}^2 =  \sigma \left\| u_h-g_D \right\|_f^2
+@f]
+The only difference is that we use $\sigma = \gamma/h_f$ instead of $\gamma^2/h_f$ for the jump terms of $u$ (the first term in $\eta^2_{f_i}$ and $\eta_{f_b}^2$).
+
+In each cell $K$, we compute
+@f[
+\eta_{c}^2 = h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2,
+@f]
+@f[
+\eta_{f}^2 = \sum_{f\in \partial K}\lbrace \sigma \left\| \jump{u_h}  \right\|_f^2   +  h_f \left\|  \jump{\nu \nabla u_h} \cdot \mathbf n  \right\|_f^2 \rbrace,
+@f]
+@f[
+\eta_{b}^2 = \sum_{f\in \partial K \cap \partial \Omega}  \sigma \left\| (u_h -g_D)  \right\|_f^2.
+@f]
+Then the error estimate square per cell is
+@f[
+\eta_{local}^2 =\eta_{c}^2+0.5\eta_{f}^2+\eta_{b}^2.
+@f]
+Note that we compute $\eta_{local}^2$ instead of $\eta_{local}^2$ to simplify the implementation.
+The error estimate square per cell is store in a global vector, whose $L_1$ norm is equal to $\eta^2$.
diff --git a/examples/step-74/doc/kind b/examples/step-74/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-74/doc/results.dox b/examples/step-74/doc/results.dox
new file mode 100644 (file)
index 0000000..9ed00fd
--- /dev/null
@@ -0,0 +1,195 @@
+<h1>Results</h1>
+
+The output of this program consist of the console output and
+solutions in vtu format.
+
+Convergence rate for the smooth case with polynomial degree 3:
+@code
+DEAL::Cycle 0
+DEAL::Number of active cells:       16
+DEAL::Number of degrees of freedom: 256
+   Solving system...
+DEAL::Writing solution to <sol_Q3-00.vtu>
+   Error in the L2 norm       :     0.00193285
+   Error in the H1 seminorm   :     0.106087
+   Error in the enery norm    :     0.150625
+
+DEAL::Cycle 1
+DEAL::Number of active cells:       64
+DEAL::Number of degrees of freedom: 1024
+   Solving system...
+DEAL::Writing solution to <sol_Q3-01.vtu>
+   Error in the L2 norm       :     9.60497e-05
+   Error in the H1 seminorm   :     0.0089954
+   Error in the enery norm    :     0.0113265
+
+DEAL::Cycle 2
+DEAL::Number of active cells:       256
+DEAL::Number of degrees of freedom: 4096
+   Solving system...
+DEAL::Writing solution to <sol_Q3-02.vtu>
+   Error in the L2 norm       :     5.60638e-06
+   Error in the H1 seminorm   :     0.000901791
+   Error in the enery norm    :     0.000973568
+
+DEAL::Cycle 3
+DEAL::Number of active cells:       1024
+DEAL::Number of degrees of freedom: 16384
+   Solving system...
+DEAL::Writing solution to <sol_Q3-03.vtu>
+   Error in the L2 norm       :     3.48399e-07
+   Error in the H1 seminorm   :     0.000107055
+   Error in the enery norm    :     0.000108757
+
+DEAL::Cycle 4
+DEAL::Number of active cells:       4096
+DEAL::Number of degrees of freedom: 65536
+   Solving system...
+DEAL::Writing solution to <sol_Q3-04.vtu>
+   Error in the L2 norm       :     2.17928e-08
+   Error in the H1 seminorm   :     1.32657e-05
+   Error in the enery norm    :     1.33065e-05
+
+DEAL::Cycle 5
+DEAL::Number of active cells:       16384
+DEAL::Number of degrees of freedom: 262144
+   Solving system...
+DEAL::Writing solution to <sol_Q3-05.vtu>
+   Error in the L2 norm       :     1.36274e-09
+   Error in the H1 seminorm   :     1.65576e-06
+   Error in the enery norm    :     1.65683e-06
+
+degree = 3
+| cycle | cells | dofs   | L2        | L2...red.rate.log2 | H1        | H1...red.rate.log2 | Energy    |
+| 0     | 16    | 256    | 1.933e-03 | -                  | 1.061e-01 | -                  | 1.506e-01 |
+| 1     | 64    | 1024   | 9.605e-05 | 4.33               | 8.995e-03 | 3.56               | 1.133e-02 |
+| 2     | 256   | 4096   | 5.606e-06 | 4.10               | 9.018e-04 | 3.32               | 9.736e-04 |
+| 3     | 1024  | 16384  | 3.484e-07 | 4.01               | 1.071e-04 | 3.07               | 1.088e-04 |
+| 4     | 4096  | 65536  | 2.179e-08 | 4.00               | 1.327e-05 | 3.01               | 1.331e-05 |
+| 5     | 16384 | 262144 | 1.363e-09 | 4.00               | 1.656e-06 | 3.00               | 1.657e-06 |
+@endcode
+Thearetically, for polynomial degree $p$, the order of convergece in $L_2$
+norm and $H_1$ seminorm should be $p+1$ and $p$, respectively. Our numerical
+results are in good agreement with the theories.
+
+Errors of the singular case with polynomial degree 3.
+@code
+DEAL::Cycle 0
+DEAL::Number of active cells:       12
+DEAL::Number of degrees of freedom: 192
+   Solving system...
+DEAL::Writing solution to <sol_Q3-00.vtu>
+   Error in the L2 norm       :     0.00278279
+   Error in the H1 seminorm   :     0.0748987
+   Error in the enery norm    :     0.106218
+
+DEAL::Cycle 1
+DEAL::Number of active cells:       15
+DEAL::Number of degrees of freedom: 240
+   Solving system...
+DEAL::Writing solution to <sol_Q3-01.vtu>
+   Error in the L2 norm       :     0.00179741
+   Error in the H1 seminorm   :     0.0568531
+   Error in the enery norm    :     0.0815378
+
+DEAL::Cycle 2
+DEAL::Number of active cells:       18
+DEAL::Number of degrees of freedom: 288
+   Solving system...
+DEAL::Writing solution to <sol_Q3-02.vtu>
+   Error in the L2 norm       :     0.00171775
+   Error in the H1 seminorm   :     0.0598664
+   Error in the enery norm    :     0.0850871
+
+DEAL::Cycle 3
+DEAL::Number of active cells:       21
+DEAL::Number of degrees of freedom: 336
+   Solving system...
+DEAL::Writing solution to <sol_Q3-03.vtu>
+   Error in the L2 norm       :     0.000939329
+   Error in the H1 seminorm   :     0.0470787
+   Error in the enery norm    :     0.0668118
+
+DEAL::Cycle 4
+DEAL::Number of active cells:       27
+DEAL::Number of degrees of freedom: 432
+   Solving system...
+DEAL::Writing solution to <sol_Q3-04.vtu>
+   Error in the L2 norm       :     0.000339722
+   Error in the H1 seminorm   :     0.0273542
+   Error in the enery norm    :     0.0394557
+
+DEAL::Cycle 5
+DEAL::Number of active cells:       33
+DEAL::Number of degrees of freedom: 528
+   Solving system...
+DEAL::Writing solution to <sol_Q3-05.vtu>
+   Error in the L2 norm       :     0.000217819
+   Error in the H1 seminorm   :     0.0226131
+   Error in the enery norm    :     0.0324468
+
+DEAL::Cycle 6
+DEAL::Number of active cells:       42
+DEAL::Number of degrees of freedom: 672
+   Solving system...
+DEAL::Writing solution to <sol_Q3-06.vtu>
+   Error in the L2 norm       :     9.95869e-05
+   Error in the H1 seminorm   :     0.0143484
+   Error in the enery norm    :     0.0205855
+
+DEAL::Cycle 7
+DEAL::Number of active cells:       54
+DEAL::Number of degrees of freedom: 864
+   Solving system...
+DEAL::Writing solution to <sol_Q3-07.vtu>
+   Error in the L2 norm       :     6.45972e-05
+   Error in the H1 seminorm   :     0.00914956
+   Error in the enery norm    :     0.0131191
+
+DEAL::Cycle 8
+DEAL::Number of active cells:       69
+DEAL::Number of degrees of freedom: 1104
+   Solving system...
+DEAL::Writing solution to <sol_Q3-08.vtu>
+   Error in the L2 norm       :     4.16999e-05
+   Error in the H1 seminorm   :     0.00584172
+   Error in the enery norm    :     0.0083634
+
+DEAL::Cycle 9
+DEAL::Number of active cells:       87
+DEAL::Number of degrees of freedom: 1392
+   Solving system...
+DEAL::Writing solution to <sol_Q3-09.vtu>
+   Error in the L2 norm       :     2.7367e-05
+   Error in the H1 seminorm   :     0.00372588
+   Error in the enery norm    :     0.00531894
+
+degree = 3
+| cycle | cells | dofs | L2        | H1        | Energy    | Estimator |
+| 0     | 12    | 192  | 2.783e-03 | 7.490e-02 | 1.062e-01 | 3.455e-01 |
+| 1     | 15    | 240  | 1.797e-03 | 5.685e-02 | 8.154e-02 | 3.155e-01 |
+| 2     | 18    | 288  | 1.718e-03 | 5.987e-02 | 8.509e-02 | 2.582e-01 |
+| 3     | 21    | 336  | 9.393e-04 | 4.708e-02 | 6.681e-02 | 2.174e-01 |
+| 4     | 27    | 432  | 3.397e-04 | 2.735e-02 | 3.946e-02 | 2.014e-01 |
+| 5     | 33    | 528  | 2.178e-04 | 2.261e-02 | 3.245e-02 | 1.277e-01 |
+| 6     | 42    | 672  | 9.959e-05 | 1.435e-02 | 2.059e-02 | 8.352e-02 |
+| 7     | 54    | 864  | 6.460e-05 | 9.150e-03 | 1.312e-02 | 5.609e-02 |
+| 8     | 69    | 1104 | 4.170e-05 | 5.842e-03 | 8.363e-03 | 3.807e-02 |
+| 9     | 87    | 1392 | 2.737e-05 | 3.726e-03 | 5.319e-03 | 2.590e-02 |
+@endcode
+
+The following figure provides a log-log plot of the errors versus
+the number of degrees of freedom. Let $n$ be the number of degrees of
+freedom, then $h$ is of order $1/\sqrt{n}$. Combining the theoretical
+results in the previous case,
+we see that the error in $L_2$ norm is of order $O(n^{-\frac{p+1}{2}})$
+and in $H_1$ seminorm is $O(n^{-\frac{p}{2}})$. From the figure, we see
+that the SIPG with adaptive mesh refinement produces desirable resutls
+that match theoretical ones.
+
+In addition, we observe that the error estimator decreases
+in almost the same rate as the errors in the energy norm and $H_1$ seminorm,
+and one order lower than $L_2$ error, which illustrates
+its ability to predict regions with large errors.
+
+<img width="600px" src="https://www.dealii.org/images/steps/developer/step-74-log-log-plot.png" alt="">
diff --git a/examples/step-74/doc/tooltip b/examples/step-74/doc/tooltip
new file mode 100644 (file)
index 0000000..321d3af
--- /dev/null
@@ -0,0 +1 @@
+Symmetric interior penalty Galerkin for Poisson's equation. 
\ No newline at end of file
diff --git a/examples/step-74/step-74.cc b/examples/step-74/step-74.cc
new file mode 100644 (file)
index 0000000..a458e71
--- /dev/null
@@ -0,0 +1,1053 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2015 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 at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Timo Heister and Jiaqi Zhang, Clemson University, 2020
+ */
+
+// The first few files have already been covered in previous examples and will
+// thus not be further commented on:
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/fe/fe_values.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/numerics/data_out.h>
+#include <deal.II/fe/mapping_q1.h>
+// Here the discontinuous finite elements and FEInterfaceValues are defined.
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_interface_values.h>
+
+#include <deal.II/numerics/derivative_approximation.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/base/convergence_table.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+namespace Step74
+{
+  using namespace dealii;
+
+  // @sect3{Equation data}
+  // Here we define two test cases: convergence_rate for a smooth function
+  // and l_singularity for the Functions::LSingularityFunction.
+  enum class Test_Case
+  {
+    convergence_rate,
+    l_singularity
+  };
+
+  // A smooth solution for the convergence test.
+  template <int dim>
+  class SmoothSolution : public Function<dim>
+  {
+  public:
+    SmoothSolution()
+      : Function<dim>()
+    {}
+    virtual void value_list(const std::vector<Point<dim>> &points,
+                            std::vector<double> &          values,
+                            const unsigned int component = 0) const override;
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> & point,
+             const unsigned int component = 0) const override;
+  };
+
+  template <int dim>
+  void SmoothSolution<dim>::value_list(const std::vector<Point<dim>> &points,
+                                       std::vector<double> &          values,
+                                       const unsigned int /*component*/) const
+  {
+    using numbers::PI;
+    for (unsigned int i = 0; i < values.size(); ++i)
+      values[i] =
+        std::sin(2. * PI * points[i][0]) * std::sin(2. * PI * points[i][1]);
+  }
+
+  template <int dim>
+  Tensor<1, dim>
+  SmoothSolution<dim>::gradient(const Point<dim> &point,
+                                const unsigned int /*component*/) const
+  {
+    Tensor<1, dim> return_value;
+    using numbers::PI;
+    return_value[0] =
+      2. * PI * std::cos(2. * PI * point[0]) * std::sin(2. * PI * point[1]);
+    return_value[1] =
+      2. * PI * std::sin(2. * PI * point[0]) * std::cos(2. * PI * point[1]);
+    return return_value;
+  }
+
+  // The corresponding right-hand side of the smooth function.
+  template <int dim>
+  class SmoothRightHandSide : public Function<dim>
+  {
+  public:
+    SmoothRightHandSide()
+      : Function<dim>()
+    {}
+    virtual void value_list(const std::vector<Point<dim>> &points,
+                            std::vector<double> &          values,
+                            const unsigned int /*component*/) const override
+    {
+      using numbers::PI;
+      for (unsigned int i = 0; i < values.size(); ++i)
+        values[i] = 8. * PI * PI * std::sin(2. * PI * points[i][0]) *
+                    std::sin(2. * PI * points[i][1]);
+    }
+  };
+
+  // The right-hand side corresponds to the function
+  // Functions::LSingularityFunction.
+  template <int dim>
+  class SingularRightHandSide : public Function<dim>
+  {
+  public:
+    SingularRightHandSide()
+      : Function<dim>()
+    {}
+    virtual void value_list(const std::vector<Point<dim>> &points,
+                            std::vector<double> &          values,
+                            const unsigned int /*component*/) const override
+    {
+      for (unsigned int i = 0; i < values.size(); ++i)
+        // We assume that the diffusion coefficient $\nu$ = 1.
+        values[i] = -ref.laplacian(points[i]);
+    }
+
+  private:
+    Functions::LSingularityFunction ref;
+  };
+
+  // @sect3{Auxiliary functions}
+  // The following two auxiliary functions are used to compute
+  // jump terms for $u_h$ and $\nabla u_h$ on the
+  // interface, respectively.
+  template <int dim>
+  void get_function_jump(const FEInterfaceValues<dim> &fe_iv,
+                         const Vector<double> &        solution,
+                         std::vector<double> &         jump)
+  {
+    const unsigned      n_q = fe_iv.n_quadrature_points;
+    std::vector<double> face_values[2];
+    jump.resize(n_q);
+    for (unsigned i = 0; i < 2; ++i)
+      {
+        face_values[i].resize(n_q);
+        fe_iv.get_fe_face_values(i).get_function_values(solution,
+                                                        face_values[i]);
+      }
+    for (unsigned int q = 0; q < n_q; ++q)
+      jump[q] = face_values[0][q] - face_values[1][q];
+  }
+
+  template <int dim>
+  void get_function_gradient_jump(const FEInterfaceValues<dim> &fe_iv,
+                                  const Vector<double> &        solution,
+                                  std::vector<Tensor<1, dim>> & gradient_jump)
+  {
+    const unsigned              n_q = fe_iv.n_quadrature_points;
+    std::vector<Tensor<1, dim>> face_gradients[2];
+    gradient_jump.resize(n_q);
+    for (unsigned i = 0; i < 2; ++i)
+      {
+        face_gradients[i].resize(n_q);
+        fe_iv.get_fe_face_values(i).get_function_gradients(solution,
+                                                           face_gradients[i]);
+      }
+    for (unsigned int q = 0; q < n_q; ++q)
+      gradient_jump[q] = face_gradients[0][q] - face_gradients[1][q];
+  }
+
+  // This function computes the penalty $\sigma$.
+  double compute_penalty(const unsigned int fe_degree,
+                         const double       cell_extend_left,
+                         const double       cell_extend_right)
+  {
+    const double degree = std::max(1., static_cast<double>(fe_degree));
+    return degree * (degree + 1.) * 0.5 *
+           (1. / cell_extend_left + 1. / cell_extend_right);
+  }
+
+
+  // @sect3{The CopyData}
+  // Here we define Copy objects for the MeshWorker::mesh_loop(),
+  // which is essentially the same as step-12. Note that the
+  // Scratch object is not defined here because we use
+  // MeshWorker::ScratchData<dim> instead.
+  struct CopyDataFace
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> joint_dof_indices;
+    double                               values[2];
+    unsigned int                         cell_indices[2];
+  };
+
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<CopyDataFace>            face_data;
+    double                               value;
+    unsigned int                         cell_index;
+    template <class Iterator>
+    void reinit(const Iterator &cell, unsigned int dofs_per_cell)
+    {
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
+      local_dof_indices.resize(dofs_per_cell);
+      cell->get_dof_indices(local_dof_indices);
+    }
+  };
+
+  // @sect3{The SIPGLaplace class}
+  // After this preparations, we proceed with the main class of this program
+  // called SIPGLaplace. Major differences will only come up in the
+  // implementation of the assemble functions, since use FEInterfaceValues to
+  // assemble face terms.
+  template <int dim>
+  class SIPGLaplace
+  {
+  public:
+    SIPGLaplace(const Test_Case &test_case);
+    void run();
+
+  private:
+    void setup_system();
+    void assemble_system();
+    void solve();
+    void refine_grid();
+    void output_results(const unsigned int cycle) const;
+
+    void   compute_errors();
+    void   compute_error_estimate();
+    double compute_energy_norm();
+
+    Triangulation<dim>   triangulation;
+    const MappingQ1<dim> mapping;
+
+    using ScratchData = MeshWorker::ScratchData<dim>;
+
+    FE_DGQ<dim>     fe;
+    DoFHandler<dim> dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    // Vectors to store error estimator square and energy norm square per cell.
+    Vector<double> estimated_error_square_per_cell;
+    Vector<double> energy_norm_square_per_cell;
+
+    // Print convergence rate and errors on the screen.
+    ConvergenceTable convergence_table;
+
+    // Diffusion coefficient $\nu$ is set to 1.
+    const double diffusion_coefficient = 1.;
+
+    const Test_Case test_case;
+
+    // Pointers that point to the correct classes of solution and right-hand
+    // side according to test_case.
+    std::unique_ptr<Function<dim>> exact_solution;
+    std::unique_ptr<Function<dim>> rhs_function;
+  };
+
+  // The constructor here reads the test case as an input and then determines
+  // the correct solution and right-hand side classes. The 3 in the constructor
+  // call of fe is the polynomial degree.
+  template <int dim>
+  SIPGLaplace<dim>::SIPGLaplace(const Test_Case &test_case)
+    : mapping()
+    , fe(3)
+    , dof_handler(triangulation)
+    , test_case(test_case)
+  {
+    if (test_case == Test_Case::convergence_rate)
+      {
+        exact_solution = std::make_unique<SmoothSolution<dim>>();
+        rhs_function   = std::make_unique<SmoothRightHandSide<dim>>();
+      }
+
+    else if (test_case == Test_Case::l_singularity)
+      {
+        exact_solution = std::make_unique<Functions::LSingularityFunction>();
+        rhs_function   = std::make_unique<SingularRightHandSide<dim>>();
+      }
+    else
+      AssertThrow(false, ExcNotImplemented());
+  }
+
+  template <int dim>
+  void SIPGLaplace<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+    sparsity_pattern.copy_from(dsp);
+
+    system_matrix.reinit(sparsity_pattern);
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+  }
+
+  // sect3{The assemble_system function}
+  // The assemble function here is similar to that in step-12.
+  // Different from assembling by hand, we just need to focus
+  // on assembling on each cell, each boundary face, and each
+  // interior face. The loops over cells and faces are handled
+  // automatically by MeshWorker::mesh_loop().
+  template <int dim>
+  void SIPGLaplace<dim>::assemble_system()
+  {
+    typedef decltype(dof_handler.begin_active()) Iterator;
+
+    // This function assembles the cell integrals.
+    auto cell_worker = [&](const Iterator &cell,
+                           ScratchData &   scratch_data,
+                           CopyData &      copy_data) {
+      const FEValues<dim> &fe_v          = scratch_data.reinit(cell);
+      const unsigned int   dofs_per_cell = fe_v.dofs_per_cell;
+      copy_data.reinit(cell, dofs_per_cell);
+
+      const auto &       q_points    = scratch_data.get_quadrature_points();
+      const unsigned int n_q_points  = q_points.size();
+      const std::vector<double> &JxW = scratch_data.get_JxW_values();
+
+      std::vector<double> rhs(n_q_points);
+      rhs_function->value_list(q_points, rhs);
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
+              copy_data.cell_matrix(i, j) +=
+                // nu \nabla u \nabla v
+                diffusion_coefficient * fe_v.shape_grad(i, point) *
+                fe_v.shape_grad(j, point) * JxW[point];
+
+            copy_data.cell_rhs(i) +=
+              rhs[point] * fe_v.shape_value(i, point) * JxW[point];
+          }
+    };
+
+    // This function assembles face integrals on the boundary.
+    auto boundary_worker = [&](const Iterator &    cell,
+                               const unsigned int &face_no,
+                               ScratchData &       scratch_data,
+                               CopyData &          copy_data) {
+      const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+
+      const auto &       q_points      = scratch_data.get_quadrature_points();
+      const unsigned int n_q_points    = q_points.size();
+      const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
+
+      const std::vector<double> &        JxW = scratch_data.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals =
+        scratch_data.get_normal_vectors();
+
+      std::vector<double> g(n_q_points);
+      exact_solution->value_list(q_points, g);
+
+
+      const double extent1 = cell->extent_in_direction(
+        GeometryInfo<dim>::unit_normal_direction[face_no]);
+      const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              copy_data.cell_matrix(i, j) +=
+                (
+                  // - nu (\nabla u . n) v
+                  -diffusion_coefficient *
+                    (fe_fv.shape_grad(j, point) * normals[point]) *
+                    fe_fv.shape_value(i, point)
+
+                  // - nu u (\nabla v . n)
+                  - diffusion_coefficient * fe_fv.shape_value(j, point) *
+                      (fe_fv.shape_grad(i, point) * normals[point])
+
+                  // + nu * penalty u v
+                  +
+                  diffusion_coefficient * penalty *
+                    fe_fv.shape_value(j, point) * fe_fv.shape_value(i, point)) *
+                JxW[point];
+
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            copy_data.cell_rhs(i) +=
+              (
+                // -nu g (\nabla v . n)
+                -diffusion_coefficient * g[point] *
+                  (fe_fv.shape_grad(i, point) * normals[point])
+
+                // +nu penalty g v
+                + diffusion_coefficient * penalty * g[point] *
+                    fe_fv.shape_value(i, point)) *
+              JxW[point];
+        }
+    };
+
+    // This function assembles face integrals on interior faces.
+    // To reinitialize FEInterfaceValues, we need to pass cells,
+    // face and subface indices (for adaptive refinement)
+    // to the reinit() function of FEInterfaceValues.
+    auto face_worker = [&](const Iterator &    cell,
+                           const unsigned int &f,
+                           const unsigned int &sf,
+                           const Iterator &    ncell,
+                           const unsigned int &nf,
+                           const unsigned int &nsf,
+                           ScratchData &       scratch_data,
+                           CopyData &          copy_data) {
+      const FEInterfaceValues<dim> &fe_iv =
+        scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+
+      const auto &       q_points   = fe_iv.get_quadrature_points();
+      const unsigned int n_q_points = q_points.size();
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &     copy_data_face = copy_data.face_data.back();
+      const unsigned int n_dofs_face    = fe_iv.n_current_interface_dofs();
+      copy_data_face.joint_dof_indices  = fe_iv.get_interface_dof_indices();
+      copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
+
+      const std::vector<double> &        JxW     = fe_iv.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+      const double extent1 =
+        cell->extent_in_direction(GeometryInfo<dim>::unit_normal_direction[f]);
+      const double extent2 = ncell->extent_in_direction(
+        GeometryInfo<dim>::unit_normal_direction[nf]);
+      const double penalty = compute_penalty(fe.get_degree(), extent1, extent2);
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          for (unsigned int i = 0; i < n_dofs_face; ++i)
+            for (unsigned int j = 0; j < n_dofs_face; ++j)
+              copy_data_face.cell_matrix(i, j) +=
+                (
+                  // - nu {\nabla u}.n [v] (consistency)
+                  -diffusion_coefficient *
+                    (fe_iv.average_gradient(j, point) * normals[point]) *
+                    fe_iv.jump(i, point)
+
+                  // - nu [u] {\nabla v}.n  (symmetry) // NIPG: use +
+                  - diffusion_coefficient * fe_iv.jump(j, point) *
+                      (fe_iv.average_gradient(i, point) * normals[point])
+
+                  // nu sigma [u] [v] (penalty)
+                  + diffusion_coefficient * penalty * fe_iv.jump(j, point) *
+                      fe_iv.jump(i, point)
+
+                    ) *
+                JxW[point];
+        }
+    };
+
+    // The following lambda function will copy data to
+    // the global matrix and right-hand side.
+    // Though there are no hanging node constraints in DG discretization,
+    // we define an empty AffineConstraints oject that
+    // allows us to use copy_local_to_global functionality.
+    AffineConstraints<double> constraints;
+    constraints.close();
+    auto copier = [&](const CopyData &c) {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             system_rhs);
+
+      // Copy data from interior face assembly to the global matrix.
+      for (auto &cdf : c.face_data)
+        {
+          const unsigned int joint_dofs_per_face = cdf.joint_dof_indices.size();
+          for (unsigned int i = 0; i < joint_dofs_per_face; ++i)
+            for (unsigned int k = 0; k < joint_dofs_per_face; ++k)
+              system_matrix.add(cdf.joint_dof_indices[i],
+                                cdf.joint_dof_indices[k],
+                                cdf.cell_matrix(i, k));
+        }
+    };
+
+    // Here we define ScratchData and CopyData objects,
+    // and pass them together with the lambda functions
+    // above to MeshWorker::mesh_loop. In addition, we
+    // need to specify that we want to assemble interior faces once.
+    const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+    QGauss<dim>        quadrature(n_gauss_points);
+    QGauss<dim - 1>    face_quadrature(n_gauss_points);
+
+    UpdateFlags cell_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+    UpdateFlags face_flags = update_values | update_gradients |
+                             update_quadrature_points | update_normal_vectors |
+                             update_JxW_values;
+
+    ScratchData scratch_data(
+      mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
+    CopyData cd;
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          cd,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_boundary_faces |
+                            MeshWorker::assemble_own_interior_faces_once,
+                          boundary_worker,
+                          face_worker);
+  }
+
+  template <int dim>
+  void SIPGLaplace<dim>::solve()
+  {
+    std::cout << "   Solving system..." << std::endl;
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+    A_direct.vmult(solution, system_rhs);
+  }
+
+  template <int dim>
+  void SIPGLaplace<dim>::output_results(const unsigned int cycle) const
+  {
+    std::string filename = "sol_Q" +
+                           Utilities::int_to_string(fe.get_degree(), 1) + "-" +
+                           Utilities::int_to_string(cycle, 2) + ".vtu";
+    std::cout << "Writing solution to <" << filename << ">" << std::endl;
+    std::ofstream output(filename);
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
+    data_out.build_patches();
+    data_out.write_vtu(output);
+  }
+
+  // The assembly of the error estimator here is quite similar to
+  // that of the global matrix and right-had side.
+  template <int dim>
+  void SIPGLaplace<dim>::compute_error_estimate()
+  {
+    typedef decltype(dof_handler.begin_active()) Iterator;
+    estimated_error_square_per_cell.reinit(triangulation.n_active_cells());
+
+    // Assemble cell residual $h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2$.
+    auto cell_worker = [&](const Iterator &cell,
+                           ScratchData &   scratch_data,
+                           CopyData &      copy_data) {
+      const FEValues<dim> &fe_v = scratch_data.reinit(cell);
+
+      copy_data.cell_index = cell->active_cell_index();
+
+      const auto &               q_points   = fe_v.get_quadrature_points();
+      const unsigned int         n_q_points = q_points.size();
+      const std::vector<double> &JxW        = fe_v.get_JxW_values();
+
+      std::vector<Tensor<2, dim>> hessians(n_q_points);
+      fe_v.get_function_hessians(solution, hessians);
+
+      std::vector<double> rhs(n_q_points);
+      rhs_function->value_list(q_points, rhs);
+
+      const double hk                   = cell->diameter();
+      double       residual_norm_square = 0;
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          const double residual =
+            rhs[point] + diffusion_coefficient * trace(hessians[point]);
+          residual_norm_square += residual * residual * JxW[point];
+        }
+      copy_data.value = hk * hk * residual_norm_square;
+    };
+
+    // Assemble boundary terms $\sum_{f\in \partial K \cap \partial \Omega}
+    // \sigma \left\| [  u_h-g_D ]  \right\|_f^2  $.
+    auto boundary_worker = [&](const Iterator &    cell,
+                               const unsigned int &face_no,
+                               ScratchData &       scratch_data,
+                               CopyData &          copy_data) {
+      const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+
+      const auto &   q_points   = fe_fv.get_quadrature_points();
+      const unsigned n_q_points = q_points.size();
+
+      const std::vector<double> &JxW = fe_fv.get_JxW_values();
+
+      std::vector<double> g(n_q_points);
+      exact_solution->value_list(q_points, g);
+
+      std::vector<double> sol_u(n_q_points);
+      fe_fv.get_function_values(solution, sol_u);
+
+      const double extent1 = cell->extent_in_direction(
+        GeometryInfo<dim>::unit_normal_direction[face_no]);
+      const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+
+      double difference_norm_square = 0.;
+      for (unsigned int point = 0; point < q_points.size(); ++point)
+        {
+          const double diff = (g[point] - sol_u[point]);
+          difference_norm_square += diff * diff * JxW[point];
+        }
+      copy_data.value += penalty * difference_norm_square;
+    };
+
+    // Assemble interior face terms $\sum_{f\in \partial K}\lbrace \sigma
+    // \left\| [u_h]  \right\|_f^2   +  h_f \left\|  [\nu \nabla u_h \cdot
+    // \mathbf n ] \right\|_f^2 \rbrace$.
+    auto face_worker = [&](const Iterator &    cell,
+                           const unsigned int &f,
+                           const unsigned int &sf,
+                           const Iterator &    ncell,
+                           const unsigned int &nf,
+                           const unsigned int &nsf,
+                           ScratchData &       scratch_data,
+                           CopyData &          copy_data) {
+      const FEInterfaceValues<dim> &fe_iv =
+        scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+      copy_data_face.cell_indices[0] = cell->active_cell_index();
+      copy_data_face.cell_indices[1] = ncell->active_cell_index();
+
+      const std::vector<double> &        JxW     = fe_iv.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+      const auto &       q_points   = fe_iv.get_quadrature_points();
+      const unsigned int n_q_points = q_points.size();
+
+      std::vector<double> jump(n_q_points);
+      get_function_jump(fe_iv, solution, jump);
+
+      std::vector<Tensor<1, dim>> grad_jump(n_q_points);
+      get_function_gradient_jump(fe_iv, solution, grad_jump);
+
+      const double h = cell->face(f)->diameter();
+
+      const double extent1 =
+        cell->extent_in_direction(GeometryInfo<dim>::unit_normal_direction[f]);
+      const double extent2 = ncell->extent_in_direction(
+        GeometryInfo<dim>::unit_normal_direction[nf]);
+      const double penalty = compute_penalty(fe.get_degree(), extent1, extent2);
+
+      double flux_jump_square = 0;
+      double u_jump_square    = 0;
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          u_jump_square += jump[point] * jump[point] * JxW[point];
+          const double flux_jump = grad_jump[point] * normals[point];
+          flux_jump_square +=
+            diffusion_coefficient * flux_jump * flux_jump * JxW[point];
+        }
+      copy_data_face.values[0] =
+        0.5 * h * (flux_jump_square + penalty * u_jump_square);
+      copy_data_face.values[1] = copy_data_face.values[0];
+    };
+
+    auto copier = [&](const CopyData &copy_data) {
+      if (copy_data.cell_index != numbers::invalid_unsigned_int)
+        estimated_error_square_per_cell[copy_data.cell_index] +=
+          copy_data.value;
+      for (auto &cdf : copy_data.face_data)
+        for (unsigned int j = 0; j < 2; ++j)
+          estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
+    };
+
+    const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+    QGauss<dim>        quadrature(n_gauss_points);
+    QGauss<dim - 1>    face_quadrature(n_gauss_points);
+
+    UpdateFlags cell_flags =
+      update_hessians | update_quadrature_points | update_JxW_values;
+    UpdateFlags face_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values |
+                             update_normal_vectors;
+
+    ScratchData scratch_data(
+      mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
+
+    CopyData cd;
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          cd,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_own_interior_faces_once |
+                            MeshWorker::assemble_boundary_faces,
+                          boundary_worker,
+                          face_worker);
+  }
+
+  // Here we compute the error in the energy norm, which
+  // is similar to the assembling of the error estimator.
+  template <int dim>
+  double SIPGLaplace<dim>::compute_energy_norm()
+  {
+    typedef decltype(dof_handler.begin_active()) Iterator;
+    energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
+
+    auto cell_worker = [&](const Iterator &cell,
+                           ScratchData &   scratch_data,
+                           CopyData &      copy_data) {
+      const FEValues<dim> &fe_v = scratch_data.reinit(cell);
+
+      copy_data.cell_index = cell->active_cell_index();
+
+      const auto &               q_points   = fe_v.get_quadrature_points();
+      const unsigned int         n_q_points = q_points.size();
+      const std::vector<double> &JxW        = fe_v.get_JxW_values();
+
+      std::vector<Tensor<1, dim>> grad_u(n_q_points);
+      fe_v.get_function_gradients(solution, grad_u);
+
+      std::vector<Tensor<1, dim>> grad_exact(n_q_points);
+      exact_solution->gradient_list(q_points, grad_exact);
+
+      double norm_square = 0;
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          norm_square +=
+            (grad_u[point] - grad_exact[point]).norm_square() * JxW[point];
+        }
+      copy_data.value = norm_square;
+    };
+
+    auto boundary_worker = [&](const Iterator &    cell,
+                               const unsigned int &face_no,
+                               ScratchData &       scratch_data,
+                               CopyData &          copy_data) {
+      const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+
+      const auto &   q_points   = fe_fv.get_quadrature_points();
+      const unsigned n_q_points = q_points.size();
+
+      const std::vector<double> &JxW = fe_fv.get_JxW_values();
+
+      std::vector<double> g(n_q_points);
+      exact_solution->value_list(q_points, g);
+
+      std::vector<double> sol_u(n_q_points);
+      fe_fv.get_function_values(solution, sol_u);
+
+      const double extent1 = cell->extent_in_direction(
+        GeometryInfo<dim>::unit_normal_direction[face_no]);
+      const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+
+      double difference_norm_square = 0.;
+      for (unsigned int point = 0; point < q_points.size(); ++point)
+        {
+          const double diff = (g[point] - sol_u[point]);
+          difference_norm_square += diff * diff * JxW[point];
+        }
+      copy_data.value += penalty * difference_norm_square;
+    };
+
+    auto face_worker = [&](const Iterator &    cell,
+                           const unsigned int &f,
+                           const unsigned int &sf,
+                           const Iterator &    ncell,
+                           const unsigned int &nf,
+                           const unsigned int &nsf,
+                           ScratchData &       scratch_data,
+                           CopyData &          copy_data) {
+      const FEInterfaceValues<dim> &fe_iv =
+        scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+      copy_data_face.cell_indices[0] = cell->active_cell_index();
+      copy_data_face.cell_indices[1] = ncell->active_cell_index();
+
+      const std::vector<double> &JxW = fe_iv.get_JxW_values();
+
+      const auto &       q_points   = fe_iv.get_quadrature_points();
+      const unsigned int n_q_points = q_points.size();
+
+      std::vector<double> jump(n_q_points);
+      get_function_jump(fe_iv, solution, jump);
+
+      const double extent1 =
+        cell->extent_in_direction(GeometryInfo<dim>::unit_normal_direction[f]);
+      const double extent2 = ncell->extent_in_direction(
+        GeometryInfo<dim>::unit_normal_direction[nf]);
+      const double penalty = compute_penalty(fe.get_degree(), extent1, extent2);
+
+      double u_jump_square = 0;
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          u_jump_square += jump[point] * jump[point] * JxW[point];
+        }
+      copy_data_face.values[0] = 0.5 * penalty * u_jump_square;
+      copy_data_face.values[1] = copy_data_face.values[0];
+    };
+
+    auto copier = [&](const CopyData &copy_data) {
+      if (copy_data.cell_index != numbers::invalid_unsigned_int)
+        energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
+      for (auto &cdf : copy_data.face_data)
+        for (unsigned int j = 0; j < 2; ++j)
+          energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
+    };
+
+    const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+    QGauss<dim>        quadrature(n_gauss_points);
+    QGauss<dim - 1>    face_quadrature(n_gauss_points);
+
+    UpdateFlags cell_flags =
+      update_gradients | update_quadrature_points | update_JxW_values;
+    UpdateFlags face_flags =
+      update_values | update_quadrature_points | update_JxW_values;
+
+    ScratchData scratch_data(
+      mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
+
+    CopyData cd;
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          cd,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_own_interior_faces_once |
+                            MeshWorker::assemble_boundary_faces,
+                          boundary_worker,
+                          face_worker);
+    const double energy_error =
+      std::sqrt(energy_norm_square_per_cell.l1_norm());
+    return energy_error;
+  }
+
+  template <int dim>
+  void SIPGLaplace<dim>::refine_grid()
+  {
+    const double refinement_fraction = 0.1;
+
+    GridRefinement::refine_and_coarsen_fixed_number(
+      triangulation, estimated_error_square_per_cell, refinement_fraction, 0.);
+
+    triangulation.execute_coarsening_and_refinement();
+  }
+
+  // We compute three errors in $L_2$ norm, $H_1$ seminorm, and the energy norm,
+  // respectively.
+  template <int dim>
+  void SIPGLaplace<dim>::compute_errors()
+  {
+    double L2_error, H1_error;
+
+    {
+      Vector<float> difference_per_cell(triangulation.n_active_cells());
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        *(exact_solution.get()),
+                                        difference_per_cell,
+                                        QGauss<dim>(fe.degree + 2),
+                                        VectorTools::L2_norm);
+
+      L2_error = VectorTools::compute_global_error(triangulation,
+                                                   difference_per_cell,
+                                                   VectorTools::L2_norm);
+    }
+
+    {
+      Vector<float> difference_per_cell(triangulation.n_active_cells());
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        *(exact_solution.get()),
+                                        difference_per_cell,
+                                        QGauss<dim>(fe.degree + 2),
+                                        VectorTools::H1_seminorm);
+
+      H1_error = VectorTools::compute_global_error(triangulation,
+                                                   difference_per_cell,
+                                                   VectorTools::H1_seminorm);
+    }
+
+    convergence_table.add_value("L2", L2_error);
+    convergence_table.add_value("H1", H1_error);
+    const double energy_error = compute_energy_norm();
+    convergence_table.add_value("Energy", energy_error);
+
+    std::cout << "   Error in the L2 norm       :     " << L2_error << std::endl
+              << "   Error in the H1 seminorm   :     " << H1_error << std::endl
+              << "   Error in the energy norm    :     " << energy_error
+              << std::endl;
+  }
+
+  template <int dim>
+  void SIPGLaplace<dim>::run()
+  {
+    unsigned int max_cycle = test_case == Test_Case::convergence_rate ? 6 : 10;
+    for (unsigned int cycle = 0; cycle < max_cycle; ++cycle)
+      {
+        std::cout << "Cycle " << cycle << std::endl;
+
+        switch (test_case)
+          {
+            case Test_Case::convergence_rate:
+              {
+                if (cycle == 0)
+                  {
+                    GridGenerator::hyper_cube(triangulation);
+
+                    triangulation.refine_global(2);
+                  }
+                else
+                  {
+                    triangulation.refine_global(1);
+                  }
+                break;
+              }
+            case Test_Case::l_singularity:
+              {
+                if (cycle == 0)
+                  {
+                    GridGenerator::hyper_L(triangulation);
+                    triangulation.refine_global(2);
+                  }
+                else
+                  {
+                    refine_grid();
+                  }
+              }
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+        std::cout << "Number of active cells:       "
+                  << triangulation.n_active_cells() << std::endl;
+        setup_system();
+
+        std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+                  << std::endl;
+
+        assemble_system();
+        solve();
+        output_results(cycle);
+        {
+          convergence_table.add_value("cycle", cycle);
+          convergence_table.add_value("cells", triangulation.n_active_cells());
+          convergence_table.add_value("dofs", dof_handler.n_dofs());
+        }
+        compute_errors();
+
+        if (test_case == Test_Case::l_singularity)
+          {
+            compute_error_estimate();
+            convergence_table.add_value(
+              "Estimator",
+              std::sqrt(estimated_error_square_per_cell.l1_norm()));
+          }
+        std::cout << std::endl;
+      }
+    {
+      convergence_table.set_precision("L2", 3);
+      convergence_table.set_precision("H1", 3);
+      convergence_table.set_precision("Energy", 3);
+
+      convergence_table.set_scientific("L2", true);
+      convergence_table.set_scientific("H1", true);
+      convergence_table.set_scientific("Energy", true);
+
+      if (test_case == Test_Case::l_singularity)
+        {
+          convergence_table.set_precision("Estimator", 3);
+          convergence_table.set_scientific("Estimator", true);
+        }
+      if (test_case == Test_Case::convergence_rate)
+        {
+          convergence_table.evaluate_convergence_rates(
+            "L2", ConvergenceTable::reduction_rate_log2);
+          convergence_table.evaluate_convergence_rates(
+            "H1", ConvergenceTable::reduction_rate_log2);
+        }
+
+      std::cout << "degree = " << fe.get_degree() << std::endl;
+      convergence_table.write_text(
+        std::cout, TableHandler::TextOutputFormat::org_mode_table);
+    }
+  }
+} // namespace Step74
+
+
+// The following <code>main</code> function is similar to previous examples as
+// well, and need not be commented on.
+int main()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Step74;
+      Test_Case      test_case = Test_Case::l_singularity;
+      SIPGLaplace<2> problem(test_case);
+      problem.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+
+  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.