1
\label{sec-transmission}
2
This chapter is related to the so-called transmission problem.
6
% 1d: maillages non-uniformes aussi
7
% 2d: uniforme ok, non-unif non !
8
% mais super-convergence ?
11
% quadrature et P1d-P2 ?
12
% => voir whalbin ! et book E.F. zienkiewicz
13
% proj P1-C0 pi_h(gh) du gradient P1d gh : superconvergence...?
14
% estimateur d'erreur gh-pi_h(gh) directionnel ?
17
\cindex{Dirichlet boundary condition}
19
\cindex{Neumann boundary condition}
21
\cindex{transmission problem}
22
\pbindex{transmission}
24
\section*{Formulation}
25
% -------------------------------------
27
Let us consider the diffusion problem
28
with a non-constant diffusion coefficient $\eta$
29
in a domain bounded $\Omega \subset R^N$,
32
{\it $(P)$: find $u$ defined in $\Omega$ such that:}
34
-{\rm div}(\eta \bnabla u) &=& 1 \ {\rm in}\ \Omega \\
35
u &=& 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm right} \\
36
\Frac{\partial u}{\partial n} &=& 0 \ {\rm on}\ \Gamma_{\rm top} \cup \Gamma_{\rm bottom} \mbox{ when } N \geq 2
38
\Frac{\partial u}{\partial n} &=& 0 \ {\rm on}\ \Gamma_{\rm front} \cup \Gamma_{\rm back} \mbox{ when } N = 3
40
We consider here very important special case when $\eta$ is
41
picewise constant. Let:
43
\eta (x) = \left\{ \begin{array}{ll}
44
\varepsilon & \mbox{ when } x \in \Omega_{\rm west} \\
45
1 & \mbox{ when } x \in \Omega_{\rm east}
48
where $(\Omega_{\rm west}, \Omega_{east})$ is a partition of
49
$\Omega$. This problem is known as the {\bf transmission} problem,
50
since the solution and the flux are continuous on the interface:
52
u_{\Omega_{\rm west}} &=& u_{\Omega_{\rm east}}
53
\mbox { on } \Gamma_0 \\
54
\varepsilon \Frac{\partial u_{/\Omega_{\rm west}} }{\partial n}
56
\Frac{\partial u_{\Omega_{\rm east}} }{\partial n}
57
\mbox { on } \Gamma_0 \\
61
\Gamma_0 = \partial \Omega_{\rm west} \cap \partial \Omega_{\rm east}
64
It expresses the transmission of the quantity $u$ and its flux
65
accross the interface $\Gamma_0$
66
between two regions that have different diffusion properties:
67
We go back to the standard problem when $\varepsilon=1$.
69
The variational formulation of this problem expresses:
71
{\it $(VF)$: find $u\in V$ such that:}
76
where the bilinear forms $a(.,.)$ and $m(.,.)$ are defined by
78
a(u,v) &=& \int_\Omega \eta \nabla u . \nabla v \, dx,
79
\ \ \forall u,v \in H^1(\Omega)
81
m(u,v) &=& \int_\Omega u v \, dx,
82
\ \ \forall u,v \in L^2(\Omega)
84
V &=& \{ v \in H^1(\Omega); \
85
v = 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm right} \}
89
The bilinear form $a(.,.)$ defines a scalar product
90
in $V$ and is related to the {\it energy} form.
91
This form is associated to the $-{\rm div}\,\eta\bnabla$ operator.
92
The $m(.,.)$ is here the classical scalar product on $L^2(\Omega)$,
93
and is related to the {\it mass} form as usual.
95
The approximation of this problem could performed by
96
a standard Lagrange $P_k$ continuous approximation.
97
We switch here to a mixed formulation that implements easily.
99
\section*{Mixed Formulation}
100
% -------------------------------------
101
\cindex{mixed formulation}
103
Let us introduce the following vector-valued field:
105
{\bf p} = \eta \, \bnabla u
107
The problem can be rewritten as:
109
{\it (M): find ${\bf p}$ and $u$ defined in $\Omega$ such that:}
110
\[ \begin{array}{rlrlrl}
111
\Frac{1}{\eta} \, {\bf p} &-& \bnabla u &=& 0 & \mbox{ in } \Omega \\
112
{\rm div}\, {\bf p} & & &=& -1 & \mbox{ in } \Omega \\
114
\mbox{ on } \Gamma_{\rm left} \cup \Gamma_{\rm right}
116
&& \Frac{\partial u}{\partial n} &=& 0 &
117
\mbox{ on } \Gamma_{\rm top} \cup \Gamma_{\rm bottom} \mbox{ when } N \geq 2
119
&& \Frac{\partial u}{\partial n} &=& 0 &
120
\mbox{ on } \Gamma_{\rm front} \cup \Gamma_{\rm back} \mbox{ when } N = 3
123
The variationnal formulation writes:
125
{\it $(MVF)$: find ${\bf p}\in (L^2(\Omega))^N$ and $u\in V$ such that:}
126
\[ \begin{array}{llllll}
127
\tilde{a}({\bf p},{\bf q}) &+& \tilde{b}({\bf q},u) &=& 0, & \forall {\bf q}\in (L^2(\Omega))^N \\
128
\tilde{b}({\bf p},v) & & &=& m(-1,v), & \forall v\in V
133
\tilde{a}({\bf p},{\bf q}) &=& \int_\Omega \Frac{1}{\eta} {\bf p}. {\bf q}\, dx,
135
\tilde{b}({\bf q},v) &=& -\int_\Omega \bnabla v . {\bf q} \, dx,
136
\ \ \forall u,v \in L^2(\Omega)
138
This problem appears as more complex than the previous one,
139
since there is two unknown instead of one.
140
Nevertheless, the ${\bf p}$ unknown could be eliminated
141
by inverting the operator associated to the $\tilde{a}$ form.
142
This can be done easily since ${\bf p}$ is approximated by
143
discontinuous piecewise polynomial functions.
145
\section*{Mixed approximation}
146
% -------------------------------------
150
{\bf Q}_h &=& \{ {\bf q} \in (L^2(\Omega))^N; \ {\bf q}_{/K} \in (P_{k-1})^N \} \\
151
V_h &=& \{ v \in V; \ v_{/K} \in P_k \}
153
The mixed finite element approximation writes:
155
{\it $(MVF)_h$: find ${\bf p}_h\in {\bf Q}_h$ and $u_h\in V_h$ such that:}
156
\[ \begin{array}{llllll}
157
\tilde{a}({\bf p}_h,{\bf q}) &+& \tilde{b}({\bf q},u_h) &=& 0, & \forall {\bf q}\in {\bf Q}_h \\
158
\tilde{b}({\bf p}_h,v) & & &=& m(-1,v), & \forall v\in V_h
161
Here the operator associated to $\tilde{a}$ is block-diagonal
162
and can be inverted at the element level.
163
The problem reduced to a standard linear system, and
164
solved by a direct method.
166
% -------------------------------------
167
\section*{File \file{transmission.cc}}
168
% -------------------------------------
174
\exindex{transmission.cc}
175
\verbatiminput{transmission.cc}
177
% -------------------------------------
179
% -------------------------------------
180
\cindex{diffusion tensor}
182
The code begins as usual.
183
In the second part, the space ${\bf Q}_h$ is defined.
184
For convenience, the diffusion coefficient $\eta$ is introduced
185
diagonal diffusion tensor operator $d$, i.e. for $N=2$:
187
d = \left( \begin{array}{cc}
191
It is first defined as an element of ${\bf Q}_h$ i.e. a vector-valued field:
195
and then converted to a diagonal form representing the
196
diffusion tensor operator:
200
Notice that, by this way, an anisotropic diffusion operator
205
cout << setprecision(numeric_limits<Float>::digits10)
206
<< catchmark("epsilon") << epsilon << endl
207
<< catchmark("u") << uh;
209
writes $\varepsilon$ and $u_h$
210
in a labeled format, suitable for post-traitement,
211
visualization and error analysis.
213
\cindex{floating point precision}%
214
The number of digits printed depends upon the \code{Float}
215
type. When \code{Float} is a \code{double}, this is roughly
216
$15$ disgits; this number depends upon the machine precision.
217
Here, we choose to print all the relevant digits: this
218
choice is suitable for the following error analysis.
221
\section*{How to run the program ?}
222
% -------------------------------------
226
Build the program as usual: \code{make transmission}.
227
Then, creates a one-dimensional geometry with two
230
mkgeo_grid -e 100 -region > line-region.geo
231
geo -gnuplot line-region.geo
233
The trivial mesh generator \code{mkgeo\_grid},
235
regions \code{east} and \code{west} when used with the
236
\code{-region} option.
237
This correspond to the situation:
241
\Omega_{\rm west} = \Omega \cap \{ x_0 < 1/2 \}
243
\Omega_{\rm east} = \Omega \cap \{ x_0 > 1/2 \}.
245
In order to avoid mistakes with the {\tt C++} style indexation,
246
we denote by $(x_0,\ldots,x_{N-1})$ the cartesian coordinate
249
Finaly, run the program and look at the solution:
251
./transmission line-region.geo > transmission1d.mfield
252
mfield transmission1d.mfield -u | field -
254
The two dimensionnal case corresponds to the commands:
256
mkgeo_grid -t 10 -region > square-region.geo
257
geo -plotmtv square-region.geo
258
./transmission square-region.geo > transmission2d.mfield
259
mfield transmission2d.mfield -u | field -
261
while the tridimensional to
263
mkgeo_grid -T 10 -region > cube-region.geo
264
./transmission cube-region.geo > transmission3d.mfield
265
mfield transmission3d.mfield -u | field -mayavi -
267
% See also apendix~\ref{sec-region} for building an unstructured
268
% mesh on a general geometry with regions.
270
% ================================
271
% \subsection*{Error analysis}
272
% ================================
274
This problem is convenient, since the
275
exact solution is known and is piecewise second order polynomial:
279
\Frac{x_0}{2\varepsilon}
280
\left( \Frac{1+3\varepsilon}{2(1+\varepsilon)} - x_0 \right)
281
& \mbox{ when } x_0 < 1/2 \\
284
\left(x_0 + \Frac{1-\varepsilon}{2(1+\varepsilon)} \right)
288
Since the change in the diffusion
289
coefficient value fits the element boundaries,
290
obtain the exact solution at the
291
vertices of the mesh for the $P_1$ approximation,
292
as shown on Fig.~\ref{fig-transmission}.
296
\input{transmission-1d-fig.tex}
298
\caption{Transmission problem: $u_h = \pi_h(u)$
299
($\varepsilon=10^{-2}$, $N=1$, $P_1$ approximation).}
300
\label{fig-transmission}
added
removed
doc/usrman/transmission.tex
