This tutorial is very similar to the two dimensional hp-FEM tutorial in hpFEM2d.cc; you should have read it to understand this tutorial.
Maxwell's Equations
The equations solved are the time-harmonic Maxwell's equations:
with the divergence conditions (assuming no charge density)
Extracting from the second equation and inserting in into the first formally yields
the so-called electric source problem. The appropriate space is
However, the curl-curl form above is not defined for electric fields in . This is resolved in the variational formulation of the electric source problem:
Find such that
where is with the prefect conductor boundary conditions , and . The associated bilinear operator of this variational formulation is not elliptic and the divergence condition is an independent constraint. Normally, Maxwell's equations are discretised using Nedelec's elements. However, there are good reasons why one would like to use standard H1-conforming FEM. This is possible using weighted regularization.
Weighted Regularization
The divergence can be forced to zero by adding
to the bilinear form. However, this works only in non-convex domains. With the weighted regularization and the new bilinear form
where and
is dense in (the space we are discretising with H1-conforming FEM). A simple choice for the weight in three dimensions is
where is the set of reentrant corners and the set of reentrant edges. Then,
See also
Martin Costabel and Monique Dauge. Weighted regularization of Maxwell equations in polyhedral domains. A rehabilitation of nodal finite elements. Numer. Math., 93(2):239-277, 2002.
Maxwell Eigenproblem
The Maxwell Eigenproblem is solved using the weighted regularization:
Find frequencies and functions such that
As said above, this tutorial is very similar in large parts to hpFEM2d.cc and the documentation below is terser (ie. only new stuff is commented).
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <memory>
#include <fstream>
#include <unistd.h>
#include "basics.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "integration.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"
The new includes files are hp3D.hh for the three dimensional hp-FEM, vectorial.hh for the classes to build vector valued FE spaces etc. from the scalar components and eigensolver.hh for the Eigensolvers needed below.
#include "hp3D.hh"
#include "vectorial.hh"
#include "eigensolver.hh"
FEM Routine
Most abbreviations are the same (including the same meaning) as already known.
New parameters to configure the Eigensolvers are kmax (how many Eigenvalues should be computed at the same time), solverType (which type of Eigensolver, there are a few available) and shift (parameter for the shift-invert method).
The following parameters are used for the weighted regularization: rho is to distinguish physical from spurious Eigenvalues by the following criterion (using the Eigenfunction ):
The perfect conductor boundary conditions can be set up by carefully choosing the Dirichlet boundary conditions of the three scalar components of the vector valued FE space.
Finally, the FE spaces can be set up. First, a vector valued FE space with three components is created. This is essentially a container (as most of the classes of vectorial) which takes care of different scalar components which can be added to it.
Then, one after the other scalar FE space is created and added to the vector valued FE space spc. Note that the scalar FE spaces have to be rebuilt before adding them to spc. This call to hp3D::Space::rebuild creates the elements.
The stiffness matrix has two parts: rot-rot and div-div. These two parts are built separately. Before calling the Eigensolver, a linear combination depending on s of the two is created. The classes hp3D::RotRot and hp3D::DivDiv both come with a setup routine to take care of some technicalities. The use of this routine is not mandatory but recommended.
The rot-rot part of the stiffness matrix is simple to compute.
The div-div part of the stiffness matrix is somewhat more involved: the singularities of the solutions have to be taken into account by the weight. The singular edes and vertices are marked, this information is later used by the div-div bilinear form.
std::cout << " div div matrix: " << divdiv << std::endl;
The mass matrix consists of three scalar mass matrices in the diagonal blocks of the mass matrix. They are put together in the canonical way to build a vector valued bilnear form.
Each call to vectorial::BilinearForm::put adds a scalar bilinear form into a block of the vector valued bilinear form. In this case, only the diagonal blocks have to be filled.
The Eigensolver shall be constructed by a fabric, ie. a class derived from eigensolver::SolverFabric. solverType contains which solver to use and it is evaluated in the following switchcase clause.
After having chosen the solver by creating the appropriate fabric, the solver can be instantiated.
A call to eigensolver::Solver::getEV returns the Eigenvalues. They are printed to screen in a list.
Postprocessing
Then, a list of the Eigenvalues together with the categorization (the formula is given above) is printed to screen. In addition, the Eigenvalue together with the rotrot energy and the divdiv energy are stored in the output data for later postprocessing (or plotting).
Compute the rotrot energy:
Compute the divdiv energy:
Evaluate the criterion and print the results.
Finally, the data is stored in output.
Refinement of the Space
Each component of the vector valued space has to be refined individually towards the singularly marked vertices and edges.
Main Program
The main program works very similar to the previous tutorials.
Command Line Arguments
The arrays to store the Eigenvalues and the divdiv and rotrot energies have to be prepared.
Results
To restrict the amount of output, we are calling the program with the following input file:
solver = ArPack(shift-invert mode(3), computed largest (algebraic) eigenvalues, sigma = 9, n = 26, tol = 1.11022302462516e-16, conv. eigenpairs = 3(max = 3), Arnoldi iter = 13(max = 1000), 34 OP*x ops., 112 B*x ops., 33 steps of re-orth.)
EV 0 = 9.90411138346484 has rotEnergy = 9.9041, divEnergy*s = 7.1322e-32, L2 norm = 1
EV 1 = 14.167902552732 has rotEnergy = 13.806, divEnergy*s = 0.36222, L2 norm = 1
EV 2 = 14.5499375050742 has rotEnergy = 13.526, divEnergy*s = 1.0238, L2 norm = 1
--
Writing gathered data to disk: ./hpFEM3d-EV.out
The mesh normally used is a thick L shaped domain in three hexahedra. The mesh looks like this: With 18726 dofs, the Eigenvalues plotted versus s: The physical Eigenvalues (green) are lined up on horizontal lines whereas the spurious Eigenvalues are on straight lines through the origin. The exact values (horizontal blue lines) are taken from Monique Dauge's Benchmax page. Below are the convergence histories for the first three Eigenvalues for selected values of s: These convergence histories show exponential convergence: straight lines in a plot.
Complete Source Code
Author
Philipp Frauenfelder, 2004
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <memory>
#include <fstream>
#include <unistd.h>
#include "basics.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "integration.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"
#include "hp3D.hh"
#include "vectorial.hh"
#include "eigensolver.hh"
// ************************************************************* fem routine **
![\[ - \varepsilon i \omega \vec E + \text{curl} \vec H = \sigma \vec E + \vec J \]](form_153_dark.png)
![\[ \mu i \omega \vec H + \text{curl} \vec E = 0 \]](form_154_dark.png)
![\[ \text{div} \mu \vec H = 0 \qquad \text{and} \qquad \text{div} \varepsilon \vec E = 0. \]](form_155_dark.png)
![\[ \text{curl}( \mu^{-1} \text{curl} \vec E) - \omega^2 (\varepsilon + \frac{\sigma}{i\omega}) \vec E = -i \omega \vec J, \]](form_157_dark.png)
![\[ H(\text{curl}; D) := \{\vec u \in L^2(D)^3 : \text{curl} \vec u \in L^2(D)^3 \}. \]](form_158_dark.png)
![\[ \int_D (\mu^{-1} \text{curl} \vec E \cdot \text{curl} \vec v - \omega^2 \tilde\varepsilon \vec E \cdot \vec v) \, d\vec x = \int_D \vec f \cdot \vec v \, d\vec x \quad \forall \vec v \in H_0(\text{curl}; D), \]](form_161_dark.png)
![\[ \int_D \text{div} \vec E \text{div} \vec v \, d\vec x \]](form_166_dark.png)
![\[ \int_D (\mu^{-1} \text{curl} \vec E \cdot \text{curl} \vec v - \omega^2 \tilde\varepsilon \vec E \cdot \vec v) \, d\vec x + \langle\text{div}\vec E, \text{div}\vec v\rangle_Y \]](form_167_dark.png)
![\[ X_n[Y] := \{\vec u \in H_0(\text{curl}; D): \text{div}\vec u \in Y\} \]](form_169_dark.png)
![\[ w(\vec x) = \text{dist}\left( \vec x, \bigcup_{c \in C} c \cup \bigcup_{e \in E} \bar e \right), \]](form_172_dark.png)
![\[ \langle\text{div} \vec u, \text{div} \vec v\rangle_Y = \int_D w^2 \text{div} \vec u \text{div} \vec v \, d\vec x. \]](form_175_dark.png)
![$0 \neq \vec E \in X_n[Y]$](form_177_dark.png)
![\[ \int_D \mu^{-1} \text{curl}\vec E \cdot \text{curl} \vec v\, d\vec x + s\langle\text{div} \vec E, \text{div} \vec v\rangle_Y = \omega^2 \int_D \varepsilon \vec E \cdot \vec v \, d\vec x \quad \forall \vec v \in X_n[Y]. \]](form_178_dark.png)
![\[ \frac{\|\text{curl} \vec E\|^2_0}{s \|\text{div} \vec E\|_Y^2}
\begin{cases}
\geq \rho & \text{physical Eigenvalue} \\
\leq \rho^{-1} & \text{spurious Eigenvalue} \\
\text{otherwise} & \text{undecided}.
\end{cases} \]](form_181_dark.png)
With 18726 dofs, the Eigenvalues plotted versus s: 
The physical Eigenvalues (green) are lined up on horizontal lines whereas the spurious Eigenvalues are on straight lines through the origin. The exact values (horizontal blue lines) are taken from Monique Dauge's Benchmax page. Below are the convergence histories for the first three Eigenvalues for selected values of s: 
These convergence histories show exponential convergence: straight lines in a ![$ \log-\sqrt[4]{.} $](form_184_dark.png)
1.9.8