From a64a6f4293acce20a1dbc77403bdc1a36a571774 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 13 Jan 2020 16:06:49 -0700 Subject: [PATCH] Address the comments by peterrum. --- examples/step-71/doc/intro.dox | 50 +++++++++++++++++++++---------- examples/step-71/step-71.cc | 55 ++++++++++++++++++++++++---------- 2 files changed, 74 insertions(+), 31 deletions(-) diff --git a/examples/step-71/doc/intro.dox b/examples/step-71/doc/intro.dox index 72943194ab..f17e8df24d 100644 --- a/examples/step-71/doc/intro.dox +++ b/examples/step-71/doc/intro.dox @@ -86,8 +86,8 @@ piecewise linear "hat functions", we would get integrals of the sort \int_0^L (\Delta \varphi_i) (\Delta \varphi_j) = \int_0^L - \frac 1h \right[\delta(x-x_{i-1}) - 2\delta(x-x_i) + \delta(x-x_{i+1})\right] - \frac 1h \right[\delta(x-x_{j-1}) - 2\delta(x-x_j) + \delta(x-x_{j+1})\right] + \frac 1h \left[\delta(x-x_{i-1}) - 2\delta(x-x_i) + \delta(x-x_{i+1})\right] + \frac 1h \left[\delta(x-x_{j-1}) - 2\delta(x-x_j) + \delta(x-x_{j+1})\right] @f} where $x_i$ is the node location at which the shape function $\varphi_i$ is defined, and $h$ is the mesh size (assumed @@ -230,8 +230,16 @@ values on each interface $e\in \mathbb{F}$ (shared by the two cells $K_{+},K_{-} \in \mathbb{T}$), we cope with this discontinuity by defining the following single-valued functions on $e$: @f{align*}{ -\Bigg[\!\Bigg[\frac{\partial^k v_h}{\partial \mathbf n^k}\Bigg]\!\Bigg]&= \frac{\partial^k v_h|_{K_+}}{\partial \mathbf n^k} \bigg |_e - \frac{\partial^k v_h|_{K_-}}{\partial \mathbf n^k} \bigg |_e,\\ -\Bigg\{\!\Bigg\{\frac{\partial^k v_h}{\partial \mathbf n^k}\Bigg\}\!\Bigg\}&=\frac{1}{2} \bigg( \frac{\partial^k v_h|_{K_+}}{\partial \mathbf n^k} \bigg |_e + \frac{\partial^k v_h|_{K_-}}{\partial \mathbf n^k} \bigg |_e \bigg ) + \jump{\frac{\partial^k v_h}{\partial \mathbf n^k}} + &= + \frac{\partial^k v_h|_{K_+}}{\partial \mathbf n^k} \bigg |_e + - \frac{\partial^k v_h|_{K_-}}{\partial \mathbf n^k} \bigg |_e, + \\ + \average{\frac{\partial^k v_h}{\partial \mathbf n^k}} + &= + \frac{1}{2} + \bigg( \frac{\partial^k v_h|_{K_+}}{\partial \mathbf n^k} \bigg |_e + + \frac{\partial^k v_h|_{K_-}}{\partial \mathbf n^k} \bigg |_e \bigg ) @f} for $k =1,2$ (i.e., for the gradient and the matrix of second derivatives), and where $\mathbf n$ denotes a unit vector normal to @@ -300,10 +308,13 @@ stable as the mesh size goes to zero, and to obtain a stable formulation that converges to the correct solution, we need to add the following terms: @f{align*}{ --\sum_{e \in \mathbb{F}} \int_{e} \bigg \{\!\bigg\{ \frac{\partial^2 v_h}{\partial \mathbf n^2}\bigg\}\!\bigg\} \bigg [\!\bigg[ \frac{\partial u_h}{\partial \mathbf n}\bigg]\!\bigg] +-\sum_{e \in \mathbb{F}} \int_{e} + \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} + \jump{\frac{\partial u_h}{\partial \mathbf n}} + \sum_{e \in \mathbb{F}} -\frac{\sigma}{h_e}\int_e \bigg[\!\bigg[\! \frac{\partial v_h}{\partial \mathbf n} \bigg]\!\bigg] -\bigg[\!\bigg[ \frac{\partial u_h}{\partial \mathbf n} \bigg]\!\bigg]. + \frac{\sigma}{h_e}\int_e + \jump{\frac{\partial v_h}{\partial \mathbf n}} + \jump{\frac{\partial u_h}{\partial \mathbf n}}. @f} Then, after making cancellations that arise, we arrive at the following $C^0$ IP formulation of the biharmonic equation: find $u_h$ such that $u_h = @@ -314,19 +325,28 @@ g$ on $\partial \Omega$ and where @f{align*}{ \mathcal{A}(v_h,u_h):=&\sum_{K \in \mathbb{T}}\int_K D^2v_h:D^2u_h \ dx --\sum_{e \in \mathbb{F}} \int_{e} \bigg \{\!\bigg\{ \frac{\partial^2 v_h}{\partial \mathbf n^2}\bigg\}\!\bigg\} \bigg [\!\bigg[ \frac{\partial u_h}{\partial \mathbf n}\bigg]\!\bigg] - ds \\ -&-\sum_{e \in \mathbb{F}} \int_{e} \bigg \{\!\bigg\{ \frac{\partial^2 v_h}{\partial \mathbf n^2}\bigg\}\!\bigg\} \bigg [\!\bigg[ \frac{\partial u_h}{\partial \mathbf n}\bigg]\!\bigg] \ ds +& + -\sum_{e \in \mathbb{F}} \int_{e} + \jump{ \frac{\partial^2 v_h}{\partial \mathbf n^2}} + \average{\frac{\partial u_h}{\partial \mathbf n}} \ ds + -\sum_{e \in \mathbb{F}} \int_{e} + \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} + \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds \\ &+ \sum_{e \in \mathbb{F}} -\frac{\gamma}{h_e}\int_e \bigg[\!\bigg[\! \frac{\partial v_h}{\partial \mathbf n} \bigg]\!\bigg] -\bigg[\!\bigg[ \frac{\partial u_h}{\partial \mathbf n} \bigg]\!\bigg] \ ds, + \frac{\gamma}{h_e} + \int_e + \jump{\frac{\partial v_h}{\partial \mathbf n}} + \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds, @f} and @f{align*}{ -\mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx +\sum_{e \in \mathbb{F}^b} \frac{\gamma}{h_e}\int_e -\frac{\partial v_h}{\partial \mathbf n} j \ ds. +\mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx ++ +\sum_{e \in \mathbb{F}, e\subset\partial\Omega} +\frac{\gamma}{h_e} +\int_e \frac{\partial v_h}{\partial \mathbf n} j \ ds. @f} Here, $\gamma$ is the penalty parameter which both weakly enforces the boundary condition @@ -367,7 +387,7 @@ energy norm, as follows: + \sum\limits_{e \in \mathbb{F} } \frac{\gamma }{h_e} \left\| - \ \!\bigg[ \!\bigg[ \frac{\partial u_h}{\partial \mathbf n}\bigg]\!\bigg] \right\|_{L^2(e)}^2. + \jump{\frac{\partial u_h}{\partial \mathbf n}} \right\|_{L^2(e)}^2. @f} In this (semi)norm, the theory in the paper mentioned above yields that we diff --git a/examples/step-71/step-71.cc b/examples/step-71/step-71.cc index 04be59b3a1..b8f2b44e15 100644 --- a/examples/step-71/step-71.cc +++ b/examples/step-71/step-71.cc @@ -397,10 +397,10 @@ namespace Step71 copy_data.cell_matrix = 0; copy_data.cell_rhs = 0; - scratch_data.fe_values.reinit(cell); - cell->get_dof_indices(copy_data.local_dof_indices); - const FEValues &fe_values = scratch_data.fe_values; + fe_values.reinit(cell); + + cell->get_dof_indices(copy_data.local_dof_indices); const ExactSolution::RightHandSide right_hand_side; @@ -512,13 +512,29 @@ namespace Step71 ncell->extent_in_direction( GeometryInfo::unit_normal_direction[nf]))); - // Finally, and as usual, we loop over the quadrature points - // and indices `i` and `j` to add up the contributions of this - // face or sub-face. These are then stored in the `copy_data.face_data` - // object created above. As for the cell worker, we pull the evalation - // of averages and jumps out of the loops if possible, introducing - // local variables that store these results. The assembly then only - // needs to use these local variables in the innermost loop. + // Finally, and as usual, we loop over the quadrature points and + // indices `i` and `j` to add up the contributions of this face + // or sub-face. These are then stored in the + // `copy_data.face_data` object created above. As for the cell + // worker, we pull the evalation of averages and jumps out of + // the loops if possible, introducing local variables that store + // these results. The assembly then only needs to use these + // local variables in the innermost loop. Regarding the concrete + // formula this code implements, recall that the interface terms + // of the bilinear form were as follows: + // @f{align*}{ + // -\sum_{e \in \mathbb{F}} \int_{e} + // \jump{ \frac{\partial^2 v_h}{\partial \mathbf n^2}} + // \average{\frac{\partial u_h}{\partial \mathbf n}} \ ds + // -\sum_{e \in \mathbb{F}} \int_{e} + // \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} + // \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds + // + \sum_{e \in \mathbb{F}} + // \frac{\gamma}{h_e} + // \int_e + // \jump{\frac{\partial v_h}{\partial \mathbf n}} + // \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds. + // @f} for (unsigned int qpoint = 0; qpoint < fe_interface_values.n_quadrature_points; ++qpoint) @@ -600,12 +616,19 @@ namespace Step71 cell->extent_in_direction( GeometryInfo::unit_normal_direction[face_no])); - // The third piece is the assembly of terms. This is now slightly more - // involved since these contains both terms for the matrix and for - // the right hand side. The latter requires us to evaluate the - // boundary conditions $j(\mathbf x)$, which in the current - // case (where we know the exact solution) we compute from - // $j(\mathbf x) = \frac{\partial u(\mathbf x)}{\partial {\mathbf n}}$: + // The third piece is the assembly of terms. This is now + // slightly more involved since these contains both terms for + // the matrix and for the right hand side. The former is exactly + // the same as for the interior faces stated above if one just + // defines the jump and average appropriately (which is what the + // FEInterfaceValues class does). The latter requires us to + // evaluate the boundary conditions $j(\mathbf x)$, which in the + // current case (where we know the exact solution) we compute + // from $j(\mathbf x) = \frac{\partial u(\mathbf x)}{\partial + // {\mathbf n}}$. The term to be added to the right hand side + // vector is then + // $\frac{\gamma}{h_e}\int_e + // \frac{\partial v_h}{\partial \mathbf n} j \ ds$. for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) { const auto &n = normals[qpoint]; -- 2.39.5