]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step-81: minor typos fixes, default parameters changed
authorManaswinee Bezbaruah <bezba004@tamu.edu>
Tue, 2 Nov 2021 04:57:28 +0000 (23:57 -0500)
committerMatthias Maier <tamiko@43-1.org>
Thu, 26 May 2022 05:16:02 +0000 (00:16 -0500)
doc/doxygen/references.bib
examples/step-81/CMakeLists.txt
examples/step-81/doc/intro.dox
examples/step-81/doc/intro.dox-template [deleted file]
examples/step-81/doc/results.dox-template [deleted file]
examples/step-81/step-81.cc

index 6f97c4e3706c9ccbafb1107ef1e1d93a580e80cf..d042c7b9849f4d5d12a6cc4d2d5fad8e8d836876 100644 (file)
   url =       {https://doi.org/10.1023/a:1020533003783}
 }
 
+% ------------------------------------
+% Step 81
+% ------------------------------------
+
+@BOOK{Schwartz1972,
+  AUTHOR={M. Schwartz},
+  TITLE={Principles of Electrodynamics},
+  Year={1972},
+  PUBLISHER={McGraw-Hill Book Company},
+  ADDRESS={New York},
+  SERIES={International Series in Pure and Applied Physics},
+}
+
+@BOOK{Monk2003,
+  AUTHOR={P. Monk},
+  TITLE={Finite Element Methods for {Maxwell's} Equations},
+  YEAR={2003},
+  PUBLISHER={Oxford University Press},
+  SERIES={Numerical Mathematics and Scientific Computation},
+}
+
+@ARTICLE{Geim2004,
+  AUTHOR={K.S. Novoselov and A.K. Geim, S.V. Morozov and D. Jiang and Y. Zhang and S.V. Dubonos and I.V. Grigorieva and A.A. Firsov},
+  TITLE={Electric Field Effect in Atomically Thin Carbon Films},
+  JOURNAL={Science},
+  VOLUME={306},
+  PAGES={666},
+  YEAR={2004}
+}
+
+@ARTICLE{Maier2017,
+    author = {M. Maier and D. Margetis and M. Luskin},
+    journal = {Journal of Computational Physics},
+    pages = {126--145},
+    title = {Dipole excitation of surface plasmon on a conducting sheet: finite element approximation and validation},
+    volume = {339},
+    year = {2017}
+}
+
 % ------------------------------------
 % Step 82
 % ------------------------------------
index c8e17bb6865d45078dcbebf4b48ca7bbdfd5ada0..d28e396a5a6dd9fa4499af6e607a754af0dce9b3 100644 (file)
@@ -23,7 +23,7 @@ SET(TARGET_SRC
 
 CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0)
 
-FIND_PACKAGE(deal.II 10.0.0
+FIND_PACKAGE(deal.II 9.2.0
   HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
   )
 IF(NOT ${deal.II_FOUND})
index a4e114dc951c007b8a88742480dd0eda21ffa55b..b9e5d19196f7c5446928775262839db0656d6b5a 100644 (file)
 <a name="Intro"></a>
 
-<h1>Introduction</h1>
-
-A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic wave,
-excited along a conducting sheet by an electric Hertzian dipole. Suitably
-rescaled time-harmonic Maxwell's equations can be used to derive a variational
-form, which in turn enables a numerical observation of these SPPs by appropriate
-curl-conforming finite elements. The conducting sheet is modeled as an idealized
-hypersuface with an effective electric conductivity, and the weak discontinuity
-for the tangential surface appears naturally in the variational form. <br />
-
-The following tutorial is a direct solver for the 2D time-harmonic Maxwell
-equations describing a scattering configuration on a lower-dimensional
-interface (with absorbing impedance boundary conditions). These variational
-equations aim to numerically simulate a SPP on an infite sheet with with constant
-isotropic conductivity embedded in two spatial dimensions. Following is a
-detailed discussion on the derivation of the variational form and the appropriate
-boundary conditions. <br />
-
-<h2> Defining the Problem </h2>
-
-<h3> Time-Harmonic Maxwell's Equations </h3>
-
-Consider an electromagnetic wave $(\mathbf{E},\mathbf{H})$ in a surface
-$\Omega\backslash\Sigma \subset \mathbb{R}^n$, where $n=2$ or $3$.
-Assume all material parameters are time-independent.
-
-The Maxwell's equations for this wave are,
+<h1> Introduction </h1>
+
+A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic
+wave, confined near a metal-air (or similar) interfaces. SPP structures on
+novel "2D" materials such as graphene, a monoatomic layer of carbon atoms
+arranged in a hexagonal lattice, typically have wavelengths much shorter
+than the wavelength of the free-space radiation. This scale separation
+makes SPPs on 2D materials a promising ingredient in the design of novel
+optical devices.
+
+In the following, we discuss a method for observing SPPs numerically by
+solving a suitable electromagnetic model based on time-harmonic Maxwell's
+equations that incorporates jump conditions on lower-dimensional material
+interfaces: The conducting sheet is modeled as an idealized hypersurface
+with an effective electric conductivity, and the weak discontinuity for the
+tangential surface appears naturally in the variational form.
+
+This tutorial presents a direct solver for the time-harmonic Maxwell
+equations for scattering configurations with lower-dimensional interfaces.
+We discuss in particular how to set up a complex-valued (time-harmonic),
+how to implement simple first-order absorbing boundary conditions and a
+more sophisticated "perfectly matched layer" for electromagnetic waves.
+
+
+<h2>Time-Harmonic Maxwell's Equations with interface conditions</h2>
+
+We start the discussion with a short derivation of the governing equations
+and some pointers to literature.
+
+
+<h3>Derivation of time-harmonic Maxwell's equations</h3>
+
+In two ($d=2$) or three ($d=3$) spatial dimensions,
+the time evolution of an electromagnetic
+wave $(\mathbf{E},\mathbf{H})$ that consists of an electric field component
+$\mathbf{E}(t,\mathbf{x})\;:\;\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}^d$
+and a magnetic field component
+$\mathbf{H}(t,\mathbf{x})\;:\;\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}^d$
+is described by
+<a href="https://en.wikipedia.org/wiki/Maxwell%27s_equations">Maxwell's
+equations</a>
+@cite Schwartz1972, @cite Monk2003 :
 @f[
 \begin{cases}
--i\omega \mathbf{H} + \nabla \times \mathbf{E} = -\mathbf{M}_a,\\
-\nabla \cdot \mathbf{H} = \frac{1}{i\omega}\nabla \cdot \mathbf{M}_a,\\
-i\omega\varepsilon\mathbf{E} + \nabla\times(\mu^{-1}\mathbf{H}) = \mathbf{J}_a,\\
-\nabla\cdot(\varepsilon\mathbf{E}) = \frac{1}{i\omega}\nabla\cdot\mathbf{J}_a.
+\frac{\partial}{\partial t} \mathbf{H} + \nabla \times \mathbf{E} = -\mathbf{M}_a,
+\\
+\nabla \cdot \mathbf{H} = \rho,
+\\
+\frac{\partial}{\partial t} (\varepsilon\mathbf{E}) - \nabla\times(\mu^{-1}\mathbf{H}) = - \mathbf{J}_a,
+\\
+\nabla\cdot(\varepsilon\mathbf{E}) = \rho_m,
 \end{cases}
 @f]
+where $\nabla\times\mathbf{F}(\mathbf{x})$ denotes the curl and
+$\nabla\cdot\mathbf{F}(\mathbf{x})$ denotes the divergence of a vector
+field $\mathbf{F}:\Omega\to\mathbb{R}^d$ and where we have set $d=2,3$. We
+have introduced two (time-independent) material parameters, the
+<a href="https://en.wikipedia.org/wiki/Permittivity">electric permittivity</a>
+$\varepsilon$
+and the
+<a href="https://en.wikipedia.org/wiki/Permeability">magnetic permeability</a>
+$\mu$. In addition, $\rho$ is the (electric) charge density and $\rho_m$ is
+a corresponding (hypothetical)
+<a href="https://en.wikipedia.org/wiki/Magnetic_monopole">magnetic monopole</a>
+density. $\mathbf{J}_a$ and $\mathbf{M}_a$ are the electric and magnetic
+flux densities. Both are related to their respective charge densities by a
+conservation equation @cite Schwartz1972 :
+@f[
+\frac{\partial}{\partial t} \rho + \nabla\cdot\mathbf{J}_a \,=\, 0,
+\qquad
+\frac{\partial}{\partial t} \rho_m + \nabla\cdot\mathbf{M}_a \,=\, 0.
+@f]
 
-Here, the positive parameter $\omega$ is the temporal angular frequency.
-$\mathbf{J}_a$ and $\mathbf{M}_a$ are time-independent externally applied
-electric-current and magnetic-current densities respectively, arising from the
-time-harmonic densities $\mathcal{J}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{J}_a(x)\}$
-and $\mathcal{M}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{M}_a(x)\}$. Moreover,
-$\varepsilon(x)$ and $\mu(x)$ are rank $2$ tensors representing the complex
-permittivity and the relative magnetic permeability of the corresponding medium.
-$\varepsilon = \varepsilon_0(x) + i\sigma(x)/\omega$, where $\varepsilon_0(x)$ is
-the dielectric permittivity and $\sigma(x)$ is the surface conductivity.<br />
-We assume some (weak) regularity of the $x$ dependent variables $(\mathbf{E}$,
-$\mathbf{H})$, $(\mathbf{J}_a,\mathbf{M}_a)$, and $(\varepsilon, \mu)$ to ensure
-well-posedness.
-
-The surface conductivity $\sigma$ gives rise to a current density, which in turn
-gives rise to a jump conditions on $\Sigma$ in the tangential component (away
-from the boundary) of the magnetic field. The tangential electric field is
-continuous. On the idealized, oriented hypersurface $\Sigma \subset \mathbb{R}^n$,
-with unit normal $\nu$ and effective surface conductivity $\sigma^{\Sigma}$, this
-is modelled as,
-
+We now make the important assumption that the material parameters
+$\varepsilon$ and $\mu$ are time-independent and that the fields
+$\mathbf{E}$ and $\mathbf{H}$, the fluxes $\mathbf{M}_a$ and
+$\mathbf{J}_a$, as well as the densities $\rho$ and $\rho_m$ are all
+<i>time-harmonic</i>, i.e., their time evolution is completely described by
+@f[
+  \mathbf{F}(\mathbf{x},t) = \text{Re}\{e^{-i\omega
+  t}\tilde{\mathbf{F}}(\mathbf{x})\},
+@f]
+where $\omega$ is the temporal angular frequency and
+$\tilde{\mathbf{F}}(\mathbf{x})$ is a corresponding complex-valued vector
+field (or density). Inserting this ansatz into Maxwell's equations,
+substituting the charge conservation equations and some minor algebra then
+yields the so-called <i>time-harmonic</i> Maxwell's equations, viz.,
 @f[
 \begin{cases}
-\nu \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma}
-= \sigma^{\Sigma}\left[(\nu\times \mathbf{E}\times \nu)\right]|_{\Sigma},\\
-\nu \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} = 0.
+-i\omega \tilde{\mathbf{H}} + \nabla \times \tilde{\mathbf{E}} =
+-\tilde{\mathbf{M}}_a,\\
+\nabla \cdot \tilde{\mathbf{H}} = \frac{1}{i\omega}\nabla \cdot
+\tilde{\mathbf{M}}_a,\\
+i\omega\varepsilon\tilde{\mathbf{E}} +
+\nabla\times(\mu^{-1}\tilde{\mathbf{H}}) = \tilde{\mathbf{J}}_a,\\
+\nabla\cdot(\varepsilon\tilde{\mathbf{E}}) =
+\frac{1}{i\omega}\nabla\cdot\tilde{\mathbf{J}}_a.
 \end{cases}
 @f]
 
-Moreover, the above equations are supplemented by the Silver-Müller radiation
-condition, if the ambient (unbounded) medium is isotropic. This amounts to the
-requirement that $\mathbf{E}, \mathbf{H}$ approach a spherical wave uniformly in
-the radial direction for points at infinity and away from the conducting sheet.
-
+For the sake of better readability we will now drop the tilde and simply
+write $\mathbf{E}(\mathbf{x})$, $\mathbf{H}(\mathbf{x})$, etc., when
+referring to the time-harmonic fields.
+
+
+<h3>Jump conditions on lower dimensional interfaces</h3>
+
+Graphene is a two-dimensional carbon allotrope with a <i>single</i> atom
+layer that is arranged in a honeycomb lattice @cite Geim2004. Due to its
+atomic thickness it is an example of a so-called 2D material: Compared to
+the other spatial dimensions (where graphene samples can reach up to
+several centimeters) the atomistic thickness of graphene typically ranges
+around 2.5 ångstrom ($2.5\times10^{-10}\text{m}$). We will thus model
+graphene as a lower-dimensional interface $\Sigma$ imbedded into the
+computational domain $\Omega\subset\mathbb{R}^d$. More precisely, $\Sigma$
+is a two-dimensional sheet in three spatial dimensions, or a
+one-dimensional line in two spatial dimensions. The special electronic
+structure of graphene gives rise to a current density on the
+lower-dimensional interface that is modeled with an effective surface
+conductivity $\sigma^\Sigma$ obeying <a
+href="https://en.wikipedia.org/wiki/Ohm%27s_law">Ohm's Law</a>, viz,
 @f[
-\lim\limits_{|x|\to\infty} \{\mathbf{H}\times x - c^{-1}|x|\mathbf{E}\} = 0;\qquad
-\lim\limits_{|x|\to\infty} \{\mathbf{E}\times x - c^{-1}|x|\mathbf{H}\} = 0;\qquad
-x \not\in \Sigma
+  \mathbf{J}^\Sigma=\sigma^\Sigma\,\mathbf{E}_T.
 @f]
+Here, $\mathbf{J}^\Sigma$ is the surface current density, $\mathbf{E}_T$
+denotes the tangential part of the electric field $\mathbf{E}$, and
+$\sigma^\Sigma$ is an appropriately chosen surface conductivity that will
+be discussed in more detail below. The surface current density gives rise
+to a jump condition on $\Sigma$ in the tangential component of the magnetic
+field. This is best seen by visualizing <a
+href="https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law">Ampère's
+law</a>,
+
+@f[\text{TODO}@f]
+
+and then taking the limit of the upper and lower part of the line integral
+approaching the sheet. In contrast, the tangential part of the electric
+field is continuous. By fixing a unit normal $\mathbf{\nu}$ on the hypersurface
+$\Sigma$ both jump conditions read,
+@f[
+\begin{cases}
+\mathbf{\nu} \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma}
+= \sigma^{\Sigma}\left[(\mathbf{\nu}\times \mathbf{E}\times \mathbf{\nu})\right]|_{\Sigma},\\
+\mathbf{\nu} \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} = 0.
+\end{cases}
+@f]
+Here, the notation $\mathbf{F}^\pm$ indicates the limit values of the field
+when approaching the interface from above or below the interface:
+$\mathbf{F}^\pm(\mathbf{x})=\lim_{\delta\to0,\delta>0}\mathbf{F}(\mathbf{x}\pm\delta\mathbf{\nu})$.
 
-In our case, we eliminate reflection from infinity by implementing a PML and
-avoid the explicit use of the last condition.
 
 <h3> Rescaling </h3>
 
@@ -236,6 +309,23 @@ Find a unique $\mathbf{E} \in X(\Omega)$ such that for all $\varphi \in X(\Omega
 A(\mathbf{E},\varphi) = F(\varphi)
 @f]
 
+
+<h3>Absorbing boundary conditions and perfectly matched layer</h3>
+
+Moreover, the above equations are supplemented by the Silver-Müller radiation
+condition, if the ambient (unbounded) medium is isotropic. This amounts to the
+requirement that $\mathbf{E}, \mathbf{H}$ approach a spherical wave uniformly in
+the radial direction for points at infinity and away from the conducting sheet.
+
+@f[
+\lim\limits_{|x|\to\infty} \{\mathbf{H}\times x - c^{-1}|x|\mathbf{E}\} = 0;\qquad
+\lim\limits_{|x|\to\infty} \{\mathbf{E}\times x - c^{-1}|x|\mathbf{H}\} = 0;\qquad
+x \not\in \Sigma
+@f]
+
+In our case, we eliminate reflection from infinity by implementing a PML and
+avoid the explicit use of the last condition.
+
 <h2> Discretization Scheme</h2>
 
 The variational form is discretized on a non-uniform quadrilateral mesh with
diff --git a/examples/step-81/doc/intro.dox-template b/examples/step-81/doc/intro.dox-template
deleted file mode 100644 (file)
index 5e67a7f..0000000
+++ /dev/null
@@ -1,252 +0,0 @@
-<a name="Intro"></a>
-
-<h1>Introduction</h1>
-
-A surface plasmon-polariton (SPP) is a slowly decaying electromagnetic wave,
-excited along a conducting sheet by an electric Hertzian dipole. Suitably
-rescaled time-harmonic Maxwell's equations can be used to derive a variational
-form, which in turn enables a numerical observation of these SPPs.
-
-The following tutorial is a direct solver for the 2D time-harmonic Maxwell
-equations describing a scattering configuration on a lower-dimensional
-interface (with absorbing impedance boundary conditions).<br>
-
-<h3> Defining the Problem </h3>
-
-<h4> Time-Harmonic Maxwell's Equations </h4>
-
-Consider an electromagnetic wave $(\mathbf{E},\mathbf{H})$ in a surface
-$\Omega\backslash\Sigma \subset \mathbb{R}^n$, where $n=2 or 3.
-Assume all material parameters are time-independent.
-The Maxwell's equations for this wave are,
-
-\f{align*}{
-\begin{cases}
--i\omega \mathbf{H} + \nabla \times \mathbf{E} = -\mathbf{M}_a,\\
-\nabla \cdot \mathbf{H} = \frac{1}{i\omega}\nabla \cdot \mathbf{M}_a,\\
-i\omega\varepsilon\mathbf{E} + \nabla\times(\mu^{-1}\mathbf{H}) = \mathbf{J}_a,\\
-\nabla\cdot(\varepsilon\mathbf{E}) = \frac{1}{i\omega}\nabla\cdot\mathbf{J}_a.
-\f}{align*}<br>
-
-Here, the positive parameter $\omega$ is the temporal angular frequency.
-$\mathbf{J}_a$ and $\mathbf{M}_a$ are time-independent externally applied
-electric-current and magnetic-current densities respectively, arising from the
-time-harmonic densities $\mathcal{J}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{J}_a(x)$
-and $\mathcal{M}_a(x,t) = \text{Re}\{e^{-i\omega t}\mathbf{M}_a(x)$. Moreover,
-$\varepsilon(x)$ and $\mu(x)$ are rank 2 tensors representing the complex permittivity
-and the relative magnetic permeability of the corresponding medium.
-$\varepsilon = \varepsilon_0(x) + i\sigma(x)/\omega$, where $\varepsilon_0(x)$ is
-the dielectric permittivity and $\sigma(x)$ is the surface conductivity.<br>
-We assume some (weak) regularity of the $x$ dependent variables $(\mathbf{E}$,
-$\mathbf{H})$, $(\mathbf{J}_a,\mathbf{M}_a)$, and $(\varepsilon, \mu)$ to ensure
-well-posedness.
-
-The surface conductivity $\sigma$ gives rise to a current density, which in turn
-gives rise to a jump conditions on $\Sigma$ in the tangential component (away
-from the boundary) of the magnetic field. The tangential electric field is
-continuous. On the idealized, oriented hypersurface $\Sigma \subset \mathbb{R}^n$,
-with unit normal $\nu$ and effective surface conductivity $\sigma^{\Sigma}$, this
-is modelled as,
-\f{align*}{
-\begin{cases}
-\nu \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma}
-= \sigma^{\Sigma}\left[(\nu\times \mathbf{E}\times \nu)\right]|_{\Sigma},\\
-\nu \times \left[\mathbf{E}^+ - \mathbf{E}^-\right]|_{\Sigma} = 0.
-\end{cases}
-\f}{align*}
-
-Moreover, the above equations are supplemented by the Silver-Müller radiation
-condition, if the ambient (unbounded) medium is isotropic. In our case, we eliminate
-reflection from infinity by implementing a PML and avoid the explicit use of this
-condition.
-
-<h4> Rescaling </h4>
-
-We will be using a rescaled version of the Maxwell's equations where:<br>
-1. Every length is rescaled by the free-space wavelength $2\pi k^{-1}
-:= 2\pi(\omega\sqrt{\varepsilon_0\mu_0})^{-1}, where $\varepsilon_0$ and $\mu_0$
-denote the vacuum dielectric permittivity and magnetic permeability, respectively.<br>
-2. $\mathbf{E}$, $\mathbf{H}$, $\mathbf{J}_a$, $\mathbf{M}_a$ are all rescaled by
-typical electric current strength $J_0$, where $J_0$ is the strength of the
-prescribed dipole source at location $a$ in the $e_i$ direction in Cartesian
-coordinates.<br>
-
-We introduce the rescaled variables $\mu_r$, $\varepsilon_r$, $sigma^{\Sigma}_r$,
-$\hat{x}$, $\hat{\mathbf{H}}$, $\hat{\mathbf{E}}$, $\hat{\mathbf{J}}_a$,
-$\hat{\mathbf{M}}_a$.
-
-\f{\begin{align*}
-\begin{cases}
-\mu_r = \frac{1}{\mu_0}\mu\\
-\varepsilon_r = \frac{1}{\varepsilon_0}\varepsilon\\
-\sigma^{Sigma}_r = \sqrt{\frac{\mu_0}{\varepsilon_0}}\sigma^{\Sigma},\\
-\hat{x} = k_0x\\
-\hat{\mathbf{H}} = \frac{k_0}{J_0}\mu^{-1}\mathbf{H}\\
-\hat{\mathbf{E}} = \frac{k_0^2}{\omega\mu_0 J_0}\mathbf{E}\\
-\hat{\mathbf{J}}_a = \frac{1}{J_0}\mathbf{J}_a\\
-\hat{\mathbf{M}}_a = \frac{k_0}{\omega\mu_0 J_0}\mathbf{M}_a
-\end{cases}
-\end{align*}}\f
-
-Accordingly, our rescaled equations are:
-\f{align*}{
-\begin{cases}
--i\mu_r \hat{\mathbf{H}} + \hat{\nabla} \times \hat{\mathbf{E}}
-= -\hat{\mathbf{M}}_a,\\
-\hat{\nabla} \cdot (\mu_r\hat{\mathbf{H}}) = \frac{1}{i}\hat{\nabla}
-\cdot \hat{\mathbf{M}}_a,\\
-i\varepsilon_r\hat{\mathbf{E}} + \nabla\times(\mu^{-1}\mathbf{H})
-= \mathbf{J}_a,\\
-\nabla\cdot(\varepsilon\mathbf{E}) = \frac{1}{i\omega}\hat{\nabla}
-\cdot\hat{\mathbf{J}}_a.
-\f}{align*}<br>
-We will omit the hat further discuss for ease of notation.
-
-<h4> Variational Form</h4>
-
-Let $\Omega \subset \mathbb{R}^n$ be a simply connected and bounded
-domain with Lipschitz-continuous and piecewise smooth boundary, $\partial\Omega$.
-Let $\Sigma$ be an oriented, Lipschitz-continuous, piecewise smooth hypersurface.
-Fix a normal field $\nu$ on $\Sigma$ and let $n$ denote the outer normal vector
-on $\partial\Omega$.<br>
-
-In order to arrive at the variational form, we will substitute $\mathbf{H}$ in
-the first equation as follows:
-\f{\begin{align*}
-\nabla \times (\mu_r^{-1}\nabla\times\mathbf{E}) - \varepsilon_r \mathbf{E}
-= i\mathbf{J}_a - \nabla\times (\mu_r^{-1}mathbf{M}_a)
-\end{align*}}\f <br>
-
-Now consider a smooth test function $\varphi$ with complex conjugate $\bar{\varphi}$.
-Multiply both sides of the above equation by $\bar{\varphi}$ and integrate by parts
-in $\Omega\backslash\Sigma$. Moreover, we use the subscript $T$ to denote the
-tangential part of the given vector i.e. $F_T = (\nu\times F)\times\nu. We arrive
-at:
-\f{\begin{align*}
-\int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot (\nabla\times\bar{\varphi})\text{d}x
-- \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\text{d}x
-- i\int_\Sigma (\sigma_r^{\Sigma}(\mathbf{E})_T) \cdot (\bar{\varphi})_T\text{do}x
-- i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\mathbf{E})_T) \cdot
-(\nabla\times(\bar{\varphi})_T)\text{d}x =
-i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\text{d}x
-- \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi}) \text{d}x.
-\end{align*}}\f <br>
-
-Define
-\f{\begin{align*}
-A(\mathbf{E},\varphi) := \int_\Omega (\mu_r^{-1}\nabla\times\mathbf{E})\cdot
-(\nabla\times\bar{\varphi})\text{d}x
-- \int_\Omega \varepsilon_r\mathbf{E} \cdot \bar{\varphi}\text{d}x
-- i\int_\Sigma (\sigma_r^{\Sigma}(\mathbf{E})_T) \cdot (\bar{\varphi})_T\text{do}x
-- i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\mathbf{E})_T) \cdot
-(\nabla\times(\bar{\varphi})_T)\text{d}x.\\
-F(\varphi) := i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi}\text{d}x
-- \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi}) \text{d}x.
-\end{align*}}\f<br>
-
-Then, our rescaled weak formulation is:<br>
-Find a unique $\mathbf{E} \in X(\Omega)$ such that for all $\varphi \in $X(\Omega)$
-\f{\begin{align*}
-A(\mathbf{E},\varphi) = F(\varphi)
-\end{align*}}\f<br>
-
-<h4> Discretization Scheme</h4>
-
-The variational form is discretized on a non-uniform quadrilateral mesh with
-higher-order, curl-conforming Nédélec elements. This way the interface with a
-weak discontinuity can be aligned with or away from the mesh, and the convergence
-rate is high. Specifically, we use second-order Nédélec elements, which under our
-conditions will have a convergence rate $\mathcal{O}(#dofs)$.
-
-Now, consider the finite element subspace $X_h(\Omega) \subset X(\Omega)$. Define
-the matrices
-\f{align*}{
-  A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
-          (\nabla\times\bar{\varphi}_j)\text{d}x
-          - \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
-          - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot
-          (\bar{\varphi}_j)_T\text{do}x
-          - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T)
-          \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x,
-\f}
-\f{align}{
-  F_i = i\int_\Omega \mathbf{J}_a \cdot \bar{\varphi_i}\text{d}x
-      - \int_\Omega \mu_r^{-1}\mathbf{M}_a \cdot (\nabla \times \bar{\varphi_i})
-      \text{d}x.
-\f} <br>
-Then under the assumption of a sufficiently refined initial mesh
-the discretized variational problem is: Find a $\varphi_j \in \X_h(\Omega)$ such
-that for all $\varphi_i \in X_h(\Omega)$:
-\f{\begin{align*}
-A_{ij} = F_i
-\end{align*}}\f <br>
-
-Using a skeleton similar to step-4, we have constructed a Maxwell class and we
-have used complex-valued FENedelec elements to solve our equations. <br>
-The material coefficients such as $\mu_r^{-1}$, $J_a$, etc. are instantiated
-using the Parameters class, which is also documented below.
-
-<h4> Perfectly Matched Layer </h4>
-The SPP amplitude is negatively effected by the absorbing boundary condition and
-this causes the solution image to be distorted. In order to reduce the resonance
-and distortion in our solutions, we are implementing a Perfectly Matched Layer
-(PML) in the scattering configuration. <br>
-
-The concept of a PML was pioneered by Bérenger and it is is an indispensable tool
-for truncating unbounded domains for wave equations and often used in the
-numerical approximation of scattering problems. It is essentially a thin layer with
-modified material parameters placed near the boundary such that all outgoing
-electromagnetic waves decay exponentially with no “artificial” reflection due to
-truncation of the domain. <br>
-
-Our PML is essentially a concentric circle with modified material coefficients
-($\varepsilon_r, \mu_r, \sigma$). It is located in a small region near the boundary
-$\partial\Omega$ and the transformation of the material coordinates is chosen to
-be a function of the radial distance $\rho$ from the origin $e_r$. The normal field
-$\nu$ of $\Sigma$ is orthogonal to the radial direction $e_r$, which makes
-$\mathbf{J}_a \equiv 0$ and $\mathbf{M}_a \equiv 0$ within the PML. <br>
-
-\\TODO: Insert image of the PML
-
-Introduce a change of coordinates
-\f{\begin{align*}
-x \to \bar{x} =
-\begin{cases}
-x + ie_r\int\limits_\rho^r s(\tau)\text{d}\tau,\;\;\;\;\;\;\; r\ge\rho
-x\;\;\;\;\;\;\;\;\;\otherwise
-\end{cases}
-\end{align*}}\f<br>
-
-where $r = e_r \cdot x$ and $s(\tau)$ is an appropriately chosen, nonnegative
-scaling function.<br>
-
-We introduce the following $2\times2$ matrices
-\f{\begin{align*}
-A = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2},\frac{1}{d\bar{d}}\right)T_{e_xe_r},\\
-B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\\
-C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right)T_{e_xe_r}.\\
-\end{align*}}\f <br>
-
-where
-\f{\begin{align*}
-d = 1 + is(r),\\
-\bar{d} = 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau.
-\end{align*}}\f<br>
-
-and $T_{e_xe_r}$ is the rotation matrix that rotates $e_r$ onto $e_x$.<br>
-
-Thus, after applying the rescaling we get the following modified parameters
-\f{\begin{align*}
-\bar{\mu}_r^{-1} = \frac{\mu_r^{-1}}{d},\\
-\bar{\varepsilon}_r = A^{-1} \varepsilon_r B^{-1},\\
-\bar{\sigma}^{\Sigma}_r = C^{-1} \sigma^{\Sigma}_r B^{-1}.
-\end{align*}}\f
-
-These PML transformations are implement in our PMLParameters class. After the PML
-is implemented, the electromagnetic wave essentially decays exponentially within
-the PML region near the boundary, therefore reducing reflection from the boundary
-of our finite domain. The decay function also depends on the strength of our PML,
-which can be adjusted to see more or less visible decaying in the PML region.<br>
-
-We also add an interface at $y = 0$ in our domain, in order to test the
-validity of our PML.<br>
diff --git a/examples/step-81/doc/results.dox-template b/examples/step-81/doc/results.dox-template
deleted file mode 100644 (file)
index d195483..0000000
+++ /dev/null
@@ -1,15 +0,0 @@
-<h1>Results</h1>
-
-Using the above code (without the PML) and a forcing term of a Hertzian dipole
-at the center, we have generated the following solution wave. The two complex
-plane solutions are followed by the two solutions in the real plane.
-Furthermore, a significant resonance is observes, causing distorted images
-and necessitating a PML boundary condition
-TODO: add the images
-The solution of the same problem with a PML of strength 2 of radii 8 and 10
-is shown below (in the same order or complex and real solutions). Clearly,
-the PML significantly reduces the distortion in our solution.
-TODO: add the images
-Additionally, an interface is added at y = 0, and by adjusting the surface
-conductivity value and the position of our dipole, we arrive at a standing wave
-TODO: add the images
index 12fb81427882dc626e737b6753d8419f473281db..9d6b7fe62ceed70137b16b7ba180cf4d22fbea75 100644 (file)
@@ -66,7 +66,9 @@
 #include <memory>
 
 
-// @sect3{Class Template Declartions}
+// @sect3{Class Template Declarations}
+// We begin our actual implementation by declaring all classes with their
+// data structures and methods upfront.
 
 namespace Step81
 {
@@ -82,24 +84,20 @@ namespace Step81
    // More explanation on the use and inheritance from the ParameterAcceptor
    // can be found in step-60.
   
-
-/**
-   *
-   * epsilon is the Electric Permitivitty coefficient and it is a rank 2 tensor. Depending on the material,
-   * we assign the i^th diagonal element of the tensor to the material epsilon value
-   * (one of the private epsilon_1_ or epsilon_2_ variables).
-   *
-   * mu_inv  is the inverese of the Magnetic Permiabillity coefficient and it is a complex number.
-   *
-   * sigma is the Surface Conductivity coefficient between material left and material right
-   * and it is a rank 2 tensor. It is only changed if we are at the interface between two
-   * materials. If we are at an interface, we assign the i^th diagonal element of the
-   * tensor to the private sigma_ value.
-   *
-   * J_a is the strength and orientation of the dipole. It is a rank 1 tensor that depends
-   * on the private dipole_position_, dipole_radius_, dipole_strength_, dipole_orientation_
-   * variables.
-*/
+   // epsilon is the Electric Permitivitty coefficient and it is a rank 2 tensor. Depending on the material,
+   // we assign the i^th diagonal element of the tensor to the material epsilon value
+   // (one of the private epsilon_1_ or epsilon_2_ variables).
+   //
+   // mu_inv  is the inverese of the Magnetic Permiabillity coefficient and it is a complex number.
+   
+   // sigma is the Surface Conductivity coefficient between material left and material right
+   // and it is a rank 2 tensor. It is only changed if we are at the interface between two
+   // materials. If we are at an interface, we assign the i^th diagonal element of the
+   // tensor to the private sigma_ value.
+   
+   // J_a is the strength and orientation of the dipole. It is a rank 1 tensor that depends
+   // on the private dipole_position_, dipole_radius_, dipole_strength_, dipole_orientation_
+   // variables.
 
   template <int dim>
   class Parameters : public ParameterAcceptor
@@ -233,9 +231,12 @@ namespace Step81
   }
 
    // @sect4{PerfectlyMatchedLayer Class}
-   // The PerfectlyMatchedLayer class inherits ParameterAcceptor, and it modifies our coefficients from Parameters.
-   // The radii and the strength of the PML is specified, and the coefficients will be modified using transformation
-   // matrices within the PML region. The radii and strength of the PML are editable through a .prm file
+   // The PerfectlyMatchedLayer class inherits ParameterAcceptor,
+   // and it modifies our coefficients from Parameters.
+   // The radii and the strength of the PML is specified, and the
+   // coefficients will be modified using transformation
+   // matrices within the PML region. The radii and strength of
+   // the PML are editable through a .prm file
 
   template <int dim>
   class PerfectlyMatchedLayer : public ParameterAcceptor
@@ -284,7 +285,7 @@ namespace Step81
     add_parameter("outer radius",
                   outer_radius,
                   "outer radius of the PML shell");
-    strength = 0.;
+    strength = 8.;
     add_parameter("strength", strength, "strength of the PML");
   };
 
@@ -362,7 +363,9 @@ namespace Step81
   }
 
 
-
+  // @sect4{Maxwell Class}
+  // At this point we are ready to instantiate all the major functions of
+  // the finite element program and also a list of variables.
 
   template <int dim>
   class Maxwell : public ParameterAcceptor
@@ -403,6 +406,9 @@ namespace Step81
 
 
   // @sect4{The Constructor}
+  // The Constructor simply consists specifications for the mesh
+  // and the order of the fnite elements. These are editable through
+  // the .prm file.
 
   template <int dim>
   Maxwell<dim>::Maxwell()
@@ -415,7 +421,7 @@ namespace Step81
     scaling = 20;
     add_parameter("scaling", scaling, "scale of the hypercube geometry");
 
-    refinements = 8;
+    refinements = 10;
     add_parameter("refinements",
                   refinements,
                   "number of refinements of the geometry");
@@ -427,9 +433,6 @@ namespace Step81
     add_parameter("quadrature order",
                   quadrature_order,
                   "order of the quadrature");
-
-    n_outputs = 2;
-    add_parameter("number of outputs", n_outputs, "number of output images");
   }
 
 
@@ -733,7 +736,7 @@ namespace Step81
   {
     DataOut<2> data_out;
     data_out.attach_dof_handler(dof_handler);
-    data_out.add_data_vector(solution, "solution");
+      data_out.add_data_vector(solution, {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
     data_out.build_patches();
     std::ofstream output("solution.vtk");
     data_out.write_vtk(output);
@@ -752,7 +755,7 @@ namespace Step81
 
 } // namespace Step81
 
-// The following main function just calls the class step-81(), initializes the ParameterAcceptor,
+// The following main function calls the class step-81(), initializes the ParameterAcceptor,
 // and calls the run() function.
 
 int main()

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.