]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step-81: Add a skeleton
authorManaswinee Bezbaruah <bezba004@tamu.edu>
Tue, 22 Jun 2021 18:26:33 +0000 (13:26 -0500)
committerMatthias Maier <tamiko@43-1.org>
Thu, 26 May 2022 05:16:01 +0000 (00:16 -0500)
This commit introduces a time-harmonic Maxwell solver example step.

examples/step-81/CMakeLists.txt [new file with mode: 0644]
examples/step-81/doc/builds-on [new file with mode: 0644]
examples/step-81/doc/intro.dox [new file with mode: 0644]
examples/step-81/doc/intro.dox-template [new file with mode: 0644]
examples/step-81/doc/kind [new file with mode: 0644]
examples/step-81/doc/results.dox [new file with mode: 0644]
examples/step-81/doc/results.dox-template [new file with mode: 0644]
examples/step-81/doc/tooltip [new file with mode: 0644]
examples/step-81/step-81.cc [new file with mode: 0644]

diff --git a/examples/step-81/CMakeLists.txt b/examples/step-81/CMakeLists.txt
new file mode 100644 (file)
index 0000000..c8e17bb
--- /dev/null
@@ -0,0 +1,55 @@
+##
+#  CMake script
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-81")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0)
+
+FIND_PACKAGE(deal.II 10.0.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+#
+# Are all dependencies fulfilled?
+#
+IF(NOT (DEAL_II_WITH_UMFPACK AND DEAL_II_WITH_COMPLEX_VALUES)) # keep in one line
+  MESSAGE(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_UMFPACK        = ON
+    DEAL_II_WITH_COMPLEX_VALUES = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
+    DEAL_II_WITH_UMFPACK        = ${DEAL_II_WITH_UMFPACK}
+    DEAL_II_WITH_COMPLEX_VALUES = ${DEAL_II_WITH_COMPLEX_VALUES}
+This conflicts with the requirements."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+SET(CLEAN_UP_FILES *.log *.vtu *.pvtu)
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-81/doc/builds-on b/examples/step-81/doc/builds-on
new file mode 100644 (file)
index 0000000..850b582
--- /dev/null
@@ -0,0 +1 @@
+step-8
diff --git a/examples/step-81/doc/intro.dox b/examples/step-81/doc/intro.dox
new file mode 100644 (file)
index 0000000..a4e114d
--- /dev/null
@@ -0,0 +1,337 @@
+<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,
+@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.
+\end{cases}
+@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,
+
+@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.
+\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.
+
+@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.
+
+<h3> Rescaling </h3>
+
+We will be using a rescaled version of the Maxwell's equations described above.
+The rescaling has the following key differences:<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.
+@f[
+\mathbf{J}_a = J_0 e_i\delta(x-a)
+@f]
+<br />
+
+Accordingly, our electric permittivity and magnetic permeability are rescaled by
+$\varepsilon_0$ and $\mu_0$ as follows:
+@f[
+\mu_r = \frac{1}{\mu_0}\mu,\qquad
+\varepsilon_r = \frac{1}{\varepsilon_0}\varepsilon.
+@f]
+
+We use the free space wave number $k_0 = \omega\sqrt{\varepsilon_0\mu_0}$, and
+the dipole strength, $J_0$, to arrive at the following rescaling of the vector
+fields and coordinates:
+@f[
+\begin{align*}
+\hat{x} = k_0x,\qquad
+\hat{\nabla} = \frac{1}{k_0}\nabla,\\
+\hat{\mathbf{H}} = \frac{k_0}{J_0}\mu^{-1}\mathbf{H},\qquad
+\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,\qquad
+\hat{\mathbf{M}}_a = \frac{k_0}{\omega\mu_0 J_0}\mathbf{M}_a.
+\end{align*}
+@f]
+
+Finally, the interface conductivity is rescaled as follows:
+@f[
+\sigma^{\Sigma}_r = \sqrt{\frac{\mu_0}{\varepsilon_0}}\sigma^{\Sigma}.
+@f]
+
+Accordingly, our rescaled equations are:
+@f[
+\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.
+\end{cases}
+@f]
+
+We will omit the hat in further discussion for ease of notation.
+
+<h3> Variational Statement</h3>
+
+Let $\Omega \subset \mathbb{R}^n$, $(n = 2,3)$ 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[
+\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)
+@f]
+
+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$.
+@f[
+\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
+- \int_\Sigma [\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} +
+\mu^{-1}\mathbf{M}_a)]_{\Sigma}\cdot \bar{\varphi}_T\;\text{d}o_x
+- \int_{\partial\Omega} (\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} +
+\mu^{-1}\mathbf{M}_a)) \cdot \bar{\varphi}_T\;\text{d}o_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.
+@f]
+
+We use the subscript $T$ to denote the tangential part of the given vector i.e.
+$F_T = (\nu\times F)\times\nu$ and $[\cdot]_{\Sigma}$ to denote a jump over
+$\Sigma$ i.e. $[F]_{\Sigma}(x) = \lim\limits_{s\searrow 0}(F(x+s\nu)-F(x-s\nu))$
+for $x\in \Sigma$.<br />
+
+For the computational domain $\Omega$, we introduce the absorbing boundary condition
+at $\partial\Omega$, which is obtained by using a first-order approximation of
+the Silver-Müller radiation condition, truncated at $\partial\Omega$.
+@f[
+\nu\times\mathbf{H}+\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}=0\qquad x\in\partial\Omega
+@f]
+We assume that $\mu_r^{-1}$ and $\varepsilon$ have well-defined square root. In
+our numerical computation, we combine the above absorbing boundary condition
+with a Perfectly Matched Layer (PML). <br />
+
+The jump condition can be expressed as a weak discontinuity as follows:
+@f[
+[\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)]_{\Sigma}
+= i\sigma_r^{\Sigma}\mathbf{E}_T,\qquad \text{on }\Sigma\\
+\nu \times (\mu_r^{-1}\nabla\times\mathbf{E} + \mu^{-1}\mathbf{M}_a)
+= -i\sqrt{\mu_r^{-1}\varepsilon_r}\mathbf{E}_T,\qquad \text{on }\partial\Omega.
+@f]
+
+Combining, our weak form is as follows:
+@f[
+\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{d}o_x
+- i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\mathbf{E}_T) \cdot
+(\nabla\times\bar{\varphi}_T)\;\text{d}o_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.
+@f]
+
+Assume that $\sigma_r^{\Sigma} \in L^{\infty}(\Sigma)^{2\times 2}$ is matrix-valued
+and symmetric, and has a semidefinite real and complex part. Let $\varepsilon_r$
+be a smooth scalar function with $–\text{Im}(\varepsilon_r) = 0$, or
+$\text{Im}(\varepsilon_r)\ge c > 0$ in $\Omega$. $\mu_r^{-1}$ is a smooth scalar
+such that $\sqrt{\mu_r^{-1}\varepsilon_r}$ is real valued and strictly positive
+in $\partial\Omega$. <br />
+
+$\mathbf{H}(curl;\Omega)$ is space of vector-valued, measurable and square
+integrable functions whose (distributive) curl admits a representation by a
+square integrable function. Define a Hilbert space
+@f[
+X(\Omega) = \{\varphi \in \mathbf{H}(curl;\Omega)\;\;:\;\; \varphi_T|_{\Sigma}
+\in L^2(\Sigma)^2,\;\varphi_T|_{\partial\Omega} \in L^2(\partial\Omega)^2\}
+@f]
+equipped with the norm $\|\varphi\|^2_X = \|\varphi\|^2_{L^2(\Omega)} +
+\|\nabla\times\varphi\|^2_{L^2(\Omega)} + \|\varphi_T\|^2_{L^2(\Sigma)} +
+\|\varphi_T\|^2_{L^2(\partial\Omega)}.$
+
+Define
+@f[
+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{d}o_x
+- i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\mathbf{E}_T) \cdot
+(\nabla\times\bar{\varphi}_T)\;\text{d}o_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.
+@f]
+
+Then, our rescaled weak formulation is:<br />
+Find a unique $\mathbf{E} \in X(\Omega)$ such that for all $\varphi \in X(\Omega)$
+@f[
+A(\mathbf{E},\varphi) = F(\varphi)
+@f]
+
+<h2> Discretization Scheme</h2>
+
+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}(\#\text{dofs})$. <br />
+
+Now, consider the finite element subspace $X_h(\Omega) \subset X(\Omega)$. Define
+the matrices
+@f[
+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{d}o_x
+          - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}\varphi_{i_T}
+          \cdot (\nabla\times \bar{\varphi}_{j_T})\;\text{d}o_x,
+@f]
+@f[
+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]
+
+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[
+A_{ij} = F_i
+@f]
+
+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 />
+
+<h2> Perfectly Matched Layer </h2>
+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 />
+
+<img src = "PML.png"></img>
+
+Introduce a change of coordinates
+@f[
+x \to \bar{x} =
+\begin{cases}
+x + ie_r\int\limits_\rho^r s(\tau)\text{d}\tau,\;\;\;\;\;\;\; r\ge\rho\\
+x\;\;\;\;\;\;\;\;\;\text{otherwise}
+\end{cases}
+@f]
+
+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[
+A = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2},
+\frac{1}{d\bar{d}}\right)T_{e_xe_r},\qquad
+B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\qquad
+C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right)
+T_{e_xe_r}.\qquad
+@f]
+
+where
+
+@f[
+d = 1 + is(r),\qquad
+\bar{d} = 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau.
+@f]
+
+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[
+\bar{\mu}_r^{-1} = \frac{\mu_r^{-1}}{d},\qquad
+\bar{\varepsilon}_r = A^{-1} \varepsilon_r B^{-1},\qquad
+\bar{\sigma}^{\Sigma}_r = C^{-1} \sigma^{\Sigma}_r B^{-1}.
+@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 />
diff --git a/examples/step-81/doc/intro.dox-template b/examples/step-81/doc/intro.dox-template
new file mode 100644 (file)
index 0000000..5e67a7f
--- /dev/null
@@ -0,0 +1,252 @@
+<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/kind b/examples/step-81/doc/kind
new file mode 100644 (file)
index 0000000..e62f4e7
--- /dev/null
@@ -0,0 +1 @@
+fluids
diff --git a/examples/step-81/doc/results.dox b/examples/step-81/doc/results.dox
new file mode 100644 (file)
index 0000000..a5b174b
--- /dev/null
@@ -0,0 +1,16 @@
+<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 observed, causing distorted images and necessitating
+a PML boundary condition. <br />
+TODO: add the images <br />
+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.<br />
+TODO: add the images<br />
+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.<br />
+TODO: add the images
diff --git a/examples/step-81/doc/results.dox-template b/examples/step-81/doc/results.dox-template
new file mode 100644 (file)
index 0000000..d195483
--- /dev/null
@@ -0,0 +1,15 @@
+<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
diff --git a/examples/step-81/doc/tooltip b/examples/step-81/doc/tooltip
new file mode 100644 (file)
index 0000000..ae43afa
--- /dev/null
@@ -0,0 +1 @@
+Creating a mesh. Refining it. Writing it to a file.
diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc
new file mode 100644 (file)
index 0000000..cf3debe
--- /dev/null
@@ -0,0 +1,787 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2021 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Manaswinee Bezbaruah, Matthias Maier, Texas A&M University, 2021.
+ */
+
+// @sect3{Include files}
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/parameter_acceptor.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_ilu.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+
+
+namespace Step81
+{
+  using namespace dealii;
+  using namespace std::complex_literals;
+
+
+
+  template <int dim>
+  class Parameters : public ParameterAcceptor
+  {
+  public:
+    Parameters();
+
+    using rank0_type = std::complex<double>;
+
+    using rank1_type = Tensor<1, dim, std::complex<double>>;
+
+    using rank2_type = Tensor<2, dim, rank0_type>;
+
+    using curl_type = Tensor<1, dim == 2 ? 1 : dim, rank0_type>;
+
+  public:
+    rank2_type epsilon(const Point<dim, double> &x,
+                       types::material_id        material);
+
+    std::complex<double> mu_inv(const Point<dim, double> &x,
+                                types::material_id        material);
+
+    rank2_type sigma(const dealii::Point<dim, double> &x,
+                     types::material_id                left,
+                     types::material_id                right);
+
+    rank1_type J_a(const dealii::Point<dim, double> &point,
+                   types::material_id                id);
+
+  private:
+    rank2_type           epsilon_1;
+    rank2_type           epsilon_2;
+    std::complex<double> mu_inv_1;
+    std::complex<double> mu_inv_2;
+    rank2_type           sigma_tensor;
+
+    double                 dipole_radius;
+    Point<dim>             dipole_position;
+    Tensor<1, dim, double> dipole_orientation;
+    rank0_type             dipole_strength;
+  };
+
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+
+
+  template <int dim>
+  Parameters<dim>::Parameters()
+    : ParameterAcceptor("Parameters")
+  {
+    epsilon_1[0][0] = 1.;
+    epsilon_1[1][1] = 1.;
+    add_parameter("material 1 epsilon",
+                  epsilon_1,
+                  "relative permittivity of material 1");
+
+    epsilon_2[0][0] = 1.;
+    epsilon_2[1][1] = 1.;
+    add_parameter("material 2 epsilon",
+                  epsilon_2,
+                  "relative permittivity of material 2");
+
+    mu_inv_1 = 1.;
+    add_parameter("material 1 mu_inv",
+                  mu_inv_1,
+                  "inverse of relative permeability of material 1");
+
+    mu_inv_2 = 1.;
+    add_parameter("material 2 mu_inv",
+                  mu_inv_2,
+                  "inverse of relative permeability of material 2");
+
+    sigma_tensor[0][0] = 0.001 + 0.2i;
+    sigma_tensor[1][1] = 0.001 + 0.2i;
+    add_parameter("sigma",
+                  sigma_tensor,
+                  "surface conductivity between material 1 and material 2");
+
+    dipole_radius = 0.3;
+    add_parameter("dipole radius", dipole_radius, "radius of the dipole");
+
+    dipole_position = Point<dim>(0., 0.8);
+    add_parameter("dipole position",
+                  dipole_position,
+                  "posititon of the dipole");
+
+    dipole_orientation = Tensor<1, dim, double>{{0., 1.}};
+    add_parameter("dipole orientation",
+                  dipole_orientation,
+                  "orientation of the dipole");
+
+    dipole_strength = 1.;
+    add_parameter("dipole strength", dipole_strength, "strength of the dipole");
+  }
+
+  template <int dim>
+  typename Parameters<dim>::rank2_type
+  Parameters<dim>::epsilon(const Point<dim, double> & /*x*/,
+                           types::material_id material)
+  {
+    return (material == 1 ? epsilon_1 : epsilon_2);
+  }
+
+  template <int dim>
+  std::complex<double> Parameters<dim>::mu_inv(const Point<dim, double> & /*x*/,
+                                               types::material_id material)
+  {
+    return (material == 1 ? mu_inv_1 : mu_inv_2);
+  }
+
+  template <int dim>
+  typename Parameters<dim>::rank2_type
+  Parameters<dim>::sigma(const dealii::Point<dim, double> & /*x*/,
+                         types::material_id left,
+                         types::material_id right)
+  {
+    return (left == right ? rank2_type() : sigma_tensor);
+  }
+
+  template <int dim>
+  typename Parameters<dim>::rank1_type
+  Parameters<dim>::J_a(const dealii::Point<dim, double> &point,
+                       types::material_id /*id*/)
+  {
+    rank1_type J_a;
+    const auto distance = (dipole_position - point).norm() / dipole_radius;
+    if (distance > 1.)
+      return J_a;
+    double scale = std::cos(distance * M_PI / 2.) *
+                   std::cos(distance * M_PI / 2.) / (M_PI / 2. - 2. / M_PI) /
+                   dipole_radius / dipole_radius;
+    J_a = dipole_strength * dipole_orientation * scale;
+    return J_a;
+  }
+
+
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+
+
+  template <int dim>
+  class PerfectlyMatchedLayer : public ParameterAcceptor
+  {
+  public:
+    static_assert(dim == 2, "dim == 2"); /* only works in 2D */
+
+    Parameters<dim> parameters;
+
+    using rank1_type = Tensor<1, dim, std::complex<double>>;
+
+    using rank2_type = Tensor<2, dim, std::complex<double>>;
+
+    PerfectlyMatchedLayer();
+
+    double inner_radius;
+    double outer_radius;
+    double strength;
+
+    std::complex<double> d_tensor(const Point<dim, double> point);
+
+    std::complex<double> d_bar_tensor(const Point<dim, double> point);
+
+
+    rank2_type T_exer(std::complex<double> d_1,
+                      std::complex<double> d_2,
+                      Point<dim>           point);
+
+    rank2_type a_matrix(const Point<dim, double> point);
+
+    rank2_type b_matrix(const Point<dim, double> point);
+
+    rank2_type c_matrix(const Point<dim, double> point);
+  };
+
+
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+
+
+  template <int dim>
+  PerfectlyMatchedLayer<dim>::PerfectlyMatchedLayer()
+    : ParameterAcceptor("PerfectlyMatchedLayer")
+  {
+    inner_radius = 12.;
+    add_parameter("inner radius",
+                  inner_radius,
+                  "inner radius of the PML shell");
+    outer_radius = 15.;
+    add_parameter("outer radius",
+                  outer_radius,
+                  "outer radius of the PML shell");
+    strength = 0.;
+    add_parameter("strength", strength, "strength of the PML");
+  };
+
+
+  template <int dim>
+  typename std::complex<double>
+  PerfectlyMatchedLayer<dim>::d_tensor(const Point<dim, double> point)
+  {
+    const auto   radius = point.norm();
+    const double s =
+      strength * ((radius - inner_radius) * (radius - inner_radius)) /
+      ((outer_radius - inner_radius) * (outer_radius - inner_radius));
+    return 1 + 1.0i * s;
+  }
+
+
+  template <int dim>
+  typename std::complex<double>
+  PerfectlyMatchedLayer<dim>::d_bar_tensor(const Point<dim, double> point)
+  {
+    const auto   radius = point.norm();
+    const double s_bar =
+      strength / 3. *
+      ((radius - inner_radius) * (radius - inner_radius) *
+       (radius - inner_radius)) /
+      (radius * (outer_radius - inner_radius) * (outer_radius - inner_radius));
+    return 1 + 1.0i * s_bar;
+  }
+
+
+  template <int dim>
+  typename PerfectlyMatchedLayer<dim>::rank2_type
+  PerfectlyMatchedLayer<dim>::T_exer(std::complex<double> d_1,
+                                     std::complex<double> d_2,
+                                     Point<dim>           point)
+  {
+    rank2_type result;
+    result[0][0] = point[0] * point[0] * d_1 + point[1] * point[1] * d_2;
+    result[0][1] = point[0] * point[1] * (d_1 - d_2);
+    result[1][0] = point[0] * point[1] * (d_1 - d_2);
+    result[1][1] = point[1] * point[1] * d_1 + point[0] * point[0] * d_2;
+    return result;
+  }
+
+
+  template <int dim>
+  typename PerfectlyMatchedLayer<dim>::rank2_type
+  PerfectlyMatchedLayer<dim>::a_matrix(const Point<dim, double> point)
+  {
+    const auto d     = d_tensor(point);
+    const auto d_bar = d_bar_tensor(point);
+    return invert(T_exer(d * d, d * d_bar, point)) *
+           T_exer(d * d, d * d_bar, point);
+  }
+
+
+  template <int dim>
+  typename PerfectlyMatchedLayer<dim>::rank2_type
+  PerfectlyMatchedLayer<dim>::b_matrix(const Point<dim, double> point)
+  {
+    const auto d     = d_tensor(point);
+    const auto d_bar = d_bar_tensor(point);
+    return invert(T_exer(d, d_bar, point)) * T_exer(d, d_bar, point);
+  }
+
+
+  template <int dim>
+  typename PerfectlyMatchedLayer<dim>::rank2_type
+  PerfectlyMatchedLayer<dim>::c_matrix(const Point<dim, double> point)
+  {
+    const auto d     = d_tensor(point);
+    const auto d_bar = d_bar_tensor(point);
+    return invert(T_exer(1. / d_bar, 1. / d, point)) *
+           T_exer(1. / d_bar, 1. / d, point);
+  }
+
+
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+
+
+  template <int dim>
+  class Maxwell : public ParameterAcceptor
+  {
+  public:
+    Maxwell();
+    void run();
+
+  private:
+    /* run time parameters */
+    double       scaling;
+    unsigned int refinements;
+    unsigned int fe_order;
+    unsigned int quadrature_order;
+    bool         postprocess;
+    unsigned int n_outputs;
+
+    void parse_parameters_callback();
+    void make_grid();
+    void setup_system();
+    void assemble_system();
+    void solve();
+    void output_results(bool postprocess);
+
+    Parameters<dim>            parameters;
+    PerfectlyMatchedLayer<dim> perfectly_matched_layer;
+
+    Triangulation<dim> triangulation;
+    DoFHandler<dim>    dof_handler;
+
+    std::unique_ptr<FiniteElement<dim>> fe;
+
+    AffineConstraints<double> constraints;
+    SparsityPattern           sparsity_pattern;
+    SparseMatrix<double>      system_matrix;
+    Vector<double>            solution;
+    Vector<double>            system_rhs;
+  };
+
+
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////
+
+
+  template <int dim>
+  Maxwell<dim>::Maxwell()
+    : ParameterAcceptor("Maxwell")
+    , dof_handler(triangulation)
+  {
+    ParameterAcceptor::parse_parameters_call_back.connect(
+      std::bind(&Maxwell<dim>::parse_parameters_callback, this));
+
+    scaling = 20;
+    add_parameter("scaling", scaling, "scale of the hypercube geometry");
+
+    refinements = 8;
+    add_parameter("refinements",
+                  refinements,
+                  "number of refinements of the geometry");
+
+    fe_order = 0;
+    add_parameter("fe order", fe_order, "order of the finite element space");
+
+    quadrature_order = 1;
+    add_parameter("quadrature order",
+                  quadrature_order,
+                  "order of the quadrature");
+
+    postprocess = true;
+    add_parameter("postprocess", postprocess, "solution postprocessing option");
+
+    n_outputs = 2;
+    add_parameter("number of outputs", n_outputs, "number of output images");
+  }
+
+
+  template <int dim>
+  void Maxwell<dim>::parse_parameters_callback()
+  {
+    fe = std::make_unique<FESystem<dim>>(FE_Nedelec<dim>(fe_order), 2);
+  }
+
+
+  template <int dim>
+  void Maxwell<dim>::make_grid()
+  {
+    GridGenerator::hyper_cube(triangulation, -scaling, scaling);
+    triangulation.refine_global(refinements);
+
+    for (auto &cell : triangulation.active_cell_iterators())
+      if (cell->center()[1] > 0.)
+        cell->set_material_id(1);
+      else
+        cell->set_material_id(2);
+
+    std::cout << "Number of active cells: " << triangulation.n_active_cells()
+              << std::endl;
+  }
+
+
+  template <int dim>
+  void Maxwell<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(*fe);
+    std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl;
+
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+    constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    constraints.close();
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    dsp,
+                                    constraints,
+                                    /* keep_constrained_dofs = */ true);
+    sparsity_pattern.copy_from(dsp);
+    system_matrix.reinit(sparsity_pattern);
+  }
+
+
+  template <int dim>
+  DEAL_II_ALWAYS_INLINE inline Tensor<1, dim, std::complex<double>>
+  tangential_part(const dealii::Tensor<1, dim, std::complex<double>> &tensor,
+                  const Tensor<1, dim> &                              normal)
+  {
+    auto result = tensor;
+    result[0]   = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
+    result[1]   = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
+    return result;
+  }
+
+
+  template <int dim>
+  void Maxwell<dim>::assemble_system()
+  {
+    QGauss<dim>     quadrature_formula(quadrature_order);
+    QGauss<dim - 1> face_quadrature_formula(quadrature_order);
+
+    FEValues<dim, dim>     fe_values(*fe,
+                                 quadrature_formula,
+                                 update_values | update_gradients |
+                                   update_quadrature_points |
+                                   update_JxW_values);
+    FEFaceValues<dim, dim> fe_face_values(*fe,
+                                          face_quadrature_formula,
+                                          update_values | update_gradients |
+                                            update_quadrature_points |
+                                            update_normal_vectors |
+                                            update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe->dofs_per_cell;
+
+    const unsigned int n_q_points      = quadrature_formula.size();
+    const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        fe_values.reinit(cell);
+        cell_matrix = 0.;
+        cell_rhs    = 0.;
+
+        cell->get_dof_indices(local_dof_indices);
+
+        FEValuesViews::Vector<dim> real_part(fe_values, 0);
+        FEValuesViews::Vector<dim> imag_part(fe_values, dim);
+        const auto &quadrature_points = fe_values.get_quadrature_points();
+
+        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          {
+            const Point<dim, double> &position = quadrature_points[q_point];
+            const auto                radius   = position.norm();
+            const auto                id       = cell->material_id();
+            const auto inner_radius = perfectly_matched_layer.inner_radius;
+
+            auto       mu_inv  = parameters.mu_inv(position, id);
+            auto       epsilon = parameters.epsilon(position, id);
+            const auto J_a     = parameters.J_a(position, id);
+
+            if (radius >= inner_radius)
+              {
+                auto A = perfectly_matched_layer.a_matrix(position);
+                auto B = perfectly_matched_layer.b_matrix(position);
+                auto d = perfectly_matched_layer.d_tensor(position);
+
+                mu_inv  = mu_inv / d;
+                epsilon = invert(A) * epsilon * invert(B);
+              };
+
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                const auto phi_i = real_part.value(i, q_point) -
+                                   1.0i * imag_part.value(i, q_point);
+                const auto curl_phi_i = real_part.curl(i, q_point) -
+                                        1.0i * imag_part.curl(i, q_point);
+
+                const auto rhs_value =
+                  (1.0i * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
+                cell_rhs(i) += rhs_value.real();
+
+                for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                  {
+                    const auto phi_j = real_part.value(j, q_point) +
+                                       1.0i * imag_part.value(j, q_point);
+                    const auto curl_phi_j = real_part.curl(j, q_point) +
+                                            1.0i * imag_part.curl(j, q_point);
+
+                    const auto temp =
+                      (scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
+                       scalar_product(epsilon * phi_j, phi_i)) *
+                      fe_values.JxW(q_point);
+                    cell_matrix(i, j) += temp.real();
+                  }
+              }
+          }
+
+        for (const auto &face : cell->face_iterators())
+          {
+            if (face->at_boundary())
+              {
+                fe_face_values.reinit(cell, face);
+                FEValuesViews::Vector<dim> real_part(fe_face_values, 0);
+                FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);
+
+                for (unsigned int q_point = 0; q_point < n_face_q_points;
+                     ++q_point)
+                  {
+                    const Point<dim, double> position =
+                      quadrature_points[q_point];
+                    const auto radius = position.norm();
+                    const auto id     = cell->material_id();
+                    const auto inner_radius =
+                      perfectly_matched_layer.inner_radius;
+
+
+                    auto mu_inv  = parameters.mu_inv(position, id);
+                    auto epsilon = parameters.epsilon(position, id);
+
+                    if (radius >= inner_radius)
+                      {
+                        auto A = perfectly_matched_layer.a_matrix(position);
+                        auto B = perfectly_matched_layer.b_matrix(position);
+                        auto d = perfectly_matched_layer.d_tensor(position);
+
+                        mu_inv  = mu_inv / d;
+                        epsilon = invert(A) * epsilon * invert(B);
+                      };
+
+                    const auto normal = fe_face_values.normal_vector(q_point);
+
+                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                      {
+                        const auto phi_i = real_part.value(i, q_point) -
+                                           1.0i * imag_part.value(i, q_point);
+                        const auto phi_i_T = tangential_part(phi_i, normal);
+
+                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                          {
+                            const auto phi_j =
+                              real_part.value(j, q_point) +
+                              1.0i * imag_part.value(j, q_point);
+                            const auto phi_j_T =
+                              tangential_part(phi_j, normal) *
+                              fe_face_values.JxW(q_point);
+
+                            const auto prod      = mu_inv * epsilon;
+                            const auto sqrt_prod = prod;
+
+                            const auto temp =
+                              -1.0i *
+                              scalar_product((sqrt_prod * phi_j_T), phi_i_T);
+                            cell_matrix(i, j) += temp.real();
+                          } /* j */
+                      }     /* i */
+                  }         /* q_point */
+              }
+            else
+              {
+                Assert(!face->at_boundary(), ExcMessage("oops!"));
+                const auto face_index = cell->face_iterator_to_index(face);
+
+                const auto id1 = cell->material_id();
+                const auto id2 = cell->neighbor(face_index)->material_id();
+
+                if (id1 == id2)
+                  continue; /* skip this face */
+
+                fe_face_values.reinit(cell, face);
+                FEValuesViews::Vector<dim> real_part(fe_face_values, 0);
+                FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);
+
+                for (unsigned int q_point = 0; q_point < n_face_q_points;
+                     ++q_point)
+                  {
+                    const Point<dim, double> position =
+                      quadrature_points[q_point];
+                    const auto radius = position.norm();
+                    const auto inner_radius =
+                      perfectly_matched_layer.inner_radius;
+
+                    auto sigma = parameters.sigma(position, id1, id2);
+
+                    if (radius >= inner_radius)
+                      {
+                        auto B = perfectly_matched_layer.b_matrix(position);
+                        auto C = perfectly_matched_layer.c_matrix(position);
+                        sigma  = invert(C) * sigma * invert(B);
+                        ;
+                      };
+
+                    const auto normal = fe_face_values.normal_vector(q_point);
+
+                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                      {
+                        const auto phi_i = real_part.value(i, q_point) -
+                                           1.0i * imag_part.value(i, q_point);
+                        const auto phi_i_T = tangential_part(phi_i, normal);
+
+                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                          {
+                            const auto phi_j =
+                              real_part.value(j, q_point) +
+                              1.0i * imag_part.value(j, q_point);
+                            const auto phi_j_T = tangential_part(phi_j, normal);
+
+                            const auto temp =
+                              -1.0i *
+                              scalar_product((sigma * phi_j_T), phi_i_T) *
+                              fe_face_values.JxW(q_point);
+                            cell_matrix(i, j) += temp.real();
+                          } /* j */
+                      }     /* i */
+                  }         /* q_point */
+              }
+          }
+
+        constraints.distribute_local_to_global(
+          cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+      }
+  };
+
+
+  template <int dim>
+  void Maxwell<dim>::solve()
+  {
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+    A_direct.vmult(solution, system_rhs);
+  }
+
+
+  template <int dim>
+  void Maxwell<dim>::output_results(bool postprocess)
+  {
+    DataOut<2> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "solution");
+    data_out.build_patches();
+    std::ofstream output("solution.vtk");
+    data_out.write_vtk(output);
+
+
+    if (postprocess)
+      {
+        for (unsigned int alpha = 1; alpha <= n_outputs; alpha++)
+          {
+            std::cout << "Running step:" << alpha << std::endl;
+            Vector<double> postprocessed;
+            postprocessed.reinit(solution);
+            for (unsigned int i = 0; i < dof_handler.n_dofs(); i += 2)
+              {
+                postprocessed[i] =
+                  std::cos(-2 / n_outputs * M_PI * alpha) * solution[i] +
+                  std::sin(-2 / n_outputs * M_PI * alpha) * solution[i + 1];
+              }
+            data_out.add_data_vector(postprocessed, "postprocessed");
+            data_out.build_patches();
+            const std::string filename =
+              "postprocessed-" + Utilities::int_to_string(alpha) + ".vtk";
+            std::ofstream output(filename);
+            data_out.write_vtk(output);
+            std::cout << "Done running step:" << alpha << std::endl;
+          }
+      }
+  }
+
+
+  template <int dim>
+  void Maxwell<dim>::run()
+  {
+    make_grid();
+    setup_system();
+    assemble_system();
+    solve();
+    output_results(postprocess);
+  }
+
+} // namespace Step81
+
+
+int main()
+{
+  try
+    {
+      Step81::Maxwell<2> maxwell_2d;
+      dealii::ParameterAcceptor::initialize("parameters.prm");
+      maxwell_2d.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}

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.