This commit introduces a time-harmonic Maxwell solver example step.
--- /dev/null
+##
+# 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()
--- /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 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 />
--- /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 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
--- /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
--- /dev/null
+Creating a mesh. Refining it. Writing it to a file.
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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;
+}
projects / dealii.git / commitdiff