\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
$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
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 =
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
+
\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
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;
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)
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];
projects / dealii.git / commitdiff