From 0ac288a73cc17500e12e5012134a84bbac6b0d3a Mon Sep 17 00:00:00 2001 From: bangerth Date: Sun, 28 Oct 2007 23:19:29 +0000 Subject: [PATCH] Fix the properties of all these files. git-svn-id: https://svn.dealii.org/trunk@15383 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-29/doc/intro.dox | 374 +++++++++++++------------ 1 file changed, 190 insertions(+), 184 deletions(-) diff --git a/deal.II/examples/step-29/doc/intro.dox b/deal.II/examples/step-29/doc/intro.dox index f2fc0d7231..3f18293ba8 100644 --- a/deal.II/examples/step-29/doc/intro.dox +++ b/deal.II/examples/step-29/doc/intro.dox @@ -1,184 +1,190 @@ - -

Introduction

- - -This program was contributed by Moritz Allmaras at Texas A&M -University. Some of the work on this tutorial program has been funded -by NSF under grant DMS-0604778. - - -A question that comes up frequently is how to solve problems involving complex -valued functions with deal.II. For many problems, instead of working with -complex valued finite elements directly, which are not readily available in -the library, it is often much more convenient to split complex valued -functions into their real and imaginary parts and use separate scalar finite -element fields for discretizing each one of them. Basically this amounts to -viewing a single complex valued equation as a system of two real valued -equations. This short example demonstrates how this can be implemented in -deal.II by using an FE_system object to stack two finite element -fields representing real and imaginary parts. We also revisit the -ParameterHandler class introduced in @ref step_19 "step-19", which provides a -convenient way for reading parameters from a configuration file at runtime -without the need to recompile the program code. - -

Problem setting

- -The original purpose of this program is to simulate the focussing properties -of an ultrasound wave generated by a transducer lens with variable -geometry. Recent applications in medical imaging use ultrasound waves not only -for imaging porposes, but also to excite certain local effects in a -material, like changes in optical properties, that can then be measured by -other imaging techniques. A vital ingredient for these methods is the ability -to focus the intensity of the ultrasound wave in a particular part of the -material, ideally in a point, to be able to examine the properties of the -material at that particular location. - -To derive a model for this problem, we think of ultrasound as a pressure wave -governed by the wave equation: -@f[ - \frac{\partial^2 U}{\partial t^2} - c^2 \Delta U = 0 -@f] -where $c$ is the wave speed (that for simplicity we assume to be constant), $U -= U(x,t),\;x \in \Omega,\;t\in\mathrm{R}$. The boundary -$\Gamma=\partial\Omega$ is divided into two parts $\Gamma_1$ and -$\Gamma_2=\Gamma\setminus\Gamma_1$, with $\Gamma_1$ representing the -transducer lens and $\Gamma_2$ an absorbing boundary (that is, we want to -choose boundary conditions on $\Gamma_2$ in such a way that they imitate a -larger domain). On $\Gamma_1$, the transducer generates a wave of constant -frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here): -@f[ -U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1 -@f] - -If there are no other (interior or boundary) sources, and since the only -source has frequency $\omega$, then the solution admits a separation of -variables of the form $U(x,t) = \textrm{Re}\left(u(x)\,e^{i\omega -t})\right)$. The complex-valued function $u(x)$ describes the spatial -dependency of amplitude and phase (relative to the source) of the waves of -frequency ${\omega}$, with the amplitude being the quantity that we are -interested in. By plugging this form of the solution into the wave equation, -we see that for $u$ we have -@f{eqnarray*} --\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\ -u(x) &=& 1, \qquad x\in\Gamma_1. -@f} - -For finding suitable conditions on $\Gamma_2$ that model an absorbing -boundary, consider a wave of the form $V(x,t)=e^{i(k\cdot x -\omega t)}$ with -frequency ${\omega}$ traveling in direction $k\in {\mathrm{R}^2}$. In order -for $V$ to solve the wave equation, $|k|={\frac{\omega}{c}}$ must -hold. Suppose that this wave hits the boundary in $x_0\in\Gamma_2$ at a right -angle, i.e. $n=\frac{k}{|k|}$ with $n$ denoting the outer unit normal of -$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation -@f[ -c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0. -@f] -Hence, by enforcing the boundary condition -@f[ -c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2, -@f] -waves that hit the boundary $\Gamma_2$ at a right angle will be perfectly -absorbed. On the other hand, those parts of the wave field that do not hit a -boundary at a right angle do not satisfy this condition and enforcing it as a -boundary condition will yield partial reflections, i.e. only parts of the wave -will pass through the boundary as if it wasn't here whereas the remaining -fraction of the wave will be reflected back into the domain. - -If we are willing to accept this as a suffient approximation to an absorbing boundary we finally arrive at the following problem for $u$: -@f{eqnarray*} --\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\ -c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\ -u &=& 1, \qquad x\in\Gamma_1. -@f} -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: -@f{eqnarray*} - \left.\begin{array}{ccc} - -\omega^2 v - c^2\Delta v &=& 0 \quad\\ - -\omega^2 w - c^2\Delta w &=& 0 \quad - \end{array}\right\} &\;& x\in\Omega, - \\ - \left.\begin{array}{ccc} - c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\ - c (n\cdot\nabla w) + \omega\,v &=& 0 \quad - \end{array}\right\} &\;& x\in\Gamma_2, - \\ - \left.\begin{array}{ccc} - v &=& 1 \quad\\ - w &=& 0 \quad - \end{array}\right\} &\;& x\in\Gamma_1. -@f} - -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 -@f{eqnarray*} --\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, \\ --\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. -@f} - -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 -@f[ -v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j. -@f] -Plugging into the variational form yields the equation system -@f[ -\renewcommand{\arraystretch}{2.0} -\left.\begin{array}{ccc} -\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 \\ -\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 -\end{array}\right\}\;\;\forall\; i =1,\ldots,n. -@f] -In matrix notation: -@f[ -\renewcommand{\arraystretch}{2.0} -\left( -\begin{array}{cc} --\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)} \\ -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)} -\end{array} -\right) -\left( -\begin{array}{c} -\alpha \\ \beta -\end{array} -\right) -= -\left( -\begin{array}{c} -0 \\ 0 -\end{array} -\right) -@f] -(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.) -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. -Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be -the same. However, we expect real and imaginary part of the solution to -have similar properties and will therefore indeed take $V_h=W_h$ in the -implementation, and also use the same basis functions $\phi_i = \psi_i$ for -both spaces. The reason for the notation using different symbols is just that -it allows us to distinguish between shape functions for $v$ and $w$, as this -distinction plays an important role in the implementation. - - -

The test case

- -For the computations, we will consider wave propagation in the unit square, -with ultrasound generated by a transducer lens that is shaped like a segment -of the circle with center at $(0.5, d)$ and a -radius slightly greater than $d$; this shape should lead to a focusing of the sound -wave at the center of the circle. Varying $d$ changes the "focus" of the lens -and affects the spatial distribution of the intensity of $u$, where our main -concern is how well $|u|=\sqrt{v^2+w^2}$ is.focussed. - -In the program below, we will implement the complex-valued Helmholtz equations -using the formulation with split real and imaginary parts. We will also -discuss how to generate a domain that looks like a square with a slight bulge -simulating the transducer (in the -UltrasoundProblem::make_grid() function), and how to -generate graphical output that not only contains the solution components $v$ and -$w$, but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in -UltrasoundProblem::output_results()). Finally, we use the -ParameterHandler class to easily read parameters like the focal distance $d$, -wave speed $c$, frequency $\omega$, and a number of other parameters from an -input file at run-time, rather than fixing those parameters in the source code -where we would have to re-compile every time we want to change parameters. + +This program was contributed by Moritz Allmaras at Texas A&M +University. Some of the work on this tutorial program has been funded +by NSF under grant DMS-0604778. + + +Note: In order to run this program, deal.II must be configured to use +the UMFPACK sparse direct solver. Refer to the ReadMe for instructions how to do this. + + + +

Introduction

+ + +A question that comes up frequently is how to solve problems involving complex +valued functions with deal.II. For many problems, instead of working with +complex valued finite elements directly, which are not readily available in +the library, it is often much more convenient to split complex valued +functions into their real and imaginary parts and use separate scalar finite +element fields for discretizing each one of them. Basically this amounts to +viewing a single complex valued equation as a system of two real valued +equations. This short example demonstrates how this can be implemented in +deal.II by using an FE_system object to stack two finite element +fields representing real and imaginary parts. We also revisit the +ParameterHandler class first used in @ref step_19 "step-19", which provides a +convenient way for reading parameters from a configuration file at runtime +without the need to recompile the program code. + +

Problem setting

+ +The original purpose of this program is to simulate the focussing properties +of an ultrasound wave generated by a transducer lens with variable +geometry. Recent applications in medical imaging use ultrasound waves not only +for imaging porposes, but also to excite certain local effects in a +material, like changes in optical properties, that can then be measured by +other imaging techniques. A vital ingredient for these methods is the ability +to focus the intensity of the ultrasound wave in a particular part of the +material, ideally in a point, to be able to examine the properties of the +material at that particular location. + +To derive a model for this problem, we think of ultrasound as a pressure wave +governed by the wave equation: +@f[ + \frac{\partial^2 U}{\partial t^2} - c^2 \Delta U = 0 +@f] +where $c$ is the wave speed (that for simplicity we assume to be constant), $U += U(x,t),\;x \in \Omega,\;t\in\mathrm{R}$. The boundary +$\Gamma=\partial\Omega$ is divided into two parts $\Gamma_1$ and +$\Gamma_2=\Gamma\setminus\Gamma_1$, with $\Gamma_1$ representing the +transducer lens and $\Gamma_2$ an absorbing boundary (that is, we want to +choose boundary conditions on $\Gamma_2$ in such a way that they imitate a +larger domain). On $\Gamma_1$, the transducer generates a wave of constant +frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here): +@f[ +U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1 +@f] + +If there are no other (interior or boundary) sources, and since the only +source has frequency $\omega$, then the solution admits a separation of +variables of the form $U(x,t) = \textrm{Re}\left(u(x)\,e^{i\omega +t})\right)$. The complex-valued function $u(x)$ describes the spatial +dependency of amplitude and phase (relative to the source) of the waves of +frequency ${\omega}$, with the amplitude being the quantity that we are +interested in. By plugging this form of the solution into the wave equation, +we see that for $u$ we have +@f{eqnarray*} +-\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\ +u(x) &=& 1, \qquad x\in\Gamma_1. +@f} + +For finding suitable conditions on $\Gamma_2$ that model an absorbing +boundary, consider a wave of the form $V(x,t)=e^{i(k\cdot x -\omega t)}$ with +frequency ${\omega}$ traveling in direction $k\in {\mathrm{R}^2}$. In order +for $V$ to solve the wave equation, $|k|={\frac{\omega}{c}}$ must +hold. Suppose that this wave hits the boundary in $x_0\in\Gamma_2$ at a right +angle, i.e. $n=\frac{k}{|k|}$ with $n$ denoting the outer unit normal of +$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation +@f[ +c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0. +@f] +Hence, by enforcing the boundary condition +@f[ +c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2, +@f] +waves that hit the boundary $\Gamma_2$ at a right angle will be perfectly +absorbed. On the other hand, those parts of the wave field that do not hit a +boundary at a right angle do not satisfy this condition and enforcing it as a +boundary condition will yield partial reflections, i.e. only parts of the wave +will pass through the boundary as if it wasn't here whereas the remaining +fraction of the wave will be reflected back into the domain. + +If we are willing to accept this as a suffient approximation to an absorbing boundary we finally arrive at the following problem for $u$: +@f{eqnarray*} +-\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\ +c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\ +u &=& 1, \qquad x\in\Gamma_1. +@f} +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: +@f{eqnarray*} + \left.\begin{array}{ccc} + -\omega^2 v - c^2\Delta v &=& 0 \quad\\ + -\omega^2 w - c^2\Delta w &=& 0 \quad + \end{array}\right\} &\;& x\in\Omega, + \\ + \left.\begin{array}{ccc} + c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\ + c (n\cdot\nabla w) + \omega\,v &=& 0 \quad + \end{array}\right\} &\;& x\in\Gamma_2, + \\ + \left.\begin{array}{ccc} + v &=& 1 \quad\\ + w &=& 0 \quad + \end{array}\right\} &\;& x\in\Gamma_1. +@f} + +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 +@f{eqnarray*} +-\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, \\ +-\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. +@f} + +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 +@f[ +v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j. +@f] +Plugging into the variational form yields the equation system +@f[ +\renewcommand{\arraystretch}{2.0} +\left.\begin{array}{ccc} +\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 \\ +\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 +\end{array}\right\}\;\;\forall\; i =1,\ldots,n. +@f] +In matrix notation: +@f[ +\renewcommand{\arraystretch}{2.0} +\left( +\begin{array}{cc} +-\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)} \\ +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)} +\end{array} +\right) +\left( +\begin{array}{c} +\alpha \\ \beta +\end{array} +\right) += +\left( +\begin{array}{c} +0 \\ 0 +\end{array} +\right) +@f] +(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.) +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. +Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be +the same. However, we expect real and imaginary part of the solution to +have similar properties and will therefore indeed take $V_h=W_h$ in the +implementation, and also use the same basis functions $\phi_i = \psi_i$ for +both spaces. The reason for the notation using different symbols is just that +it allows us to distinguish between shape functions for $v$ and $w$, as this +distinction plays an important role in the implementation. + + +

The test case

+ +For the computations, we will consider wave propagation in the unit square, +with ultrasound generated by a transducer lens that is shaped like a segment +of the circle with center at $(0.5, d)$ and a +radius slightly greater than $d$; this shape should lead to a focusing of the sound +wave at the center of the circle. Varying $d$ changes the "focus" of the lens +and affects the spatial distribution of the intensity of $u$, where our main +concern is how well $|u|=\sqrt{v^2+w^2}$ is.focussed. + +In the program below, we will implement the complex-valued Helmholtz equations +using the formulation with split real and imaginary parts. We will also +discuss how to generate a domain that looks like a square with a slight bulge +simulating the transducer (in the +UltrasoundProblem::make_grid() function), and how to +generate graphical output that not only contains the solution components $v$ and +$w$, but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in +UltrasoundProblem::output_results()). Finally, we use the +ParameterHandler class to easily read parameters like the focal distance $d$, +wave speed $c$, frequency $\omega$, and a number of other parameters from an +input file at run-time, rather than fixing those parameters in the source code +where we would have to re-compile every time we want to change parameters. -- 2.39.5