This program was contributed by Jake J. Harmon <[email protected]>.
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.
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:
- Solve the forward problem (when using either the DWR or Kelly refinement strategies) and the adjoint problem (when using the DWR estimator)
- Generate error contribution estimates (DWR) or error indicators (Kelly)
- Mark the top 20% elements with the largest refinement indicators by magnitude
- For those cells marked for refinement, estimate the local smoothness of the forward solution (Legendre decay rate)
- According to the Legendre decay rate, reclassify the refinement to \(p\) if necessary
- 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.
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.
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
#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/cell_data.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>
* * * struct InterferenceTaperTransform *
For parallelization (using WorkStream and Intel TBB)
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/work_stream.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>
curlcurl(const ::FEValues<2> &fe_values,
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]);
dot_term(const ::FEValues<dim> &fe_values,
const unsigned int & q_point)
{
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);
{
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();
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;
{
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();
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;
set_refinement_cycle(
const unsigned int cycle);
unsigned int refinement_cycle = 0;
std::unique_ptr<ParameterHandler> parameters;
unsigned int n_eigenpairs = 1;
unsigned int eigenpair_selection_scheme;
unsigned int max_cycles = 0;
ompi_communicator_t * mpi_communicator = PETSC_COMM_SELF;
: triangulation(&coarse_grid)
parameters->declare_entry(
"Eigenpair selection scheme",
"The type of eigenpairs to find (0 - smallest, 1 - target)");
parameters->declare_entry(
"Number of eigenvalues/eigenfunctions",
"The number of eigenvalues/eigenfunctions "
parameters->declare_entry(
"Target eigenvalue",
"The target eigenvalue (if scheme == 1)");
parameters->declare_entry(
"Cycles number",
"The number of cycles in refinement");
parameters->parse_input(prm_file);
eigenpair_selection_scheme =
parameters->get_integer(
"Eigenpair selection scheme");
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
void refine_global(const unsigned int times=1)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
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!");
parameters->get_integer(
"Number of eigenvalues/eigenfunctions");
"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)
Base<dim>::set_refinement_cycle(
const unsigned int cycle)
{
refinement_cycle = cycle;
class EigenSolver :
public virtual Base<dim>
EigenSolver(
const std::string & prm_file,
const unsigned int &minimum_degree,
const unsigned int &maximum_degree,
const unsigned int &starting_degree);
solve_problem()
override;
template <
class SolverType>
initialize_eigensolver(SolverType &eigensolver);
const std::unique_ptr<hp::FECollection<dim>> fe_collection;
std::unique_ptr<hp::QCollection<dim>> quadrature_collection;
const unsigned int max_degree, min_degree;
for the actual solution
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> eigenfunctions;
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)
for (unsigned
int degree = min_degree; degree <= max_degree; ++degree)
* * for(const auto &cell :triangulation.active_cell_iterators())
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Generate quadrature collection with sorted quadrature weights
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);
EigenSolver<dim>::get_lambda_h()
return &(*eigenvalues)[0];
EigenSolver<dim>::get_solution()
EigenSolver<dim>::convert_solution()
for (
unsigned int i = 0; i < solution.size(); ++i)
solution[i] = (*eigenfunctions)[0][i];
template <
class SolverType>
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
From the parameters class, initialize the eigensolver...
switch (this->eigenpair_selection_scheme)
eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
eigensolver.set_target_eigenvalue(this->target);
eigensolver.set_which_eigenpairs(EPS_SMALLEST_MAGNITUDE);
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
this->mpi_communicator, additional_data);
eigensolver.set_transformation(spectral_transformation);
eigensolver.set_target_eigenvalue(this->target);
EigenSolver<dim>::solve_problem()
initialize_eigensolver(eigensolver);
solve the problem
eigensolver.solve(stiffness_matrix,
for (
auto &entry : *eigenfunctions)
constraints.distribute(entry);
return solver_control.last_step();
EigenSolver<dim>::n_dofs() const
return dof_handler.n_dofs();
EigenSolver<dim>::setup_system()
dof_handler.distribute_dofs(*fe_collection);
eigenfunctions->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);
EigenSolver<dim>::assemble_system()
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.max_couplings_between_dofs());
mass_matrix.reinit(dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
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);
hp_fe_values.reinit(cell);
for (
unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
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) *
cell_mass_matrix(i, j) +=
(Operations::dot_term(fe_values, i, j, q_point)) *
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_stiffness_matrix,
constraints.distribute_local_to_global(cell_mass_matrix,
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
if (constraints.is_constrained(i))
stiffness_matrix.set(i, i, 10000.0);
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
since we have just set individual elements, we need the following
class PrimalSolver :
public EigenSolver<dim>
PrimalSolver(
const std::string & prm_file,
const unsigned int &min_degree,
const unsigned int &max_degree,
const unsigned int &starting_degree);
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,
PrimalSolver<dim>::output_solution()
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_") +
<< " 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;
data_out.build_patches();
std::ofstream output(
"eigenvectors-" +
std::to_string(this->refinement_cycle) +
".vtu");
data_out.write_vtu(output);
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)
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).
class DualSolver :
public EigenSolver<dim>
DualSolver(
const std::string & prm_file,
const unsigned int &min_degree,
const unsigned int &max_degree,
const unsigned int &starting_degree);
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,
namespace ErrorIndicators
template <
int dim,
bool report_dual>
class DualWeightedResidual :
public PrimalSolver<dim>,
public DualSolver<dim>
output_eigenvalue_data(std::ofstream &os);
output_qoi_error_estimates(std::ofstream &os);
DualWeightedResidual(
const std::string & prm_file,
const unsigned int &min_primal_degree,
const unsigned int &max_primal_degree,
const unsigned int &starting_primal_degree);
solve_problem()
override;
output_solution()
override;
get_primal_FECollection();
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
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>> &
synchronize_discretization();
return PrimalSolver<dim>::fe_collection->max_degree();
double qoi_error_estimate = 0;
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;
typename std::map<typename DoFHandler<dim>::face_iterator,
double>;
initialize_error_estimation_data();
FaceIntegrals & face_integrals);
integrate_over_regular_face(
const unsigned int & face_no,
FaceIntegrals & face_integrals);
integrate_over_irregular_face(
const unsigned int & face_no,
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,
, DualSolver<dim>(prm_file,
starting_primal_degree + 1)
initialize_error_estimation_data();
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::get_DoFHandler()
return &(PrimalSolver<dim>::dof_handler);
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()
return &*(PrimalSolver<dim>::fe_collection);
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()
return (PrimalSolver<dim>::eigenfunctions);
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>
DualWeightedResidual<dim, report_dual>::output_solution()
PrimalSolver<dim>::output_solution();
Solves the primal problem
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::solve_primal_problem()
return PrimalSolver<dim>::solve_problem();
Solves the dual problem
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::solve_dual_problem()
return DualSolver<dim>::solve_problem();
template <
int dim,
bool report_dual>
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>
DualWeightedResidual<dim, report_dual>::n_dofs() const
return PrimalSolver<dim>::n_dofs();
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::synchronize_discretization()
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);
for (; cell1 < endc1; ++cell1, ++cell2)
cell2->set_active_fe_index(cell1->active_fe_index());
template <
int dim,
bool report_dual>
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>
DualWeightedResidual<dim, report_dual>::normalize_solutions(
{
double sum_primal = 0.0, sum_dual = 0.0;
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
std::vector<Vector<double>> cell_primal_values(
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)
cell_primal_values[p] * cell_primal_values[p] * fe_values.JxW(p);
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>
DualWeightedResidual<dim, report_dual>::estimate_error(
{
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
The constraints could be grabbed directly, but this is simple
primal_hanging_node_constraints);
primal_hanging_node_constraints.close();
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()),
normalize_solutions(primal_solution, dual_solution);
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()),
Now embed this back to the space of the dual solution
embed(PrimalSolver<dim>::dof_handler,
DualSolver<dim>::dof_handler,
dual_hanging_node_constraints,
Subtract this from the full dual solution
dual_weights -= *(DualSolver<dim>::get_solution());
*(DualSolver<dim>::get_solution()) -= primal_solution;
FaceIntegrals face_integrals;
DualSolver<dim>::dof_handler.active_cell_iterators())
for (const auto &face : cell->face_iterators())
face_integrals[face] = -1e20;
DualSolver<dim>::dof_handler.active_cell_iterators())
estimate_on_one_cell(cell,
*(PrimalSolver<dim>::get_lambda_h()),
unsigned int present_cell = 0;
DualSolver<dim>::dof_handler.active_cell_iterators())
for (const auto &face : cell->face_iterators())
Assert(face_integrals.find(face) != face_integrals.
end(),
error_indicators(present_cell) -= 0.5 * face_integrals[face];
#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()),
std::cout <<
"Estimated QoI error: " << std::setprecision(20)
<< qoi_error_estimate <<
std::endl;
template <
int dim,
bool report_dual>
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,
Vector<
double> & error_indicators,
FaceIntegrals & face_integrals)
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;
if ((cell->neighbor(face_no)->has_children() ==
false) &&
(cell->neighbor(face_no)->level() == cell->level()) &&
(cell->neighbor(face_no)->index() < cell->index()))
if (cell->at_boundary(face_no) ==
false)
if (cell->neighbor(face_no)->level() < cell->level())
if (cell->face(face_no)->has_children() ==
false)
integrate_over_regular_face(
cell, face_no, primal_solution, dual_weights, face_integrals);
integrate_over_irregular_face(
cell, face_no, primal_solution, dual_weights, face_integrals);
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::integrate_over_cell(
{
cell_hp_fe_values->reinit(cell);
Grab the fe_values object
std::vector<std::vector<Tensor<2, dim, double>>> cell_hessians(
std::vector<Vector<double>> cell_primal_values(
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);
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
( (cell_hessians[p][1][1][0] -
cell_hessians[p][0][1][1]) *
(cell_dual_values[p](0)) +
(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))) *
error_indicators(cell->active_cell_index()) += sum;
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::integrate_over_regular_face(
const unsigned int & face_no,
FaceIntegrals & face_integrals)
{
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);
std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
fe_face_values_cell.n_quadrature_points,
neighbor_primal_grads(fe_face_values_cell.n_quadrature_points,
fe_face_values_cell.get_function_gradients(primal_solution,
face_hp_fe_values_neighbor->reinit(neighbor,
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,
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);
(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(),
Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
face_integrals[cell->face(face_no)] = face_integral;
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::integrate_over_irregular_face(
const unsigned int & face_no,
FaceIntegrals & face_integrals)
{
Assert(neighbor->has_children(), ExcInternalError());
const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
for (
unsigned int subface_no = 0; subface_no < face->n_children();
cell->neighbor_child_on_subface(face_no, subface_no);
Assert(neighbor_child->face(neighbor_neighbor) ==
cell->face(face_no)->child(subface_no),
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,
std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
subface_fe_values_cell.n_quadrature_points,
neighbor_primal_grads(subface_fe_values_cell.n_quadrature_points,
subface_fe_values_cell.get_function_gradients(primal_solution,
const FESubfaceValues< dim, spacedim > & get_present_fe_values() const
initialize fe_face_values_neighbor
face_hp_fe_values_neighbor->reinit(neighbor_child,
const unsigned int n_q_points =
subface_fe_values_cell.n_quadrature_points;
std::vector<Vector<double>> cell_dual_values(n_q_points,
face_fe_values_neighbor.get_function_values(dual_weights,
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);
(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;
for (
unsigned int subface_no = 0; subface_no < face->n_children();
Assert(face_integrals.find(face->child(subface_no)) !=
Assert(face_integrals[face->child(subface_no)] != -1e20,
sum += face_integrals[face->child(subface_no)];
face_integrals[face] =
sum;
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::get_global_QoI_error(
{
double scaling_factor = 0.0;
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
std::vector<Vector<double>> cell_values(fe_values.n_quadrature_points,
fe_values.get_function_values(dual_less_primal, cell_values);
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
(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);
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::embed(
{
assert(u2.size() == dof2.n_dofs() &&
"Incorrect input vector size!");
for (; cell1 < endc1; ++cell1, ++cell2)
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();
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);
Applies the constraints of the target finite element space
constraints.distribute(u2);
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::extract(
{
Maps from fe1 to fe2
assert(u2.size() == dof2.n_dofs() &&
"Incorrect input vector size!");
for (; cell1 < endc1; ++cell1, ++cell2)
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();
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);
Applies the constraints of the target finite element space
constraints.distribute(u2);
template <
int dim,
bool report_dual>
DualWeightedResidual<dim, report_dual>::output_eigenvalue_data(
{
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>
DualWeightedResidual<dim, report_dual>::output_qoi_error_estimates(
{
os << qoi_error_estimate << std::endl;
class KellyErrorIndicator :
public PrimalSolver<dim>
output_eigenvalue_data(std::ofstream &os);
output_qoi_error_estimates(std::ofstream &);
KellyErrorIndicator(
const std::string & prm_file,
const unsigned int &min_degree,
const unsigned int &max_degree,
const unsigned int &starting_degree);
solve_problem()
override;
output_solution()
override;
get_primal_FECollection();
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
get_primal_eigenfunctions();
std::unique_ptr<std::vector<double>> &
get_primal_eigenvalues();
synchronize_discretization();
return PrimalSolver<dim>::fe_collection->max_degree();
double qoi_error_estimate = 0;
prune_eigenpairs(
const double &TOL);
std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
std::vector<const double *> eigenvalue_ptrs;
std::vector<std::shared_ptr<Vector<float>>> errors;
KellyErrorIndicator<dim>::KellyErrorIndicator(
const std::string & prm_file,
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,
KellyErrorIndicator<dim>::solve_problem()
return PrimalSolver<dim>::solve_problem();
KellyErrorIndicator<dim>::get_FECollection()
return &*(PrimalSolver<dim>::fe_collection);
KellyErrorIndicator<dim>::get_primal_FECollection()
return &*(PrimalSolver<dim>::fe_collection);
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
KellyErrorIndicator<dim>::get_eigenfunctions()
return (PrimalSolver<dim>::eigenfunctions);
std::unique_ptr<std::vector<double>> &
KellyErrorIndicator<dim>::get_primal_eigenvalues()
return PrimalSolver<dim>::eigenvalues;
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
KellyErrorIndicator<dim>::get_primal_eigenfunctions()
return (PrimalSolver<dim>::eigenfunctions);
KellyErrorIndicator<dim>::get_DoFHandler()
return &(PrimalSolver<dim>::dof_handler);
KellyErrorIndicator<dim>::get_primal_DoFHandler()
return &(PrimalSolver<dim>::dof_handler);
KellyErrorIndicator<dim>::synchronize_discretization()
This function does nothing for this error indicator
KellyErrorIndicator<dim>::output_solution()
PrimalSolver<dim>::output_solution();
KellyErrorIndicator<dim>::prune_eigenpairs(
const double &TOL)
{
for (
size_t eigenpair_index = 0;
eigenpair_index < this->eigenfunctions->size();
if (count >= this->n_eigenpairs)
if (
abs((*this->eigenvalues)[eigenpair_index]) < TOL)
eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
KellyErrorIndicator<dim>::estimate_error(
Vector<double> &error_indicators)
{
std::cout <<
"Marking cells via Kelly indicator..." << std::endl;
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
deallocate the errors vector
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
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();
*this->face_quadrature_collection,
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;
KellyErrorIndicator<dim>::output_eigenvalue_data(std::ofstream &os)
{
os << (*this->get_primal_eigenvalues())[0] <<
" "
<< (this->get_primal_DoFHandler())->n_dofs() << std::endl;
KellyErrorIndicator<dim>::output_qoi_error_estimates(std::ofstream &)
{
namespace RegularityIndicators
std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
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>
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)
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,
smoothness_v(smoothness_indicators.size());
for (
unsigned int i = 0; i < smoothness_indicators.size(); ++i)
smoothness_indicators[i] =
std::min(smoothness_u[i], smoothness_v[i]);
template <
class VectorType>
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
const unsigned int & index_of_goal,
std::vector<double> & smoothness_indicators);
LegendreIndicator<dim>::attach_FE_info_and_initialize(
{
template <
class VectorType>
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],
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
class Refiner :
public ErrorIndicator,
public RegularityIndicator
Refiner(
const std::string & prm_file,
const unsigned int &min_degree,
const unsigned int &max_degree,
const unsigned int &starting_degree);
execute_refinement(
const double &smoothness_threshold_fraction);
output_solution()
override;
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,
const unsigned int &min_degree,
const unsigned int &max_degree,
const unsigned int &starting_degree)
: Base<dim>(prm_file, coarse_grid)
, ErrorIndicator(prm_file,
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);
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
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)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
For generating samples of the curl of the electric field
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>
Refiner<dim, ErrorIndicator, RegularityIndicator>::output_solution()
CurlPostprocessor<dim> curl_u;
auto & output_dof = *(ErrorIndicator::get_primal_DoFHandler());
for (
const auto &cell : output_dof.active_cell_iterators())
fe_degrees(cell->active_cell_index()) =
(*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
data_out.add_data_vector(fe_degrees,
"fe_degree");
data_out.add_data_vector(estimated_error_per_cell,
"error");
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) = 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],
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()
data_out.build_patches();
std::to_string(this->refinement_cycle) + +".vtu");
data_out.write_vtu(output);
template <
int dim, class ErrorIndicator, class RegularityIndicator>
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;
<<
"(This may take a while if the max expansion order is set too high)"
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...
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
const double threshold_smoothness = smoothness_threshold_fraction;
unsigned int num_refined = 0, num_coarsened = 0;
if (ErrorIndicator::PrimalSolver::min_degree !=
ErrorIndicator::PrimalSolver::max_degree)
ErrorIndicator::get_DoFHandler()->active_cell_iterators())
if (cell->refine_flag_set())
if (cell->coarsen_flag_set())
if (cell->refine_flag_set() &&
smoothness_indicators[cell->active_cell_index()] >
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()] <
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();
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();
main(
int argc,
char **argv)
{
using namespace Refinement;
using namespace ErrorIndicators;
using namespace RegularityIndicators;
ExcMessage(
"This program can only be run in serial, use ./maxwell-hp"));
Structures::create_L_waveguide(triangulation_DWR, 2.0);
Structures::create_L_waveguide(triangulation_Kelly, 2.0);
Refiner<2, KellyErrorIndicator<2>, LegendreIndicator<2>> problem_Kelly(
Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
problem_DWR(
"maxwell-hp.prm",
* * int main(int argc, char **argv)
#define AssertThrow(cond, exc)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
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 <<
"Exception on processing: " << std::endl
<< exc.what() << std::endl
<<
"Aborting!" << std::endl
<<
"----------------------------------------------------"
<<
"----------------------------------------------------"
std::cerr <<
"Unknown exception!" << std::endl
<<
"Aborting!" << std::endl
<<
"----------------------------------------------------"
std::cout << std::endl <<
" Job done." << std::endl;
1.9.8