]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added tutorial program step-82
authorDiane Guignard <dguignar@uottawa.ca>
Mon, 28 Jun 2021 18:41:51 +0000 (14:41 -0400)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 29 Sep 2021 11:37:09 +0000 (13:37 +0200)
doc/doxygen/references.bib
examples/step-82/CMakeLists.txt [new file with mode: 0644]
examples/step-82/doc/builds-on [new file with mode: 0644]
examples/step-82/doc/intro.dox [new file with mode: 0644]
examples/step-82/doc/kind [new file with mode: 0644]
examples/step-82/doc/results.dox [new file with mode: 0644]
examples/step-82/doc/tooltip [new file with mode: 0644]
examples/step-82/step-82.cc [new file with mode: 0644]

index 2a69a20f419aaf732cf94ee506a3fbee4e02a5c0..435dc4bc899ac3b5e2c30ba46738964e33e3c4f2 100644 (file)
   url =       {https://doi.org/10.1023/a:1020533003783}
 }
 
+% ------------------------------------
+% Step 82
+% ------------------------------------
+
+@article{Pryer2014,
+author = {T. Pryer},
+title = {Discontinuous {G}alerkin methods for the p-biharmonic equation from a discrete variational perspective},
+journal = {Electronic Transactions of Numerical Analysis},
+volume = {41},
+pages = {328--349},
+year = {2014}
+}
+
+@article{Bonito2021,
+author = {A. Bonito and R.H. Nochetto and D. Ntogkas},
+title = {{DG} Approach to Large Bending Plate Deformations with Isometry Constraint},
+journal = {Mathematical Models and Methods in Applied Sciences},
+volume = {31},
+number = {1},
+pages = {133--175},
+year = {2021}
+}
+
+@book{DiPietro2011,
+author = {D.A. Di Pietro and A. Ern},
+title = {Mathematical Aspects of Discontinuous {G}alerkin Methods},
+series = {Math{\'e}matiques et Applications},
+publisher={Springer Berlin Heidelberg},
+year={2011},
+}
 
 % ------------------------------------
 % References used elsewhere
diff --git a/examples/step-82/CMakeLists.txt b/examples/step-82/CMakeLists.txt
new file mode 100644 (file)
index 0000000..5558f0a
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-82 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-82")
+
+# 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 3.1.0)
+
+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-82/doc/builds-on b/examples/step-82/doc/builds-on
new file mode 100644 (file)
index 0000000..e7977b2
--- /dev/null
@@ -0,0 +1 @@
+step-47
diff --git a/examples/step-82/doc/intro.dox b/examples/step-82/doc/intro.dox
new file mode 100644 (file)
index 0000000..25f9a9f
--- /dev/null
@@ -0,0 +1,265 @@
+<i>
+  This program was contributed by Andrea Bonito (Texas A&M University) and Diane Guignard (University of Ottawa).
+  
+  This material is based upon work supported by the National Science Foundation under Grant No. DMS-1817691.
+  Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.
+</i>
+
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+<h3>Problem statment</h3>
+
+In this example, we consider the <i>local discontinuous Galerkin</i> (LDG) method for approximating the solution to the bi-Laplacian problem
+@f{align*}{
+\Delta^2 u & = f \quad \mbox{in } \Omega, \\
+\nabla u & = \mathbf{0} \quad \mbox{on } \partial\Omega, \\ 
+u & = 0 \quad \mbox{on } \partial\Omega,
+@f}
+where $\Omega\subset\mathbb{R}^d$ $(d=2,3)$ is an open bounded Lipschitz domain and $f\in L^2(\Omega)$. The weak formulation reads: find $u\in H_0^2(\Omega)$ such that
+@f[
+\int_{\Omega}D^2u:D^2v = \int_{\Omega}fv \qquad \forall \, v\in H_0^2(\Omega),
+@f]
+where $D^2v$ denotes the Hessian of $v$ and $H_0^2(\Omega):=\{v\in H^2(\Omega): \,\, v=0 \mbox{ and } \nabla v=\mathbf{0} \,\, \mbox{ on } \partial\Omega\}$. Using so-called lifting operators as well as the Nitsche approach to impose the homogeneous Dirichlet boundary conditions, the LDG approximation of this problem consists of replacing the Hessians by discrete Hessians (see below) and adding penalty terms involving properly scaled jump terms.
+In particular, the versatility of the method described below is of particular interest for nonlinear problems or problems with intricate weak formulations for which the design of discrete algorithms is challenging.  
+
+<h3>Discretization</h3>
+<h4>Finite element spaces</h4>
+For $h>0$, let $\mathcal{T}_h$ be a partition of $\Omega$ into quadrilateral (hexahedral if $d=3$) elements $K$ of diameter $h_{K}\leq h$ and let $\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^b$ denote the set of (interior and boundary) faces. We restrict the discussion to conforming subdivisions to avoid technicalities already addressed in previous tutorials. The diameter of $e \in \mathcal{E}_h$ is denoted$h_e$. For any integer $k\ge 2$, we introduce the (broken) finite element space
+@f[
+\mathbb{V}_h:=\left\{v_h\in L^2(\Omega): \,\, v_h|_K\circ F_{K}\in\mathbb{Q}_k \quad \forall \, K \in\mathcal{T}_h \right\},
+@f]
+where $F_{K}$ is the map from the reference element $\widehat{K}$ (unit square/cube) to the physical element $K$. For $v_h\in\mathbb{V}_h$, the piecewise differential operators are denoted with a subscript $h$, for instance $\nabla_h v_h|_K=\nabla(v_h|_K)$ and $D_h^2 v_h=\nabla_h\nabla_h v_h$. For $e\in\mathcal{E}_h$, we assign a normal $\mathbf{n}_e$. The choice of normal is irrelevant except that when $e$ is a boundary face, $\mathbf{n}_e$ is the normal pointing outward $\Omega$.
+
+<h4>Jumps, averages, and discrete reconstruction of differential operators</h4>
+The piecewise differential operators do not have enough information to be accurate approximations of their continuous counterparts.
+They are missing inter-elements information. 
+
+This leads to the introductions of jump and average operators:
+@f[
+\jump{v_h}|_e := 
+\left\{\begin{array}{ll}
+v_h|_{K_1}-v_h|_{K_2} & e\in\mathcal{E}_h^0 \\
+v_h|_{K_1} & e\in\mathcal{E}_h^b
+\end{array}\right. \quad \mbox{and} \quad \average{v_h}|_e := 
+\left\{\begin{array}{ll}
+\frac{1}{2}(v_h|_{K_1}+v_h|_{K_2}) & e\in\mathcal{E}_h^0 \\
+v_h|_{K_1} & e\in\mathcal{E}_h^b,
+\end{array}\right. 
+@f]
+respectively, where $K_1$ and $K_2$ are the two elements adjacent to $e$ so that $\mathbf{n}_e$ points from $K_1$ to $K_2$ (with obvious modification when $e$ is a boundary edge). 
+
+With these notations, we are now in position to define the discrete/reconstructed Hessian $H_h(v_h)\in\left[L^2(\Omega)\right]^{d\times d}$ of $v_h\in\mathbb{V}_h$. We first consider two <i>local lifting operators</i> $r_e:[L^2(e)]^d\rightarrow[\mathbb{V}_h]^{d\times d}$ and $b_e:L^2(e)\rightarrow[\mathbb{V}_h]^{d\times d}$ defined for $e\in\mathcal{E}_h$ by, respectively,
+@f[ 
+r_e(\boldsymbol{\phi}) \in [\mathbb{V}_h]^{d\times d}: \,
+\int_{\Omega}r_e(\boldsymbol{\phi}):\tau_h = \int_e\average{\tau_h}\mathbf{n}_e\cdot\boldsymbol{\phi} \qquad \forall \, \tau_h\in [\mathbb{V}_h]^{d\times d}
+@f]
+and
+@f[
+b_e(\phi) \in [\mathbb{V}_h]^{d\times d}: \,
+\int_{\Omega} b_e(\phi):\tau_h = \int_e\average{{\rm div}\, \tau_h}\cdot\mathbf{n}_e\phi \qquad \forall \, \tau_h\in [\mathbb{V}_h]^{d\times d}.
+@f]
+We have ${\rm supp}\,(r_e(\boldsymbol{\phi}))={\rm supp}\,(b_e(\phi))=\omega_e$, where $\omega_e$ denotes the patch of elements having $e$ has part of their boundaries.
+
+The discrete Hessian operator $H_h:\mathbb{V}_h\rightarrow\left[L^2(\Omega)\right]^{2\times 2}$ is then given by
+@f[
+H_h(v_h) := D_h^2 v_h -R_h(\jump{\nabla_h v_h})+B_h(\jump{v_h}) := D_h^2 v_h - \sum_{e\in\mathcal{E}_h}r_e(\jump{\nabla_h v_h})+\sum_{e\in\mathcal{E}_h}b_e(\jump{v_h}).
+@f]
+
+@note
+In general, the polynomial degree of the finite element space for the two lifting terms do not need to be the same as the one used for the approximate solution. A different polynomial degree for each lifting term can also be considered.
+
+Note that other differential operators (e.g., gradient or divergence) can be reconstructed in a similar fashion, see for instance @cite DiPietro2011.
+
+<h4>Motivation for the lifting operators</h4>
+
+The discrete Hessian is designed such that it weakly converges to the continuous Hessian, see the note in the next section for a precise statement. As already mentioned above, the broken Hessian is not a suitable candidate as it contains no information about inter-element jumps. We provide here an informal discussion motivating the definition of the two lifting operators and we refer to @cite Pryer2014 and @cite Bonito2021 for more details (although the definitions are slightly different unless the mesh is affine). The goal is then to construct a discrete operator $H_h$ such that for all $\tau\in [C_0^{\infty}(\Omega)]^{d\times d}$ we have
+@f[
+\int_{\Omega}H_h(v_h):\tau\longrightarrow \int_{\Omega}D^2v:\tau \qquad \mbox{as } \,\, h\rightarrow 0
+@f]
+for any sequence $\{v_h\}_{h>0}$ in $\mathbb{V}_h$ such that $v_h\rightarrow v$ in $L^2(\Omega)$ as $h\rightarrow 0$ for some $v\in H^2(\Omega)$. Let $\tau\in [C_0^{\infty}(\Omega)]^{d\times d}$. Integrating by parts twice we get
+@f[
+\int_{\Omega}D^2v:\tau = -\int_{\Omega}\nabla v\cdot \mbox{div}(\tau) = \int_{\Omega}v \mbox{ div}(\mbox{div}(\tau))
+@f]
+while
+@f[
+\int_{\Omega}v_h \mbox{ div}(\mbox{div}(\tau)) \longrightarrow \int_{\Omega}v \mbox{ div}(\mbox{div}(\tau)) \qquad \mbox{as } \,\, h\rightarrow 0.
+@f]
+Now, we integrate two times by parts the left term, taking into account that $v_h$ is not necessarily continuous across interior faces. For any $K\in\mathcal{T}_h$ we have
+@f[
+\int_K v_h \mbox{ div}(\mbox{div}(\tau)) = -\int_K \nabla v_h\cdot \mbox{div}(\tau) + \int_{\partial K} v_h \mbox{ div}(\tau)\cdot \mathbf{n}_K =\int_K D^2v_h:\tau - \int_{\partial K}\nabla v_h\cdot (\tau\mathbf{n}_K) + \int_{\partial K} v_h \mbox{ div}(\tau)\cdot \mathbf{n}_K,
+@f]
+where $\mathbf{n}_K$ denotes the outward unit normal to $K$. Then, summing over the elements $K\in\mathcal{T}_h$ and using that $\tau$ is smooth, we obtain
+@f[
+\int_{\Omega} v_h \mbox{ div}(\mbox{div}(\tau)) = \int_{\Omega} D_h^2v_h:\tau - \sum_{e\in\mathcal{E}_h}\int_e\jump{\nabla_h v_h}\cdot \average{\tau}\mathbf{n}_e + \sum_{e\in\mathcal{E}_h}\int_e v_h \average{\mbox{div}(\tau)}\cdot \mathbf{n}_e
+@f]
+which reveals the motivation for the definition of the two lifting operators: if $\tau$ was an admissible test function, then the right-hand side would be equal to $\int_{\Omega}H_h(v_h):\tau$ and we would have shown the desired (weak) convergence. Actually, if we add and subtract $\tau_h$, the Lagrange interpolant of $\tau$ in $[\mathbb{V}_h\cap H_0^1(\Omega)]^{d\times d}$, we can show that the right-hand side is indeed equal to $\int_{\Omega}H_h(v_h):\tau$ up to terms that tends to zero as $h\rightarrow 0$ under appropriate assumptions on $v_h$.
+
+It is worth mentioning that defining $H_h$ without the lifting operators $r_e$ and $b_e$ for $e\in\mathcal{E}_h^b$ would not affect the weak convergence property (the integrals over boundary faces are zero since $\tau$ is compactly supported in $\Omega$). However, they are included in $H_h$ to ensure that the solution of the discrete problem introduced in the next section satisfies the homogeneous Dirichlet boundary conditions in the limit $h\rightarrow 0$.
+
+
+<h4>LDG approximations</h4>
+The proposed LDG approximation of the bi-Laplacian problem reads: find $u_h\in\mathbb{V}_h$ such that
+@f[
+A_h(u_h,v_h):= a_h(u_h,v_h)+j_h(u_h,v_h) = F_h(v_h) \qquad \forall \, v_h\in\mathbb{V}_h,
+@f]
+where
+@f{align*}{
+       a_h(u_h,v_h) & := \int_{\Omega}H_h(u_h):H_h(v_h) \\
+       j_h(u_h,v_h) & := \gamma_1\sum_{e\in\mathcal{E}_h}h_e^{-1}\int_e\jump{\nabla_h u_h}\cdot\jump{\nabla_h v_h}+\gamma_0\sum_{e\in\mathcal{E}_h}h_e^{-3}\int_e\jump{u_h}\jump{v_h} \\
+       F_h(v_h) & := \int_{\Omega}fv_h.
+@f}
+Here, $\gamma_0,\gamma_1>0$ are penalty parameters.
+
+Let $\{\varphi_i\}_{i=1}^{N_h}$ be the standard basis functions that generate $\mathbb{V}_h$. We can then express the solution as $u_h=\sum_{j=1}^{N_h}U_j\varphi_j$ and the problem reads: find $\boldsymbol{U}=(U_j)_{j=1}^{N_h}\in\mathbb{R}^{N_h}$ such that
+@f[
+A\boldsymbol{U} = \boldsymbol{F},
+@f]
+where $A=(A_{ij})_{i,j=1}^{N_h}\in\mathbb{R}^{N_h\times N_h}$ and $\boldsymbol{F}=(F_i)_{i=1}^{N_h}\in\mathbb{R}^{N_h}$ are defined by
+@f[
+A_{ij}:=A_h(\varphi_j,\varphi_i) \quad \text{and} \quad F_i:=F_h(\varphi_i), \qquad 1\leq i,j \leq N_h.
+@f]
+
+@note
+The sparsity pattern associated with the above LDG method is slightly larger than that of, e.g., the symmetric interior penalty discontinuous Galerkin (SIPG) method. However, we have the following interesting properties:
+<ol>
+  <li>The bilinear form $A_h(\cdot,\cdot)$ is coercive with respect to the DG $H^2$ norm
+  @f[
+  \|v_h\|_{H_h^2(\Omega)}^2:=\|D_h^2v_h\|_{L^2(\Omega)}^2+\sum_{e\in\mathcal{E}_h}h_e^{-1}\|\jump{\nabla_h v_h}\|_{L^2(e)}^2+\sum_{e\in\mathcal{E}_h}h_e^{-3}\|\jump{v_h}\|_{L^2(e)}^2
+  @f]
+for any choice of penalty parameters $\gamma_0,\gamma_1>0$. In other words, the stability of the method is ensured for any positive parameters. This is in contrast with interior penalty methods for which they need to be large enough.
+  <li>If $\{v_h\}_{h>0}\subset \mathbb{V}_h$ is a sequence uniformly bounded in the $\|\cdot\|_{H_h^2(\Omega)}$ norm such that $v_h\rightarrow v$ in $L^2(\Omega)$ as $h\rightarrow 0$ for some $v\in H^2(\Omega)$, then the discrete Hessian $H_h(v_h)$ weakly converges to $D^2v$ in $[L^2(\Omega)]^{2\times 2}$ as $h\rightarrow 0$. Note that the uniform boundedness assumption implies that the limit $v$ belongs to $H_0^2(\Omega)$.</li>
+  <li>The use of a reconstructed operator simplifies the design of the numerical algorithm. In particular, no integration by parts is needed to derive the discrete problem. This strategy of replacing differential operators by appropriate discrete counter-parts can be applied to nonlinear and more general problems, for instance variational problems without a readily accessible strong formulation.
+  </li>
+</ol>
+
+As in step-47, we could consider $C^0$ finite element approximations by replacing <code>FE_DGQ<dim></code> by <code>FE_Q<dim></code> (and include the appropriate header file <code>deal.II/fe/fe_q.h</code>) in the program below. In this case, the jump of the basis functions across any interior face is zero, and thus $b_e([\varphi_i])=\mathbf{0}$ for all $e\in\mathcal{E}_h^0$, and could be dropped to save computational time. While an overkill for the bi-Laplacian problem, the flexibility of fully discontinuous methods combined with reconstructed differential operators is advantageous for nonlinear problems.
+
+
+<h3>Implementation</h3>
+
+As customary, we assemble the matrix $A$ and the right-hand side $\boldsymbol{F}$ by looping over the elements $K\in\mathcal{T}_h$. Since we are using discontinuous finite elements, the support of each $\varphi_i$ is only one element $K\in\mathcal{T}_h$. However, due to the lifting operators, the support of $H_h(\varphi_i)$ is $K$ plus all the neighbors of $K$ (recall that for $e\in \mathcal{E}_h$, the support of the lifting operators $r_e$ and $b_e$ is $\omega_e$). Therefore, when integrating over a cell $K_c$, we need to consider the following interactions (case $d=2$)
+
+<table align="center">
+  <tr>
+    <td><img src="https://www.dealii.org/images/steps/developer/step-82.schema_cells.png" width="250px" height="250px" alt=""></td>
+    <td>&nbsp;</td>
+    <td>
+      <table  align="center">
+        <tr>
+          <td align="left">dofs $K_c$ $\leftrightarrow$ dofs $K_c$</td>
+          <td>&nbsp;</td>
+          <td align="left">(stored in <code>stiffness_matrix_cc</code>)</td>
+        </tr>
+        <tr>
+          <td align="left">dofs $K_c$ $\leftrightarrow$ dofs $K_{n_k}$</td>
+          <td>&nbsp;</td>
+          <td align="left">(stored in <code>stiffness_matrix_cn</code> and <code>stiffness_matrix_cn</code>)</td>
+        </tr>
+        <tr>
+          <td align="left">dofs $K_{n_k}$ $\leftrightarrow$ dofs $K_{n_k}$</td>
+          <td>&nbsp;</td>
+          <td align="left">(stored in <code>stiffness_matrix_nn</code>)</td>
+        </tr>
+        <tr>
+          <td align="left">dofs $K_{n_k}$ $\leftrightarrow$ dofs $K_{n_l}$, $k\ne l$</td>
+          <td>&nbsp;</td>
+          <td align="left">(stored in <code>stiffness_matrix_n1n2</code> and <code>stiffness_matrix_n2n1</code>)</td>
+        </tr>
+      </table>         
+    </td>
+  </tr>
+</table>
+
+namely we need to compute the discrete Hessian of all the basis functions with support on $K_c$ as well as all the basis functions with support on the neighboring cells of $K_c$. This is done in the function <code>compute_discrete_hessians</code>. A cell $K_c$ can have fewer than four neighbors (six when $d=3$) when at least one face $e\subset\partial K_c$ belongs to $\mathcal{E}_h^b$. It can also have more neighbors when hanging nodes are present. To simplify the presentation we do not discuss the later. 
+
+Due to the local support of the basis functions, many of the terms of the discrete Hessian are zero. For any basis function $\varphi^c$ with support on $K_c$ we have $r_e(\jump{\nabla_h\varphi^c})\not\equiv 0$ only if $e\subset\partial K_c$, and similarly for $b_e(\jump{\varphi^c})$. Therefore, the discrete Hessian of $\varphi^c$ reduces to
+@f[
+H_h(\varphi^c)=D_h^2\varphi^c-\sum_{e\subset\partial K}r_e(\jump{\nabla_h \varphi^c})+\sum_{e\subset\partial K}b_e(\jump{\varphi^c}).
+@f]
+Furthermore, since we integrate on $K_c$, we only need to evaluate the discrete Hessian at quadrature points $x_q$ that belong to $K_c$, namely $H_h(\varphi^c)(x_q)$. We store this information in
+@f[
+{\rm compute\_discrete\_hessians[i][q]}, \qquad 0\leq {\rm i} < {\rm n\_dofs}, \,\, 0\leq {\rm q} < {\rm n\_q\_points},
+@f]
+where <code>n_dofs = fe_values.dofs_per_cell</code> is the number of degrees of freedom per cell and <code>n_q_points = quad.size()</code> is the number of quadrature points on $K_c$. For any basis function $\varphi^n$ with support on a neighboring cell, the discrete Hessian $H_h(\varphi^n)$ evaluated on $K_c$ contains only the two lifting terms since $\varphi^n|_{K}\equiv 0$. Moreover, only the lifting over the common face $e$ is nonzero on $K_c$, namely for all $x_q\in K_c$
+@f[
+H_h(\varphi^n)(x_q)=-r_e(\jump{\nabla_h\varphi^n})(x_q)+b_e(\jump{\varphi^n})(x_q).
+@f]
+This information is stored in 
+@f[
+{\rm compute\_discrete\_hessians\_neigh[face\_no][i][q]}, \qquad 0\leq {\rm face\_no} < {\rm n\_faces}, \,\, 0\leq {\rm i} < {\rm n\_dofs}, \,\, 0\leq {\rm q} < {\rm n\_q\_points},
+@f]
+where <code>n_dofs</code> and <code>n_q_points</code> are as above, and <code>n_faces = GeometryInfo<dim>::faces_per_cell</code> is the number of faces of $K_c$. As we shall see in the next section, we will only need to solve <i>half</i> of the local problems for the lifting terms.
+
+@note
+The variable <code>discrete_hessians_neigh</code> is of size <code>n_faces x n_dofs x n_q_points</code>. However, we only need to consider the interior faces, namely we do not need to fill <code>discrete_hessians_neigh[face_no][i][q]</code> whenever <code>face_no</code> corresponds to a boundary face. We could then save a little bit of storage by considering $0\leq {\rm face\_no} < {\rm n\_faces}$ with <code>n_faces</code> the actual number of neighboring cells, i.e., not counting the boundary faces. By doing so, we could also avoid testing if a face lies on the boundary in the assembly of the matrix.
+
+
+<h4>Computation of the lifting terms</h4>
+
+We now describe the computation of the lifting operators $r_e$ and $b_e$ on each cell $K_c$.
+We focus on $b_e$ for an interior face $e\in\mathcal{E}_h^0$, but the other cases are treated similarly.
+
+We have $e=\partial K_c\cap \partial K_n$ for some neighbor $K_n$ of $K_c$. For a basis function $\varphi\in\mathbb{V}_h$ with support on $K_c$ or $K_n$ (for the other basis functions we have $b_e(\jump{\varphi})\equiv 0$), we write $b_e(\jump{\varphi})\in[\mathbb{V}_h]^{d\times d}$ as
+@f[
+b_e(\jump{\varphi})=\sum_{n=1}^{N_c+N_n}B_n\psi_n,
+@f]
+where $\{\psi_n\}_{n=1}^{N_c}$ and $\{\psi_n\}_{n=N_c+1}^{N_c+N_n}$ are the basis functions of $[\mathbb{V}_h]^{d\times d}$ which have support on $K_c$ and $K_n$, respectively. The coefficients $\boldsymbol{B}=(B_n)_{n=1}^{N_c+N_n}\in\mathbb{R}^{N_c+N_c}$ of the lifting operator $b_e$ are obtain upon solving the linear system
+@f[
+M\boldsymbol{B}=\boldsymbol{G},
+@f]
+where the components of the (local) mass matrix and the right-hand side are given respectively by
+@f[
+M_{mn}:=\int_{\Omega}\psi_n:\psi_m \quad \mbox{and} \quad G_m:=\int_e\average{{\rm div}\, \psi_n}\cdot \mathbf{n}_e\jump{\varphi}, \qquad 1\leq m,n \leq N_c+N_n.
+@f]
+Note that this system has the decoupled form
+@f[
+\left[\begin{array}{cc}
+M_c & \mathbf{0} \\
+\mathbf{0} & M_n
+\end{array}\right]\left[\begin{array}{c}
+\boldsymbol{B}_c \\ \boldsymbol{B}_n
+\end{array}\right]=\left[\begin{array}{c}
+\boldsymbol{G}_c \\ \boldsymbol{G}_n
+\end{array}\right]
+@f]
+with $M_c\in\mathbb{R}^{N_c\times N_c}$, $M_n\in\mathbb{R}^{N_n\times N_n}$, $\boldsymbol{B}_c,\boldsymbol{G}_c\in\mathbb{R}^{N_c}$, and $\boldsymbol{B}_n,\boldsymbol{G}_n\in\mathbb{R}^{N_n}$.
+
+In fact, since we evaluate the discrete Hessians at quadrature points $x_q\in K_c$ and $\psi_n|_{K_c}\equiv 0$ for $n=N_c+1,\ldots,N_c+N_n$, we have
+@f[
+b_e(\jump{\varphi})(x_q)=\sum_{n=1}^{N_c+N_n}B_n\psi_n(x_q)=\sum_{n=1}^{N_c}B_n\psi_n(x_q).
+@f]
+As a consequence, only the coefficients $B_n$, $n=1,\ldots,N_c$, are needed.
+
+To compute the components $G_m$, $m=1,\ldots,N_c$, we take advantage of the relation
+@f[
+\mathbf{n}_e\jump{\varphi}=\mathbf{n}_{K_c}\varphi|_{K_c}+\mathbf{n}_{K_n}\varphi|_{K_n},
+@f]
+where $\mathbf{n}_{K_c}$ (resp. $\mathbf{n}_{K_n}$) denotes the outward unit normal to $K_c$ (resp. $K_n$). Therefore, if $\varphi=\varphi^c$, namely $\varphi$ has support on the current cell $K_c$, then
+@f[
+G_m=\int_e\average{{\rm div}\, \psi_m}\cdot\mathbf{n}_e\jump{\varphi^c}=\frac{1}{2}\int_e{\rm div}\, \psi_m\cdot\mathbf{n}_{K_c}\varphi^c,
+@f]
+while if $\varphi=\varphi^n$, namely $\varphi$ has support on the neighbhoring cell $K_n$, then
+@f[
+G_m=\int_e\average{{\rm div}\, \psi_m}\cdot\mathbf{n}_e\jump{\varphi^n}=\frac{1}{2}\int_e{\rm div}\, \psi_m\cdot\mathbf{n}_{K_n}\varphi^n.
+@f]
+The factor $\frac{1}{2}$ comes from the average operator as $e$ is assumed to be an interior face.
+
+<h3>Test case</h3>
+
+The performances of the numerical algorithm will be assessed using a manufactured solution $u:(0,1)^d\rightarrow\mathbb{R}$ given by
+@f[
+u(x,y)=x^2(1-x)^2y^2(1-y)^2
+@f]
+if $d=2$, while if $d=3$ we take
+@f[
+u(x,y,z)=x^2(1-x)^2y^2(1-y)^2z^2(1-z)^2.
+@f]
+
+For different values of $h$, we will report the distortion $u-u_h$ measured in the discrete $H^2$ metric (defined above but extended to piecewise $H^2$ functions), the discrete $H^1$ metric
+@f[
+\|v\|_{H_h^1(\Omega)}^2 := \|\nabla_h v\|_{L^2(\Omega)}^2+\sum_{e\in\mathcal{E}_h}h_e^{-1}\|\jump{v}\|_{L^2(e)}^2, \quad v\in \prod_{K\in\mathcal{T}_h}H^1(K),
+@f]
+as well as the $L^2$ metric.
+
+
diff --git a/examples/step-82/doc/kind b/examples/step-82/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-82/doc/results.dox b/examples/step-82/doc/results.dox
new file mode 100644 (file)
index 0000000..8aee1cf
--- /dev/null
@@ -0,0 +1,144 @@
+<h1>Results</h1>
+
+
+When running the program, the sparsity pattern is written to an svg file, the solution is written to a vtk file, and some results are printed to the console. With the current setup, the output should read
+
+@code
+
+Number of active cells: 64
+Number of degrees of freedom: 576
+Assembling the system.............
+Done. 
+DG H2 norm of the error: 0.0151063
+DG H1 norm of the error: 0.000399747
+   L2 norm of the error: 5.33856e-05
+
+@endcode
+
+This corresponds to the bi-Laplacian problem with the manufactured solution mentioned above for $d=2$, 3 refinements of the mesh, degree $k=2$, and $\gamma_0=\gamma_1=1$ for the penalty coefficients. By changing the number of refinements, we get the following results:
+
+<table align="center" class="doxtable">
+  <tr>
+    <th>n_ref</th>
+    <th>n_cells</th>
+    <th>n_dofs</th>
+    <th>error H2 </th>
+    <th>rate</th>
+    <th>error H1</th>
+    <th>rate</th>
+    <th>error L2</th>
+    <th>rate</th>
+  </tr>
+  <tr>
+    <td align="center">1</td>
+    <td align="right">4</td>
+    <td align="right">36</td>
+    <td align="center">5.651e-02</td>
+    <td align="center">--</td>
+    <td align="center">3.366e-03</td>
+    <td align="center">--</td>
+    <td align="center">3.473e-04</td>
+    <td align="center">--</td>
+  </tr>
+  <tr>
+    <td align="center">2</td>
+    <td align="right">16</td>
+    <td align="right">144</td>
+    <td align="center">3.095e-02</td>
+    <td align="center">0.87</td>
+    <td align="center">1.284e-03</td>
+    <td align="center">1.39</td>
+    <td align="center">1.369e-04</td>
+    <td align="center">1.34</td>
+  </tr>
+  <tr>
+    <td align="center">3</td>
+    <td align="right">64</td>
+    <td align="right">576</td>
+    <td align="center">1.511e-02</td>
+    <td align="center">1.03</td>
+    <td align="center">3.997e-04</td>
+    <td align="center">1.68</td>
+    <td align="center">5.339e-05</td>
+    <td align="center">1.36</td>
+  </tr>
+  <tr>
+    <td align="center">4</td>
+    <td align="right">256</td>
+    <td align="right">2304</td>
+    <td align="center">7.353e-03</td>
+    <td align="center">1.04</td>
+    <td align="center">1.129e-04</td>
+    <td align="center">1.82</td>
+    <td align="center">1.691e-05</td>
+    <td align="center">1.66</td>
+  </tr>
+  <tr>
+    <td align="center">5</td>
+    <td align="right">1024</td>
+    <td align="right">9216</td>
+    <td align="center">3.609e-03</td>
+    <td align="center">1.03</td>
+    <td align="center">3.024e-05</td>
+    <td align="center">1.90</td>
+    <td align="center">4.789e-06</td>
+    <td align="center">1.82</td>
+  </tr>
+  <tr>
+    <td align="center">6</td>
+    <td align="right">4096</td>
+    <td align="right">36864</td>
+    <td align="center">1.785e-03</td>
+    <td align="center">1.02</td>
+    <td align="center">7.850e-06</td>
+    <td align="center">1.95</td>
+    <td align="center">1.277e-06</td>
+    <td align="center">1.91</td>
+  </tr>
+  <tr>
+    <td align="center">7</td>
+    <td align="right">16384</td>
+    <td align="right">147456</td>
+    <td align="center">8.870e-04</td>
+    <td align="center">1.01</td>
+    <td align="center">2.001e-06</td>
+    <td align="center">1.97</td>
+    <td align="center">3.298e-07</td>
+    <td align="center">1.95</td>
+  </tr>
+</table>
+
+
+<h3>Possible extensions</h3>
+
+The code can be easily adapted to deal with the following cases:
+
+<ol>
+  <li>Non-homogeneous Dirichlet boundary conditions on (part of) the boundary $\partial \Omega$ of $\Omega$.</li>
+  <li>Hanging-nodes (proceed as in step-14 to not visit a sub-face twice when computing the lifting terms in <code>compute_discrete_hessians</code> and the penalty terms in <code>assemble_matrix</code>).</li>
+  <li>LDG method for the Poisson problem (use the discrete gradient consisting of the broken gradient and the lifting of the jump of $u_h$).</li>
+</ol>
+
+We give below additional details for the first point.
+
+<h4>Non-homogeneous Dirichlet boundary conditions</h4>
+If we prescribe non-homogeneous Dirichlet conditions, say
+@f[
+\nabla u=\mathbf{g} \quad \mbox{and} \quad u=g \qquad \mbox{on } \partial \Omega, 
+@f]
+then the right-hand side $\boldsymbol{F}$ of the linear system needs to be modified as follows
+@f[
+F_i:=\int_{\Omega}f\varphi_i-\sum_{e\in\mathcal{E}_h^b}\int_{\Omega}r_e(\mathbf{g}):H_h(\varphi_i)+\sum_{e\in\mathcal{E}_h^b}\int_{\Omega}b_e(g):H_h(\varphi_i)+\gamma_1\sum_{e\in\mathcal{E}_h^b}h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0\sum_{e\in\mathcal{E}_h^b}h_e^{-3}\int_e g\varphi_i, \qquad 1\leq i \leq N_h.
+@f]
+Note that for any given index $i$, many of the terms are zero. Indeed, for $e\in \mathcal{E}_h^b$ we have ${\rm supp}\,(r_e(\mathbf{g}))={\rm supp}\,(b_e(g))=K$, where $K$ is the element for which $e\subset\partial K$. Therefore, the liftings $r_e(\mathbf{g})$ and $b_e(g)$ contributes to $F_i$ only if $\varphi_i$ has support on $K$ or a neighbor of $K$. In other words, when integrating on a cell $K$, we need to add
+@f[
+\int_{K}f\varphi_i+\sum_{e\in\mathcal{E}_h^b, e\subset\partial K}\left[-\int_{K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{K}b_e(g):H_h(\varphi_i)+\gamma_1h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0h_e^{-3}\int_e g\varphi_i\right]
+@f]
+to $F_i$ for the indices $i$ such that $\varphi_i$ has support on $K$ and
+@f[
+\sum_{e\in\mathcal{E}_h^b, e\subset\partial K}\left[-\int_{K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{K}b_e(g):H_h(\varphi_i)\right]
+@f]
+to $F_i$ for the indices $i$ such that $\varphi_i$ has support on a neighbor of $K$.
+
+@note
+Note that we can easily consider the case where Dirichlet boundary conditions are imposed only on a subset $\emptyset\neq\Gamma_D\subset\partial \Omega$. In this case, we simply need to replace $\mathcal{E}_h^b$ by $\mathcal{E}_h^D\subset\mathcal{E}_h^b$ consisting of the faces belonging to $\Gamma_D$. This also affects the matrix $A$ (simply replace $\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^b$ by $\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^D$).
diff --git a/examples/step-82/doc/tooltip b/examples/step-82/doc/tooltip
new file mode 100644 (file)
index 0000000..1f4fa0a
--- /dev/null
@@ -0,0 +1 @@
+Local discontinuous Galerkin method for the bi-Laplacian problem
diff --git a/examples/step-82/step-82.cc b/examples/step-82/step-82.cc
new file mode 100644 (file)
index 0000000..6e84beb
--- /dev/null
@@ -0,0 +1,1203 @@
+
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2021 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Authors: Andrea Bonito and Diane Guignard, 2021.
+ */
+
+// @sect3{Include files}
+
+// All the include files have already been discussed in previous tutorials.
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+// The following three header files are for the solvers.
+// The linear system is solved using a direct method
+// while the conjugate gradient method is used to solve
+// the local problems for the lifting terms.
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <fstream>
+#include <iostream>
+
+
+namespace Step82
+{
+
+  using namespace dealii;
+
+  // @sect3{The <code>BiLaplacianLDGLift</code> class template}
+
+  // The main class of this program is similar to that of step-3
+  // or step-20, as well as many other tutorial programs. The key
+  // function here is <code>discrete_hessians</code> which compute
+  // the discrete Hessians needed for the assembly of the matrix $A$.
+  template <int dim>
+  class BiLaplacianLDGLift
+  {
+  public:
+
+    BiLaplacianLDGLift(const unsigned int fe_degree, double penalty_jump_grad, double penalty_jump_val);
+
+    void run();
+
+  private:
+
+    void make_grid();
+    void setup_system();
+    void assemble_system();
+    void assemble_matrix();
+    void assemble_rhs();
+
+    void solve();
+  
+    void compute_errors();
+    void output_results() const;
+
+    // As indicated by its name, the function <code>assemble_local_matrix</code>
+    // is used for the assembly of the (local) mass matrix used to compute the
+    // two lifting terms (see the matrix $\boldsymbol{M}_c$ introduced in
+    // the introduction when describing the computation of $b_e$). The function
+    // <code>compute_discrete_hessians</code> computes the required discrete Hessians:
+    // the discrete Hessians of the basis functions with support on the current
+    // <code>cell</code> (stored in the output variable <code>discrete_hessians</code>)
+    // and the basis functions with support on a neighbor of the current <code>cell</code>
+    // (stored in the output variable <code>discrete_hessians_neigh</code>).
+    // More precisely, <code>discrete_hessians[i][q_point]</code> stores
+    // $H_h(\varphi_i)(x_q)$, where $\varphi_i$ is a basis function with support
+    // on cell, while <code>discrete_hessians_neigh[face_no][i][q_point]</code> stores
+    // $H_h(\varphi_i)(x_q)$, where $\varphi_i$ is a basis function of the neighboring
+    // cell adjacent to the face <code>face=cell->face(face_no)</code>.
+    void assemble_local_matrix(const FEValues<dim> &fe_values_lift, const unsigned int n_q_points, FullMatrix<double> &local_matrix);
+
+    void compute_discrete_hessians(const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                   const typename DoFHandler<dim>::active_cell_iterator &cell_lift,
+                                   std::vector<std::vector<Tensor<2,dim>>> &discrete_hessians,
+                                   std::vector<std::vector<std::vector<Tensor<2,dim>>>> &discrete_hessians_neigh);
+
+    Triangulation<dim> triangulation;
+
+    FE_DGQ<dim>     fe;
+    DoFHandler<dim> dof_handler;
+
+    // We also need variables for the finite element space
+    // $[\mathbb{V}_h]^{d\times d}$ used for the two lifting
+    // operators.
+    FESystem<dim>   fe_lift;
+    DoFHandler<dim> dof_handler_lift;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> matrix;
+    Vector<double>       rhs;
+    Vector<double>       solution;
+
+    // Finaly, the last two variables correspond to the penalty coefficients
+    // $\gamma_1$ and $\gamma_0$ for the jump of $\nabla_hu_h$ and $u_h$,
+    // respectively.
+    double penalty_jump_grad;
+    double penalty_jump_val;
+
+  };
+
+
+
+  // @sect3{Equation data}
+
+  // This class implement the right-hand side $f=\Delta^2 u$ corresponding to
+  // the manufactured solution $u$ defined in the introduction.
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+  public:
+    RightHandSide () : Function<dim>() {}
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const override;
+  };
+
+  template <int dim>
+  double RightHandSide<dim>::value (const Point<dim> &p,
+                                    const unsigned int /*component*/) const
+  {
+    double return_value = 0.0;
+
+    if (dim==2){
+
+      return_value = 24.0*std::pow(p(1)*(1.0-p(1)),2)+
+                     +24.0*std::pow(p(0)*(1.0-p(0)),2)
+                     +2.0*(2.0-12.0*p(0)+12.0*p(0)*p(0))*(2.0-12.0*p(1)+12.0*p(1)*p(1));
+
+    } else if (dim==3){
+
+      return_value = 24.0*std::pow(p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2)
+                     +24.0*std::pow(p(0)*(1.0-p(0))*p(2)*(1.0-p(2)),2)
+                     +24.0*std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2)
+                     +2.0*(2.0-12.0*p(0)+12.0*p(0)*p(0))*(2.0-12.0*p(1)+12.0*p(1)*p(1))*std::pow(p(2)*(1.0-p(2)),2)
+                     +2.0*(2.0-12.0*p(0)+12.0*p(0)*p(0))*(2.0-12.0*p(2)+12.0*p(2)*p(2))*std::pow(p(1)*(1.0-p(1)),2)
+                     +2.0*(2.0-12.0*p(1)+12.0*p(1)*p(1))*(2.0-12.0*p(2)+12.0*p(2)*p(2))*std::pow(p(0)*(1.0-p(0)),2);
+
+    }
+
+    return return_value;
+  }
+
+
+
+  // This class implement the manufactured (exact) solution $u$. To compute the
+  // errors, we need the value of $u$ as well as its gradient and its Hessian. 
+  template <int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    ExactSolution () : Function<dim>() {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const override;
+
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                    const unsigned int  component = 0) const override;
+    
+    virtual SymmetricTensor<2,dim> hessian (const Point<dim>   &p,
+                                            const unsigned int  component = 0) const override;  
+  };
+
+
+
+  template <int dim>
+  double ExactSolution<dim>::value (const Point<dim> &p,
+                                    const unsigned int /*component*/) const
+  {
+    double return_value = 0.0;
+    
+    if (dim==2){
+      return_value = std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2);
+    } else if (dim==3){
+      return_value = std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2);
+    }
+    
+    return return_value;
+  }
+
+
+
+  template <int dim>
+  Tensor<1,dim> ExactSolution<dim>::gradient (const Point<dim> &p,
+                                              const unsigned int /*component*/) const
+  {
+    Tensor<1,dim> return_gradient;
+    return_gradient = 0.0;
+
+    if (dim==2){
+      return_gradient[0] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * std::pow(p(1)*(1.0-p(1)),2);
+      return_gradient[1] = (2.0*p(1)-6.0*std::pow(p(1),2)+4.0*std::pow(p(1),3)) * std::pow(p(0)*(1.0-p(0)),2);
+    } else if (dim==3){
+      return_gradient[0] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * std::pow(p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2);
+      return_gradient[1] = (2.0*p(1)-6.0*std::pow(p(1),2)+4.0*std::pow(p(1),3)) * std::pow(p(0)*(1.0-p(0))*p(2)*(1.0-p(2)),2);
+      return_gradient[2] = (2.0*p(2)-6.0*std::pow(p(2),2)+4.0*std::pow(p(2),3)) * std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2);
+    }
+    
+    return return_gradient;
+  }
+
+
+
+  template <int dim>
+  SymmetricTensor<2,dim> ExactSolution<dim>::hessian (const Point<dim> &p,
+                                                      const unsigned int /*component*/) const
+  {
+    SymmetricTensor<2,dim> return_hessian;
+    return_hessian = 0.0;
+    
+    if (dim==2){
+      return_hessian[0][0] = (2.0-12.0*p(0)+12.0*p(0)*p(0)) * std::pow(p(1)*(1.0-p(1)),2);
+      return_hessian[0][1] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * (2.0*p(1)-6.0*std::pow(p(1),2)
+                             +4.0*std::pow(p(1),3));
+      return_hessian[1][1] = (2.0-12.0*p(1)+12.0*p(1)*p(1)) * std::pow(p(0)*(1.0-p(0)),2);
+    } else if (dim==3){
+      return_hessian[0][0] = (2.0-12.0*p(0)+12.0*p(0)*p(0)) * std::pow(p(1)*(1.0-p(1))*p(2)*(1.0-p(2)),2);
+      return_hessian[0][1] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * (2.0*p(1)-6.0*std::pow(p(1),2)
+                             +4.0*std::pow(p(1),3)) * std::pow(p(2)*(1.0-p(2)),2);
+      return_hessian[0][2] = (2.0*p(0)-6.0*std::pow(p(0),2)+4.0*std::pow(p(0),3)) * (2.0*p(2)-6.0*std::pow(p(2),2)
+                             +4.0*std::pow(p(2),3)) * std::pow(p(1)*(1.0-p(1)),2);
+      return_hessian[1][1] = (2.0-12.0*p(1)+12.0*p(1)*p(1)) * std::pow(p(0)*(1.0-p(0))*p(2)*(1.0-p(2)),2);
+      return_hessian[1][2] = (2.0*p(1)-6.0*std::pow(p(1),2)+4.0*std::pow(p(1),3)) * (2.0*p(2)-6.0*std::pow(p(2),2)
+                             +4.0*std::pow(p(2),3)) * std::pow(p(0)*(1.0-p(0)),2);
+      return_hessian[2][2] = (2.0-12.0*p(2)+12.0*p(2)*p(2)) * std::pow(p(0)*(1.0-p(0))*p(1)*(1.0-p(1)),2);
+    }
+
+    return return_hessian;
+  }
+
+
+
+  // @sect3{Implementation of the <code>BiLaplacianLDGLift</code> class}
+
+  // @sect4{BiLaplacianLDGLift::BiLaplacianLDGLift}
+
+  // In the constructor, we set the polynomial degree of the two finite element
+  // spaces, we associate the corresponding DoF handlers to the triangulation,
+  // and we set the two penalty coefficients.
+  template <int dim>
+  BiLaplacianLDGLift<dim>::BiLaplacianLDGLift (const unsigned int fe_degree,double penalty_jump_grad, double penalty_jump_val):
+    fe(fe_degree),
+    dof_handler(triangulation),
+    fe_lift(FE_DGQ<dim>(fe_degree),dim*dim),
+    dof_handler_lift(triangulation),
+    penalty_jump_grad(penalty_jump_grad),
+    penalty_jump_val(penalty_jump_val)
+  {}
+
+
+
+  // @sect4{BiLaplacianLDGLift::make_grid}
+
+  // To build a mesh for $\Omega=(0,1)^d$, we simply call the function
+  // <code>GridGenerator::hyper_cube</code> and then refine it using 
+  // <code>refine_global</code>. The number of refinements is hard-coded
+  // here.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::make_grid()
+  {
+    std::cout << "Building the mesh............." << std::endl;
+
+    GridGenerator::hyper_cube(triangulation,0.0,1.0);
+
+    triangulation.refine_global(3);
+  
+    std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl;
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::setup_system}
+
+  // In the following function, we set up the degrees of freedom, the sparsity pattern,
+  // the size of the matrix $A$, and the size of the solution and right-hand side vectors
+  // $\boldsymbol{U}$ and $\boldsymbol{F}$. For the sparsity pattern, we cannot directly
+  // use the function <code>DoFTools::make_flux_sparsity_pattern</code> (as we would do for
+  // instance for the SIPG method) because we need to take into account the interactions
+  // of a neighboring cell with another neighboring cell as described in the introduction.
+  // The extended sparsity pattern is build by iterating over all the active cells. For
+  // the current cell, we collect all its degrees of freedom as well as the degrees of
+  // freedom of all its neighboring cells, and then couple everything with everything.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+    dof_handler_lift.distribute_dofs(fe_lift);
+    
+    std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl;
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(),dof_handler.n_dofs());
+
+    const auto dofs_per_cell = fe.dofs_per_cell;
+
+    for (const auto cell : dof_handler.active_cell_iterators()){
+
+      std::vector<types::global_dof_index> dofs(dofs_per_cell);
+      cell->get_dof_indices(dofs);
+
+      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f){
+        if (!cell->face(f)->at_boundary()){
+          const auto neighbor_cell = cell->neighbor(f);
+
+          std::vector<types::global_dof_index> tmp(dofs_per_cell);
+          neighbor_cell->get_dof_indices(tmp);
+
+          dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp));
+        }
+      }
+
+      for (const auto i : dofs){
+        for (const auto j : dofs){
+          dsp.add(i, j);
+          dsp.add(j, i);
+        }
+      }
+    }
+
+    sparsity_pattern.copy_from(dsp);
+
+    std::ofstream out("sparsity_pattern.svg");
+    sparsity_pattern.print_svg(out);
+  
+    matrix.reinit(sparsity_pattern);
+    rhs.reinit(dof_handler.n_dofs());
+
+    solution.reinit(dof_handler.n_dofs());
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::assemble_system}
+
+  // This function simply call the two functions responsible
+  // for the assembly of the matrix and the right-hand side.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::assemble_system()
+  {
+    std::cout << "Assembling the system............." << std::endl;
+
+    assemble_matrix();
+    assemble_rhs();
+  
+    std::cout << "Done. " << std::endl;
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::assemble_matrix}
+
+  // This function assemble the matrix $A$ whose entries are defined
+  // by $A_{ij}=A_h(\varphi_j,\varphi_i)$ which involves the product of
+  // discrete Hessians and the penalty terms.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::assemble_matrix()
+  {
+    matrix = 0;
+
+    QGauss<dim> quad(fe.degree+1);
+    QGauss<dim-1> quad_face(fe.degree+1);
+
+    const unsigned int   n_q_points = quad.size();
+    const unsigned int   n_q_points_face = quad_face.size();
+
+    FEValues<dim> fe_values (fe, quad, update_hessians |
+                                       update_JxW_values);
+
+    FEFaceValues<dim> fe_face (fe, quad_face, update_values |
+                                              update_gradients |
+                                              update_normal_vectors);
+
+    FEFaceValues<dim> fe_face_neighbor (fe, quad_face, update_values |
+                                                       update_gradients |
+                                                       update_normal_vectors);
+
+    const unsigned int   n_dofs = fe_values.dofs_per_cell;
+
+    std::vector<types::global_dof_index> local_dof_indices(n_dofs),
+                                         local_dof_indices_neighbor (n_dofs),
+                                         local_dof_indices_neighbor_2 (n_dofs);
+
+    // As indicated in the introduction, the following matrices are used for
+    // the contributions of the products of the discrete Hessians.
+    FullMatrix<double> stiffness_matrix_cc (n_dofs,n_dofs); // interactions cell / cell
+    FullMatrix<double> stiffness_matrix_cn (n_dofs,n_dofs); // interactions cell / neighboor
+    FullMatrix<double> stiffness_matrix_nc (n_dofs,n_dofs); // interactions neighboor / cell
+    FullMatrix<double> stiffness_matrix_nn (n_dofs,n_dofs); // interactions neighboor / neighboor
+    FullMatrix<double> stiffness_matrix_n1n2 (n_dofs,n_dofs); // interactions neighboor_1 / neighboor_2
+    FullMatrix<double> stiffness_matrix_n2n1 (n_dofs,n_dofs); // interactions neighboor_2 / neighboor_1
+
+    // The following matrices are used for the contributions of the two
+    // penalty terms.
+    FullMatrix<double>   ip_matrix_cc (n_dofs,n_dofs); // interactions cell / cell
+    FullMatrix<double>   ip_matrix_cn (n_dofs,n_dofs); // interactions cell / neighboor
+    FullMatrix<double>   ip_matrix_nc (n_dofs,n_dofs); // interactions neighboor / cell
+    FullMatrix<double>   ip_matrix_nn (n_dofs,n_dofs); // interactions neighboor / neighboor
+
+    std::vector<std::vector<Tensor<2,dim>>> discrete_hessians (n_dofs, std::vector<Tensor<2,dim>>(n_q_points) );
+    std::vector<std::vector<std::vector<Tensor<2,dim>>>> discrete_hessians_neigh (GeometryInfo<dim>::faces_per_cell, discrete_hessians);
+
+    Tensor<2,dim> H_i,H_j;
+    Tensor<2,dim> H_i_neigh,H_j_neigh;
+    Tensor<2,dim> H_i_neigh2,H_j_neigh2;
+
+    double mesh_inv,mesh3_inv;
+    bool at_boundary,at_boundary_2;
+    unsigned int face_no_neighbor = 0;
+
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+  
+    typename DoFHandler<dim>::active_cell_iterator neighbor_cell,neighbor_cell_2;
+
+    typename DoFHandler<dim>::active_cell_iterator cell_lift = dof_handler_lift.begin_active();
+
+    for (; cell != endc; ++cell, ++cell_lift){
+
+      fe_values.reinit(cell);
+      cell->get_dof_indices (local_dof_indices); 
+
+      // We now compute all the discrete Hessians that are not vanishing
+      // on the current cell, i.e., the discrete Hessian of all the basis
+      // functions with support on the current cell or on one of its neighbors.
+      compute_discrete_hessians(cell,cell_lift,
+                                discrete_hessians,discrete_hessians_neigh);
+
+      // First, we compute and add the interactions of the degrees of freedom
+      // of the current cell.
+      stiffness_matrix_cc = 0;
+      for (unsigned int q=0; q<n_q_points; ++q){
+        const double dx = fe_values.JxW(q);
+        
+        for (unsigned int i=0; i<n_dofs; ++i){
+          for (unsigned int j=0; j<n_dofs; ++j){
+
+            H_i = discrete_hessians[i][q];
+            H_j = discrete_hessians[j][q];
+            
+            stiffness_matrix_cc(i,j) += dx * scalar_product(H_j,H_i);
+
+          } 
+        }
+      }
+          
+      for (unsigned int i=0; i<n_dofs; ++i){
+        for (unsigned int j=0; j<n_dofs; ++j){
+          matrix(local_dof_indices[i],local_dof_indices[j]) += stiffness_matrix_cc(i,j);
+        }
+      }
+
+      // Next, we compute and add the interactions of the degrees of freedom of the current
+      // cell with those of its neighbors. Note that the interactions of the degrees of
+      // freedom of a neighbor with those of the same neighbor are included here.
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
+        const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+                
+        at_boundary = face->at_boundary();
+
+       if (!at_boundary){ // nothing to be done if boundary face (the liftings of the Dirichlet BCs are accounted for in the assembly of the RHS; in fact, nothing to be done in this program since we prescribe homogeneous BCs)
+
+          neighbor_cell =cell->neighbor(face_no);
+          neighbor_cell->get_dof_indices (local_dof_indices_neighbor);
+
+          stiffness_matrix_cn=0;
+          stiffness_matrix_nc=0;
+          stiffness_matrix_nn=0;
+          for (unsigned int q=0; q<n_q_points; ++q){
+            const double dx = fe_values.JxW(q);
+        
+            for (unsigned int i=0; i<n_dofs; ++i){
+              for (unsigned int j=0; j<n_dofs; ++j){
+
+                H_i = discrete_hessians[i][q];
+                H_j = discrete_hessians[j][q];
+              
+                H_i_neigh = discrete_hessians_neigh[face_no][i][q];
+                H_j_neigh = discrete_hessians_neigh[face_no][j][q];
+              
+                stiffness_matrix_cn(i,j) += dx * scalar_product(H_j_neigh,H_i);
+                stiffness_matrix_nc(i,j) += dx * scalar_product(H_j,H_i_neigh);
+                stiffness_matrix_nn(i,j) += dx * scalar_product(H_j_neigh,H_i_neigh);
+
+             } 
+           }
+          }
+          
+          for (unsigned int i=0; i<n_dofs; ++i){
+            for (unsigned int j=0; j<n_dofs; ++j){
+              matrix(local_dof_indices[i],local_dof_indices_neighbor[j]) += stiffness_matrix_cn(i,j);
+              matrix(local_dof_indices_neighbor[i],local_dof_indices[j]) += stiffness_matrix_nc(i,j);
+              matrix(local_dof_indices_neighbor[i],local_dof_indices_neighbor[j]) += stiffness_matrix_nn(i,j);
+            }
+          }
+
+        } // boundary check
+      } // for face
+
+      // We now compute and add the interactions of the degrees of freedom of a
+      // neighboring cells with those of another neighboring cell (this is where we
+      // need the extended sparsity pattern).
+      for (unsigned int face_no=0; face_no < GeometryInfo<dim>::faces_per_cell-1; ++face_no){
+        const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+                
+        at_boundary = face->at_boundary();
+
+       if (!at_boundary){ // nothing to be done if boundary face (the liftings of the Dirichlet BCs are accounted for in the assembly of the RHS; in fact, nothing to be done in this program since we prescribe homogeneous BCs)
+
+
+          for (unsigned int face_no_2=face_no+1; face_no_2 < GeometryInfo<dim>::faces_per_cell; ++face_no_2){
+            const typename DoFHandler<dim>::face_iterator face_2=cell->face(face_no_2);
+            at_boundary_2 = face_2->at_boundary();
+
+           if (!at_boundary_2){
+
+             neighbor_cell = cell->neighbor(face_no);
+             neighbor_cell->get_dof_indices (local_dof_indices_neighbor);
+             neighbor_cell_2 = cell->neighbor(face_no_2);
+             neighbor_cell_2->get_dof_indices (local_dof_indices_neighbor_2);
+
+              stiffness_matrix_n1n2=0;
+              stiffness_matrix_n2n1=0;
+
+              for (unsigned int q=0; q<n_q_points; ++q){
+                const double dx = fe_values.JxW(q);
+        
+                for (unsigned int i=0; i<n_dofs; ++i){
+                  for (unsigned int j=0; j<n_dofs; ++j){
+
+                    H_i_neigh = discrete_hessians_neigh[face_no][i][q];
+                    H_j_neigh = discrete_hessians_neigh[face_no][j][q];
+                    
+                    H_i_neigh2 = discrete_hessians_neigh[face_no_2][i][q];
+                    H_j_neigh2 = discrete_hessians_neigh[face_no_2][j][q];
+                    
+                    stiffness_matrix_n1n2(i,j) += dx * scalar_product(H_j_neigh2,H_i_neigh);
+                    stiffness_matrix_n2n1(i,j) += dx * scalar_product(H_j_neigh,H_i_neigh2);
+
+                 } 
+               }
+              }
+          
+              for (unsigned int i=0; i<n_dofs; ++i){
+                for (unsigned int j=0; j<n_dofs; ++j){
+                  matrix(local_dof_indices_neighbor[i],local_dof_indices_neighbor_2[j]) += stiffness_matrix_n1n2(i,j);
+                  matrix(local_dof_indices_neighbor_2[i],local_dof_indices_neighbor[j]) += stiffness_matrix_n2n1(i,j);
+                }
+              }
+            } // boundary check face_2
+          } // for face_2
+        } // boundary check face_1
+      } // for face_1
+
+
+      // Finally, we compute and add the two penalty terms.
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
+        const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+
+        mesh_inv = 1.0/face->diameter(); // h_e^{-1}
+        mesh3_inv = 1.0/std::pow(face->diameter(),3); // Ä¥_e^{-3}
+         
+        fe_face.reinit(cell,face_no);
+
+        ip_matrix_cc = 0; // filled in any case (boundary or interior face)
+        
+        at_boundary = face->at_boundary();
+        if (at_boundary){
+
+          for (unsigned int q=0; q<n_q_points_face; ++q){
+            const double dx = fe_face.JxW(q);
+
+            for (unsigned int i=0; i<n_dofs; ++i){
+              for (unsigned int j=0; j<n_dofs; ++j){
+                ip_matrix_cc(i,j) += penalty_jump_grad * mesh_inv * dx * fe_face.shape_grad(j,q) * fe_face.shape_grad(i,q);
+                ip_matrix_cc(i,j) += penalty_jump_val * mesh3_inv * dx * fe_face.shape_value(j,q) * fe_face.shape_value(i,q);
+              }
+            }
+          }
+
+        } else{ // interior face
+   
+          neighbor_cell =cell->neighbor(face_no);        
+          face_no_neighbor = cell->neighbor_of_neighbor (face_no);
+
+          if(neighbor_cell->id().operator<(cell->id())){ //we need to have a global way to compare the cells in order to not calculate the same jump term twice
+            continue; // skip this face (already considered)
+         } else{
+
+            fe_face_neighbor.reinit(neighbor_cell,face_no_neighbor);
+            neighbor_cell->get_dof_indices (local_dof_indices_neighbor);
+                        
+            ip_matrix_cn = 0;
+            ip_matrix_nc = 0;
+            ip_matrix_nn = 0;
+
+            for (unsigned int q=0; q<n_q_points_face; ++q){
+              const double dx = fe_face.JxW(q);
+
+              for (unsigned int i=0; i<n_dofs; ++i){
+                for (unsigned int j=0; j<n_dofs; ++j){
+                  ip_matrix_cc(i,j) += penalty_jump_grad * mesh_inv * dx * fe_face.shape_grad(j,q) * fe_face.shape_grad(i,q);
+                  ip_matrix_cc(i,j) += penalty_jump_val * mesh3_inv * dx * fe_face.shape_value(j,q) * fe_face.shape_value(i,q);
+                  
+                  ip_matrix_cn(i,j) -= penalty_jump_grad * mesh_inv * dx * fe_face_neighbor.shape_grad(j,q) * fe_face.shape_grad(i,q);
+                  ip_matrix_cn(i,j) -= penalty_jump_val * mesh3_inv * dx * fe_face_neighbor.shape_value(j,q) * fe_face.shape_value(i,q);
+                  
+                  ip_matrix_nc(i,j) -= penalty_jump_grad * mesh_inv * dx * fe_face.shape_grad(j,q) * fe_face_neighbor.shape_grad(i,q);
+                  ip_matrix_nc(i,j) -= penalty_jump_val * mesh3_inv * dx * fe_face.shape_value(j,q) * fe_face_neighbor.shape_value(i,q);
+                  
+                  ip_matrix_nn(i,j) += penalty_jump_grad * mesh_inv * dx * fe_face_neighbor.shape_grad(j,q) * fe_face_neighbor.shape_grad(i,q);
+                  ip_matrix_nn(i,j) += penalty_jump_val * mesh3_inv * dx * fe_face_neighbor.shape_value(j,q) * fe_face_neighbor.shape_value(i,q);
+                }
+              }
+            }
+          } // face not visited yet
+
+        } // boundary check
+
+        for (unsigned int i=0; i<n_dofs; ++i){
+          for (unsigned int j=0; j<n_dofs; ++j){
+            matrix(local_dof_indices[i],local_dof_indices[j]) += ip_matrix_cc(i,j);
+          }
+        }
+
+        if (!at_boundary){
+
+          for (unsigned int i=0; i<n_dofs; ++i){
+            for (unsigned int j=0; j<n_dofs; ++j){
+              matrix(local_dof_indices[i],local_dof_indices_neighbor[j]) += ip_matrix_cn(i,j);
+              matrix(local_dof_indices_neighbor[i],local_dof_indices[j]) += ip_matrix_nc(i,j);
+              matrix(local_dof_indices_neighbor[i],local_dof_indices_neighbor[j]) += ip_matrix_nn(i,j);
+            }
+          }
+        }
+
+      } // for face
+    } // for cell
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::assemble_rhs}
+
+  // This function assemble the right-hand side of the system. Since we consider
+  // homogeneous Dirichlet boundary conditions, imposed weakly in the bilinear
+  // form using the Nitsche approach, it only involves the contribution of the
+  // forcing term $\int_{\Omega}fv_h$.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::assemble_rhs()
+  {
+    rhs = 0;
+    QGauss<dim> quad(fe.degree+1);  
+    FEValues<dim> fe_values(fe, quad, update_values |
+                                      update_quadrature_points |
+                                      update_JxW_values);
+
+    const unsigned int n_dofs = fe_values.dofs_per_cell;
+    const unsigned int n_quad_pts = quad.size();
+
+    const RightHandSide<dim> right_hand_side;
+
+    Vector<double> local_rhs(n_dofs);
+    std::vector<types::global_dof_index> local_dof_indices(n_dofs);
+
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+
+    for (; cell!=endc; ++cell){
+
+     fe_values.reinit(cell);
+      cell->get_dof_indices(local_dof_indices);
+
+      local_rhs = 0;
+      for (unsigned int q=0; q<n_quad_pts; ++q){
+        const double dx = fe_values.JxW(q);
+
+        for (unsigned int i=0; i<n_dofs; ++i){
+          local_rhs(i) += dx * right_hand_side.value(fe_values.quadrature_point(q)) * fe_values.shape_value(i,q);
+        }
+      }
+
+      for (unsigned int i=0; i<n_dofs; ++i){
+        rhs(local_dof_indices[i]) += local_rhs(i);
+      }
+    }
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::solve}
+
+  // To solve the linear system $A\boldsymbol{U}=\boldsymbol{F}$,
+  // we proceed as in step-74 and use a direct method. We could
+  // as well use an iterative method, for instance the conjugate
+  // gradient method as in step-3.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::solve()
+  {
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(matrix);
+    A_direct.vmult(solution, rhs);
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::compute_errors}
+
+  // This function computes the discrete $H^2$, $H^1$ and $L^2$ norms of
+  // the error $u-u_h$, where $u$ is the exact solution and $u_h$ is
+  // the approximate solution. See the introduction for the definition
+  // of the norms.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::compute_errors()
+  {  
+
+    double error_H2 = 0; // sqrt( ||D_h^2(u-u_h)||_{L^2(Omega)}^2 + ||h^{-1/2}[grad_h(u-u_h)]||_{L^2(Sigma)}^2 + ||h^{-3/2}[u-u_h]||_{L^2(Sigma)}^2 )
+    double error_H1 = 0; // sqrt( ||grad_h(u-u_h)||_{L^2(Omega)}^2 + ||h^{-1/2}[u-u_h]||_{L^2(Sigma)}^2 )
+    double error_L2 = 0; // ||u-u_h||_{L^2(Omega)}
+    
+    QGauss<dim> quad(fe.degree+1);
+    QGauss<dim-1> quad_face(fe.degree+1);
+    
+    FEValues<dim> fe_values (fe, quad, update_values |
+                                       update_gradients |
+                                       update_hessians |
+                                       update_quadrature_points |
+                                       update_JxW_values);
+    
+    FEFaceValues<dim> fe_face (fe, quad_face, update_values |
+                                              update_gradients |
+                                              update_quadrature_points |
+                                              update_JxW_values);
+    
+    FEFaceValues<dim> fe_face_neighbor (fe, quad_face, update_values |
+                                                       update_gradients);
+    
+    const unsigned int   n_q_points      = quad.size();
+    const unsigned int   n_q_points_face = quad_face.size();
+    
+    // We introduce some variables for the exact solution
+    const ExactSolution<dim> u_exact;
+    double u_exact_q;
+    Tensor<1,dim> u_exact_grad_q;
+    
+    // and for the approximate solution
+    std::vector<double>        solution_values_cell(n_q_points);
+    std::vector<Tensor<1,dim>> solution_gradients_cell(n_q_points);
+    std::vector<Tensor<2,dim>> solution_hessians_cell(n_q_points);
+
+    std::vector<double>        solution_values(n_q_points_face);
+    std::vector<double>        solution_values_neigh(n_q_points_face);    
+    std::vector<Tensor<1,dim>> solution_gradients(n_q_points_face);
+    std::vector<Tensor<1,dim>> solution_gradients_neigh(n_q_points_face);
+    
+    double mesh_inv;
+    double mesh3_inv;
+    bool at_boundary;
+    
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    
+    typename DoFHandler<dim>::active_cell_iterator neighbor_cell;
+    unsigned int face_no_neighbor = 0;
+    
+    for (; cell!=endc; ++cell){
+        
+      fe_values.reinit (cell);
+        
+      fe_values.get_function_values(solution,solution_values_cell);
+      fe_values.get_function_gradients(solution,solution_gradients_cell);
+      fe_values.get_function_hessians(solution,solution_hessians_cell);
+                
+      // We first add the <i>bulk</i> terms.          
+      for (unsigned int q=0; q<n_q_points; ++q){
+        const double dx = fe_values.JxW(q);
+            
+        error_H2 += dx * (u_exact.hessian(fe_values.quadrature_point(q))-solution_hessians_cell[q]).norm_square();
+        error_H1 += dx * (u_exact.gradient(fe_values.quadrature_point(q))-solution_gradients_cell[q]).norm_square();
+        error_L2 += dx * std::pow(u_exact.value(fe_values.quadrature_point(q))-solution_values_cell[q],2);
+      } // for quad
+        
+      // We then add the face contributions.
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
+        const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+            
+        mesh_inv = 1.0/face->diameter(); // h^{-1}
+        mesh3_inv = 1.0/std::pow(face->diameter(),3); // h^{-3}
+            
+        fe_face.reinit(cell,face_no);
+            
+        fe_face.get_function_values(solution,solution_values);
+        fe_face.get_function_gradients(solution,solution_gradients);
+            
+        at_boundary = face->at_boundary();
+        if (at_boundary){
+                
+          for (unsigned int q=0; q<n_q_points_face; ++q){
+            const double dx = fe_face.JxW(q);
+            u_exact_q = u_exact.value(fe_face.quadrature_point(q));
+            u_exact_grad_q = u_exact.gradient(fe_face.quadrature_point(q));
+                    
+            error_H2 += dx * mesh_inv * (u_exact_grad_q-solution_gradients[q]).norm_square();
+            error_H2 += dx * mesh3_inv * std::pow(u_exact_q-solution_values[q],2);
+            error_H1 += dx * mesh_inv * std::pow(u_exact_q-solution_values[q],2);
+          }
+                
+        } else{ // interior face
+                
+          neighbor_cell =cell->neighbor(face_no);
+          face_no_neighbor = cell->neighbor_of_neighbor (face_no);
+                
+          if(neighbor_cell->id().operator<(cell->id())){ // we need to have a global way to compare the cells in order to not calculate the same jump term twice
+            continue; // skip this face (already considered)
+          } else{
+                    
+            fe_face_neighbor.reinit(neighbor_cell,face_no_neighbor);
+                    
+            fe_face.get_function_values(solution,solution_values);
+            fe_face_neighbor.get_function_values(solution,solution_values_neigh);
+            fe_face.get_function_gradients(solution,solution_gradients);
+            fe_face_neighbor.get_function_gradients(solution,solution_gradients_neigh);
+                    
+            for (unsigned int q=0; q<n_q_points_face; ++q){
+              const double dx = fe_face.JxW(q);
+                        
+              // To compute the jump term, we use the fact that $\jump{u}=0$ and
+              // $\jump{\nabla u}=\mathbf{0}$ since $u\in H^2(\Omega)$.
+              error_H2 += dx * mesh_inv * (solution_gradients_neigh[q]-solution_gradients[q]).norm_square();
+              error_H2 += dx * mesh3_inv * std::pow(solution_values_neigh[q]-solution_values[q],2);
+              error_H1 += dx * mesh_inv * std::pow(solution_values_neigh[q]-solution_values[q],2);
+            }
+          } // face not visited yet
+                
+        } // boundary check
+            
+      } // for face
+        
+    } // for cell
+    
+    error_H2 = std::sqrt(error_H2);
+    error_H1 = std::sqrt(error_H1);
+    error_L2 = std::sqrt(error_L2);
+
+    std::cout << "DG H2 norm of the error: " << error_H2 << std::endl;
+    std::cout << "DG H1 norm of the error: " << error_H1 << std::endl;
+    std::cout << "   L2 norm of the error: " << error_L2 << std::endl;     
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::output_results}
+
+  // This function, which write the solution to a vtk file,
+  // is copied from step-3.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::output_results() const
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "solution"); 
+    data_out.build_patches();
+
+    std::ofstream output ("solution.vtk");
+    data_out.write_vtk (output);
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::assemble_local_matrix}
+
+  // As already mentioned above, this function is used to assemble
+  // the (local) mass matrices needed for the computations of the
+  // lifting terms. We reiterate that only the basis functions with
+  // support on the current cell are accounting for.
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::assemble_local_matrix(const FEValues<dim> &fe_values_lift, const unsigned int n_q_points, FullMatrix<double> &local_matrix)
+  {
+    const FEValuesExtractors::Tensor<2> tau_ext(0);
+
+    const unsigned int n_dofs = fe_values_lift.dofs_per_cell;
+
+    local_matrix = 0;
+    for (unsigned int q=0; q<n_q_points; ++q){
+
+      const double dx = fe_values_lift.JxW(q);
+
+      for (unsigned int m=0; m<n_dofs; ++m){
+        for (unsigned int n=0; n<n_dofs; ++n){
+          local_matrix(m,n) += dx * scalar_product(fe_values_lift[tau_ext].value(n,q),fe_values_lift[tau_ext].value(m,q));
+        }
+      }
+    }
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::compute_discrete_hessians}
+
+  // This function is the main novelty of this program. It computes the discrete
+  // Hessian $H_h(\varphi)$ for all the basis functions $\varphi$ of $\mathbb{V}_h$
+  // supported on the current cell and those supported on a neighboring cell. The
+  // first two arguments are inputs indicating the current cell (both refer to the
+  // same cell but are attached to different DoF Handlers), while the last two
+  // arguments are output variables that are filled-in in this function. 
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::compute_discrete_hessians(const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                                          const typename DoFHandler<dim>::active_cell_iterator &cell_lift,
+                                                          std::vector<std::vector<Tensor<2,dim>>> &discrete_hessians,
+                                                          std::vector<std::vector<std::vector<Tensor<2,dim>>>> &discrete_hessians_neigh)
+  {
+    QGauss<dim> quad(fe.degree+1);
+    QGauss<dim-1> quad_face(fe.degree+1);
+
+    const unsigned int   n_q_points = quad.size();
+    const unsigned int   n_q_points_face = quad_face.size();
+
+    // The information we need from the basis functions of
+    // $\mathbb{V}_h$: <code>fe_values</code> is needed to add
+    // the broken Hessian part of the discrete Hessian, while
+    // <code>fe_face</code> and <code>fe_face_neighbor</code>
+    // are used to compute the right-hand sides for the local
+    // problems.
+    FEValues<dim> fe_values (fe, quad, update_hessians |
+                                       update_JxW_values);
+
+    FEFaceValues<dim> fe_face (fe, quad_face, update_values |
+                                              update_gradients |
+                                              update_normal_vectors);
+
+    FEFaceValues<dim> fe_face_neighbor (fe, quad_face, update_values |
+                                                       update_gradients |
+                                                       update_normal_vectors);
+
+    const unsigned int n_dofs = fe_values.dofs_per_cell;
+
+    typename DoFHandler<2,dim>::active_cell_iterator neighbor_cell;
+    unsigned int face_no_neighbor = 0;
+
+    // The information needed from the basis functions
+    // of the finite element space for the lifting terms:
+    // <code>fe_values_lift</code> is used for the (local)
+    // mass matrix (see $\boldsymbol{M}_c$ in the introduction),
+    // while <code>fe_face_lift</code> is used to compute the
+    // right-hand sides (see $\boldsymbol{G}_c$ for $b_e$).
+    FEValues<dim> fe_values_lift (fe_lift, quad, update_values |
+                                                 update_JxW_values);
+
+    FEFaceValues<dim> fe_face_lift (fe_lift, quad_face, update_values |
+                                                        update_gradients |
+                                                        update_JxW_values);
+
+    const FEValuesExtractors::Tensor<2> tau_ext(0);
+    
+    const unsigned int n_dofs_lift = fe_values_lift.dofs_per_cell;
+    FullMatrix<double> local_matrix_lift (n_dofs_lift,n_dofs_lift);
+
+    Vector<double> local_rhs_re(n_dofs_lift), local_rhs_be(n_dofs_lift),
+                   coeffs_re(n_dofs_lift), coeffs_be(n_dofs_lift),
+                   coeffs_tmp(n_dofs_lift);
+     
+    SolverControl solver_control(1000, 1e-12);
+    SolverCG<> solver(solver_control);
+    
+    bool at_boundary;    
+    double factor_avg; // 0.5 for interior faces, 1.0 for boundary faces
+    
+    fe_values.reinit(cell);
+    fe_values_lift.reinit(cell_lift);  
+    // We start by assembling the (local) mass matrix used for the computation
+    // of the lifting terms $r_e$ and $b_e$.
+    assemble_local_matrix(fe_values_lift,n_q_points,local_matrix_lift);
+
+    for (unsigned int i=0; i<n_dofs; ++i){
+      for (unsigned int q=0; q<n_q_points; ++q){
+        discrete_hessians[i][q]=0;
+
+        for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
+          discrete_hessians_neigh[face_no][i][q]=0;
+        }
+
+      }
+    }
+
+    // In this loop, we compute the discrete Hessian at each quadrature point $x_q$
+    // of <code>cell</code> for each basis function supported on <code>cell</code>,
+    // namely we fill-in the variable <code>discrete_hessians[i][q]</code>.
+    // For the lifting terms, we need to add the contribution of all the faces of
+    // <code>cell</code>.
+    for (unsigned int i=0; i<n_dofs; ++i){
+
+      coeffs_re=0; coeffs_be=0;
+
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
+        const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+            
+        at_boundary = face->at_boundary();
+        
+        // Recall that by convention, the average of a function accross a boundary
+        // face $e$ reduces to the trace of the function on the only element
+        // adjacent to $e$, namely there is no factor $\frac{1}{2}$. We distinguish
+        // between the two cases (the current face lies in the interior or on
+        // the boundary of the domain) using the variable <code>factor_avg</code>.
+        factor_avg = 0.5;
+        if (at_boundary){
+          factor_avg = 1.0;
+        }
+        
+        fe_face.reinit(cell,face_no);
+        fe_face_lift.reinit(cell_lift,face_no);
+        
+        local_rhs_re=0;
+        for (unsigned int q=0; q<n_q_points_face; ++q){
+          const double dx = fe_face_lift.JxW(q);
+          const Tensor<1,dim> normal = fe_face.normal_vector(q); // same as fe_face_lift.normal_vector(q)
+          
+          for (unsigned int m=0; m<n_dofs_lift; ++m){
+            local_rhs_re(m) += factor_avg * dx * (fe_face_lift[tau_ext].value(m,q)*normal) * fe_face.shape_grad(i,q);
+          }
+        }
+
+        // Here, <code>local_rhs_be(m)</code> corresponds to $G_m$ introduced in the
+        // comments about the implementation of the lifting $b_e$ in the case
+        // $\varphi=\varphi^c$.
+        local_rhs_be=0;
+        for (unsigned int q=0; q<n_q_points_face; ++q){
+          const double dx = fe_face_lift.JxW(q);
+          const Tensor<1,dim> normal = fe_face.normal_vector(q); // same as fe_face_lift.normal_vector(q)
+
+          for (unsigned int m=0; m<n_dofs_lift; ++m){
+            local_rhs_be(m) += factor_avg * dx * fe_face_lift[tau_ext].divergence(m,q)*normal * fe_face.shape_value(i,q);
+          }
+        }
+        
+        coeffs_tmp=0;
+        solver.solve(local_matrix_lift,coeffs_tmp,local_rhs_re,PreconditionIdentity());
+        coeffs_re += coeffs_tmp;
+        
+        coeffs_tmp=0;
+        solver.solve(local_matrix_lift,coeffs_tmp,local_rhs_be,PreconditionIdentity());
+        coeffs_be += coeffs_tmp;
+
+      } // for face
+
+      for (unsigned int q=0; q<n_q_points; ++q){
+        discrete_hessians[i][q] += fe_values.shape_hessian(i,q);
+
+       for (unsigned int m=0; m<n_dofs_lift; ++m){
+         discrete_hessians[i][q] -= coeffs_re[m]*fe_values_lift[tau_ext].value(m,q);
+        }
+
+       for (unsigned int m=0; m<n_dofs_lift; ++m){
+         discrete_hessians[i][q] += coeffs_be[m]*fe_values_lift[tau_ext].value(m,q);
+        }
+      }
+    } // for dof i
+
+
+
+    // In this loop, we compute the discrete Hessian at each quadrature point
+    // $x_q$ of <code>cell</code> for each basis function supported on a neighboring
+    // <code>neighbor_cell</code> of <code>cell</code>, namely we fill-in the
+    // variable <code>discrete_hessians_neigh[face_no][i][q]</code>.
+    // For the lifting terms, we only need to add the contribution of the
+    // face adjecent to <code>cell</code> and <code>neighbor_cell</code>.    
+    for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no){
+
+      const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+      
+      at_boundary = face->at_boundary();
+
+      if (!at_boundary){ // for non-homogeneous Dirichlet BCs, we would need to compute the lifting of the prescribed BC (see Section Possible Extensions for more details)
+
+        neighbor_cell =cell->neighbor(face_no);                    
+        face_no_neighbor = cell->neighbor_of_neighbor (face_no);
+        fe_face_neighbor.reinit(neighbor_cell,face_no_neighbor);
+
+        for (unsigned int i=0; i<n_dofs; ++i){
+               
+          coeffs_re=0; coeffs_be=0;
+          
+          fe_face_lift.reinit(cell_lift,face_no);
+          
+          local_rhs_re=0;
+          for (unsigned int q=0; q<n_q_points_face; ++q){
+            const double dx = fe_face_lift.JxW(q);
+            const Tensor<1,dim> normal = fe_face_neighbor.normal_vector(q);
+
+            for (unsigned int m=0; m<n_dofs_lift; ++m){
+              local_rhs_re(m) += 0.5 * dx * (fe_face_lift[tau_ext].value(m,q)*normal) * fe_face_neighbor.shape_grad(i,q);
+            }
+          }
+
+          // Here, <code>local_rhs_be(m)</code> corresponds to $G_m$ introduced in
+          // the comments about the implementation of the lifting $b_e$ in the case
+          // $\varphi=\varphi^n$.
+          local_rhs_be=0;
+          for (unsigned int q=0; q<n_q_points_face; ++q){
+            const double dx = fe_face_lift.JxW(q);
+            const Tensor<1,dim> normal = fe_face_neighbor.normal_vector(q);  
+
+            for (unsigned int m=0; m<n_dofs_lift; ++m){
+              local_rhs_be(m) += 0.5 * dx * fe_face_lift[tau_ext].divergence(m,q)*normal * fe_face_neighbor.shape_value(i,q);
+            }
+          }
+          
+          solver.solve(local_matrix_lift,coeffs_re,local_rhs_re,PreconditionIdentity());
+          solver.solve(local_matrix_lift,coeffs_be,local_rhs_be,PreconditionIdentity());
+
+          for (unsigned int q=0; q<n_q_points; ++q){
+
+            for (unsigned int m=0; m<n_dofs_lift; ++m){
+              discrete_hessians_neigh[face_no][i][q] -= coeffs_re[m]*fe_values_lift[tau_ext].value(m,q);
+            }
+            
+            for (unsigned int m=0; m<n_dofs_lift; ++m){
+              discrete_hessians_neigh[face_no][i][q] += coeffs_be[m]*fe_values_lift[tau_ext].value(m,q);
+            }
+
+          }
+
+        } // for dof i
+      } // boundary check
+    } // for face
+  }
+
+
+
+  // @sect4{BiLaplacianLDGLift::run}
+  template <int dim>
+  void BiLaplacianLDGLift<dim>::run()
+  {
+    make_grid();
+    
+    setup_system();
+    assemble_system();
+
+    solve();
+
+    compute_errors();
+    output_results();
+  }
+
+}  // namespace Step82
+
+
+
+// @sect3{The <code>main</code> function}
+
+// The is the <code>main</code> function. We define here the polynomial degree
+// for the two finite element spaces (for the solution and the two liftings) and
+// the two penalty coefficients. We can also change the dimension to run the
+// code in 3D. 
+int main()
+{
+  int degree=2;    // FE degree for u_h and the two lifting terms
+
+  double penalty_grad=1.0; // penalty coefficient for the jump of the gradients
+  double penalty_val=1.0; // penalty coefficient for the jump of the values
+
+  Step82::BiLaplacianLDGLift<2> problem(degree,penalty_grad,penalty_val);
+  
+  problem.run();
+  
+  return 0;
+}

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.