]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Address the comments by peterrum.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 13 Jan 2020 23:06:49 +0000 (16:06 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Tue, 14 Jan 2020 00:52:45 +0000 (17:52 -0700)
examples/step-71/doc/intro.dox
examples/step-71/step-71.cc

index 72943194ab5d2608931107778d8b953de4c75549..f17e8df24dc45f4ba5b84d4001aa26d31bb3c0aa 100644 (file)
@@ -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
index 04be59b3a122d38afa1c097a27f351b805237d5b..b8f2b44e153e5a8aff973e6476f0183a7889d4fb 100644 (file)
@@ -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<dim> &fe_values = scratch_data.fe_values;
+      fe_values.reinit(cell);
+
+      cell->get_dof_indices(copy_data.local_dof_indices);
 
       const ExactSolution::RightHandSide<dim> right_hand_side;
 
@@ -512,13 +512,29 @@ namespace Step71
                   ncell->extent_in_direction(
                     GeometryInfo<dim>::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<dim>::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];

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.