From: Jiaqi Zhang Date: Wed, 12 Aug 2020 00:13:40 +0000 (-0400) Subject: add step-74 X-Git-Tag: v9.3.0-rc1~687^2~20 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c0d1ebad736caf31060d93b2743a6634cd6d817b;p=dealii.git add step-74 --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 5ed6044eb3..8daefc8652 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -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 index 0000000000..6d166823cd --- /dev/null +++ b/examples/step-74/CMakeLists.txt @@ -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 index 0000000000..20c3064395 --- /dev/null +++ b/examples/step-74/doc/builds-on @@ -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 index 0000000000..7bd4844c96 --- /dev/null +++ b/examples/step-74/doc/intro.dox @@ -0,0 +1,125 @@ +
+ + +This program was contributed by Timo Heister and Jiaqi Zhang. +
+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. +
+ + + +

Symmetric interior penalty Galerkin (SIPG) method for Poisson's equation

+ +

Overview

+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. +
    +
  1. The SIPG method for Poisson's equation, which has been used in step-39 and step-59. +
  2. Assembling of face terms using FEInterfaceValues and the system matrix using MeshWorker::mesh_loop(), which is similar to step-12. +
  3. Adaptive mesh refinement using an error estimator. +
  4. Two test cases: convergence test for a smooth function and adaptive mesh refinement test for a singular solution. +
+ +

The equation

+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_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_F + \biggr\}. +@f} + +

The penalty parameter

+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)$. + +

Posteriori error estimator

+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 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-74/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-74/doc/results.dox b/examples/step-74/doc/results.dox new file mode 100644 index 0000000000..9ed00fd82b --- /dev/null +++ b/examples/step-74/doc/results.dox @@ -0,0 +1,195 @@ +

Results

+ +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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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. + + diff --git a/examples/step-74/doc/tooltip b/examples/step-74/doc/tooltip new file mode 100644 index 0000000000..321d3af5f4 --- /dev/null +++ b/examples/step-74/doc/tooltip @@ -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 index 0000000000..a458e7159c --- /dev/null +++ b/examples/step-74/step-74.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// Here the discontinuous finite elements and FEInterfaceValues are defined. +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +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 + class SmoothSolution : public Function + { + public: + SmoothSolution() + : Function() + {} + virtual void value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; + virtual Tensor<1, dim> + gradient(const Point & point, + const unsigned int component = 0) const override; + }; + + template + void SmoothSolution::value_list(const std::vector> &points, + std::vector & 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 + Tensor<1, dim> + SmoothSolution::gradient(const Point &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 + class SmoothRightHandSide : public Function + { + public: + SmoothRightHandSide() + : Function() + {} + virtual void value_list(const std::vector> &points, + std::vector & 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 + class SingularRightHandSide : public Function + { + public: + SingularRightHandSide() + : Function() + {} + virtual void value_list(const std::vector> &points, + std::vector & 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 + void get_function_jump(const FEInterfaceValues &fe_iv, + const Vector & solution, + std::vector & jump) + { + const unsigned n_q = fe_iv.n_quadrature_points; + std::vector 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 + void get_function_gradient_jump(const FEInterfaceValues &fe_iv, + const Vector & solution, + std::vector> & gradient_jump) + { + const unsigned n_q = fe_iv.n_quadrature_points; + std::vector> 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(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 instead. + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + double values[2]; + unsigned int cell_indices[2]; + }; + + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + double value; + unsigned int cell_index; + template + 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 + 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 triangulation; + const MappingQ1 mapping; + + using ScratchData = MeshWorker::ScratchData; + + FE_DGQ fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + // Vectors to store error estimator square and energy norm square per cell. + Vector estimated_error_square_per_cell; + Vector 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> exact_solution; + std::unique_ptr> 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 + SIPGLaplace::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>(); + rhs_function = std::make_unique>(); + } + + else if (test_case == Test_Case::l_singularity) + { + exact_solution = std::make_unique(); + rhs_function = std::make_unique>(); + } + else + AssertThrow(false, ExcNotImplemented()); + } + + template + void SIPGLaplace::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 + void SIPGLaplace::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 &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 &JxW = scratch_data.get_JxW_values(); + + std::vector 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 &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 & JxW = scratch_data.get_JxW_values(); + const std::vector> &normals = + scratch_data.get_normal_vectors(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + + const double extent1 = cell->extent_in_direction( + GeometryInfo::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 &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 & JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + const double extent1 = + cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); + const double extent2 = ncell->extent_in_direction( + GeometryInfo::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 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 quadrature(n_gauss_points); + QGauss 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 + void SIPGLaplace::solve() + { + std::cout << " Solving system..." << std::endl; + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); + } + + template + void SIPGLaplace::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 data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u", DataOut::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 + void SIPGLaplace::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 &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 &JxW = fe_v.get_JxW_values(); + + std::vector> hessians(n_q_points); + fe_v.get_function_hessians(solution, hessians); + + std::vector 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 &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 &JxW = fe_fv.get_JxW_values(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); + + const double extent1 = cell->extent_in_direction( + GeometryInfo::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 &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_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 & JxW = fe_iv.get_JxW_values(); + const std::vector> &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 jump(n_q_points); + get_function_jump(fe_iv, solution, jump); + + std::vector> 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::unit_normal_direction[f]); + const double extent2 = ncell->extent_in_direction( + GeometryInfo::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 ©_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 quadrature(n_gauss_points); + QGauss 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 + double SIPGLaplace::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 &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 &JxW = fe_v.get_JxW_values(); + + std::vector> grad_u(n_q_points); + fe_v.get_function_gradients(solution, grad_u); + + std::vector> 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 &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 &JxW = fe_fv.get_JxW_values(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); + + const double extent1 = cell->extent_in_direction( + GeometryInfo::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 &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_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 &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 jump(n_q_points); + get_function_jump(fe_iv, solution, jump); + + const double extent1 = + cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); + const double extent2 = ncell->extent_in_direction( + GeometryInfo::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 ©_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 quadrature(n_gauss_points); + QGauss 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 + void SIPGLaplace::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 + void SIPGLaplace::compute_errors() + { + double L2_error, H1_error; + + { + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + *(exact_solution.get()), + difference_per_cell, + QGauss(fe.degree + 2), + VectorTools::L2_norm); + + L2_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + } + + { + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + *(exact_solution.get()), + difference_per_cell, + QGauss(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 + void SIPGLaplace::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 main 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; +}