]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Updated intro for step-29
authorallmaras <allmaras@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 10 Sep 2007 03:32:23 +0000 (03:32 +0000)
committerallmaras <allmaras@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 10 Sep 2007 03:32:23 +0000 (03:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@15185 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-29/doc/intro.dox

index 1a5015de61d64b6259b8e63b6241703afe5c29ce..0360423334055bfa6ceb2ac642ce0e14ac4a580c 100644 (file)
@@ -1,6 +1,185 @@
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-This program was contributed by Moritz Allmaras at Texas A&amp;M
-University. Some of the work on this tutorial program has been funded
-by NSF under grant DMS-0604778.
+<a name="Intro"></a>\r
+<h1>Introduction</h1>\r
+\r
+<i>\r
+This program was contributed by Moritz Allmaras at Texas A&amp;M\r
+University. Some of the work on this tutorial program has been funded\r
+by NSF under grant DMS-0604778.\r
+</i>\r
+\r
+A question that comes up frequently is how to solve problems involving complex\r
+valued functions with deal.II. For many problems, instead of working with\r
+complex valued finite elements directly, which are not readily available in\r
+the library, it is often much more convenient to split complex valued\r
+functions into their real and imaginary parts and use separate scalar finite\r
+element fields for discretizing each one of them. Basically this amounts to\r
+viewing a single complex valued equation as a system of two real valued\r
+equations. This short example demonstrates how this can be implemented in\r
+deal.II by using an <code>FE_system</code> object to stack two finite element\r
+fields representing real and imaginary parts. We will also revisit the\r
+ParameterHandler class introduced in @ref step_19 "step-19", which provides a\r
+convenient way for reading parameters from a configuration file at runtime\r
+without the need to recompile the program code.  \r
+\r
+<h3>Problem setting</h3>\r
+\r
+The original purpose of this program is to simulate the focussing properties\r
+of an ultrasound wave generated by a transducer lens with variable\r
+geometry. Recent applications in medical imaging use ultrasound waves not only\r
+for imaging porposes, but also to excite certain local effects in a\r
+material, like changes in optical properties, that can then be measured by\r
+other imaging techniques. A vital ingredient for these methods is the ability\r
+to focus the intensity of the ultrasound wave in a particular part of the\r
+material, ideally in a point, to be able to examine the properties of the\r
+material at that particular location.  \r
+\r
+To derive a model for this problem, we think of ultrasound as a pressure wave\r
+governed by the wave equation:  \r
+@f[\r
+       \frac{\partial^2 U}{\partial t^2}       -       c^2 \Delta U = 0\r
+@f]\r
+where $c$ is the wave speed (that for simplicity we assume to be constant), $U\r
+= U(x,t),\;x \in \Omega,\;t\in\mathrm{R}$. The boundary\r
+$\Gamma=\partial\Omega$ is divided into two parts $\Gamma_1$ and\r
+$\Gamma_2=\Gamma\setminus\Gamma_1$, with $\Gamma_1$ representing the\r
+transducer lens and $\Gamma_2$ an absorbing boundary (that is, we want to\r
+choose boundary conditions on $\Gamma_2$ in such a way that they imitate a\r
+larger domain). On $\Gamma_1$, the transducer generates a wave of constant\r
+frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here):  \r
+@f[\r
+U(x,t) = \sin{\omega t}, \qquad x\in \Gamma_1\r
+@f]\r
+\r
+If there are no other (interior or boundary) sources, and since the only\r
+source has frequency $\omega$, then the solution admits a separation of\r
+variables of the form $U(x,t) = \textrm{Re}\left(u(x)\,e^{i\omega\r
+t})\right)$. The complex-valued function $u(x)$ describes the spatial\r
+dependency of amplitude and phase (relative to the source) of the waves of\r
+frequency ${\omega}$, with the amplitude being the quantity that we are\r
+interested in. By plugging this form of the solution into the wave equation,\r
+we see that for $u$ we have\r
+@f{eqnarray*}\r
+-\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\\r
+u(x) &=& 1,  \qquad x\in\Gamma_1. \r
+@f}\r
+\r
+For finding suitable conditions on $\Gamma_2$ that model an absorbing\r
+boundary, consider a wave of the form $V(x,t)=e^{i(k\cdot x -\omega t)}$ with\r
+frequency ${\omega}$ traveling in direction $k\in {\mathrm{R}^2}$. In order\r
+for $V$ to solve the wave equation, $|k|={\frac{\omega}{c}}$ must\r
+hold. Suppose that this wave hits the boundary in $x_0\in\Gamma_2$ at a right\r
+angle, i.e. $n=\frac{k}{|k|}$ with $n$ denoting the outer unit normal of\r
+$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation \r
+@f[\r
+c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0.\r
+@f]\r
+Hence, by enforcing the boundary condition \r
+@f[\r
+c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2,\r
+@f]\r
+waves that hit the boundary $\Gamma_2$ at a right angle will be perfectly\r
+absorbed. On the other hand, those parts of the wave field that do not hit a\r
+boundary at a right angle do not satisfy this condition and enforcing it as a\r
+boundary condition will yield partial reflections, i.e. only parts of the wave\r
+will pass through the boundary as if it wasn't here whereas the remaining\r
+fraction of the wave will be reflected back into the domain.\r
+\r
+If we are willing to accept this as a suffient approximation to an absorbing boundary we finally arrive at the following problem for $u$: \r
+@f{eqnarray*}\r
+-\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\\r
+c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\\r
+u &=& 1,  \qquad x\in\Gamma_1. \r
+@f}\r
+This is a Helmholtz equation (similar to the one in @ref step_7 "step-7", but this time with ''the bad sign'') with Dirichlet data on $\Gamma_1$ and mixed boundary conditions on $\Gamma_2$. Because of the condition on $\Gamma_2$, we cannot just treat the equations for real and imaginary parts of $u$ separately. What we can do however is to view the PDE for $u$ as a system of two PDEs for the real and imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system:\r
+@f{eqnarray*}\r
+  \left.\begin{array}{ccc}\r
+    -\omega^2 v - c^2\Delta v &=& 0 \quad\\ \r
+    -\omega^2 w - c^2\Delta w &=& 0 \quad\r
+  \end{array}\right\} &\;& x\in\Omega,\r
+       \\\r
+  \left.\begin{array}{ccc}\r
+    c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\ \r
+    c (n\cdot\nabla w) + \omega\,v &=& 0 \quad\r
+  \end{array}\right\} &\;& x\in\Gamma_2,\r
+       \\\r
+       \left.\begin{array}{ccc}\r
+    v &=& 1 \quad\\ \r
+    w &=& 0 \quad\r
+  \end{array}\right\} &\;& x\in\Gamma_1.\r
+@f}\r
+\r
+For test functions $\phi,\psi$ with $\phi|_{\Gamma_1}=\psi|_{\Gamma_1}=0$, after the usual multiplication, integration over $\Omega$ and applying integration by parts, we get the weak formulation\r
+@f{eqnarray*}\r
+-\omega^2 \langle \phi, v \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \phi, \nabla v \rangle_{\mathrm{L}^2(\Omega)} - c \omega \langle \phi, w \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0, \\\r
+-\omega^2 \langle \psi, w \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \psi, \nabla w \rangle_{\mathrm{L}^2(\Omega)} + c \omega \langle \psi, v \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0.\r
+@f}\r
+\r
+We choose finite element spaces $V_h$ and $W_h$ with bases $\{\phi_j\}_{j=1}^n, \{\psi_j\}_{j=1}^n$ and look for approximate solutions \r
+@f[\r
+v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j. \r
+@f]\r
+Plugging into the variational form yields the equation system\r
+@f[\r
+\renewcommand{\arraystretch}{2.0}\r
+\left.\begin{array}{ccc}\r
+\sum_{j=1}^n \left(-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\alpha_j - \left(c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\beta_j &=& 0 \\ \r
+\sum_{j=1}^n \left(-\omega^2 \langle \psi_i, \psi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \psi_i, \nabla \psi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\beta_j + \left(c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\alpha_j &=& 0\r
+\end{array}\right\}\;\;\forall\; i =1,\ldots,n.\r
+@f]\r
+In matrix notation: \r
+@f[\r
+\renewcommand{\arraystretch}{2.0}\r
+\left(\r
+\begin{array}{cc}\r
+-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)} & -c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)} \\\r
+c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)} & -\omega^2 \langle \psi_{i}, \psi_j \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \psi_{i}, \nabla \psi_j  \rangle_{\mathrm{L}^2(\Omega)}\r
+\end{array}\r
+\right)\r
+\left(\r
+\begin{array}{c}\r
+\alpha \\ \beta\r
+\end{array}\r
+\right)\r
+=\r
+\left(\r
+\begin{array}{c}\r
+0 \\ 0\r
+\end{array}\r
+\right)\r
+@f]\r
+(One should not be fooled by the right hand side being zero here, that is because we haven't included the Dirichlet boundary data yet.) \r
+Because of the alternating sign in the off-diagonal blocks, we can already see that this system is non-symmetric, in fact it is even indefinite.\r
+Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be\r
+the same. However, we expect the real and imaginary parts of the solution to\r
+have similar properties and will therefore indeed choose $V_h=W_h$. We choose\r
+the current notation using different symbols merely to be able to distinguish\r
+between test functions for $v$ and $w$, as this distinction plays an important\r
+role in the implementation.  \r
+\r
+\r
+<h3>The test case</h3>\r
+\r
+For the computations, we will consider wave propagation in a rectangular\r
+area. Sound is generated by a transducer that we model using a part of the\r
+boundary shaped like a segment of the circle with center at $(0.5, d)$ and a\r
+radius slightly greater than $d$; this shape leads to a focusing of the sound\r
+wave to the center of the circle. Consequently, we take $\Omega$ to be the\r
+unit square $[0,1]^2$ with the piece $[0.4,0.6]\times\{0\}$ removed and\r
+replaced by the transducer lens (see the Results section below for a graphical\r
+depiction of this domain). We let $d$ vary to change the "focus" of the lens\r
+and see how that affects the spatial distribution of the amplitude of $u$,\r
+i.e. we will consider how well focussed $|u|=\sqrt{v^2+w^2}$ is.\r
+\r
+In the program below, we will implement the complex-valued Helmholtz equations\r
+using the formulation with split real and imaginary parts. We will also\r
+discuss how to generate a domain that looks like a square with a slight bulge\r
+simulating the transducer (in the\r
+<code>UltrasoundProblem<dim>::make_grid()</code> function), and how to\r
+generate graphical output that not only contains the solution components $v,w$\r
+but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in\r
+<code>UltrasoundProblem<dim>::output_results()</code>). Finally, we will\r
+explain how to use the ParameterHandler class to easily handle the situation\r
+where we want to prescribe the focal distance $d$, wave speed $c$, frequency\r
+$\omega$, and a number of other parameters, in an input file that is read at\r
+run-time, rather than in the source code where we would have to re-compile the\r
+file every time we want to change parameters.\r

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.