]> https://gitweb.dealii.org/ - dealii.git/commitdiff
update
authorJiaqi Zhang <zjiaqi@vt.edu>
Mon, 9 Nov 2020 15:24:47 +0000 (10:24 -0500)
committerJiaqi Zhang <zjiaqi@vt.edu>
Mon, 30 Nov 2020 17:06:41 +0000 (12:06 -0500)
doc/doxygen/references.bib
examples/step-74/doc/intro.dox
examples/step-74/doc/results.dox
examples/step-74/doc/tooltip
examples/step-74/step-74.cc

index 9421a8205c0da6ff49761c3253a5ed2d9dc6be9b..d920803a5e73a368950b06a2cbc0b14af798a4a6 100644 (file)
@@ -830,7 +830,7 @@ year = {2008},
 
 @book{di2011mathematical,
   title={Mathematical aspects of discontinuous Galerkin methods},
-  author={Di Pietro and Daniele Antonio and Ern, Alexandre},
+  author={Di Pietro, Daniele Antonio and Ern, Alexandre},
   volume={69},
   year={2011},
   publisher={Springer Science \& Business Media}
index 54a00110b4ed46fd639b0cc63fa63148077ed1ba..90965c30cffeda6725bcbfd72dbed4849b1a8101 100644 (file)
@@ -28,7 +28,7 @@ This tutorial includes the following topics.
 <h3>The equation</h3>
 In this example, we consider Poisson's equation
 @f[
-- \nu  \Delta u = f  \qquad   \mbox{in } \Omega,
+- \nabla \cdot \left( \nu  \nabla u\right) = f  \qquad   \mbox{in } \Omega,
 @f]
 subject to the boundary condition
 @f[
@@ -38,8 +38,8 @@ 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$
+We denote the mesh by ${\mathbb T}_h$, and $K\in{\mathbb T}_h$ is a mesh cell.
+The sets of interior and boundary faces are denoted by ${\mathbb F}^i_h$ and ${\mathbb 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
@@ -55,24 +55,24 @@ $\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_{K\in {\mathbb T}_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^i} \left\{
+    \left< \jump{v_h}, \nu\average{ \nabla u_h} \cdot  \mathbf n \right>_F
+   +\left<\average{ \nabla v_h }\cdot \mathbf n,\nu\jump{u_h}\right>_F
+   -\left<\jump{v_h},\nu \sigma \jump{u_h} \right>_F
+  \right\}
   \\
-  - \sum_{F \in F_h^b} \biggl\{
-    \bigl<v_h, \nu  \nabla u_h\cdot \mathbf n \bigr>_F
-  + \bigl< \nabla v_h \cdot \mathbf n , \nu u_h\bigr>_F
-  - \bigl< v_h,\nu \sigma u_h\bigr>_F
-  \biggr\}
+  - \sum_{F \in F_h^b} \left\{
+    \left<v_h, \nu  \nabla u_h\cdot \mathbf n \right>_F
+  + \left< \nabla v_h \cdot \mathbf n , \nu u_h\right>_F
+  - \left< v_h,\nu \sigma u_h\right>_F
+  \right\}
   \\
   = (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\}.
+  - \sum_{F \in F_h^b} \left\{
+    \left< \nabla v_h \cdot \mathbf n, \nu g_D\right>_F - \left<v_h,\nu \sigma g_D\right>_F
+  \right\}.
 @f}
 
 <h3>The penalty parameter</h3>
@@ -85,9 +85,9 @@ One can just pick a large constant, while other options could be the multiples o
 we follow step-39 and use $\gamma = p(p+1)$.
 
 <h3>A posteriori error estimator</h3>
-In this example, we use the error estimator by Karakashian and Pascal @cite karakashian2003posteriori
+In this example, we use the error estimator by Karakashian and Pascal @cite karakashian2003posteriori with a slight modification
 @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}
+\eta^2 = \sum_{K \in {\mathbb T}_h} \eta^2_{K} +  \sum_{f_i \in {\mathbb F}^i_h}  \eta^2_{f_i} + \sum_{f_b \in F^i_b}\eta^2_{f_b}
 @f]
 where
 @f[
@@ -99,7 +99,7 @@ where
 @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$).
+Here we use $\sigma = \gamma/h_f$ instead of $\gamma^2/h_f$ for the jump terms of $u_h$ (the first term in $\eta^2_{f_i}$ and $\eta_{f_b}^2$).
 
 In each cell $K$, we compute
 @f[
@@ -117,3 +117,19 @@ Then the error estimate square per cell is
 @f]
 Note that we compute $\eta_{local}^2$ instead of $\eta_{local}$ to simplify the implementation.
 The error estimate square per cell is stored in a global vector, whose $L_1$ norm is equal to $\eta^2$.
+
+<h3>The test case</h3>
+In the first test problem, we run a convergence test using a smooth manufactured solution with $\nu =1$ in 2D
+@f[
+u=\sin(2\pi x)\sin(2\pi y), (x,y)\in (0,1)\times (0,1),
+@f]
+correspondingly,
+@f[
+u=0, \qquad \mbox{on } \partial \Omega.
+@f]
+and $f= 8\pi^2 u$. We compute errors against the manufactured solution and evaluate the convergence rate.
+
+In the second test, we choose Functions::LSingularityFunction on a L-shaped domain (GridGenerator::hyper_L) in 2D.
+The solution is given in the polar coordinates by $u(r,\phi) = r^{\frac{2}{3}}\sin \left(\frac{2}{3}\phi \right)$,
+which has a singularity at the origin. An error estimator is constructed to detect the region with large errors,
+according to which the mesh is refined adaptively.
index 9896f504e4e411375bf16a8f8021e492e874faa0..b14ec3acfd96a021164bde9c515c7bef80280d66 100644 (file)
 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:
+In the first test case, when you run the program, the screen output should look like the following:
 @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 |
+Cycle 0
+  Number of active cells       : 16
+  Number of degrees of freedom : 256
+  Error in the L2 norm         : 0.00193285
+  Error in the H1 seminorm     : 0.106087
+  Error in the energy norm     : 0.150625
+
+Cycle 1
+  Number of active cells       : 64
+  Number of degrees of freedom : 1024
+  Error in the L2 norm         : 9.60497e-05
+  Error in the H1 seminorm     : 0.0089954
+  Error in the energy norm     : 0.0113265
+
+Cycle 2
+.
+.
+.
 @endcode
+Convergence rate for the smooth case with polynomial degree 3:
+<table align="center" class="doxtable">
+  <tr>
+    <th>cycle</th>
+    <th>n_cellss</th>
+    <th>n_dofs</th>
+    <th>L2 </th>
+    <th>rate</th>
+    <th>H1</th>
+    <th>rate</th>
+    <th>Energy</th>
+  </tr>
+  <tr>
+    <td align="center">0</td>
+    <td align="right">16</td>
+    <td align="right">256</td>
+    <td align="center">1.933e-03</td>
+    <td>&nbsp;</td>
+    <td align="center">1.061e-01</td>
+    <td>&nbsp;</td>
+    <td align="center">1.506e-01</td>
+  </tr>
+  <tr>
+    <td align="center">1</td>
+    <td align="right">64</td>
+    <td align="right">1024</td>
+    <td align="center">9.605e-05</td>
+    <td align="center">4.33</td>
+    <td align="center">8.995e-03</td>
+    <td align="center">3.56</td>
+    <td align="center">1.133e-02</td>
+  </tr>
+  <tr>
+    <td align="center">2</td>
+    <td align="right">256</td>
+    <td align="right">4096</td>
+    <td align="center">5.606e-06</td>
+    <td align="center">4.10</td>
+    <td align="center">9.018e-04</td>
+    <td align="center">3.32</td>
+    <td align="center">9.736e-04</td>
+  </tr>
+  <tr>
+    <td align="center">3</td>
+    <td align="right">1024</td>
+    <td align="right">16384</td>
+    <td align="center">3.484e-07</td>
+    <td align="center">4.01</td>
+    <td align="center">1.071e-04</td>
+    <td align="center">3.07</td>
+    <td align="center">1.088e-04</td>
+  </tr>
+  <tr>
+    <td align="center">4</td>
+    <td align="right">4096</td>
+    <td align="right">65536</td>
+    <td align="center">2.179e-08</td>
+    <td align="center">4.00</td>
+    <td align="center">1.327e-05</td>
+    <td align="center">3.01</td>
+    <td align="center">1.331e-05</td>
+  </tr>
+  <tr>
+    <td align="center">5</td>
+    <td align="right">16384</td>
+    <td align="right">262144</td>
+    <td align="center">1.363e-09</td>
+    <td align="center">4.00</td>
+    <td align="center">1.656e-06</td>
+    <td align="center">3.00</td>
+    <td align="center">1.657e-06</td>
+  </tr>
+</table>
+
 Theoretically, for polynomial degree $p$, the order of convergence in $L_2$
 norm and $H_1$ seminorm should be $p+1$ and $p$, respectively. Our numerical
 results are in good agreement with theory.
 
-Errors of the singular case with polynomial degree 3.
+In the second test case, when you run the program, the screen output should look like the following:
 @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 |
+Cycle 0
+  Number of active cells       : 192
+  Number of degrees of freedom : 3072
+  Error in the L2 norm         : 0.000323585
+  Error in the H1 seminorm     : 0.0296202
+  Error in the energy norm     : 0.0420478
+  Estimated error              : 0.136067
+
+Cycle 1
+  Number of active cells       : 249
+  Number of degrees of freedom : 3984
+  Error in the L2 norm         : 0.000114739
+  Error in the H1 seminorm     : 0.0186571
+  Error in the energy norm     : 0.0264879
+  Estimated error              : 0.0857186
+
+Cycle 2
+.
+.
+.
 @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}$ in 2d. Combining the theoretical
+freedom, then $h$ is of order $1/\sqrt{n}$ in 2D. 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
@@ -193,3 +141,9 @@ 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="">
+
+While this tutorial is focused on the implementation, the step-59 tutorial program achieves an efficient
+large-scale solver in terms of computing time with matrix-free solution techniques.
+Note that the step-59 tutorial does not work with meshes containing hanging nodes at this moment,
+because the multigrid interface matrices are not as easily determined,
+but that is merely the lack of some interfaces in deal.II, nothing fundamental.
\ No newline at end of file
index 321d3af5f4814fc0988ffdb04b01d9b2ef92e444..fbb04bcfc1236dcab580cf2e3362f2e574e70ce6 100644 (file)
@@ -1 +1 @@
-Symmetric interior penalty Galerkin for Poisson's equation. 
\ No newline at end of file
+Symmetric interior penalty Galerkin for Poisson's equation. 
index 04080a4bc95c91900ff9808ae6e043665a82e0e0..9d42c971aedd30c2b722276198cd724cb3eab456 100644 (file)
@@ -1,6 +1,6 @@
 /* ---------------------------------------------------------------------
  *
- * Copyright (C) 2009 - 2015 by the deal.II authors
+ * Copyright (C) 2020 by the deal.II authors
  *
  * This file is part of the deal.II library.
  *
@@ -40,7 +40,6 @@
 #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>
@@ -115,15 +114,21 @@ namespace Step74
     {}
     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]);
-    }
+                            const unsigned int /*component*/) const override;
   };
 
+  template <int dim>
+  void
+  SmoothRightHandSide<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] = 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>
@@ -135,17 +140,23 @@ namespace Step74
     {}
     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]);
-    }
+                            const unsigned int /*component*/) const override;
 
   private:
     Functions::LSingularityFunction ref;
   };
 
+  template <int dim>
+  void
+  SingularRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
+                                         std::vector<double> &          values,
+                                         const unsigned int /*component*/) const
+  {
+    for (unsigned int i = 0; i < values.size(); ++i)
+      // We assume that the diffusion coefficient $\nu$ = 1.
+      values[i] = -ref.laplacian(points[i]);
+  }
+
   // @sect3{Auxiliary functions}
   // The following two auxiliary functions are used to compute
   // jump terms for $u_h$ and $\nabla u_h$ on the
@@ -155,8 +166,8 @@ namespace Step74
                          const Vector<double> &        solution,
                          std::vector<double> &         jump)
   {
-    const unsigned      n_q = fe_iv.n_quadrature_points;
-    std::vector<double> face_values[2];
+    const unsigned                     n_q = fe_iv.n_quadrature_points;
+    std::array<std::vector<double>, 2> face_values;
     jump.resize(n_q);
     for (unsigned i = 0; i < 2; ++i)
       {
@@ -206,8 +217,8 @@ namespace Step74
   {
     FullMatrix<double>                   cell_matrix;
     std::vector<types::global_dof_index> joint_dof_indices;
-    double                               values[2];
-    unsigned int                         cell_indices[2];
+    std::array<double, 2>                values;
+    std::array<unsigned int, 2>          cell_indices;
   };
 
   struct CopyData
@@ -252,6 +263,9 @@ namespace Step74
     double compute_energy_norm();
 
     Triangulation<dim>   triangulation;
+    const unsigned       degree;
+    QGauss<dim>          quadrature;
+    QGauss<dim - 1>      face_quadrature;
     const MappingQ1<dim> mapping;
 
     using ScratchData = MeshWorker::ScratchData<dim>;
@@ -283,12 +297,14 @@ namespace Step74
   };
 
   // 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.
+  // the correct solution and right-hand side classes.
   template <int dim>
   SIPGLaplace<dim>::SIPGLaplace(const Test_Case &test_case)
-    : mapping()
-    , fe(3)
+    : degree(3)
+    , quadrature(degree + 1)
+    , face_quadrature(degree + 1)
+    , mapping()
+    , fe(degree)
     , dof_handler(triangulation)
     , test_case(test_case)
   {
@@ -329,44 +345,41 @@ namespace Step74
   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) +=
-                diffusion_coefficient *     // nu
-                fe_v.shape_grad(i, point) * // grad v_h
-                fe_v.shape_grad(j, point) * // grad u_h
-                JxW[point];                 // dx
-
-            copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
-                                     rhs[point] *                 // f
-                                     JxW[point];                  // dx
-          }
-    };
+    const auto cell_worker =
+      [&](const auto &cell, auto &scratch_data, auto &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) +=
+                  diffusion_coefficient *     // nu
+                  fe_v.shape_grad(i, point) * // grad v_h
+                  fe_v.shape_grad(j, point) * // grad u_h
+                  JxW[point];                 // dx
+
+              copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
+                                       rhs[point] *                 // f
+                                       JxW[point];                  // dx
+            }
+      };
 
     // 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 auto boundary_worker = [&](const auto &        cell,
+                                     const unsigned int &face_no,
+                                     auto &              scratch_data,
+                                     auto &              copy_data) {
       const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
 
       const auto &       q_points      = scratch_data.get_quadrature_points();
@@ -382,7 +395,7 @@ namespace Step74
 
       const double extent1 = cell->extent_in_direction(
         GeometryInfo<dim>::unit_normal_direction[face_no]);
-      const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+      const double penalty = compute_penalty(degree, extent1, extent1);
 
       for (unsigned int point = 0; point < n_q_points; ++point)
         {
@@ -426,14 +439,14 @@ namespace Step74
     // 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 auto face_worker = [&](const auto &        cell,
+                                 const unsigned int &f,
+                                 const unsigned int &sf,
+                                 const auto &        ncell,
+                                 const unsigned int &nf,
+                                 const unsigned int &nsf,
+                                 auto &              scratch_data,
+                                 auto &              copy_data) {
       const FEInterfaceValues<dim> &fe_iv =
         scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
 
@@ -453,7 +466,7 @@ namespace Step74
         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);
+      const double penalty = compute_penalty(degree, extent1, extent2);
 
       for (unsigned int point = 0; point < n_q_points; ++point)
         {
@@ -483,10 +496,10 @@ namespace Step74
     // 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.
+    // allows us to use distribute_local_to_global functionality.
     AffineConstraints<double> constraints;
     constraints.close();
-    auto copier = [&](const CopyData &c) {
+    const auto copier = [&](const auto &c) {
       constraints.distribute_local_to_global(c.cell_matrix,
                                              c.cell_rhs,
                                              c.local_dof_indices,
@@ -496,12 +509,9 @@ namespace Step74
       // 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));
+          constraints.distribute_local_to_global(cdf.cell_matrix,
+                                                 cdf.joint_dof_indices,
+                                                 system_matrix);
         }
     };
 
@@ -509,9 +519,6 @@ namespace Step74
     // 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;
@@ -538,7 +545,6 @@ namespace Step74
   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);
@@ -547,10 +553,8 @@ namespace Step74
   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) + "-" +
+    std::string filename = "sol_Q" + Utilities::int_to_string(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;
@@ -565,45 +569,43 @@ namespace Step74
   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);
+    const auto cell_worker =
+      [&](const auto &cell, auto &scratch_data, auto &copy_data) {
+        const FEValues<dim> &fe_v = scratch_data.reinit(cell);
 
-      copy_data.cell_index = cell->active_cell_index();
+        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();
+        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<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);
+        std::vector<double> rhs(n_q_points);
+        rhs_function->value_list(q_points, rhs);
 
-      const double hk                   = cell->diameter();
-      double       residual_norm_square = 0;
+        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;
-    };
+        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 auto boundary_worker = [&](const auto &        cell,
+                                     const unsigned int &face_no,
+                                     auto &              scratch_data,
+                                     auto &              copy_data) {
       const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
 
       const auto &   q_points   = fe_fv.get_quadrature_points();
@@ -619,7 +621,7 @@ namespace Step74
 
       const double extent1 = cell->extent_in_direction(
         GeometryInfo<dim>::unit_normal_direction[face_no]);
-      const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+      const double penalty = compute_penalty(degree, extent1, extent1);
 
       double difference_norm_square = 0.;
       for (unsigned int point = 0; point < q_points.size(); ++point)
@@ -633,14 +635,14 @@ namespace Step74
     // 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 auto face_worker = [&](const auto &        cell,
+                                 const unsigned int &f,
+                                 const unsigned int &sf,
+                                 const auto &        ncell,
+                                 const unsigned int &nf,
+                                 const unsigned int &nsf,
+                                 auto &              scratch_data,
+                                 auto &              copy_data) {
       const FEInterfaceValues<dim> &fe_iv =
         scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
 
@@ -668,7 +670,7 @@ namespace Step74
         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);
+      const double penalty = compute_penalty(degree, extent1, extent2);
 
       double flux_jump_square = 0;
       double u_jump_square    = 0;
@@ -684,7 +686,7 @@ namespace Step74
       copy_data_face.values[1] = copy_data_face.values[0];
     };
 
-    auto copier = [&](const CopyData &copy_data) {
+    const auto copier = [&](const auto &copy_data) {
       if (copy_data.cell_index != numbers::invalid_unsigned_int)
         estimated_error_square_per_cell[copy_data.cell_index] +=
           copy_data.value;
@@ -693,10 +695,6 @@ namespace Step74
           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 |
@@ -725,39 +723,37 @@ namespace Step74
   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);
+    const auto cell_worker =
+      [&](const auto &cell, auto &scratch_data, auto &copy_data) {
+        const FEValues<dim> &fe_v = scratch_data.reinit(cell);
 
-      copy_data.cell_index = cell->active_cell_index();
+        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();
+        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_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);
+        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;
-    };
+        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 auto boundary_worker = [&](const auto &        cell,
+                                     const unsigned int &face_no,
+                                     auto &              scratch_data,
+                                     auto &              copy_data) {
       const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
 
       const auto &   q_points   = fe_fv.get_quadrature_points();
@@ -773,7 +769,7 @@ namespace Step74
 
       const double extent1 = cell->extent_in_direction(
         GeometryInfo<dim>::unit_normal_direction[face_no]);
-      const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+      const double penalty = compute_penalty(degree, extent1, extent1);
 
       double difference_norm_square = 0.;
       for (unsigned int point = 0; point < q_points.size(); ++point)
@@ -784,14 +780,14 @@ namespace Step74
       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 auto face_worker = [&](const auto &        cell,
+                                 const unsigned int &f,
+                                 const unsigned int &sf,
+                                 const auto &        ncell,
+                                 const unsigned int &nf,
+                                 const unsigned int &nsf,
+                                 auto &              scratch_data,
+                                 auto &              copy_data) {
       const FEInterfaceValues<dim> &fe_iv =
         scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
 
@@ -813,7 +809,7 @@ namespace Step74
         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);
+      const double penalty = compute_penalty(degree, extent1, extent2);
 
       double u_jump_square = 0;
       for (unsigned int point = 0; point < n_q_points; ++point)
@@ -824,7 +820,7 @@ namespace Step74
       copy_data_face.values[1] = copy_data_face.values[0];
     };
 
-    auto copier = [&](const CopyData &copy_data) {
+    const auto copier = [&](const auto &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)
@@ -832,17 +828,17 @@ namespace Step74
           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);
+    ScratchData scratch_data(mapping,
+                             fe,
+                             QGauss<dim>(fe.degree + 2),
+                             cell_flags,
+                             QGauss<dim - 1>(fe.degree + 2),
+                             face_flags);
 
     CopyData cd;
     MeshWorker::mesh_loop(dof_handler.begin_active(),
@@ -914,16 +910,16 @@ namespace Step74
     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::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;
+    unsigned int max_cycle = test_case == Test_Case::convergence_rate ? 6 : 20;
     for (unsigned int cycle = 0; cycle < max_cycle; ++cycle)
       {
         std::cout << "Cycle " << cycle << std::endl;
@@ -949,7 +945,7 @@ namespace Step74
                 if (cycle == 0)
                   {
                     GridGenerator::hyper_L(triangulation);
-                    triangulation.refine_global(2);
+                    triangulation.refine_global(3);
                   }
                 else
                   {
@@ -962,11 +958,11 @@ namespace Step74
                 Assert(false, ExcNotImplemented());
               }
           }
-        std::cout << "Number of active cells:       "
+        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::cout << "  Number of degrees of freedom : " << dof_handler.n_dofs()
                   << std::endl;
 
         assemble_system();
@@ -982,6 +978,10 @@ namespace Step74
         if (test_case == Test_Case::l_singularity)
           {
             compute_error_estimate();
+            std::cout << "  Estimated error              : "
+                      << std::sqrt(estimated_error_square_per_cell.l1_norm())
+                      << std::endl;
+
             convergence_table.add_value(
               "Estimator",
               std::sqrt(estimated_error_square_per_cell.l1_norm()));
@@ -1009,8 +1009,7 @@ namespace Step74
           convergence_table.evaluate_convergence_rates(
             "H1", ConvergenceTable::reduction_rate_log2);
         }
-
-      std::cout << "degree = " << fe.get_degree() << std::endl;
+      std::cout << "degree = " << degree << std::endl;
       convergence_table.write_text(
         std::cout, TableHandler::TextOutputFormat::org_mode_table);
     }

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.