deal.II version GIT relicensing-2657-g43ecde61a7 2025-02-18 02:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The 'Goal-Oriented hp-Adaptivity for the Maxwell Eigenvalue Problem' code gallery program

This program was contributed by Jake J. Harmon <jake.harmon@ieee.org>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of Readme.md

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}\),

\begin{align*} \nabla\times(\mu_r^{-1}\nabla\times\mathbf{E})-k_0^2\varepsilon_r\mathbf{E} = 0 \text{ in } \Omega, \end{align*}

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 \text{ 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:

\begin{align*} \text{Find } U_{hp}=\left\{ \mathbf{u}_{hp},\,\lambda_{hp}\right\}\in V_{hp}\times \mathbb{R}_{>0} \text{ such that} \end{align*}

\begin{align*} a(\textbf{u}_{hp},\,\boldsymbol{\phi}_{hp}) = \lambda_{hp} m(\textbf{u}_{hp},\,\boldsymbol{\phi}_{hp}) \quad \forall\boldsymbol{\phi}\in V_{hp}, \end{align*}

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\). The finite dimensional subspace \(V_{hp}\) will 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.,

\begin{align*} e_{\lambda_{hp}} := \lambda-\lambda_{hp}. \end{align*}

Some comments on the discretization of Maxwell's equations

In the proper solution of variational problem, \(V_{hp}\) is not arbitrary, but should instead a subspace of \(H(\mathrm{curl};\,\Omega)\) or, with the boundary conditions indicated above, \(H_{0}(\mathrm{curl};\,\Omega)\), where

\begin{align*} 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\} \end{align*}

and

\begin{align*} H_{0}(\mathrm{curl};\,\Omega) = \left\{ \textbf{u}\in H(\mathrm{curl};\,\Omega) \, \mathrm{s.t.} \, \hat{\textbf{n}}\times\textbf{u} = 0 \text{ on @f$\partial\Omega@f$}\right\}. \end{align*}

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 is 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}\)

\begin{align*} \nabla\cdot\textbf{E} = 0 \end{align*}

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.

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).

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:

\begin{align*} 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}) \end{align*}

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

\begin{align*} \langle \textbf{u},\,\textbf{u}\rangle = 1, \end{align*}

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

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 (when using either the DWR or Kelly refinement strategies) and the adjoint problem (when using the DWR estimator)
  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

Problem geometry

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

Problem geometry

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.

Annotated version of maxwell-hp.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2022 by Jake J. Harmon
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  *
  * 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>
 
  #if DEAL_II_VERSION_GTE(9, 7, 0)
  #include <deal.II/grid/cell_data.h>
  #endif
  #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>
 
#define DEAL_II_VERSION_GTE(major, minor, subminor)
Definition config.h:349

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
  {
 
  double
  curlcurl(const ::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]);
  }
 
 
  template <int dim>
  inline double
  dot_term(const ::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
 
  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
 
  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:
  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;
  };
 
 
  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",
  "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");
 
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14889
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

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;
  }
 
 
  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();
 
  get_solution();
 
  void
  convert_solution();
 
  private:
  PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
  };
 
 
  template <int dim>
  EigenSolver<dim>::EigenSolver(const std::string & prm_file,
  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));
Definition hp.h:117
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)

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;
  cell1 = dof_handler.begin_active(),
  endc1 = dof_handler.end();
  for (; cell1 < endc1; ++cell1)
  {
  cell1->set_active_fe_index(start_diff);
  }
  }
  }
 
 
  template <int dim>
  double *
  EigenSolver<dim>::get_lambda_h()
  {
  return &(*eigenvalues)[0];
  }
 
 
  template <int dim>
  EigenSolver<dim>::get_solution()
  {
  return &solution;
  }
 
 
  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];
  }
 
 
  template <int dim>
  template <class SolverType>
  void
  EigenSolver<dim>::initialize_eigensolver(SolverType &eigensolver)
  {
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
typename ActiveSelector::active_cell_iterator active_cell_iterator
std::size_t size
Definition mpi.cc:739

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

  shift_scalar);
  SLEPcWrappers::TransformationShiftInvert spectral_transformation(
  this->mpi_communicator, additional_data);
 
  eigensolver.set_transformation(spectral_transformation);
  eigensolver.set_target_eigenvalue(this->target);
  }
 
 
  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,
  *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();
  }
 
 
  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);
  }
  }
 
 
  template <int dim>
  void
  EigenSolver<dim>::assemble_system()
  {
  hp::FEValues<dim> hp_fe_values(*fe_collection,
  *quadrature_collection,
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)

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)
  {
const FEValues< dim, spacedim > & get_present_fe_values() const

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);
 
  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);
  }
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57

since we have just set individual elements, we need the following

  stiffness_matrix.compress(VectorOperation::insert);
  mass_matrix.compress(VectorOperation::insert);
  }
 
 
  template <int dim>
  class PrimalSolver : public EigenSolver<dim>
  {
  public:
  PrimalSolver(const std::string & prm_file,
  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,
  const unsigned int &min_degree,
  const unsigned int &max_degree,
  const unsigned int &starting_degree)
  : Base<dim>(prm_file, triangulation)
  , EigenSolver<dim>(prm_file,
  min_degree,
  max_degree,
  starting_degree)
  {}
 
 
  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_") +
 
  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();
  }
 
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474

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,
  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,
  const unsigned int &min_degree,
  const unsigned int &max_degree,
  const unsigned int &starting_degree)
  : Base<dim>(prm_file, triangulation)
  , EigenSolver<dim>(prm_file,
  min_degree,
  max_degree,
  starting_degree)
  {}
 
  } // namespace Maxwell
 
  namespace ErrorIndicators
  {
  using namespace Maxwell;
 
 
  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,
  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);
 
  get_DoFHandler();
 
  get_primal_DoFHandler();
 
  get_dual_DoFHandler();
 
  get_FECollection();
 
  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);
  };
 
 
  template <int dim, bool report_dual>
  DualWeightedResidual<dim, report_dual>::DualWeightedResidual(
  const std::string & prm_file,
  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,
  min_primal_degree,
  max_primal_degree,
  starting_primal_degree)
  , DualSolver<dim>(prm_file,
  min_primal_degree + 1,
  max_primal_degree + 1,
  starting_primal_degree + 1)
  {
  initialize_error_estimation_data();
  }
 
 
  template <int dim, bool report_dual>
  DualWeightedResidual<dim, report_dual>::get_DoFHandler()
  {
  if (!report_dual)
  return &(PrimalSolver<dim>::dof_handler);
  else
  return &(DualSolver<dim>::dof_handler);
  }
 
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)

See above function, but to specifically output the primal DoFHandler...

  template <int dim, bool report_dual>
  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>
  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>
  DualWeightedResidual<dim, report_dual>::get_primal_FECollection()
  {
  return &*(PrimalSolver<dim>::fe_collection);
  }
 
  template <int dim, bool report_dual>
  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();
  }
 
 
  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();
  }
 
 
  template <int dim, bool report_dual>
  unsigned int
  DualWeightedResidual<dim, report_dual>::n_dofs() const
  {
  return PrimalSolver<dim>::n_dofs();
  }
 
 
  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());
  }
  }
 
 
  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,
  face_hp_fe_values = std::make_unique<hp::FEFaceValues<dim>>(
  *DualSolver<dim>::fe_collection,
  *DualSolver<dim>::face_quadrature_collection,
  face_hp_fe_values_neighbor = std::make_unique<hp::FEFaceValues<dim>>(
  *DualSolver<dim>::fe_collection,
  *DualSolver<dim>::face_quadrature_collection,
  subface_hp_fe_values = std::make_unique<hp::FESubfaceValues<dim>>(
  *DualSolver<dim>::fe_collection,
  *DualSolver<dim>::face_quadrature_collection,
  }
 
 
  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);
 
@ update_hessians
Second derivatives of shape functions.
@ update_normal_vectors
Normal vectors.

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);
  }
 
 
  template <int dim, bool report_dual>
  void
  DualWeightedResidual<dim, report_dual>::estimate_error(
  Vector<double> &error_indicators)
  {
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

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;
  }
 
#define Assert(cond, exc)

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;
  }
 
 
 
  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);
  }
  }
 
 
  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;
  }
 
 
  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;
  }
 
 
  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());
const FEFaceValues< dim, spacedim > & get_present_fe_values() const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
@ valid
Iterator points to a valid object.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)

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);
const FESubfaceValues< dim, spacedim > & get_present_fe_values() const

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);
T sum(const T &t, const MPI_Comm mpi_communicator)

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,
  {
  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,
  {

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;
  }
 
 
  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;
 
  get_FECollection();
 
  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();
 
  get_DoFHandler();
 
  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);
 
  #if DEAL_II_VERSION_GTE(9, 6, 0)
  std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
  #else
  std::vector<const PETScWrappers::MPI::Vector *> eigenfunction_ptrs;
  #endif
  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>
  KellyErrorIndicator<dim>::get_FECollection()
  {
  return &*(PrimalSolver<dim>::fe_collection);
  }
 
  template <int 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>
  KellyErrorIndicator<dim>::get_DoFHandler()
  {
  return &(PrimalSolver<dim>::dof_handler);
  }
 
  template <int 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);
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)

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();
  }
 
  #if DEAL_II_VERSION_GTE(9, 6, 0)
  const auto solution_view = make_array_view(eigenfunction_ptrs);
  auto error_view = make_array_view(estimated_error_per_cell);
  KellyErrorEstimator<dim>::estimate(this->dof_handler,
  *this->face_quadrature_collection,
  {},
  solution_view,
  error_view);
  #else
  *this->face_quadrature_collection,
  {},
  eigenfunction_ptrs,
  estimated_error_per_cell);
  #endif
 
  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
 
 
  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>>(
  legendre_v = std::make_unique<FESeries::Legendre<2>>(
 
  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)
  {
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
virtual size_type size() const override
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)

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());
 
  *dof_handler,
  eigenfunction,
  smoothness_u);
 
  *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]);
  }
  }
  };
 
 
  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(
  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
 
 
  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);
  }
 
void coefficient_decay(FESeries::Legendre< dim, spacedim > &fe_legendre, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const VectorTools::NormType regression_strategy=VectorTools::Linfty_norm, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
double legendre(unsigned int l, double x)
Definition cmath.h:65
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)

For generating samples of the curl of the electric field

  template <int dim>
  class CurlPostprocessor : public DataPostprocessorScalar<dim>
  {
  public:
  CurlPostprocessor()
  {}
 
  virtual void
  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];
  }
  }
  };
 
 
  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_") +
  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);
  }
 
 
 
  template <int dim, class ErrorIndicator, class RegularityIndicator>
  void
  Refiner<dim, ErrorIndicator, RegularityIndicator>::execute_refinement(
  const double &smoothness_threshold_fraction)
  {
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
#define AssertDimension(dim1, dim2)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)

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);
 
 
  *this->triangulation, estimated_error_per_cell, 1. / 5., 0.000);
 
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())

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 &&
  static_cast<unsigned int>(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 (static_cast<unsigned int>(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);
 
 
  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,
  /*Minimum Degree*/ 2,
  /*Maximum Degree*/ 5,
  /*Starting Degree*/ 2);
 
  Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
  problem_DWR("maxwell-hp.prm",
  triangulation_DWR,
  /*Minimum Degree*/ 2,
  /*Maximum Degree*/ 5,
  /*Starting Degree*/ 2);
 
#define AssertThrow(cond, exc)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:97

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;
  }