Add program for goal-oriented Maxwell eigenvalue hp-refinement

diff --git a/Maxwell-Eigenvalue-hp-Refinement/CMakeLists.txt b/Maxwell-Eigenvalue-hp-Refinement/CMakeLists.txt
new file mode 100644 (file)
index 0000000..546cf14
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/CMakeLists.txt
@@ -0,0 +1,42 @@
+##
+#  CMake script for the maxwell-hp program:
+##
+
+SET(TARGET "maxwell-hp")
+
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.4 QUIET
+  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()
+
+IF(NOT DEAL_II_WITH_PETSC OR NOT DEAL_II_WITH_SLEPC OR DEAL_II_PETSC_WITH_COMPLEX)
+  MESSAGE(FATAL_ERROR "
+Error! This program requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_PETSC = ON
+    DEAL_II_WITH_SLEPC = ON
+    DEAL_II_PETSC_WITH_COMPLEX = OFF
+    DEAL_II_WITH_GSL = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options
+    DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
+    DEAL_II_WITH_SLEPC = ${DEAL_II_WITH_SLEPC}
+    DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX}
+    DEAL_II_WITH_GSL = ${DEAL_II_WITH_GSL}
+which conflict with the requirements."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()

diff --git a/Maxwell-Eigenvalue-hp-Refinement/Readme.md b/Maxwell-Eigenvalue-hp-Refinement/Readme.md
new file mode 100644 (file)
index 0000000..b738f07
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/Readme.md
@@ -0,0 +1,156 @@
+Readme file for Maxwell-Eigenvallue-hp-Refinement
+============================
+
+@note The implementation of this program is in part based on [1].
+
+Motivation for project
+----------------------
+
+From the source free Maxwell equations in differential form, we may find the following eigenvalue problem involving the electric field $\mathbf{E}$,
+$$ \nabla\times(\mu_r^{-1}\nabla\times\mathbf{E})-k_0^2\varepsilon_r\mathbf{E} = 0 \textrm{ in } \Omega, $$
+ where $\mu_r$ and $\varepsilon_r$ denote, respectively, the relative permeability and permitivity of the medium (which we assume to be homogeneous), and $k_0$ signifies the free space wavenumber, for some $\Omega \subset \mathbb{R}^d, \, d = 2,3.$ Finding (approximate) solutions of this eigenvalue problem poses a number of challenges computationally; for those interested, we refer to the excellent thesis of S. Zaglmayr [2].
+
+ In the remainder of this project, we assume $d=2$, though the methodology is largely unaffected by this choice. We further assume perfect electrical conductor (PEC) boundary conditions: $\hat{\textbf{n}}\times\textbf{E}=0 \textrm{ on }\partial\Omega$, $\hat{\textbf{n}}$ being the outward normal vector.
+
+ In the standard way, we consider *weak* solutions by solving the variational form of the eigenvalue problem, which, in the 2-D case, is found to be the following after Galerkin testing:
+ $$ \textrm{Find } U_{hp}=\left\{ \mathbf{u}_{hp},\,\lambda_{hp}\right\}\in V_{hp}\times \mathbb{R}_{>0} \textrm{ such that} $$
+ $$ a(\textbf{u}_{hp},\,\boldsymbol{\phi}_{hp}) = \lambda_{hp} m(\textbf{u}_{hp},\,\boldsymbol{\phi}_{hp}) \quad \forall\boldsymbol{\phi}\in V_{hp}, $$
+ with $a(\textbf{u}_{hp},\,\boldsymbol{\phi}_{hp}) = \langle \nabla_t\times\textbf{u}_{hp},\,\nabla_t\times\boldsymbol{\phi}_{hp}\rangle$ (note: $\nabla_t$ represents the transversal gradient operator and $\langle \cdot ,\, \cdot \rangle$ represents the $L^2$ inner-product), and $m(\textbf{u}_{hp},\,\boldsymbol{\phi}_{hp}) = \langle \textbf{u}_{hp},\,\boldsymbol{\phi}_{hp} \rangle,$ for a finite dimensional subspace $V_{hp}$ to be further specified below along with its infinite dimensional analog $V$ associated with an exact solution $U=\left\{ \textbf{u},\, \lambda \right\}$.
+
+
+
+For this problem, we consider a single family of quantities of interest (QoIs), namely the approximation error of the approximate eigenvalue $\lambda_{hp}$, i.e.,
+$$ e_{\lambda_{hp}} := \lambda-\lambda_{hp}. $$
+
+Some comments on the discretization of Maxwell's equations
+----------------------------------------------------------
+In the proper solution of variational problem, $V_{hp}$ is not arbitrary, but is instead a subspace of $H(\mathrm{curl};\,\Omega)$ or, with the boundary conditions indicated above, $H_{0}(\mathrm{curl};\,\Omega)$, where
+$$ H(\mathrm{curl};\,\Omega) = \left\{\textbf{u}\in \left[ L_2(\Omega)\right]^d \, \mathrm{s.t.} \, \nabla\times\textbf{u}\in \left[ L_2(\Omega)\right]^{2d-3} \right\}$$
+and 
+$$
+H_{0}(\mathrm{curl};\,\Omega) = \left\{ \textbf{u}\in H(\mathrm{curl};\,\Omega) \, \mathrm{s.t.} \, \hat{\textbf{n}}\times\textbf{u} = 0 \right\}.
+$$
+
+The convergence of the discrete problem (e.g., as $h\rightarrow 0$) to the continuous one may be proved via a discrete compactness property [3]. It also possible for other choices of finite element spaces to *converge*, just not necessarily to the *correct* solution, which was the case for discretizations (which typically treated the components of $\textbf{u}$ as belonging to $H^1$) of the Maxwell PDE prior to the work of Nédélec.
+
+It should be noted that while the choice of the appropriate finite element space eliminates most spurious solutions to the generalized eigenvalue problem above, not all are absent. Specifically, those which do not satisfy the source free divergence condition on the electric field $\textbf{E}$ $$\nabla\cdot\textbf{E} = 0$$ still cluster around $\lambda = 0$. For $\lambda > 0$, corresponding to physical eigenpairs, the divergence condition is automatically satisfied.
+
+A variety of techniques may be leveraged to enforce the divergence condition explicitly, e.g., mixed finite element formulations. However (as noted above), we instead restate the variational problem, restricting our eigenvalues to $\lambda\in\mathbb{R}_{>0}$, which implies some additional mechanics in how we solve the eigenvalue problem. Solving for interior eigenvalues, as is required for this approach, is less than ideal; however, formalizing the divergence condition in the discretization is left as an extension of this project by the reader.
+
+On the adaptivity
+-----------------
+While adaptive refinement based on pure error indicators may perform (surprisingly) well, in many cases we are interested in stricter interpretations of the mesh adaption problem that rely on error estimates and error contribution estimates. When the relationship between the accuracy of the goal-functional and the error indicator weakens (if it ever existed), refinement instructions may yield useless changes, even, as  we will see, to the extent of converging to the wrong value.
+
+Galerkin projection, being a global process, is *not* interpolation. The question, then, is how to relate the accuracy of our QoI, which may formally be a global quantity (e.g., the radar cross section of a target from a certain look angle) or a local one (e.g., value of the solution at a point), on the approximation quality throughout the discretization. Solving a secondary *global* problem, the adjoint (or dual) problem, provides exactly the mechanism we need via the adjoint solution, a generalized Green's function. Specifically, we study a goal-oriented error estimate based on the Dual Weighted Residual (DWR), similar to that  studied in Step-14.
+
+<strong>For this section of project, we strongly recommended the following readings: [4] (W. Bangerth and R. Rannacher; 2003), Chapter 3, 7; [5] (V. Heuveline and R. Rannacher; 2001).</strong>
+
+Several assumptions as outlined in [5], some of which may likely not be satisfied, should be stated first. Firstly, the algorithm assumes that the eigenvalues are simple. Secondly, for approximating $\lambda$, $\lambda_h$ should be closer to $\lambda$ than any other eigenvalue of the approximate problem. This outlines a necessarily vague condition on the coarseness/fineness of the starting discretization, provided the natural assumption that $\lambda_{h}\rightarrow\lambda$ as $h\rightarrow 0$ holds.
+
+With this in mind, the error may be expressed in DWR-form:
+$$ e_{\lambda_{hp}}(1-\sigma_{hp}) = a(\textbf{u}_{hp},\, \textbf{u} - \boldsymbol{\psi}_{hp}) - \lambda_{hp}m(\textbf{u}_{hp},\, \textbf{u}-\boldsymbol{\psi}_{hp})$$
+for arbitrary $\boldsymbol{\psi}_{hp}\in V_{hp}$, with $\sigma_{hp} = \frac12m(\textbf{u}-\textbf{u}_{hp},\,\textbf{u}-\textbf{u}_{hp})$.
+
+It is important to note that the expression above assumes normalization such that
+$$ \langle \textbf{u},\,\textbf{u}\rangle = 1, $$
+and likewise for $\textbf{u}_{hp}$, though naturally, given that our solutions are eigenfunctions, this normalization is not a unique choice.
+
+That $\boldsymbol{\psi}_{hp}$ is arbitrary *does not* imply that any choice is equally good (i.e., informative) for generating refinement indicators. Certainly the most obvious and easiest choice of $\boldsymbol{\psi_{hp}} = 0$ yields the same global error, yet in the accumulation of error *contributions* retains irrelevant "contributions" that hamper adaptivity.
+
+As a result we extract from $\textbf{u}$ the portions that belong to $V_{hp}$ by some interpolation or projection. According to the hierarchical nature of the `FENedelec` cell we can simply retain the coefficients associated with higher-order shape functions and discard the rest.
+
+In any event, the global error estimate (in DWR-form) is accumulated in a cellwise fashion. Based on the continuity requirements of the Nédélec cell (continuous in the tangential direction), not every DoF is constrained to a single cell, motivating some "sharing" of contributions. One way is through integration-by-parts, producing a cell and boundary residual. The sharing is conducted by averaging the boundary term for one cell with its neighbor. In total, the error is computed over every cell $K$ as
+$$\begin{align} e_{\lambda_{hp}}(1-\sigma_{hp}) = &\sum_K {\langle \nabla\times\nabla\times\textbf{u}_{hp},\,\textbf{u}-\boldsymbol\psi_{hp}\rangle}_K \\ &- \frac12 \left[ {\langle \hat{\textbf{n}}\times(\nabla\times\textbf{u}_{hp}),\,\textbf{u}-\boldsymbol\psi_{hp}\rangle}_{\partial K}\right. \\ &\,\,\,- \left. {\langle \hat{\textbf{n}}\times(\nabla\times\textbf{u}_{hp}),\,\textbf{u}-\boldsymbol\psi_{hp}\rangle}_{\partial K^{'}} \right] \\ &- \lambda_{hp} m_K(\textbf{u}_{hp}, \textbf{u}-\boldsymbol\psi_{hp}),\end{align}$$
+where $K^{'}$ denotes the neighbor cells to $K$.
+
+Given the tangential continuity of $\textbf{u}-\boldsymbol{\psi}_{hp}$, its evaluation may be performed once for each edge/face.
+
+The above expressions assume access to the exact solution $\textbf{u}$, which is not satisfiable in general (and if it were, why would we need refinement?) and therefore a substitution is necessary. Unfortunately, Galerkin orthogonality precludes taking $\textbf{u}_{hp}$. Instead, we replace $\textbf{u}$ by $\textbf{u}_{hp^+}\in V_{hp^+}$, $V_{hp^+}\supset V_{hp},$ where $\textbf{u}_{hp^+}$ belongs to an enriched finite element space generated by increasing the local expansions orders by 1 throughout the mesh.
+
+For $h$- or $p$-refinement, we now have sufficient data. For $hp$-refinement, we either need further post-processing on the error contributions or an additional mechanism to distinguish between the profitability of $h$ versus $p$. Since the theoretical conditions for exponential convergence depend on the solution regularity, *smoothness indication* suggests a promising vehicle for the $hp$-decision.
+
+A prevalent choice in the literature is the examination of the decay rates of Legendre expansions. Fourier expansions are similarly employed, though the non-polynomial integrands are slightly less convenient for numerical integration. The theory and methodology has been developed over the years, e.g., [6]-[8], culminating in the excellent work of M. Fehling [9], which is now included in the deal.II library. As such, we leave out further description except for the following: as the solution is vectorial, we apply smoothness indication to each component. As we consider only *isotropic* refinements, the smallest decay is propagated forward for the $hp$-decision, though directional refinements would benefit from retaining all decay rates.
+
+As a final note for those that may extend this project to other problems or applications, the $hp$-decision herein is really only feasible for affine mappings between the reference cell and its images as the Legendre integrals are otherwise too burdensome. 
+
+Numerical Results
+-----------------
+We finally arrive at the heart of the project.
+
+As we are paying a price in computation and implementation by studying a goal-oriented approach, a comparison with error indication is warranted. The reference method is based on the Kelly error indicator (applied to each component of the solution), which we also refer to as the "Jump" approach. 
+
+Except for the refinement indicators, the approaches are identical: the smoothness estimation is the same, the number of elements refined each iteration (20%) is the same, and the starting discretizations are the same. Since the DWR strategy utilizes the higher order solution, we use its approximation (and the higher number of DoFs required) for the comparison.
+
+The starting discretization is shown below, consisting of three unit squares arranged in an 'L'-shape. This benchmark was originally proposed by M. Dauge in [10], which includes reference values for the first five eigenvalues.
+
+![Problem geometry](./doc/model.PNG)
+
+The first nine eigenfunctions for this problem group equally into three classes: singular (unbounded electric field at the reentrant corner), sharp (nonsmooth but bounded at the reentrant corner), and globally smooth. The globally smooth eigenpairs may be resolved relatively easily and sufficiently accurately to determine that the eigenvalues are multiples of $\pi^2$. As such, we focus exclusively on the first two classes of eigenfunctions.
+
+The adaption pipeline is summarized as follows:
+1. Solve the forward problem (DWR/Kelly) and the adjoint problem (DWR)
+2. Generate error contribution estimates (DWR) or error indicators (Kelly)
+3. Mark the top 20% elements with the largest refinement indicators by magnitude
+4. For those cells marked for refinement, estimate the local smoothness of the forward solution (Legendre decay rate)
+5. According to the Legendre decay rate, reclassify the refinement to $p$ if necessary
+6. Apply refinement instructions and return to 1.
+
+Of course, this process is easily applied to alternative estimators, indicators, and selection criteria.
+
+Approximate values of the three "sharp" eigenvalues are listed below:
+- $\lambda_2 = 3.53403136678$
+- $\lambda_5  = 11.3894793979$
+- $\lambda_9 = 23.3443719571$
+
+Approximate values of the three "singular" eigenvalues are listed below:
+- $\lambda_1 = 1.47562182397$
+- $\lambda_6 = 12.5723873200$
+- $\lambda_8 = 21.4247335393$
+
+We first examine application of the two approaches to the "sharp" eigenpairs. In the first figure below, we see that the DWR approach significantly outperforms  that of the jump-based Kelly approach. In the second figure, we illustrate the close agreement between the estimated relative error and the actual relative error, which applies exclusively to the DWR-based method. 
+
+The $\log$-cuberoot scaling of the axes is based on the theoretical convergence rates for the eigenvalues as described in [11], such that a linear trend indicates exponential convergence.
+
+![Problem geometry](./doc/convergence_sharp.png)
+
+![Problem geometry](./doc/estimated_error_sharp.png)
+
+We now examine the same procedure for the  "singular" eigenpairs. Here we see that while the DWR-based approach is robust and achieves consistent exponential convergence with respect to the number of degrees of freedom, the jump-based method fails to  converge to the benchmark values. Even before saturation of accuracy by the jump-based method, the DWR-method is several orders of magnitude more accurate for the same cost.
+
+![Problem geometry](./doc/convergence_singular.png)
+
+![Problem geometry](./doc/estimated_error_singular.png)
+
+Once again, as shown directly above, the error estimate for the approximate eigenvalue provides close agreement with the actual error.
+
+We encourage the reader to study these examples further, including the auxiliary tasks of selecting and tuning an appropriate generalized eigenvalue problem solver, etc, which are beyond the scope of this discussion. Furthermore, the selection of which cells to refine (namely, the top 20%) is a fair approach for this comparison, but is evidently suboptimal especially with access to the extremely high-quality DWR-based refinement indicators, and so we invite the reader to investigate alternative marking mechanisms in the context of $hp$-refinement.
+
+To run the code
+---------------
+
+After running `cmake` and compiling via `make`, you can run the
+executable by either executing `make run` or using `./maxwell-hp`
+on the command line. Executing either `make debug` or `make release` 
+swaps between "debug mode" and "release mode," respectively.
+
+The parameters file "maxwell-hp.prm" may be modified to refine for different eigenvalues (by specifying the "set target eigenvalue" option).
+
+The end results of the program are `.vtu` files for each iteration of the refinement procedure. All files are named according to the strategy employed (Kelly or DWR). The program also generates a text file with the eigenvalues and their cost (in terms of the number  of DoFs to attain said eigenvalues). Specifically, the text file has the following structure:
+- The first column houses the eigenvalue from the lower-order discretization
+- The second column houses the number of DoFs for that discretization
+- In the case of the DWR-based method, the third and fourth columns match the first two, except they are for the higher-order discretization
+
+Finally, for the DWR-based method, "error_estimates.txt" contains the estimated eigenvalue error (the signed sum of every error contribution estimate) at each cycle of the refinement procedure.
+
+## References 
+* [1] J. J. Harmon and B. M. Notaroš, "Adaptive hp-Refinement for 2-D Maxwell Eigenvalue Problems: Method and Benchmarks," *IEEE Transactions on Antennas and Propagation*, vol. 70, no. 6, pp. 4663-4673, June 2022.
+* [2] S. Zaglmayr, “High order finite element methods for electromagnetic field computation,” Ph.D. dissertation, Institute for Numerical Mathematics, Johannes Kepler University Linz, Linz, Austria, 2006.
+* [3] D. Boffi, M. Costabel, M. Dauge, and L. F. Demkowicz, “Discrete compactness for the hp version of rectangular edge finite elements,” *SIAM Journal on Numerical Analysis*, vol. 44, no. 3, pp. 979–1004, 2006.
+* [4] W. Bangerth and R. Rannacher, *Adaptive Finite Element Methods for Differential Equations*. Birkhauser Basel, 2003.
+* [5] V. Heuveline and R. Rannacher, “A posteriori error control for finite element approximations of elliptic eigenvalue problems,” *J. Adv. Comp. Math*, vol. 15, pp. 107–138, 2001.
+* [6] C. Mavriplis, “Adaptive mesh strategies for the spectral element method,” *Computer Methods in Applied Mechanics and Engineering*, vol. 116, no. 1, pp. 77–86, 1994.
+* [7] P. Houston, B. Senior, and E. Süli, “Sobolev regularity estimation for hp-adaptive finite element methods,” *Numerical Mathematics and Advanced Applications*, pp. 619–644, 2003.
+* [8] T. Eibner and J. M. Melenk, “An adaptive strategy for hp-FEM based on testing for analyticity,” *Comput Mesh*, vol. 39, no. 39, pp. 575–595, 2007.
+* [9] M. Fehling, “Algorithms for massively parallel generic hp-adaptive finite element methods,” Ph.D. dissertation, Univ. Wuppertal, 2020.
+* [10] M. Dauge, “Benchmark computations for Maxwell equations for the approximation of highly singular solutions,” https://perso.univ-rennes1.fr/monique.dauge/benchmax.html, 2004.
+* [11] J. Coyle and P. D. Ledger, “Evidence of exponential convergence in the computation of Maxwell eigenvalues,” *Computer Methods in Applied Mechanics and Engineering*, vol. 194, pp. 587–604, 02 2005.

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/author b/Maxwell-Eigenvalue-hp-Refinement/doc/author
new file mode 100644 (file)
index 0000000..cd09b39
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/doc/author
@@ -0,0 +1 @@
+Jake J. Harmon <jake.harmon@ieee.org>

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/builds-on b/Maxwell-Eigenvalue-hp-Refinement/doc/builds-on
new file mode 100644 (file)
index 0000000..38b2f59
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/doc/builds-on
@@ -0,0 +1 @@
+step-14 step-27 step-36

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/convergence_sharp.png b/Maxwell-Eigenvalue-hp-Refinement/doc/convergence_sharp.png
new file mode 100644 (file)
index 0000000..6bf55e3
Binary files /dev/null and b/Maxwell-Eigenvalue-hp-Refinement/doc/convergence_sharp.png differ

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/convergence_singular.png b/Maxwell-Eigenvalue-hp-Refinement/doc/convergence_singular.png
new file mode 100644 (file)
index 0000000..9a9b6b5
Binary files /dev/null and b/Maxwell-Eigenvalue-hp-Refinement/doc/convergence_singular.png differ

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/dependencies b/Maxwell-Eigenvalue-hp-Refinement/doc/dependencies
new file mode 100644 (file)
index 0000000..78c3d98
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/doc/dependencies
@@ -0,0 +1 @@
+DEAL_II_WITH_PETSC DEAL_II_WITH_SLEPC

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/entry-name b/Maxwell-Eigenvalue-hp-Refinement/doc/entry-name
new file mode 100644 (file)
index 0000000..598e3cd
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/doc/entry-name
@@ -0,0 +1 @@
+Goal-Oriented hp-Adaptivity for the Maxwell Eigenvalue Problem

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/estimated_error_sharp.png b/Maxwell-Eigenvalue-hp-Refinement/doc/estimated_error_sharp.png
new file mode 100644 (file)
index 0000000..625641b
Binary files /dev/null and b/Maxwell-Eigenvalue-hp-Refinement/doc/estimated_error_sharp.png differ

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/estimated_error_singular.png b/Maxwell-Eigenvalue-hp-Refinement/doc/estimated_error_singular.png
new file mode 100644 (file)
index 0000000..8519c19
Binary files /dev/null and b/Maxwell-Eigenvalue-hp-Refinement/doc/estimated_error_singular.png differ

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/model.PNG b/Maxwell-Eigenvalue-hp-Refinement/doc/model.PNG
new file mode 100644 (file)
index 0000000..0ca5ddd
Binary files /dev/null and b/Maxwell-Eigenvalue-hp-Refinement/doc/model.PNG differ

diff --git a/Maxwell-Eigenvalue-hp-Refinement/doc/tooltip b/Maxwell-Eigenvalue-hp-Refinement/doc/tooltip
new file mode 100644 (file)
index 0000000..80e9e52
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/doc/tooltip
@@ -0,0 +1 @@
+An implementation of goal-oriented hp-refinement for the Maxwell eigenvalue problem, with a focus on computing eigenvalues accurately

diff --git a/Maxwell-Eigenvalue-hp-Refinement/maxwell-hp.cc b/Maxwell-Eigenvalue-hp-Refinement/maxwell-hp.cc
new file mode 100644 (file)
index 0000000..338f39b
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/maxwell-hp.cc
@@ -0,0 +1,2271 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2022 by the deal.II authors and Jake J. Harmon
+ *
+ * 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: Jake J. Harmon, 2022
+ *
+ */
+
+
+
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/slepc_solver.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+// For parallelization (using WorkStream and Intel TBB)
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/work_stream.h>
+
+#include "petscpc.h"
+
+// For Error Estimation/Indication and Smoothness Indication
+#include <deal.II/fe/fe_tools.h>
+
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/smoothness_estimator.h>
+// For refinement
+#include <deal.II/grid/grid_refinement.h>
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+namespace Operations
+{
+  /**
+  Computes the curl-curl term needed to fill the stiffness matrix (specific to
+  2-D)
+  */
+  double
+  curlcurl(const dealii::FEValues<2> &fe_values,
+           const unsigned int &       i,
+           const unsigned int &       j,
+           const unsigned int &       q_point)
+  {
+    auto gradu1_x1x2 = fe_values.shape_grad_component(i, q_point, 0);
+    auto gradu2_x1x2 = fe_values.shape_grad_component(i, q_point, 1);
+
+    auto gradv1_x1x2 = fe_values.shape_grad_component(j, q_point, 0);
+    auto gradv2_x1x2 = fe_values.shape_grad_component(j, q_point, 1);
+    return (gradu2_x1x2[0] - gradu1_x1x2[1]) *
+           (gradv2_x1x2[0] - gradv1_x1x2[1]);
+  }
+
+  /**
+  Computes the dot product of the shape functions needed to fill the mass matrix
+  */
+  template <int dim>
+  inline double
+  dot_term(const dealii::FEValues<dim> &fe_values,
+           const unsigned int &         i,
+           const unsigned int &         j,
+           const unsigned int &         q_point)
+  {
+    double output = 0.0;
+    for (unsigned int comp = 0; comp < dim; ++comp)
+      {
+        output += fe_values.shape_value_component(i, q_point, comp) *
+                  fe_values.shape_value_component(j, q_point, comp);
+      }
+    return output;
+  }
+} // namespace Operations
+/**
+The Structures namespace includes the necessary functions for constructing two
+examples problems, the so-called "Standard waveguide" (width/height = 2), and
+the L-domain waveguide
+*/
+namespace Structures
+{
+  using namespace dealii;
+
+  void
+  create_L_waveguide(Triangulation<2> &triangulation, const double &scaling)
+  {
+    const unsigned int dim = 2;
+
+    const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
+                                            {scaling * 0.5, scaling * 0.0},
+                                            {scaling * 0.0, scaling * 0.5},
+                                            {scaling * 0.5, scaling * 0.5},
+                                            {scaling * 0.0, scaling * 1.0},
+                                            {scaling * 0.5, scaling * 1.0},
+                                            {scaling * 1.0, scaling * 0.5},
+                                            {scaling * 1.0, scaling * 1.0}};
+
+    const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
+      cell_vertices = {{{0, 1, 2, 3}}, {{2, 3, 4, 5}}, {{3, 6, 5, 7}}};
+    const unsigned int         n_cells = cell_vertices.size();
+    std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
+    for (unsigned int i = 0; i < n_cells; ++i)
+      {
+        for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      }
+    triangulation.create_triangulation(vertices, cells, SubCellData());
+    triangulation.refine_global(1);
+  }
+
+
+  void
+  create_standard_waveguide(Triangulation<2> &triangulation,
+                            const double &    scaling)
+  {
+    const unsigned int dim = 2;
+
+    const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
+                                            {scaling * 0.6, scaling * 0.0},
+                                            {scaling * 0.0, scaling * 0.3},
+                                            {scaling * 0.6, scaling * 0.3}};
+
+    const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
+                               cell_vertices = {{{0, 1, 2, 3}}};
+    const unsigned int         n_cells       = cell_vertices.size();
+    std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
+    for (unsigned int i = 0; i < n_cells; ++i)
+      {
+        for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      }
+    triangulation.create_triangulation(vertices, cells, SubCellData());
+    triangulation.refine_global(0);
+  }
+} // namespace Structures
+/**
+The Maxwell namespace includes all of the classes for solving the Maxwell
+eigenvalue problem
+*/
+namespace Maxwell
+{
+  using namespace dealii;
+
+  /*
+  The "Base" class provides the universal functionality of any eigensolver,
+  namely the parameters for the problem, an underlying triangulation, and
+  functionality for setting the refinement cycle and to output the solution.
+
+  In this case, and for any future class, the use of raw pointers (as opposed to
+  "smart" pointers) indicates a lack of ownership. Specifically, the
+  triangulation raw pointer is pointing to a triangulation that is owned (and
+  created) elsewhere.
+  */
+  template <int dim>
+  class Base
+  {
+  public:
+    Base(const std::string &prm_file, Triangulation<dim> &coarse_grid);
+
+    virtual unsigned int
+    solve_problem() = 0; // Implemented by a derived class
+    virtual void
+    set_refinement_cycle(const unsigned int cycle);
+
+    virtual void
+    output_solution() = 0; // Implemented by a derived class
+
+
+  protected:
+    Triangulation<dim> *              triangulation;
+    unsigned int                      refinement_cycle = 0;
+    std::unique_ptr<ParameterHandler> parameters;
+    unsigned int                      n_eigenpairs = 1;
+    double                            target       = 0.0;
+    unsigned int                      eigenpair_selection_scheme;
+    unsigned int                      max_cycles       = 0;
+    ompi_communicator_t *             mpi_communicator = PETSC_COMM_SELF;
+  };
+
+  /**
+  Reads in the parameters file and the triangulation
+  */
+  template <int dim>
+  Base<dim>::Base(const std::string &prm_file, Triangulation<dim> &coarse_grid)
+    : triangulation(&coarse_grid)
+    , parameters(std::make_unique<ParameterHandler>())
+  {
+    parameters->declare_entry(
+      "Eigenpair selection scheme",
+      "1",
+      Patterns::Integer(0, 1),
+      "The type of eigenpairs to find (0 - smallest, 1 - target)");
+    parameters->declare_entry("Number of eigenvalues/eigenfunctions",
+                              "1",
+                              Patterns::Integer(0, 100),
+                              "The number of eigenvalues/eigenfunctions "
+                              "to be computed.");
+    parameters->declare_entry("Target eigenvalue",
+                              "1",
+                              Patterns::Anything(),
+                              "The target eigenvalue (if scheme == 1)");
+
+    parameters->declare_entry("Cycles number",
+                              "1",
+                              Patterns::Integer(0, 1500),
+                              "The number of cycles in refinement");
+    parameters->parse_input(prm_file);
+
+    eigenpair_selection_scheme =
+      parameters->get_integer("Eigenpair selection scheme");
+
+    // The project currently only supports selection by a target eigenvalue.
+    // Furthermore, only one eigenpair can be computed at a time.
+    assert(eigenpair_selection_scheme == 1 &&
+           "Selection by a target is the only currently supported option!");
+    n_eigenpairs =
+      parameters->get_integer("Number of eigenvalues/eigenfunctions");
+    assert(
+      n_eigenpairs == 1 &&
+      "Only the computation of a single eigenpair is currently supported!");
+
+    target     = parameters->get_double("Target eigenvalue");
+    max_cycles = parameters->get_integer("Cycles number");
+    if (eigenpair_selection_scheme == 1)
+      n_eigenpairs = 1;
+  }
+
+  template <int dim>
+  void
+  Base<dim>::set_refinement_cycle(const unsigned int cycle)
+  {
+    refinement_cycle = cycle;
+  }
+
+  /**
+  Provides the central solver (derived from the base class). Virtual inheritance
+  is crucial to eliminate compiler ambiguity in the case of the
+  DualWeightedResidual.
+  */
+  template <int dim>
+  class EigenSolver : public virtual Base<dim>
+  {
+  public:
+    EigenSolver(const std::string & prm_file,
+                Triangulation<dim> &coarse_grid,
+                const unsigned int &minimum_degree,
+                const unsigned int &maximum_degree,
+                const unsigned int &starting_degree);
+
+    virtual unsigned int
+    solve_problem() override;
+
+    virtual unsigned int
+    n_dofs() const;
+
+    template <class SolverType>
+    void
+    initialize_eigensolver(SolverType &eigensolver);
+
+    virtual void
+    setup_system();
+
+    virtual void
+    assemble_system();
+
+  protected:
+    const std::unique_ptr<hp::FECollection<dim>> fe_collection;
+    std::unique_ptr<hp::QCollection<dim>>        quadrature_collection;
+    std::unique_ptr<hp::QCollection<dim - 1>>    face_quadrature_collection;
+    DoFHandler<dim>                              dof_handler;
+    const unsigned int                           max_degree, min_degree;
+    // for the actual solution
+    std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> eigenfunctions;
+    std::unique_ptr<std::vector<double>>                     eigenvalues;
+    Vector<double>                                           solution;
+
+    double *
+    get_lambda_h();
+
+    Vector<double> *
+    get_solution();
+
+    void
+    convert_solution();
+
+  private:
+    AffineConstraints<double>   constraints;
+    PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
+  };
+
+  /**
+  Typical constructor. Executes the constructor for the base class and creates
+  the unique pointers for the fe_collection, quadrature_collection, etc.
+  */
+  template <int dim>
+  EigenSolver<dim>::EigenSolver(const std::string & prm_file,
+                                Triangulation<dim> &triangulation,
+                                const unsigned int &minimum_degree,
+                                const unsigned int &maximum_degree,
+                                const unsigned int &starting_degree)
+    : Base<dim>(prm_file, triangulation)
+    , fe_collection(std::make_unique<hp::FECollection<dim>>())
+    , quadrature_collection(std::make_unique<hp::QCollection<dim>>())
+    , face_quadrature_collection(std::make_unique<hp::QCollection<dim - 1>>())
+    , dof_handler(triangulation)
+    , max_degree(maximum_degree)
+    , min_degree(minimum_degree)
+    , eigenfunctions(
+        std::make_unique<std::vector<PETScWrappers::MPI::Vector>>())
+    , eigenvalues(std::make_unique<std::vector<double>>())
+  {
+    for (unsigned int degree = min_degree; degree <= max_degree; ++degree)
+      {
+        fe_collection->push_back(FE_Nedelec<dim>(degree - 1));
+        // Generate quadrature collection with sorted quadrature weights
+        const QGauss<dim>  quadrature(degree + 1);
+        const QSorted<dim> sorted_quadrature(quadrature);
+        quadrature_collection->push_back(sorted_quadrature);
+
+        const QGauss<dim - 1>  face_quadrature(degree + 1);
+        const QSorted<dim - 1> sorted_face_quadrature(face_quadrature);
+        face_quadrature_collection->push_back(sorted_face_quadrature);
+      }
+    // adjust the discretization
+    if (starting_degree > min_degree && starting_degree <= max_degree)
+      {
+        const unsigned int start_diff = starting_degree - min_degree;
+        typename DoFHandler<dim>::active_cell_iterator
+          cell1 = dof_handler.begin_active(),
+          endc1 = dof_handler.end();
+        for (; cell1 < endc1; ++cell1)
+          {
+            cell1->set_active_fe_index(start_diff);
+          }
+      }
+  }
+
+  /**
+  Returns the (first) eigenvalue.
+  TODO: Generalize to arbitrary, valid eigenvalues
+  */
+  template <int dim>
+  double *
+  EigenSolver<dim>::get_lambda_h()
+  {
+    return &(*eigenvalues)[0];
+  }
+
+  /**
+  Returns the (first) eigenvector.
+  TODO: Generalize to arbitrary
+  */
+  template <int dim>
+  Vector<double> *
+  EigenSolver<dim>::get_solution()
+  {
+    return &solution;
+  }
+
+  /**
+  Temporary helper function for copying the solution vector
+  */
+  template <int dim>
+  void
+  EigenSolver<dim>::convert_solution()
+  {
+    solution.reinit((*eigenfunctions)[0].size());
+    for (unsigned int i = 0; i < solution.size(); ++i)
+      solution[i] = (*eigenfunctions)[0][i];
+  }
+
+  /**
+  Initializes the eigensolver according the selection scheme in the parameters
+  file. Additionally, applies the necessary problem type (GHEP) and introduces
+  the Shift-and-Invert spectrum transformation based on the specified target
+  value
+  */
+  template <int dim>
+  template <class SolverType>
+  void
+  EigenSolver<dim>::initialize_eigensolver(SolverType &eigensolver)
+  {
+    // From the parameters class, initialize the eigensolver...
+    switch (this->eigenpair_selection_scheme)
+      {
+        case 1:
+          eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
+          // eigensolver.set_target_eigenvalue(this->target);
+          break;
+        default:
+          eigensolver.set_which_eigenpairs(EPS_SMALLEST_MAGNITUDE);
+
+          break;
+      }
+    eigensolver.set_problem_type(EPS_GHEP);
+    // apply a Shift-Invert spectrum transformation
+
+    double shift_scalar = this->parameters->get_double("Target eigenvalue");
+    //         //For the shift-and-invert transformation
+    SLEPcWrappers::TransformationShiftInvert::AdditionalData additional_data(
+      shift_scalar);
+    SLEPcWrappers::TransformationShiftInvert spectral_transformation(
+      this->mpi_communicator, additional_data);
+
+    eigensolver.set_transformation(spectral_transformation);
+    eigensolver.set_target_eigenvalue(this->target);
+  }
+
+  /**
+  Solves the eigenvalue problem and applies the constraints to the
+  eigenfunctions
+  */
+  template <int dim>
+  unsigned int
+  EigenSolver<dim>::solve_problem()
+  {
+    setup_system();
+    assemble_system();
+
+    SolverControl                    solver_control(dof_handler.n_dofs() * 10,
+                                 5.0e-8,
+                                 false,
+                                 false);
+    SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control,
+                                                 this->mpi_communicator);
+
+    initialize_eigensolver(eigensolver);
+
+    // solve the problem
+    eigensolver.solve(stiffness_matrix,
+                      mass_matrix,
+                      *eigenvalues,
+                      *eigenfunctions,
+                      eigenfunctions->size());
+    for (auto &entry : *eigenfunctions)
+      {
+        constraints.distribute(entry);
+      }
+    convert_solution();
+
+    return solver_control.last_step();
+  }
+
+  template <int dim>
+  unsigned int
+  EigenSolver<dim>::n_dofs() const
+  {
+    return dof_handler.n_dofs();
+  }
+
+  /**
+  Distributes the degrees of freedom and makes the necessary hanging_node
+  constraints, which includes the constraints for non-uniform $p$, and for the
+  Dirichlet boundary.
+  */
+  template <int dim>
+  void
+  EigenSolver<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(*fe_collection);
+    constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    DoFTools::make_zero_boundary_constraints(dof_handler, constraints);
+    constraints.close();
+
+    eigenfunctions->resize(this->n_eigenpairs);
+    eigenvalues->resize(this->n_eigenpairs);
+
+    IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs();
+
+    for (auto &entry : *eigenfunctions)
+      {
+        entry.reinit(eigenfunction_index_set, MPI_COMM_WORLD);
+      }
+  }
+
+  /**
+  Fills the mass and stiffness matrices
+  */
+  template <int dim>
+  void
+  EigenSolver<dim>::assemble_system()
+  {
+    hp::FEValues<dim> hp_fe_values(*fe_collection,
+                                   *quadrature_collection,
+                                   update_values | update_gradients |
+                                     update_quadrature_points |
+                                     update_JxW_values);
+    // Prep the system matrices for the solution
+    stiffness_matrix.reinit(dof_handler.n_dofs(),
+                            dof_handler.n_dofs(),
+                            dof_handler.max_couplings_between_dofs());
+    mass_matrix.reinit(dof_handler.n_dofs(),
+                       dof_handler.n_dofs(),
+                       dof_handler.max_couplings_between_dofs());
+
+    FullMatrix<double> cell_stiffness_matrix, cell_mass_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+        cell_stiffness_matrix.reinit(dofs_per_cell, dofs_per_cell);
+        cell_stiffness_matrix = 0;
+
+        cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
+        cell_mass_matrix = 0;
+
+        hp_fe_values.reinit(cell);
+
+        const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+
+        for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+             ++q_point)
+          {
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                  {
+                    // Note that (in general) the Nedelec element is not
+                    // primitive, namely that the shape functions are vectorial
+                    // with components in more than one direction
+
+                    cell_stiffness_matrix(i, j) +=
+                      Operations::curlcurl(fe_values, i, j, q_point) *
+                      fe_values.JxW(q_point);
+
+                    cell_mass_matrix(i, j) +=
+                      (Operations::dot_term(fe_values, i, j, q_point)) *
+                      fe_values.JxW(q_point);
+                  }
+              }
+            local_dof_indices.resize(dofs_per_cell);
+            cell->get_dof_indices(local_dof_indices);
+          }
+
+        constraints.distribute_local_to_global(cell_stiffness_matrix,
+                                               local_dof_indices,
+                                               stiffness_matrix);
+        constraints.distribute_local_to_global(cell_mass_matrix,
+                                               local_dof_indices,
+                                               mass_matrix);
+      }
+    stiffness_matrix.compress(VectorOperation::add);
+    mass_matrix.compress(VectorOperation::add);
+
+    for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
+      if (constraints.is_constrained(i))
+        {
+          stiffness_matrix.set(i, i, 10000.0);
+          mass_matrix.set(i, i, 1);
+        }
+    // since we have just set individual elements, we need the following
+    stiffness_matrix.compress(VectorOperation::insert);
+    mass_matrix.compress(VectorOperation::insert);
+  }
+
+  /**
+  The main PrimalSolver, which is derived from the EigenSolver class. Provides a
+  limited amount of additional functionality
+  */
+  template <int dim>
+  class PrimalSolver : public EigenSolver<dim>
+  {
+  public:
+    PrimalSolver(const std::string & prm_file,
+                 Triangulation<dim> &triangulation,
+                 const unsigned int &min_degree,
+                 const unsigned int &max_degree,
+                 const unsigned int &starting_degree);
+
+    virtual void
+    output_solution()
+      override; // Implements the output solution of the base class...
+    virtual unsigned int
+    n_dofs() const override;
+  };
+
+  template <int dim>
+  PrimalSolver<dim>::PrimalSolver(const std::string & prm_file,
+                                  Triangulation<dim> &triangulation,
+                                  const unsigned int &min_degree,
+                                  const unsigned int &max_degree,
+                                  const unsigned int &starting_degree)
+    : Base<dim>(prm_file, triangulation)
+    , EigenSolver<dim>(prm_file,
+                       triangulation,
+                       min_degree,
+                       max_degree,
+                       starting_degree)
+  {}
+
+  /**
+  Outputs the first eigenpair (based on the target eigenvalue)
+  TODO: Generalize to multiple eigenpairs
+  */
+  template <int dim>
+  void
+  PrimalSolver<dim>::output_solution()
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(this->dof_handler);
+    Vector<double> fe_degrees(this->triangulation->n_active_cells());
+    for (const auto &cell : this->dof_handler.active_cell_iterators())
+      fe_degrees(cell->active_cell_index()) =
+        (*this->fe_collection)[cell->active_fe_index()].degree;
+    data_out.add_data_vector(fe_degrees, "fe_degree");
+    data_out.add_data_vector((*this->eigenfunctions)[0],
+                             std::string("eigenfunction_no_") +
+                               Utilities::int_to_string(0));
+
+    std::cout << "Eigenvalue: " << (*this->eigenvalues)[0]
+              << " NDoFs: " << this->dof_handler.n_dofs() << std::endl;
+    std::ofstream eigenvalues_out(
+      "eigenvalues-" + std::to_string(this->refinement_cycle) + ".txt");
+
+    eigenvalues_out << std::setprecision(20) << (*this->eigenvalues)[0] << " "
+                    << this->dof_handler.n_dofs() << std::endl;
+
+    eigenvalues_out.close();
+
+
+    data_out.build_patches();
+    std::ofstream output("eigenvectors-" +
+                         std::to_string(this->refinement_cycle) + ".vtu");
+    data_out.write_vtu(output);
+  }
+
+  template <int dim>
+  unsigned int
+  PrimalSolver<dim>::n_dofs() const
+  {
+    return EigenSolver<dim>::n_dofs();
+  }
+
+  // Note, that at least for the demonstrated problem (i.e., a Hermitian problem
+  // and eigenvalue QoI), the dual problem is identical to the primal problem;
+  // however, it is convenient to separate them in this manner (e.g., for
+  // considering functionals of the eigenfunction).
+  template <int dim>
+  class DualSolver : public EigenSolver<dim>
+  {
+  public:
+    DualSolver(const std::string & prm_file,
+               Triangulation<dim> &triangulation,
+               const unsigned int &min_degree,
+               const unsigned int &max_degree,
+               const unsigned int &starting_degree);
+  };
+
+  template <int dim>
+  DualSolver<dim>::DualSolver(const std::string & prm_file,
+                              Triangulation<dim> &triangulation,
+                              const unsigned int &min_degree,
+                              const unsigned int &max_degree,
+                              const unsigned int &starting_degree)
+    : Base<dim>(prm_file, triangulation)
+    , EigenSolver<dim>(prm_file,
+                       triangulation,
+                       min_degree,
+                       max_degree,
+                       starting_degree)
+  {}
+
+} // namespace Maxwell
+/**
+The second major namespace, which includes all the classes for error
+estimation and error indication.
+*/
+namespace ErrorIndicators
+{
+  using namespace Maxwell;
+
+  /**
+  The DualWeightedResidual is derived from the PrimalSolver and DualSolver. In
+  this case, the DualSolver is taken with a finite element space with shape
+  functions of one polynomial degree higher.
+  */
+  template <int dim, bool report_dual>
+  class DualWeightedResidual : public PrimalSolver<dim>, public DualSolver<dim>
+  {
+  public:
+    void
+    output_eigenvalue_data(std::ofstream &os);
+    void
+    output_qoi_error_estimates(std::ofstream &os);
+
+    std::string
+    name() const
+    {
+      return "DWR";
+    }
+    DualWeightedResidual(const std::string & prm_file,
+                         Triangulation<dim> &triangulation,
+                         const unsigned int &min_primal_degree,
+                         const unsigned int &max_primal_degree,
+                         const unsigned int &starting_primal_degree);
+
+    virtual unsigned int
+    solve_problem() override;
+
+    virtual void
+    output_solution() override;
+
+    virtual unsigned int
+    n_dofs() const override;
+
+    void
+    estimate_error(Vector<double> &error_indicators);
+
+    DoFHandler<dim> *
+    get_DoFHandler();
+
+    DoFHandler<dim> *
+    get_primal_DoFHandler();
+
+    DoFHandler<dim> *
+    get_dual_DoFHandler();
+
+    hp::FECollection<dim> *
+    get_FECollection();
+
+    hp::FECollection<dim> *
+    get_primal_FECollection();
+
+    std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+    get_eigenfunctions();
+
+    std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+    get_primal_eigenfunctions();
+
+    std::unique_ptr<std::vector<double>> &
+    get_primal_eigenvalues();
+
+    std::unique_ptr<std::vector<double>> &
+    get_dual_eigenvalues();
+
+    void
+    synchronize_discretization();
+
+    unsigned int
+    get_max_degree()
+    {
+      return PrimalSolver<dim>::fe_collection->max_degree();
+    }
+    double qoi_error_estimate = 0;
+
+  private:
+    void
+    embed(const DoFHandler<dim> &          dof1,
+          const DoFHandler<dim> &          dof2,
+          const AffineConstraints<double> &constraints,
+          const Vector<double> &           solution,
+          Vector<double> &                 u2);
+
+    void
+    extract(const DoFHandler<dim> &          dof1,
+            const DoFHandler<dim> &          dof2,
+            const AffineConstraints<double> &constraints,
+            const Vector<double> &           solution,
+            Vector<double> &                 u2);
+
+
+
+    /*The following FEValues objects are unique_ptrs to 1) avoid default
+    constructors for these objects, and 2) automate memory management*/
+    std::unique_ptr<hp::FEValues<dim>>        cell_hp_fe_values;
+    std::unique_ptr<hp::FEFaceValues<dim>>    face_hp_fe_values;
+    std::unique_ptr<hp::FEFaceValues<dim>>    face_hp_fe_values_neighbor;
+    std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values;
+
+    std::unique_ptr<hp::FEValues<dim>>     cell_hp_fe_values_forward;
+    std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_forward;
+    std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor_forward;
+    std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values_forward;
+    using FaceIntegrals =
+      typename std::map<typename DoFHandler<dim>::face_iterator, double>;
+
+    unsigned int
+    solve_primal_problem();
+
+    unsigned int
+    solve_dual_problem();
+
+    void
+    normalize_solutions(Vector<double> &primal_solution,
+                        Vector<double> &dual_weights);
+
+    double
+    get_global_QoI_error(Vector<double> &dual_solution,
+                         Vector<double> &error_indicators);
+
+    void
+    initialize_error_estimation_data();
+
+    void
+    estimate_on_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      const Vector<double> &                                primal_solution,
+      const Vector<double> &                                dual_weights,
+      const double &                                        lambda_h,
+      Vector<double> &                                      error_indicators,
+      FaceIntegrals &                                       face_integrals);
+
+    void
+    integrate_over_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      const Vector<double> &                                primal_solution,
+      const Vector<double> &                                dual_weights,
+      const double &                                        lambda_h,
+      Vector<double> &                                      error_indicators);
+
+    void
+    integrate_over_regular_face(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      const unsigned int &                                  face_no,
+      const Vector<double> &                                primal_solution,
+      const Vector<double> &                                dual_weights,
+      FaceIntegrals &                                       face_integrals);
+
+    void
+    integrate_over_irregular_face(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      const unsigned int &                                  face_no,
+      const Vector<double> &                                primal_solution,
+      const Vector<double> &                                dual_weights,
+      FaceIntegrals &                                       face_integrals);
+  };
+
+  /**
+  Basic constructor, also initializes the unique pointers for evaluating the
+  cell and edge residuals in the error estimate.
+  */
+  template <int dim, bool report_dual>
+  DualWeightedResidual<dim, report_dual>::DualWeightedResidual(
+    const std::string & prm_file,
+    Triangulation<dim> &triangulation,
+    const unsigned int &min_primal_degree,
+    const unsigned int &max_primal_degree,
+    const unsigned int &starting_primal_degree)
+    : Base<dim>(prm_file, triangulation)
+    , PrimalSolver<dim>(prm_file,
+                        triangulation,
+                        min_primal_degree,
+                        max_primal_degree,
+                        starting_primal_degree)
+    , DualSolver<dim>(prm_file,
+                      triangulation,
+                      min_primal_degree + 1,
+                      max_primal_degree + 1,
+                      starting_primal_degree + 1)
+  {
+    initialize_error_estimation_data();
+  }
+
+  /**
+  If we are "reporting" the dual solution (e.g., for the purposes of smoothness
+  estimation), we must decide which dof_handler to provide.
+  */
+  template <int dim, bool report_dual>
+  DoFHandler<dim> *
+  DualWeightedResidual<dim, report_dual>::get_DoFHandler()
+  {
+    if (!report_dual)
+      return &(PrimalSolver<dim>::dof_handler);
+    else
+      return &(DualSolver<dim>::dof_handler);
+  }
+
+  // See above function, but to specifically output the primal DoFHandler...
+  template <int dim, bool report_dual>
+  DoFHandler<dim> *
+  DualWeightedResidual<dim, report_dual>::get_primal_DoFHandler()
+  {
+    return &(PrimalSolver<dim>::dof_handler);
+  }
+
+  // See above function, but for the FECollection
+  template <int dim, bool report_dual>
+  hp::FECollection<dim> *
+  DualWeightedResidual<dim, report_dual>::get_FECollection()
+  {
+    if (!report_dual)
+      return &*(PrimalSolver<dim>::fe_collection);
+    else
+      return &*(DualSolver<dim>::fe_collection);
+  }
+
+  // See above function, but for the primal FECollection
+  template <int dim, bool report_dual>
+  hp::FECollection<dim> *
+  DualWeightedResidual<dim, report_dual>::get_primal_FECollection()
+  {
+    return &*(PrimalSolver<dim>::fe_collection);
+  }
+
+  template <int dim, bool report_dual>
+  DoFHandler<dim> *
+  DualWeightedResidual<dim, report_dual>::get_dual_DoFHandler()
+  {
+    return &(DualSolver<dim>::dof_handler);
+  }
+
+  //
+  template <int dim, bool report_dual>
+  std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+  DualWeightedResidual<dim, report_dual>::get_eigenfunctions()
+  {
+    if (!report_dual)
+      return (PrimalSolver<dim>::eigenfunctions);
+    else
+      return (DualSolver<dim>::eigenfunctions);
+  }
+
+  //
+  template <int dim, bool report_dual>
+  std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+  DualWeightedResidual<dim, report_dual>::get_primal_eigenfunctions()
+  {
+    return (PrimalSolver<dim>::eigenfunctions);
+  }
+
+  //
+  template <int dim, bool report_dual>
+  std::unique_ptr<std::vector<double>> &
+  DualWeightedResidual<dim, report_dual>::get_primal_eigenvalues()
+  {
+    return PrimalSolver<dim>::eigenvalues;
+  }
+
+  //
+  template <int dim, bool report_dual>
+  std::unique_ptr<std::vector<double>> &
+  DualWeightedResidual<dim, report_dual>::get_dual_eigenvalues()
+  {
+    return DualSolver<dim>::eigenvalues;
+  }
+
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::output_solution()
+  {
+    PrimalSolver<dim>::output_solution();
+  }
+
+  // Solves the primal problem
+  template <int dim, bool report_dual>
+  unsigned int
+  DualWeightedResidual<dim, report_dual>::solve_primal_problem()
+  {
+    return PrimalSolver<dim>::solve_problem();
+  }
+
+  // Solves the dual problem
+  template <int dim, bool report_dual>
+  unsigned int
+  DualWeightedResidual<dim, report_dual>::solve_dual_problem()
+  {
+    return DualSolver<dim>::solve_problem();
+  }
+
+  /**
+  Provides the publicly accessible solve_problem function,
+  which solves both the primal and dual problems
+  */
+  template <int dim, bool report_dual>
+  unsigned int
+  DualWeightedResidual<dim, report_dual>::solve_problem()
+  {
+    DualWeightedResidual<dim, report_dual>::solve_primal_problem();
+    return DualWeightedResidual<dim, report_dual>::solve_dual_problem();
+  }
+
+  /**
+  Returns the number of dofs that the primal solver requires
+  */
+  template <int dim, bool report_dual>
+  unsigned int
+  DualWeightedResidual<dim, report_dual>::n_dofs() const
+  {
+    return PrimalSolver<dim>::n_dofs();
+  }
+
+  /**
+  This function synchronizes the expansion orders. When working with two finite
+  element spaces, we must apply the modifications made to one dof_handler
+  (according to the boolean report_dual) we must also make the same
+  modifications to the other finite element space.
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::synchronize_discretization()
+  {
+    /*Note: No additional checks need to be made ensuring that these operations
+       are legal as these checks are made prior to entering this function (i.e.,
+       if the primal attains a degree N,
+        then, by construction, a degree of N+1 must be permissible for the
+       dual)*/
+    DoFHandler<dim> *dof1 = &(PrimalSolver<dim>::dof_handler);
+    DoFHandler<dim> *dof2 = &(DualSolver<dim>::dof_handler);
+
+    if (report_dual)
+      {
+        // In this case, we have modified the polynomial orders for the dual;
+        // need to update the primal
+        dof1 = &(DualSolver<dim>::dof_handler);
+        dof2 = &(PrimalSolver<dim>::dof_handler);
+      }
+    typename DoFHandler<dim>::active_cell_iterator cell1 = dof1->begin_active(),
+                                                   endc1 = dof1->end();
+    typename DoFHandler<dim>::active_cell_iterator cell2 = dof2->begin_active();
+    for (; cell1 < endc1; ++cell1, ++cell2)
+      {
+        cell2->set_active_fe_index(cell1->active_fe_index());
+      }
+  }
+
+  /**
+  Initializes the unique pointers which contain the necessary fe_values objects
+  for computing the cell and edge residuals
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::initialize_error_estimation_data()
+  {
+    // initialize the cell fe_values...
+    cell_hp_fe_values = std::make_unique<hp::FEValues<dim>>(
+      *DualSolver<dim>::fe_collection,
+      *DualSolver<dim>::quadrature_collection,
+      update_values | update_hessians | update_quadrature_points |
+        update_JxW_values);
+    face_hp_fe_values = std::make_unique<hp::FEFaceValues<dim>>(
+      *DualSolver<dim>::fe_collection,
+      *DualSolver<dim>::face_quadrature_collection,
+      update_values | update_gradients | update_JxW_values |
+        update_normal_vectors);
+    face_hp_fe_values_neighbor = std::make_unique<hp::FEFaceValues<dim>>(
+      *DualSolver<dim>::fe_collection,
+      *DualSolver<dim>::face_quadrature_collection,
+      update_values | update_gradients | update_JxW_values |
+        update_normal_vectors);
+    subface_hp_fe_values = std::make_unique<hp::FESubfaceValues<dim>>(
+      *DualSolver<dim>::fe_collection,
+      *DualSolver<dim>::face_quadrature_collection,
+      update_gradients);
+  }
+
+  /**
+  Since any scalar multiple of an eigenvector is also an eigenvector, we must
+  choose some normalization strategy. For convenience in the QoI expression, the
+  L2 norm is taken in this case.
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::normalize_solutions(
+    Vector<double> &primal_solution,
+    Vector<double> &dual_weights)
+  {
+    double sum_primal = 0.0, sum_dual = 0.0;
+    for (const auto &cell :
+         DualSolver<dim>::dof_handler.active_cell_iterators())
+      {
+        cell_hp_fe_values->reinit(cell);
+
+        // grab the fe_values object
+        const FEValues<dim> &fe_values =
+          cell_hp_fe_values->get_present_fe_values();
+
+        std::vector<Vector<double>> cell_primal_values(
+          fe_values.n_quadrature_points, Vector<double>(dim)),
+          cell_dual_values(fe_values.n_quadrature_points, Vector<double>(dim));
+        fe_values.get_function_values(primal_solution, cell_primal_values);
+        fe_values.get_function_values(dual_weights, cell_dual_values);
+
+
+        for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
+          {
+            sum_primal +=
+              cell_primal_values[p] * cell_primal_values[p] * fe_values.JxW(p);
+            sum_dual +=
+              cell_dual_values[p] * cell_dual_values[p] * fe_values.JxW(p);
+          }
+      }
+
+    primal_solution /= sqrt(sum_primal);
+    dual_weights /= sqrt(sum_dual);
+  }
+
+  /**
+  Serves as the main control function for estimating all of the error
+  contribution estimates
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::estimate_error(
+    Vector<double> &error_indicators)
+  {
+    // The constraints could be grabbed directly, but this is simple
+    AffineConstraints<double> primal_hanging_node_constraints;
+    DoFTools::make_hanging_node_constraints(PrimalSolver<dim>::dof_handler,
+                                            primal_hanging_node_constraints);
+    primal_hanging_node_constraints.close();
+
+    AffineConstraints<double> dual_hanging_node_constraints;
+    DoFTools::make_hanging_node_constraints(DualSolver<dim>::dof_handler,
+                                            dual_hanging_node_constraints);
+    dual_hanging_node_constraints.close();
+
+    // First map the primal solution to the space of the dual solution
+    // This allows us to use just one set of FEValues objects (rather than one
+    // set for the primal, one for dual)
+
+    Vector<double> primal_solution(DualSolver<dim>::dof_handler.n_dofs());
+
+    embed(PrimalSolver<dim>::dof_handler,
+          DualSolver<dim>::dof_handler,
+          dual_hanging_node_constraints,
+          *(PrimalSolver<dim>::get_solution()),
+          primal_solution);
+
+    Vector<double> &dual_solution = *(DualSolver<dim>::get_solution());
+
+    normalize_solutions(primal_solution, dual_solution);
+
+    Vector<double> dual_weights(DualSolver<dim>::dof_handler.n_dofs()),
+      dual_weights_interm(PrimalSolver<dim>::dof_handler.n_dofs());
+
+    // First extract the dual solution to the space of the primal
+    extract(DualSolver<dim>::dof_handler,
+            PrimalSolver<dim>::dof_handler,
+            primal_hanging_node_constraints,
+            *(DualSolver<dim>::get_solution()),
+            dual_weights_interm);
+
+    // Now embed this back to the space of the dual solution
+    embed(PrimalSolver<dim>::dof_handler,
+          DualSolver<dim>::dof_handler,
+          dual_hanging_node_constraints,
+          dual_weights_interm,
+          dual_weights);
+
+
+    // Subtract this from the full dual solution
+    dual_weights -= *(DualSolver<dim>::get_solution());
+    dual_weights *= -1.0;
+
+    *(DualSolver<dim>::get_solution()) -= primal_solution;
+
+    FaceIntegrals face_integrals;
+    for (const auto &cell :
+         DualSolver<dim>::dof_handler.active_cell_iterators())
+      for (const auto &face : cell->face_iterators())
+        face_integrals[face] = -1e20;
+
+
+    for (const auto &cell :
+         DualSolver<dim>::dof_handler.active_cell_iterators())
+      {
+        estimate_on_one_cell(cell,
+                             primal_solution,
+                             dual_weights,
+                             *(PrimalSolver<dim>::get_lambda_h()),
+                             error_indicators,
+                             face_integrals);
+      }
+    unsigned int present_cell = 0;
+    for (const auto &cell :
+         DualSolver<dim>::dof_handler.active_cell_iterators())
+      {
+        for (const auto &face : cell->face_iterators())
+          {
+            Assert(face_integrals.find(face) != face_integrals.end(),
+                   ExcInternalError());
+            error_indicators(present_cell) -= 0.5 * face_integrals[face];
+          }
+        ++present_cell;
+      }
+
+    // Now, with the error indicators computed, let us produce the
+    // estimate of the QoI error
+    this->qoi_error_estimate =
+      this->get_global_QoI_error(*(DualSolver<dim>::get_solution()),
+                                 error_indicators);
+    std::cout << "Estimated QoI error: " << std::setprecision(20)
+              << qoi_error_estimate << std::endl;
+  }
+
+
+  /**
+  Accumulates the error contribution estimates for one cell
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::estimate_on_one_cell(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    const Vector<double> &                                primal_solution,
+    const Vector<double> &                                dual_weights,
+    const double &                                        lambda_h,
+    Vector<double> &                                      error_indicators,
+    FaceIntegrals &                                       face_integrals)
+  {
+    integrate_over_cell(
+      cell, primal_solution, dual_weights, lambda_h, error_indicators);
+    for (unsigned int face_no : GeometryInfo<dim>::face_indices())
+      {
+        if (cell->face(face_no)->at_boundary())
+          {
+            face_integrals[cell->face(face_no)] = 0.0;
+            continue;
+          }
+        if ((cell->neighbor(face_no)->has_children() == false) &&
+            (cell->neighbor(face_no)->level() == cell->level()) &&
+            (cell->neighbor(face_no)->index() < cell->index()))
+          continue;
+        if (cell->at_boundary(face_no) == false)
+          if (cell->neighbor(face_no)->level() < cell->level())
+            continue;
+        if (cell->face(face_no)->has_children() == false)
+          integrate_over_regular_face(
+            cell, face_no, primal_solution, dual_weights, face_integrals);
+        else
+          integrate_over_irregular_face(
+            cell, face_no, primal_solution, dual_weights, face_integrals);
+      }
+  }
+
+  /**
+  Computes the cell residual
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::integrate_over_cell(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    const Vector<double> &                                primal_solution,
+    const Vector<double> &                                dual_weights,
+    const double &                                        lambda_h,
+    Vector<double> &                                      error_indicators)
+  {
+    cell_hp_fe_values->reinit(cell);
+    // Grab the fe_values object
+    const FEValues<dim> &fe_values = cell_hp_fe_values->get_present_fe_values();
+    std::vector<std::vector<Tensor<2, dim, double>>> cell_hessians(
+      fe_values.n_quadrature_points, std::vector<Tensor<2, dim, double>>(dim));
+    std::vector<Vector<double>> cell_primal_values(
+      fe_values.n_quadrature_points, Vector<double>(dim)),
+      cell_dual_values(fe_values.n_quadrature_points, Vector<double>(dim));
+    fe_values.get_function_values(primal_solution, cell_primal_values);
+    fe_values.get_function_hessians(primal_solution, cell_hessians);
+    fe_values.get_function_values(dual_weights, cell_dual_values);
+
+
+
+    double sum = 0.0;
+    for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
+      {
+        sum +=
+          (/*x-component*/ (cell_hessians[p][1][1][0] -
+                            cell_hessians[p][0][1][1]) *
+             (cell_dual_values[p](0)) +
+           /*y-component*/
+           (cell_hessians[p][0][0][1] - cell_hessians[p][1][0][0]) *
+             (cell_dual_values[p](1)) -
+           lambda_h * (cell_primal_values[p](0) * cell_dual_values[p](0) +
+                       cell_primal_values[p](1) * cell_dual_values[p](1))) *
+          fe_values.JxW(p);
+      }
+
+    error_indicators(cell->active_cell_index()) += sum;
+  }
+
+  /**
+  Computes the edge residual when there are no hanging nodes
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::integrate_over_regular_face(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    const unsigned int &                                  face_no,
+    const Vector<double> &                                primal_solution,
+    const Vector<double> &                                dual_weights,
+    FaceIntegrals &                                       face_integrals)
+  {
+    Assert(cell->neighbor(face_no).state() == IteratorState::valid,
+           ExcInternalError());
+    const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
+    const auto         neighbor          = cell->neighbor(face_no);
+
+    const unsigned int quadrature_index =
+      std::max(cell->active_fe_index(), neighbor->active_fe_index());
+    face_hp_fe_values->reinit(cell, face_no, quadrature_index);
+    const FEFaceValues<dim> &fe_face_values_cell =
+      face_hp_fe_values->get_present_fe_values();
+    std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
+      fe_face_values_cell.n_quadrature_points,
+      std::vector<Tensor<1, dim, double>>(dim)),
+      neighbor_primal_grads(fe_face_values_cell.n_quadrature_points,
+                            std::vector<Tensor<1, dim, double>>(dim));
+    fe_face_values_cell.get_function_gradients(primal_solution,
+                                               cell_primal_grads);
+
+    face_hp_fe_values_neighbor->reinit(neighbor,
+                                       neighbor_neighbor,
+                                       quadrature_index);
+    const FEFaceValues<dim> &fe_face_values_cell_neighbor =
+      face_hp_fe_values_neighbor->get_present_fe_values();
+    fe_face_values_cell_neighbor.get_function_gradients(primal_solution,
+                                                        neighbor_primal_grads);
+    const unsigned int n_q_points    = fe_face_values_cell.n_quadrature_points;
+    double             face_integral = 0.0;
+    std::vector<Vector<double>> cell_dual_values(n_q_points,
+                                                 Vector<double>(dim));
+    fe_face_values_cell.get_function_values(dual_weights, cell_dual_values);
+    for (unsigned int p = 0; p < n_q_points; ++p)
+      {
+        auto face_normal = fe_face_values_cell.normal_vector(p);
+
+        face_integral +=
+          (cell_primal_grads[p][1][0] - cell_primal_grads[p][0][1] -
+           neighbor_primal_grads[p][1][0] + neighbor_primal_grads[p][0][1]) *
+          (cell_dual_values[p][0] * face_normal[1] -
+           cell_dual_values[p][1] * face_normal[0]) *
+          fe_face_values_cell.JxW(p);
+      }
+    Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
+           ExcInternalError());
+    Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
+    face_integrals[cell->face(face_no)] = face_integral;
+  }
+
+  /**
+  Computes the residual when there are hanging nodes
+  */
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::integrate_over_irregular_face(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    const unsigned int &                                  face_no,
+    const Vector<double> &                                primal_solution,
+    const Vector<double> &                                dual_weights,
+    FaceIntegrals &                                       face_integrals)
+  {
+    const typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+    const typename DoFHandler<dim>::cell_iterator neighbor =
+      cell->neighbor(face_no);
+
+    Assert(neighbor.state() == IteratorState::valid, ExcInternalError());
+    Assert(neighbor->has_children(), ExcInternalError());
+    (void)neighbor;
+    const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
+    for (unsigned int subface_no = 0; subface_no < face->n_children();
+         ++subface_no)
+      {
+        const typename DoFHandler<dim>::active_cell_iterator neighbor_child =
+          cell->neighbor_child_on_subface(face_no, subface_no);
+        Assert(neighbor_child->face(neighbor_neighbor) ==
+                 cell->face(face_no)->child(subface_no),
+               ExcInternalError());
+        const unsigned int quadrature_index =
+          std::max(cell->active_fe_index(), neighbor_child->active_fe_index());
+        // initialize fe_subface values_cell
+        subface_hp_fe_values->reinit(cell,
+                                     face_no,
+                                     subface_no,
+                                     quadrature_index);
+        const FESubfaceValues<dim> &subface_fe_values_cell =
+          subface_hp_fe_values->get_present_fe_values();
+        std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
+          subface_fe_values_cell.n_quadrature_points,
+          std::vector<Tensor<1, dim, double>>(dim)),
+          neighbor_primal_grads(subface_fe_values_cell.n_quadrature_points,
+                                std::vector<Tensor<1, dim, double>>(dim));
+        subface_fe_values_cell.get_function_gradients(primal_solution,
+                                                      cell_primal_grads);
+        // initialize fe_face_values_neighbor
+        face_hp_fe_values_neighbor->reinit(neighbor_child,
+                                           neighbor_neighbor,
+                                           quadrature_index);
+        const FEFaceValues<dim> &face_fe_values_neighbor =
+          face_hp_fe_values_neighbor->get_present_fe_values();
+        face_fe_values_neighbor.get_function_gradients(primal_solution,
+                                                       neighbor_primal_grads);
+        const unsigned int n_q_points =
+          subface_fe_values_cell.n_quadrature_points;
+        std::vector<Vector<double>> cell_dual_values(n_q_points,
+                                                     Vector<double>(dim));
+        face_fe_values_neighbor.get_function_values(dual_weights,
+                                                    cell_dual_values);
+
+        double face_integral = 0.0;
+
+        for (unsigned int p = 0; p < n_q_points; ++p)
+          {
+            auto face_normal = face_fe_values_neighbor.normal_vector(p);
+            face_integral +=
+              (cell_primal_grads[p][0][1] - cell_primal_grads[p][1][0] +
+               neighbor_primal_grads[p][1][0] -
+               neighbor_primal_grads[p][0][1]) *
+              (cell_dual_values[p][0] * face_normal[1] -
+               cell_dual_values[p][1] * face_normal[0]) *
+              face_fe_values_neighbor.JxW(p);
+          }
+        face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
+      }
+    double sum = 0.0;
+    for (unsigned int subface_no = 0; subface_no < face->n_children();
+         ++subface_no)
+      {
+        Assert(face_integrals.find(face->child(subface_no)) !=
+                 face_integrals.end(),
+               ExcInternalError());
+        Assert(face_integrals[face->child(subface_no)] != -1e20,
+               ExcInternalError());
+        sum += face_integrals[face->child(subface_no)];
+      }
+    face_integrals[face] = sum;
+  }
+
+  template <int dim, bool report_dual>
+  double
+  DualWeightedResidual<dim, report_dual>::get_global_QoI_error(
+    Vector<double> &dual_solution,
+    Vector<double> &error_indicators)
+  {
+    auto dual_less_primal =
+      dual_solution; // Note: We have already extracted the primal solution...
+
+
+    double scaling_factor = 0.0;
+    for (const auto &cell :
+         DualSolver<dim>::dof_handler.active_cell_iterators())
+      {
+        cell_hp_fe_values->reinit(cell);
+        // grab the fe_values object
+        const FEValues<dim> &fe_values =
+          cell_hp_fe_values->get_present_fe_values();
+
+        std::vector<Vector<double>> cell_values(fe_values.n_quadrature_points,
+                                                Vector<double>(dim));
+        fe_values.get_function_values(dual_less_primal, cell_values);
+
+        for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
+          {
+            scaling_factor +=
+              (cell_values[p] * cell_values[p]) * fe_values.JxW(p);
+          }
+      }
+    double global_QoI_error = 0.0;
+    for (const auto &indicator : error_indicators)
+      {
+        global_QoI_error += indicator;
+      }
+
+    global_QoI_error /= (1 - 0.5 * scaling_factor);
+    return global_QoI_error;
+  }
+
+
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::embed(
+    const DoFHandler<dim> &          dof1,
+    const DoFHandler<dim> &          dof2,
+    const AffineConstraints<double> &constraints,
+    const Vector<double> &           solution,
+    Vector<double> &                 u2)
+  {
+    assert(u2.size() == dof2.n_dofs() && "Incorrect input vector size!");
+
+    u2 = 0.0;
+
+    typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
+                                                   endc1 = dof1.end();
+    typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
+
+    for (; cell1 < endc1; ++cell1, ++cell2)
+      {
+        const auto &fe1 =
+          dynamic_cast<const FE_Nedelec<dim> &>(cell1->get_fe());
+        const auto &fe2 =
+          dynamic_cast<const FE_Nedelec<dim> &>(cell2->get_fe());
+
+        assert(fe1.degree < fe2.degree && "Incorrect usage of embed!");
+
+        // Get the embedding_dofs
+
+
+        std::vector<unsigned int> embedding_dofs =
+          fe2.get_embedding_dofs(fe1.degree);
+        const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
+
+
+        Vector<double> local_dof_values_1;
+        Vector<double> local_dof_values_2(dofs_per_cell2);
+
+        local_dof_values_1.reinit(fe1.dofs_per_cell);
+        cell1->get_dof_values(solution, local_dof_values_1);
+
+        for (unsigned int i = 0; i < local_dof_values_1.size(); ++i)
+          local_dof_values_2[embedding_dofs[i]] = local_dof_values_1[i];
+
+        // Now set this changes to the global vector
+        cell2->set_dof_values(local_dof_values_2, u2);
+      }
+
+    u2.compress(VectorOperation::insert);
+    // Applies the constraints of the target finite element space
+    constraints.distribute(u2);
+  }
+
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::extract(
+    const DoFHandler<dim> &          dof1,
+    const DoFHandler<dim> &          dof2,
+    const AffineConstraints<double> &constraints,
+    const Vector<double> &           solution,
+    Vector<double> &                 u2)
+  {
+    // Maps from fe1 to fe2
+    assert(u2.size() == dof2.n_dofs() && "Incorrect input vector size!");
+
+    u2 = 0.0;
+
+    typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
+                                                   endc1 = dof1.end();
+    typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
+
+    for (; cell1 < endc1; ++cell1, ++cell2)
+      {
+        const auto &fe1 =
+          dynamic_cast<const FE_Nedelec<dim> &>(cell1->get_fe());
+        const auto &fe2 =
+          dynamic_cast<const FE_Nedelec<dim> &>(cell2->get_fe());
+
+        assert(fe1.degree > fe2.degree && "Incorrect usage of extract!");
+
+        // Get the embedding_dofs
+        std::vector<unsigned int> embedding_dofs =
+          fe1.get_embedding_dofs(fe2.degree);
+        const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
+
+
+        Vector<double> local_dof_values_1;
+        Vector<double> local_dof_values_2(dofs_per_cell2);
+
+        local_dof_values_1.reinit(fe1.dofs_per_cell);
+        cell1->get_dof_values(solution, local_dof_values_1);
+
+        for (unsigned int i = 0; i < local_dof_values_2.size(); ++i)
+          local_dof_values_2[i] = local_dof_values_1[embedding_dofs[i]];
+
+        // Now set this changes to the global vector
+        cell2->set_dof_values(local_dof_values_2, u2);
+      }
+
+    u2.compress(VectorOperation::insert);
+    // Applies the constraints of the target finite element space
+    constraints.distribute(u2);
+  }
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::output_eigenvalue_data(
+    std::ofstream &os)
+  {
+    os << (*this->get_primal_eigenvalues())[0] << " "
+       << (this->get_primal_DoFHandler())->n_dofs() << " "
+       << (*this->get_dual_eigenvalues())[0] << " "
+       << (this->get_dual_DoFHandler())->n_dofs() << std::endl;
+  }
+  template <int dim, bool report_dual>
+  void
+  DualWeightedResidual<dim, report_dual>::output_qoi_error_estimates(
+    std::ofstream &os)
+  {
+    os << qoi_error_estimate << std::endl;
+  }
+
+  /**
+  Provides a secondary error estimator, based on the Kelly error indicator.
+  Requires only the primal solver.
+  */
+  template <int dim>
+  class KellyErrorIndicator : public PrimalSolver<dim>
+  {
+  public:
+    std::string
+    name() const
+    {
+      return "Kelly";
+    }
+    void
+    output_eigenvalue_data(std::ofstream &os);
+    void
+    output_qoi_error_estimates(std::ofstream &);
+    KellyErrorIndicator(const std::string & prm_file,
+                        Triangulation<dim> &coarse_grid,
+                        const unsigned int &min_degree,
+                        const unsigned int &max_degree,
+                        const unsigned int &starting_degree);
+
+    virtual unsigned int
+    solve_problem() override;
+
+    virtual void
+    output_solution() override;
+
+    hp::FECollection<dim> *
+    get_FECollection();
+
+    hp::FECollection<dim> *
+    get_primal_FECollection();
+
+    std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+    get_eigenfunctions();
+
+    std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+    get_primal_eigenfunctions();
+
+    std::unique_ptr<std::vector<double>> &
+    get_primal_eigenvalues();
+
+
+    void
+    synchronize_discretization();
+
+    DoFHandler<dim> *
+    get_DoFHandler();
+
+    DoFHandler<dim> *
+    get_primal_DoFHandler();
+
+    unsigned int
+    get_max_degree()
+    {
+      return PrimalSolver<dim>::fe_collection->max_degree();
+    }
+    double qoi_error_estimate = 0;
+
+  protected:
+    void
+    estimate_error(Vector<double> &error_indicators);
+
+  private:
+    void
+    prune_eigenpairs(const double &TOL);
+
+    std::vector<const PETScWrappers::MPI::Vector *> eigenfunction_ptrs;
+    std::vector<const double *>                     eigenvalue_ptrs;
+
+    std::vector<std::shared_ptr<Vector<float>>> errors;
+  };
+
+  template <int dim>
+  KellyErrorIndicator<dim>::KellyErrorIndicator(
+    const std::string & prm_file,
+    Triangulation<dim> &coarse_grid,
+    const unsigned int &min_degree,
+    const unsigned int &max_degree,
+    const unsigned int &starting_degree)
+    : Base<dim>(prm_file, coarse_grid)
+    , PrimalSolver<dim>(prm_file,
+                        coarse_grid,
+                        min_degree,
+                        max_degree,
+                        starting_degree)
+  {}
+
+  template <int dim>
+  unsigned int
+  KellyErrorIndicator<dim>::solve_problem()
+  {
+    return PrimalSolver<dim>::solve_problem();
+  }
+
+  template <int dim>
+  hp::FECollection<dim> *
+  KellyErrorIndicator<dim>::get_FECollection()
+  {
+    return &*(PrimalSolver<dim>::fe_collection);
+  }
+
+  template <int dim>
+  hp::FECollection<dim> *
+  KellyErrorIndicator<dim>::get_primal_FECollection()
+  {
+    return &*(PrimalSolver<dim>::fe_collection);
+  }
+
+  template <int dim>
+  std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+  KellyErrorIndicator<dim>::get_eigenfunctions()
+  {
+    return (PrimalSolver<dim>::eigenfunctions);
+  }
+
+  template <int dim>
+  std::unique_ptr<std::vector<double>> &
+  KellyErrorIndicator<dim>::get_primal_eigenvalues()
+  {
+    return PrimalSolver<dim>::eigenvalues;
+  }
+
+  template <int dim>
+  std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
+  KellyErrorIndicator<dim>::get_primal_eigenfunctions()
+  {
+    return (PrimalSolver<dim>::eigenfunctions);
+  }
+
+  template <int dim>
+  DoFHandler<dim> *
+  KellyErrorIndicator<dim>::get_DoFHandler()
+  {
+    return &(PrimalSolver<dim>::dof_handler);
+  }
+
+  template <int dim>
+  DoFHandler<dim> *
+  KellyErrorIndicator<dim>::get_primal_DoFHandler()
+  {
+    return &(PrimalSolver<dim>::dof_handler);
+  }
+
+  template <int dim>
+  void
+  KellyErrorIndicator<dim>::synchronize_discretization()
+  {
+    // This function does nothing for this error indicator
+    return;
+  }
+
+  template <int dim>
+  void
+  KellyErrorIndicator<dim>::output_solution()
+  {
+    PrimalSolver<dim>::output_solution();
+  }
+
+  template <int dim>
+  void
+  KellyErrorIndicator<dim>::prune_eigenpairs(const double &TOL)
+  {
+    unsigned int count = 0;
+    for (size_t eigenpair_index = 0;
+         eigenpair_index < this->eigenfunctions->size();
+         ++eigenpair_index)
+      {
+        if (count >= this->n_eigenpairs)
+          break;
+        if (abs((*this->eigenvalues)[eigenpair_index]) < TOL)
+          continue;
+
+        eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
+        eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
+      }
+  }
+
+  template <int dim>
+  void
+  KellyErrorIndicator<dim>::estimate_error(Vector<double> &error_indicators)
+  {
+    std::cout << "Marking cells via Kelly indicator..." << std::endl;
+    prune_eigenpairs(1e-9);
+    // deallocate the errors vector
+    errors.clear();
+    for (size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
+      {
+        errors.emplace_back(
+          new Vector<float>(this->triangulation->n_active_cells()));
+      }
+    std::vector<Vector<float> *> estimated_error_per_cell(
+      eigenfunction_ptrs.size());
+    for (size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
+      {
+        estimated_error_per_cell[i] = errors[i].get();
+      }
+
+    KellyErrorEstimator<dim>::estimate(this->dof_handler,
+                                       *this->face_quadrature_collection,
+                                       {},
+                                       eigenfunction_ptrs,
+                                       estimated_error_per_cell);
+
+    for (auto &error_vec : errors)
+      {
+        auto normalized_vec = *error_vec;
+        normalized_vec /= normalized_vec.l1_norm();
+
+        for (unsigned int i = 0; i < error_indicators.size(); ++i)
+          error_indicators(i) += double(normalized_vec(i));
+      }
+    std::cout << "...Done!" << std::endl;
+  }
+  template <int dim>
+  void
+  KellyErrorIndicator<dim>::output_eigenvalue_data(std::ofstream &os)
+  {
+    os << (*this->get_primal_eigenvalues())[0] << " "
+       << (this->get_primal_DoFHandler())->n_dofs() << std::endl;
+  }
+  template <int dim>
+  void
+  KellyErrorIndicator<dim>::output_qoi_error_estimates(std::ofstream &)
+  {
+    return;
+  }
+
+} // namespace ErrorIndicators
+
+/**
+Includes all of the classes needed for smoothness estimation.
+*/
+namespace RegularityIndicators
+{
+  using namespace dealii;
+
+  /* For the Legendre smoothness indicator*/
+  /* Adapted from M. Fehling's smoothness_estimator.cc*/
+  template <int dim>
+  class LegendreInfo
+  {};
+
+  template <>
+  class LegendreInfo<2>
+  {
+  public:
+    std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
+
+    hp::FECollection<2> *fe_collection = nullptr;
+    DoFHandler<2> *      dof_handler   = nullptr;
+
+    void
+    initialization()
+    {
+      assert(fe_collection != nullptr && dof_handler != nullptr &&
+             "A valid FECollection and DoFHandler must be accessible!");
+
+      legendre_u = std::make_unique<FESeries::Legendre<2>>(
+        SmoothnessEstimator::Legendre::default_fe_series(*fe_collection, 0));
+      legendre_v = std::make_unique<FESeries::Legendre<2>>(
+        SmoothnessEstimator::Legendre::default_fe_series(*fe_collection, 1));
+
+      legendre_u->precalculate_all_transformation_matrices();
+      legendre_v->precalculate_all_transformation_matrices();
+    }
+
+    template <class VectorType>
+    void
+    compute_coefficient_decay(const VectorType &   eigenfunction,
+                              std::vector<double> &smoothness_indicators)
+    {
+      // Compute the coefficients for the u and v components of the solution
+      // separately,
+      Vector<float> smoothness_u(smoothness_indicators.size()),
+        smoothness_v(smoothness_indicators.size());
+
+      SmoothnessEstimator::Legendre::coefficient_decay(*legendre_u,
+                                                       *dof_handler,
+                                                       eigenfunction,
+                                                       smoothness_u);
+
+      SmoothnessEstimator::Legendre::coefficient_decay(*legendre_v,
+                                                       *dof_handler,
+                                                       eigenfunction,
+                                                       smoothness_v);
+
+      for (unsigned int i = 0; i < smoothness_indicators.size(); ++i)
+        {
+          smoothness_indicators[i] = std::min(smoothness_u[i], smoothness_v[i]);
+        }
+    }
+  };
+
+  /**
+  Implements the LegendreIndicator for use with the Refiner
+  */
+  template <int dim>
+  class LegendreIndicator
+  {
+  public:
+    void
+    attach_FE_info_and_initialize(hp::FECollection<dim> *fe_ptr,
+                                  DoFHandler<dim> *      dof_ptr);
+
+  protected:
+    template <class VectorType>
+    void
+    estimate_smoothness(
+      const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
+      const unsigned int &                            index_of_goal,
+      std::vector<double> &                           smoothness_indicators);
+
+  private:
+    LegendreInfo<dim> legendre;
+  };
+
+  template <int dim>
+  void
+  LegendreIndicator<dim>::attach_FE_info_and_initialize(
+    hp::FECollection<dim> *fe_ptr,
+    DoFHandler<dim> *      dof_ptr)
+  {
+    legendre.fe_collection = fe_ptr;
+    legendre.dof_handler   = dof_ptr;
+    this->legendre.initialization();
+  }
+
+  template <int dim>
+  template <class VectorType>
+  void
+  LegendreIndicator<dim>::estimate_smoothness(
+    const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
+    const unsigned int &                            index_of_goal,
+    std::vector<double> &                           smoothness_indicators)
+  {
+    this->legendre.compute_coefficient_decay((*eigenfunctions)[index_of_goal],
+                                             smoothness_indicators);
+  }
+} // namespace RegularityIndicators
+
+/**
+The final namespace, which combines the error estimation/indication and
+smoothness estimation functionality to conduct refinement.
+*/
+namespace Refinement
+{
+  using namespace dealii;
+  using namespace Maxwell;
+
+  template <int dim, class ErrorIndicator, class RegularityIndicator>
+  class Refiner : public ErrorIndicator, public RegularityIndicator
+  {
+  public:
+    Refiner(const std::string & prm_file,
+            Triangulation<dim> &coarse_grid,
+            const unsigned int &min_degree,
+            const unsigned int &max_degree,
+            const unsigned int &starting_degree);
+
+    void
+    execute_refinement(const double &smoothness_threshold_fraction);
+
+    virtual void
+    output_solution() override;
+
+  private:
+    Vector<double>      estimated_error_per_cell;
+    std::vector<double> smoothness_indicators;
+    std::ofstream       eigenvalues_out;
+    std::ofstream       error_estimate_out;
+  };
+
+  template <int dim, class ErrorIndicator, class RegularityIndicator>
+  Refiner<dim, ErrorIndicator, RegularityIndicator>::Refiner(
+    const std::string & prm_file,
+    Triangulation<dim> &coarse_grid,
+    const unsigned int &min_degree,
+    const unsigned int &max_degree,
+    const unsigned int &starting_degree)
+    : Base<dim>(prm_file, coarse_grid)
+    , ErrorIndicator(prm_file,
+                     coarse_grid,
+                     min_degree,
+                     max_degree,
+                     starting_degree)
+    , RegularityIndicator()
+  {
+    if (ErrorIndicator::name() == "DWR")
+      {
+        error_estimate_out.open("error_estimate.txt");
+        error_estimate_out << std::setprecision(20);
+      }
+
+    eigenvalues_out.open("eigenvalues_" + ErrorIndicator::name() + "_out.txt");
+    eigenvalues_out << std::setprecision(20);
+  }
+
+  // For generating samples of the curl of the electric field
+  template <int dim>
+  class CurlPostprocessor : public DataPostprocessorScalar<dim>
+  {
+  public:
+    CurlPostprocessor()
+      : DataPostprocessorScalar<dim>("Curl", update_gradients)
+    {}
+
+    virtual void
+    evaluate_vector_field(
+      const DataPostprocessorInputs::Vector<dim> &input_data,
+      std::vector<Vector<double>> &computed_quantities) const override
+    {
+      AssertDimension(input_data.solution_gradients.size(),
+                      computed_quantities.size());
+      for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
+        {
+          computed_quantities[p](0) = input_data.solution_gradients[p][1][0] -
+                                      input_data.solution_gradients[p][0][1];
+        }
+    }
+  };
+
+  /**
+  Overrides the output_solution function in order to include the error
+  estimation and smoothness estimation information. In the case of outputting
+  the eigenfunction, the PrimalSolver result is taken.
+
+  TODO: Extend to multiple eigenpairs
+  */
+  template <int dim, class ErrorIndicator, class RegularityIndicator>
+  void
+  Refiner<dim, ErrorIndicator, RegularityIndicator>::output_solution()
+  {
+    CurlPostprocessor<dim> curl_u;
+
+    DataOut<dim> data_out;
+    auto &       output_dof = *(ErrorIndicator::get_primal_DoFHandler());
+    data_out.attach_dof_handler(output_dof);
+    Vector<double> fe_degrees(this->triangulation->n_active_cells());
+    for (const auto &cell : output_dof.active_cell_iterators())
+      fe_degrees(cell->active_cell_index()) =
+        (*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
+          .degree;
+    data_out.add_data_vector(fe_degrees, "fe_degree");
+    //
+    data_out.add_data_vector(estimated_error_per_cell, "error");
+    Vector<double> smoothness_out(this->triangulation->n_active_cells());
+    for (const auto &cell : output_dof.active_cell_iterators())
+      {
+        auto i = cell->active_cell_index();
+        if (!cell->refine_flag_set() && !cell->coarsen_flag_set())
+          smoothness_out(i) = -1;
+        else
+          smoothness_out(i) = smoothness_indicators[i];
+      }
+    data_out.add_data_vector(smoothness_out, "smoothness");
+    data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
+                             std::string("eigenfunction_no_") +
+                               Utilities::int_to_string(0));
+    data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
+                             curl_u);
+
+    ErrorIndicator::output_eigenvalue_data(eigenvalues_out);
+    ErrorIndicator::output_qoi_error_estimates(error_estimate_out);
+
+    std::cout << "Number of DoFs: " << (this->get_primal_DoFHandler())->n_dofs()
+              << std::endl;
+
+
+    data_out.build_patches();
+    std::ofstream output("eigenvectors-" + ErrorIndicator::name() + "-" +
+                         std::to_string(this->refinement_cycle) + +".vtu");
+    data_out.write_vtu(output);
+  }
+
+
+  /**
+  Solves the problem (provided by the ErrorIndicator) and estimates the
+  smoothness. For cells marked for refinement, if the smoothness_threshold
+  is exceeded, $p$-refinement is chosen, otherwise $h$-refinement is chosen.
+  */
+  template <int dim, class ErrorIndicator, class RegularityIndicator>
+  void
+  Refiner<dim, ErrorIndicator, RegularityIndicator>::execute_refinement(
+    const double &smoothness_threshold_fraction)
+  {
+    // First initialize the RegularityIndicator...
+    // Depending on the limits set, this may take a while
+    std::cout << "Initializing RegularityIndicator..." << std::endl;
+    std::cout
+      << "(This may take a while if the max expansion order is set too high)"
+      << std::endl;
+    RegularityIndicator::attach_FE_info_and_initialize(
+      ErrorIndicator::get_FECollection(), ErrorIndicator::get_DoFHandler());
+    std::cout << "Done!" << std::endl << "Starting Refinement..." << std::endl;
+
+    for (unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
+      {
+        this->set_refinement_cycle(cycle);
+        std::cout << "Cycle: " << cycle << std::endl;
+        ErrorIndicator::solve_problem();
+        this->estimated_error_per_cell.reinit(
+          this->triangulation->n_active_cells());
+
+        ErrorIndicator::estimate_error(estimated_error_per_cell);
+
+        // Depending on the source of the error estimation/indication, these
+        // values might be signed, so we address that with the following
+        for (double &error_indicator : estimated_error_per_cell)
+          error_indicator = std::abs(error_indicator);
+
+
+        GridRefinement::refine_and_coarsen_fixed_number(
+          *this->triangulation, estimated_error_per_cell, 1. / 5., 0.000);
+
+        // Now get regularity indicators
+        // For those elements which must be refined, swap to increasing $p$
+        // depending on the regularity threshold...
+
+        smoothness_indicators =
+          std::vector<double>(this->triangulation->n_active_cells(),
+                              std::numeric_limits<double>::max());
+        if (ErrorIndicator::PrimalSolver::min_degree !=
+            ErrorIndicator::PrimalSolver::max_degree)
+          RegularityIndicator::estimate_smoothness(
+            ErrorIndicator::get_eigenfunctions(), 0, smoothness_indicators);
+        // save data
+        this->output_solution();
+        const double threshold_smoothness = smoothness_threshold_fraction;
+        unsigned int num_refined = 0, num_coarsened = 0;
+        if (ErrorIndicator::PrimalSolver::min_degree !=
+            ErrorIndicator::PrimalSolver::max_degree)
+          {
+            for (const auto &cell :
+                 ErrorIndicator::get_DoFHandler()->active_cell_iterators())
+              {
+                if (cell->refine_flag_set())
+                  ++num_refined;
+                if (cell->coarsen_flag_set())
+                  ++num_coarsened;
+                if (cell->refine_flag_set() &&
+                    smoothness_indicators[cell->active_cell_index()] >
+                      threshold_smoothness &&
+                    cell->active_fe_index() + 1 <
+                      ErrorIndicator::get_FECollection()->size())
+                  {
+                    cell->clear_refine_flag();
+                    cell->set_active_fe_index(cell->active_fe_index() + 1);
+                  }
+                else if (cell->coarsen_flag_set() &&
+                         smoothness_indicators[cell->active_cell_index()] <
+                           threshold_smoothness &&
+                         cell->active_fe_index() != 0)
+                  {
+                    cell->clear_coarsen_flag();
+
+                    cell->set_active_fe_index(cell->active_fe_index() - 1);
+                  }
+                // Here we also impose a limit on how small the cells can become
+                else if (cell->refine_flag_set() && cell->diameter() < 5.0e-6)
+                  {
+                    cell->clear_refine_flag();
+                    if (cell->active_fe_index() + 1 <
+                        ErrorIndicator::get_FECollection()->size())
+                      cell->set_active_fe_index(cell->active_fe_index() + 1);
+                  }
+              }
+          }
+
+        // Check what the smallest diameter is
+        double min_diameter = std::numeric_limits<double>::max();
+        for (const auto &cell :
+             ErrorIndicator::get_DoFHandler()->active_cell_iterators())
+          if (cell->diameter() < min_diameter)
+            min_diameter = cell->diameter();
+
+        std::cout << "Min diameter: " << min_diameter << std::endl;
+
+        ErrorIndicator::synchronize_discretization();
+
+        (this->triangulation)->execute_coarsening_and_refinement();
+      }
+  }
+} // namespace Refinement
+
+int
+main(int argc, char **argv)
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Maxwell;
+      using namespace Refinement;
+      using namespace ErrorIndicators;
+      using namespace RegularityIndicators;
+
+
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+
+      AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
+                  ExcMessage(
+                    "This program can only be run in serial, use ./maxwell-hp"));
+
+      Triangulation<2> triangulation_DWR, triangulation_Kelly;
+      Structures::create_L_waveguide(triangulation_DWR, 2.0);
+      Structures::create_L_waveguide(triangulation_Kelly, 2.0);
+
+      Refiner<2, KellyErrorIndicator<2>, LegendreIndicator<2>> problem_Kelly(
+        "maxwell-hp.prm", triangulation_Kelly, 2, 12, 2);
+
+      Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
+        problem_DWR("maxwell-hp.prm", triangulation_DWR, 2, 12, 2);
+
+      // The threshold for the hp-decision: too small -> not enough
+      // $h$-refinement, too large -> not enough $p$-refinement
+      double smoothness_threshold = 0.75;
+
+      std::cout << "Executing refinement for the Kelly strategy!" << std::endl;
+      problem_Kelly.execute_refinement(smoothness_threshold);
+      std::cout << "...Done with Kelly refinement strategy!" << std::endl;
+      std::cout << "Executing refinement for the DWR strategy!" << std::endl;
+      problem_DWR.execute_refinement(smoothness_threshold);
+      std::cout << "...Done with DWR refinement strategy!" << std::endl;
+    }
+
+  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;
+    }
+
+  std::cout << std::endl << "   Job done." << std::endl;
+
+  return 0;
+}

diff --git a/Maxwell-Eigenvalue-hp-Refinement/maxwell-hp.prm b/Maxwell-Eigenvalue-hp-Refinement/maxwell-hp.prm
new file mode 100644 (file)
index 0000000..4942739
--- /dev/null
+++ b/Maxwell-Eigenvalue-hp-Refinement/maxwell-hp.prm
@@ -0,0 +1,15 @@
+# Listing of Parameters
+# ---------------------
+
+#The type of eigenpairs to find (0 - smallest, 1 - target)
+set Eigenpair selection scheme = 1
+
+#The number of eigenvalues/eigenfunctions to be computed.
+set Number of eigenvalues/eigenfunctions = 1
+
+#The target eigenvalue (if scheme == 1)
+set Target eigenvalue = 23.344371957137
+
+#The number of cycles in refinement
+set Cycles number = 20
+