<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
+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
+where $\Omega\subset\mathbb{R}^d$ $(d=2,3)$ is an open bounded Lipschitz domain and $f\in L^2(\Omega)$. This is the same problem we have already considered in step-47, but we will here take a different approach towards solving it: Rather than using continuous finite elements and enforcing the necessary continuity via penalty terms, we will use lifting operators to deal with the fact that the usual finite element spaces are not twice continuous.
+
+The weak formulation of this problem reads as follows: 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.
+where $D^2v$ denotes the Hessian of $v$ and $H_0^2(\Omega)\dealcoloneq\{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
+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 (discontinuous) 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\},
+\mathbb{V}_h\dealcoloneq\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$.
+where $F_{K}$ is the map from the reference element $\hat{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.
+They are missing inter-element information.
This leads to the introductions of jump and average operators:
@f[
-\jump{v_h}|_e :=
+\jump{v_h}|_e \dealcoloneq
\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 :=
+\end{array}\right. \quad \mbox{and} \quad \average{v_h}|_e \dealcoloneq
\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).
+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). These are the same operators that we have previously used not only in step-47, but also in other tutorials related to discontinuous Galerkin methods (e.g., step-12).
-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,
+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$; that is, something that will take the role of $D^2 v$ in the definition of the weak formulation above when moving from the continuous to the discrete formulation. 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}
+r_e\left(\boldsymbol{\phi}\right) \in [\mathbb{V}_h]^{d\times d}: \,
+\int_{\Omega} \tau_h : r_e\left(\boldsymbol{\phi}\right) = \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}.
+\int_{\Omega} \tau_h : b_e(\phi) = \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.
+We have ${\rm supp}\,(r_e\left(\boldsymbol{\phi}\right))={\rm supp}\,(b_e(\phi))=\omega_e$, where $\omega_e$ denotes the patch of (one or two) elements having $e$ as 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}).
+H_h(v_h) \dealcoloneq D_h^2 v_h -R_h(\jump{\nabla_h v_h})+B_h(\jump{v_h}) \dealcoloneq D_h^2 v_h - \sum_{e\in\mathcal{E}_h}r_e\left(\jump{\nabla_h v_h}\right)+\sum_{e\in\mathcal{E}_h}b_e\left(\jump{v_h}\right).
@f]
@note
<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
+The discrete Hessian $H_h$ is designed such that it weakly converges to the continuous Hessian $D^2$, 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]
<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,
+A_h(u_h,v_h)\dealcoloneq 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.
+ a_h(u_h,v_h) & \dealcoloneq \int_{\Omega}H_h(u_h):H_h(v_h), \\
+ j_h(u_h,v_h) & \dealcoloneq \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) & \dealcoloneq \int_{\Omega}fv_h.
@f}
Here, $\gamma_0,\gamma_1>0$ are penalty parameters.
@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.
+A_{ij}\dealcoloneq A_h(\varphi_j,\varphi_i) \quad \text{and} \quad F_i\dealcoloneq 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:
+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. This is because the lifting operators in $H_h$ extend shape functions defined on one cell to the neighboring cell where it may overlap with the lifted shape functions from a neighbor of the neighbor. 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
+ \|v_h\|_{H_h^2(\Omega)}^2\dealcoloneq\|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.
+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. (See also the discussions about penalty parameters in the step-39, step-47, and step-74 programs about penalty parameters.)
<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.
+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\left(\jump{\varphi_i}\right)=\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>
<tr>
<td align="left">dofs $K_c$ $\leftrightarrow$ dofs $K_{n_k}$</td>
<td> </td>
- <td align="left">(stored in <code>stiffness_matrix_cn</code> and <code>stiffness_matrix_cn</code>)</td>
+ <td align="left">(stored in <code>stiffness_matrix_cn</code> and <code>stiffness_matrix_nc</code>)</td>
</tr>
<tr>
<td align="left">dofs $K_{n_k}$ $\leftrightarrow$ dofs $K_{n_k}$</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.
+The last of these accounts that the lifted shape functions from one of the neighbor cells may overlap on $K_c$ with the lifted shape functions of *another* neighbor cell, as mentioned above.
+In other words, 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$ is part of the boundary of the domain. It can also have more neighbors when hanging nodes are present. To simplify the presentation we do not discuss the latter.
-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
+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\left(\jump{\nabla_h\varphi^c}\right)\not\equiv 0$ only if $e\subset\partial K_c$, and similarly for $b_e\left(\jump{\varphi^c}\right)$. 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}).
+H_h(\varphi^c)=D_h^2\varphi^c-\sum_{e\subset\partial K}r_e\left(\jump{\nabla_h \varphi^c}\right)+\sum_{e\subset\partial K}b_e\left(\jump{\varphi^c}\right).
@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$
+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, but not the term involving $D^2)h\varphi^n$, 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).
+H_h(\varphi^n)(x_q)=-r_e\left(\jump{\nabla_h\varphi^n}\right)(x_q)+b_e\left(\jump{\varphi^n}\right)(x_q).
@f]
This information is stored in
@f[
<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 now describe the computation of the lifting operators $r_e$ and $b_e$ on each cell $K_c$. This turns out to be a bit cumbersome, but it follows similar schemes as other reconstruction operators -- see, for example, the "weak Galerkin" approach on step-61 or the "hybridizable discontinuous Galerkin" method in step-51.
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
+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\left(\jump{\varphi}\right)\equiv 0$), we write $b_e\left(\jump{\varphi}\right)\in[\mathbb{V}_h]^{d\times d}$ as
@f[
-b_e(\jump{\varphi})=\sum_{n=1}^{N_c+N_n}B_n\psi_n,
+b_e\left(\jump{\varphi}\right)=\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[
@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.
+M_{mn}\dealcoloneq\int_{\Omega}\psi_n:\psi_m \quad \mbox{and} \quad G_m\dealcoloneq\int_e\average{{\rm div}\, \psi_m}\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[
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).
+b_e\left(\jump{\varphi}\right)(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.
<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
+The performance 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]
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
+For different values of $h$, we will report the error $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),
+\|v\|_{H_h^1(\Omega)}^2 \dealcoloneq \|\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.
projects / dealii.git / commitdiff