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
% ------------------------------------
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})
<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>
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
+++ /dev/null
-<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>
+++ /dev/null
-<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
#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
{
// 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
}
// @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
add_parameter("outer radius",
outer_radius,
"outer radius of the PML shell");
- strength = 0.;
+ strength = 8.;
add_parameter("strength", strength, "strength of the PML");
};
}
-
+ // @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
// @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()
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");
add_parameter("quadrature order",
quadrature_order,
"order of the quadrature");
-
- n_outputs = 2;
- add_parameter("number of outputs", n_outputs, "number of output images");
}
{
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);
} // 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()
projects / dealii.git / commitdiff