1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
\subsection{The Stokes problem}
\label{dg-sec-stokes}
\pbindex{Stokes}%
\cindex{boundary condition!Dirichlet}%
\cindex{benchmark!driven cavity flow}%
Let us consider the Stokes problem for the driven cavity
in $\Omega = ]0,1[^d$, $d = 2,3$.
The problem has been introduced in volume~1, section~\ref{sec-stokes},
page~\pageref{sec-stokes}.
\ \ $(P)$: \emph{find ${\bf u}$ and $p$, defined in $\Omega$, such that}
\[
\begin{array}{ccccl}
-\ {\bf div}(2 D({\bf u}))
&+& \bnabla p &=& {\bf f} \ {\rm in}\ \Omega, \\
-\ {\rm div}\,{\bf u} && &=& 0 \ {\rm in}\ \Omega, \\
{\bf u} && &=& {\bf g} \ {\rm on}\ \partial\Omega
\end{array}
\]
where ${\bf f}$ and ${\bf g}$ are given.
This problem is the extension to divergence free vector fields
of the elasticity problem.
The variational formulation writes:
\ \ $(VF)_h$ {\it find ${\bf u} \in {\bf V}({\bf g})$ and $p \in L^2(\Omega)$ such that:}
\begin{equation}
\begin{array}{lcccl}
a({\bf u},{\bf v}) &+& b({\bf v}, p) &=& l({\bf v}), \ \forall {\bf v}\in {\bf V}(0), \\
b({\bf u},q) & & &=& 0, \ \forall q \in L^2(\Omega)
\end{array}
\label{dg-eq-fv-stokes}
\end{equation}
where
\begin{eqnarray*}
{\bf V}({\bf g}) &=& \{ {\bf v}\in H^1(\Omega)^d; \ {\bf v}={\bf g}\ \mbox{ on }\ \partial\Omega \}
\\
a({\bf u},{\bf v})
&=&
\int_\Omega
2D({\bf u})\!:\!D({\bf v})
\, {\rm d}x
\\
b({\bf u},q)
&=&
-\int_\Omega
{\rm div}({\bf u})\, q
\, {\rm d}x
\\
l({\bf v})
&=&
\int_\Omega
{\bf f}.{\bf v}
\, {\rm d}x
\end{eqnarray*}
The discrete variational formulation writes:
\ \ $(VF)_h$ {\it find ${\bf u}_h \in {\bf X}_h$ and $p_h \in Q_h$ such that:}
\begin{equation}
\begin{array}{lcccl}
a_h({\bf u}_h,{\bf v}_h) &+& b_h({\bf v}_h, p_h) &=& l_h({\bf v}_h), \ \forall {\bf v}_h\in {\bf X}_h, \\
b_h({\bf u}_h,q_h) &-& c_h(p_h,q_h) &=& k_h(q), \ \forall q_h \in Q_h.
\end{array}
\label{dg-eq-fvh-stokes}
\end{equation}
The discontinuous finite element spaces are defined by:
\begin{eqnarray*}
{\bf X}_h &=& \{ {\bf v}_h \in L^2(\Omega)^d; {\bf v}_{h|K}\in P_k^d, \ \forall K \in \mathcal{T}_h \}
\\
Q_h &=& \{ q_h \in L^2(\Omega)^d; q_{h|K}\in P_k^d, \ \forall K \in \mathcal{T}_h \}
\end{eqnarray*}
where $k\geq 1$ is the polynomial degree.
Note that velocity and pressure are approximated by the same polynomial order.
This method was introduced by \citet{CocKanSchSch-2002}
and some recent theoretical results can be founded in \citet{PieErn-2010}.
The forms are defined for all $u,v\in H^1(\mathcal{T}_h)^d$
and $q\in L^2(\Omega)$ by
(see e.g. \citealp[p.~249]{PieErn-2012}):
\begin{eqnarray*}
a_h({\bf u},{\bf v})
&=&
\int_\Omega
2D_h({\bf u})\!:\!D_h({\bf v})
\, {\rm d}x
\\
&& + \
\sum_{S\in\mathscr{S}_h}
\int_S
\left(
\beta\varpi_s \jump{\bf u}.\jump{\bf v}
- \jump{\bf u}.\average{2D_h({\bf v}){\bf n}}
- \jump{\bf v}.\average{2D_h({\bf u}){\bf n}}
\right)
\, {\rm d}s
\\
b_h({\bf u},q)
&=&
\int_\Omega
{\bf u}.\nabla_h q
\, {\rm d}x
-
\sum_{S\in\mathscr{S}_h^{(i)}}
\int_S
\average{\bf u}.{\bf n}\,\jump{q}
\, {\rm d}s
\\
c_h(p,q)
&=&
\sum_{S\in\mathscr{S}_h^{(i)}}
\int_S
h_s \,
\jump{p}\,\jump{q}
\, {\rm d}s
\\
l_h({\bf v})
&=&
\int_\Omega
{\bf f}.{\bf v}
\,{\rm d}s
+
\int_{\partial\Omega}
{\bf g}.
\left(
\beta\varpi_s
\, {\bf v}
-
2D_h({\bf v})\,{\bf n}
\right)
\,{\rm d}s
\\
k_h(q)
&=&
\int_{\partial\Omega}
{\bf g}.{\bf n} \, q
\,{\rm d}s
\end{eqnarray*}
The stabilization form $c_h$ controls the pressure jump
across internal sides.
This stabilization term is necessary when using
equal order polynomial approximation for velocity and pressure.
The definition of the forms is grouped in a subroutine: it will
be reused later for the Navier-Stokes problem.
% ---------------------------------------
\myexamplelicense{stokes_dirichlet_dg.icc}
% ---------------------------------------
A simple test program writes:
% ---------------------------------------
\myexamplelicense{stokes_taylor_dg.cc}
% ---------------------------------------
\subsubsection*{Comments}
The data are given when $d=2$ by~\eqref{eq-taylor-benchmark}.
This choice is convenient since the exact
solution is known ${\bf u}={\bf g}$ and $p=0$.
\rawexindex{taylor.h}
The
\rawexindex{stokes_taylor_error_dg.cc}%
code \code{stokes_taylor_error_dg.cc} compute the error in $L^2$, $L^\infty$ and energy norms.
This code it is not listed here but is available in the \Rheolef\ example directory.
The computation writes:
\begin{verbatim}
make stokes_taylor_dg stokes_taylor_error_dg
mkgeo_grid -t 10 > square.geo
./stokes_taylor_dg square P1d | ./stokes_taylor_error_dg
./stokes_taylor_dg square P2d | ./stokes_taylor_error_dg
\end{verbatim}
|
