exactDtN.cc

Contents

1. Introduction
2. Variational Formulation
3. Dirichlet-to-Neumann Mapping
4. Commented Program
5. Results
6. References
7. Complete Source Code

@section intro Introduction
In this tutorial the implementation of the so called exact Dirichlet-to-Neumann map for scattering problems on a disc is shown. 

The equation solved is the 2D Helmholtz equation
\f[ 
- \nabla \cdot \alpha \nabla u - \beta u = 0
\f]
on the disc \f$ \Omega = \{\mathbf{x}=(x,y)\in\mathbb{R}^2 \mid r=\sqrt{x^2+y^2} < R\} \f$ of radius \f$ R \f$, where the coefficient functions \f$ \alpha(x) \f$ and \f$ \beta(x) \f$ depend on whether the TE or TM mode is considered. In the TE-case, i.e. \f$ u \f$ corresponds to the third component of the magnetic field, \f$ \alpha \f$ and \f$ \beta \f$ are given by
\f[
\alpha = \frac{1}{\epsilon}, \qquad \beta = \omega^2 \mu,
\f]
with the permittivity \f$ \epsilon \f$, the permeability \f$ \mu \f$ and the frequency \f$ \omega \f$. On the other hand, in the TM-case, i.e. \f$ u \f$ corresponds to the third component of the electric field, \f$ \alpha \f$ and \f$ \beta \f$ are given by
\f[
\alpha \equiv 1, \qquad \beta = \omega^2 \epsilon \mu.
\f]

In the exterior domain \f$ \Omega^{\rm ext} = \mathbb{R}^2\setminus\Omega \f$ the total field \f$ u^{\rm ext} \f$ can be split into a known incoming field \f$ u^{\rm inc} \f$ and an unknown scattered field \f$ u^{\rm sc} \f$, i.e.
\f[
u^{\rm ext} = u^{\rm inc} + u^{\rm sc}.
\f]
The coefficient functions \f$ \alpha \f$ and \f$ \beta \f$ are assumed to be constant in \f$ \Omega^{\rm ext} \f$. We will denote the constant values of the coefficient functions by \f$ \alpha_0 \f$ and \f$ \beta_0 \f$ and they are defined according to the definitions above \f$ \epsilon \equiv \epsilon_0 \f$ and \f$ \mu \equiv \mu_0 \f$ where \f$ \epsilon_0 \f$ denotes the the vacuum permittivity and \f$ \mu_0 \f$ denotes the vacuum permeability.

Finally we assume the total field \f$ u \f$ and its co-normal derivative \f$ \alpha \nabla u \cdot \mathbf{n} \f$ to be continuous over \f$ \partial\Omega \f$.

We will suppose for now that there exists an appropriate Dirichlet-to-Neumann (DtN) mapping \f$ B \f$ for the scattered field \f$ u^{sc} \f$ on the boundary \f$ \partial\Omega \f$, i.e.
\f[
\nabla u^{\rm sc} \cdot \mathbf{n} = B u^{\rm sc} 
\qquad \text{on }\partial\Omega,
\f]    
In the tutorial @ref BGT_0.cc "BGT_0.cc" this mapping was given by the BGT-0 boundary condition \f$ B = {\rm i} k_0 \f$. In the section <a href="#DtN">Dirichlet-to-Neumann Mapping</a> we will see how to construct such a linear DtN mapping that can be regarded as "exact".

@section variation Variational Formulation
Now we derive the corresponding variational formulation of the introduced problem:
find \f$ u\in H^1(\Omega) \f$ such that
\f[
\int_{\Omega}\alpha \nabla u \cdot \nabla v\;{\rm d}\mathbf{x}
- \int_{\Omega}\beta uv\;{\rm d}\mathbf{x} 
- \int_{\partial\Omega} \alpha \nabla u \cdot \mathbf{n} v \;{\rm d}s(\mathbf{x}) 
= 0 \qquad\forall v\in H^1(\Omega).
\f]
Considering its continuity over \f$ \partial\Omega \f$ the co-normal of the total field \f$ u \f$ can be written in the form 
\f[
\alpha\nabla u\cdot \mathbf{n}
= \alpha_0\nabla u^{\rm ext} \cdot \mathbf{n}
\qquad \text{on }\partial\Omega.
\f]
Using the fact that the total field in the exterior domain can be split into an incoming field and a scattered field, and incorporating the DtN mapping \f$ B \f$ gives
\f[
\alpha \nabla u \cdot \mathbf{n}
= \alpha_0 \left( B u^{\rm sc}
                 + \nabla u^{\rm inc} \cdot \mathbf{n}
  \right)
\qquad \text{on }\partial\Omega.
\f]
Now we use the linearity of the DtN mapping \f$ B \f$ and the continuity of the total field \f$ u \f$ over \f$ \partial\Omega \f$ to obtain
\f[
\alpha \nabla u \cdot \mathbf{n}
= \alpha_0 \left( B u^{\rm ext}
                 - B u^{\rm inc}
                 + \nabla u^{\rm inc} \cdot \mathbf{n}
           \right) 
= \alpha_0 \left( B u
                 - B u^{\rm inc}
                 + \nabla u^{\rm inc} \cdot \mathbf{n}
           \right) 
\qquad \text{on }\partial\Omega.
\f]
Incorporating this into the variational formulation yields
\f[
\int_{\Omega}\alpha \nabla u \cdot \nabla v\;{\rm d}\mathbf{x}
- \int_{\Omega}\beta uv\;{\rm d}\mathbf{x} 
- \alpha_0\int_{\partial\Omega} B u\,v \;{\rm d}s(\mathbf{x}) 
= \alpha_0\int_{\partial\Omega} \left(\nabla u^{\rm inc} \cdot \mathbf{n} 
  - B u^{\rm inc}\right) v \;{\rm d}s(\mathbf{x})
\qquad\forall v\in H^1(\Omega).
\f]

@section DtN Dirichlet-to-Neumann Mapping
In this section we want derive the exact DtN mapping for scattering problems on a disc. Since the total field \f$ u \f$ and the incoming field \f$ u^{\rm inc} \f$ are supposed to satisfy the Helmholtz equation in the exterior domain \f$ \Omega^{\rm ext} \f$, the scattered field \f$ u^{\rm sc} \f$ also satisfies the Helmholtz equation in \f$ \Omega^{\rm ext} \f$, i.e. we have
\f[
- \Delta u^{\rm sc} - k_0^2 u^{\rm sc} = 0
\qquad \text{in }\Omega^{\rm ext}
\f]
with the wave number \f$ k_0 = \omega\sqrt{\epsilon_0\mu_0} \f$ for both, the TE-case and the TM-case.

Rewriting this equation in polar coordinates \f$ (r,\phi) \f$ gives
\f[
-\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u^{\rm sc}}{\partial r}\right)
-\frac{1}{r^2}\frac{\partial^2 u^{\rm sc}}{\partial \phi^2}
-k_0^2 u^{\rm sc}
=0.
\f]

Using the approach 
\f$ u^{\rm sc}(r,\phi) = \sum_{n\in\mathbb{Z}}u^{\rm sc}_n(r)e^{{\rm i}n\phi} \f$, that is based on the Fourier expansion in \f$ \phi \f$, we get
\f[
\sum_{n\in\mathbb{Z}}\left[
-\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_n^{\rm sc}(r)}{\partial r}\right)e^{{\rm i}n\phi}
-\frac{1}{r^2}u_n^{\rm sc}(r)\frac{\partial^2 e^{{\rm i}n\phi}}{\partial \phi^2}
-k_0^2 u_n^{\rm sc}(r)e^{{\rm i}n\phi}
 \right] = 0,
\f]
i.e.
\f[
\sum_{n\in\mathbb{Z}}\left[
-\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_n^{\rm sc}(r)}{\partial r}\right)
+\frac{n^2}{r^2}u_n^{\rm sc}(r)
-k_0^2 u_n^{\rm sc}(r)
\right]e^{{\rm i}n\phi} = 0,
\f]
Since \f$ e^{{\rm i}n\phi} \f$, \f$ n\in\mathbb{Z} \f$, are mutually linear independent, we have
\f[
-\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_n^{\rm sc}(r)}{\partial r}\right)
+\frac{n^2}{r^2}u_n^{\rm sc}(r)
-k_0^2 u_n^{\rm sc}(r)
=0 \qquad \forall n\in\mathbb{Z}
\f]
and thus
\f[
r^2\frac{\partial^2}{\partial r^2}u_n^{\rm sc}(r)
+r\frac{\partial}{\partial r}u_n^{\rm sc}(r)
+\left(r^2 k_0^2 - n^2\right)u_n^{\rm sc}(r)
=0 \qquad \forall n\in\mathbb{Z}.
\f]
By substituting \f$ \rho = k_0 r \f$ we arrive at the Bessel differential equation
\f[
\rho^2\frac{\partial^2}{\partial \rho^2}u_n^{\rm sc}(\rho)
+\rho\frac{\partial}{\partial \rho}u_n^{\rm sc}(\rho)
+\left(\rho^2 - n^2\right)u_n^{\rm sc}(\rho)
=0 \qquad \forall n\in\mathbb{Z}.
\f]
The solutions to this equation define the Bessel functions \f$ J_n \f$ of the first kind and the Bessel functions \f$ Y_n \f$ of the second kind. 

As the Bessel differential equation is linear, any linear combination of the Bessel functions  \f$ J_n \f$ and \f$ Y_n \f$ solve the equation. In particular, we can choose the Hankel functions \f$ H^{(1)}_n \equiv J_n + {\rm i} Y_n \f$. In the limit \f$ r\rightarrow\infty \f$ the Hankel functions behave like \f$ \sqrt{\frac{2}{\pi k_0 r}}e^{{\rm i} k_0 r - n\pi/2-\pi/4} \f$ <a href="#references">[1]</a> and thus, they satisfy the Sommerfeld radiation condition
\f[
\lim_{r\rightarrow\infty}\sqrt{r}\left(\frac{\partial u^{\rm sc}}{\partial r}
- {\rm i}k_0 u^{\rm sc}\right)=0,
\f]
c.f. <a href="#references">[2]</a>.

Taking the Hankel functions \f$ H^{(1)}_n \f$ as basis for \f$ u_n^{\rm sc} \f$ we can write the scattered field in the form
\f[
u^{\rm sc}(r,\phi) = \sum_{n\in\mathbb{Z}} a_n H_n(k_0 r) e^{{\rm i} n\phi}
\f]
and its derivative
\f[
\frac{\partial}{\partial r} u^{\rm sc}(r,\phi) = \sum_{n\in\mathbb{Z}} a_n k_0 H_n'(k_0 r) e^{{\rm i} n\phi}.
\f]
In particular, we have
\f[
u^{\rm sc}(R,\phi) = \sum_{n\in\mathbb{Z}} a_n H_n(k_0 R) e^{{\rm i} n\phi}
\f]
and by multiplying with \f$ e^{-{\rm i}m\phi} \f$, \f$ m\in\mathbb{Z} \f$ and integrating over \f$ [0,2\pi] \f$ we arrive at
\f[
\int_0^{2\pi} u^{\rm sc}(R,\phi) e^{-{\rm i}m\phi} {\rm d}\phi
= 2\pi a_{m} H_{m}(k_0 R)
\f]
where we used the orthogonality of \f$ e^{{\rm i}n\phi} \f$, \f$ n\in\mathbb{Z} \f$. Thus, we obtain
\f[
a_{m} = \frac{1}{2\pi}\frac{1}{H_m(k_0 R)}\int_0^{2\pi}u(R,\phi)e^{-{\rm i}m\phi}{\rm d}\phi,
\f]
and \f$ \frac{\partial}{\partial r} u^{\rm sc}(R,\phi) \f$ can written as
\f[
\frac{\partial}{\partial r} u^{\rm sc}(R,\phi) 
=
\frac{k_0}{2\pi}\sum_{n\in\mathbb{Z}}\int_0^{2\pi}u(R,\psi)e^{-{\rm i}n\psi}{\rm d}\psi\;
\frac{H_n'(k_0 R)}{H_n(k_0 R)}e^{{\rm i}n\phi}. 
\f]

This is our sought DtN mapping which we can incorporate into the boundary integral
\f[
\int_{\partial\Omega} (Bu)v {\rm d}s(\mathbf{x})
= \frac{k_0}{2\pi R}\sum_{n\in\mathbb{Z}}\frac{H_n'(k_0 R)}{H_n(k_0 R)}
  \int_{\partial\Omega} v(\mathbf{x})e^{{\rm i}n\phi(\mathbf{x})} {\rm d}s(\mathbf{x})
  \int_{\partial\Omega} u(\mathbf{x})e^{-{\rm i}n\phi(\mathbf{x})} {\rm d}s(\mathbf{x}).
\f]


@dontinclude exactDtN.cc
@section commented Commented Program
First, system files
@until iostream
and concepts files
@skip basics
@until models/DtNmap2D
are included. With the following using directives
@skip using
@until using
@until using
we do not need to prepend <tt>concepts::</tt> to <tt>Real</tt> and <tt>Cmplx</tt> 
everytime.

<!--
@subsection class Class integrandF
Before starting with the main program we define a class <tt>integrandF</tt> derived
from <tt>concepts::Formula</tt> in order to evaluate \f$ e^{{\rm i}n\phi(\mathbf{x})} \f$
at \f$ \mathbf{x}\in\partial\Omega \f$ for different values of \f$ n\in\mathbb{Z} \f$. Since
\f$ |n| \f$ is usually very small (in the implementation below we have
\f$ -5 \leq n\leq 5 \f$) we evaluate the value of \f$ e^{{\rm i}n\phi(\mathbf{x})} \f$
by simply computing
\f$
e^{{\rm i}n\phi(\mathbf{x})} 
= \left(\cos(\phi(\mathbf{x}))+{\rm i}\sin(\phi(\mathbf{x}))\right)^n
= \left(x/r+{\rm i}\,y/r\right)^n.
\f$
@skip class
@until private:
@until os
@until }
-->

@subsection main Main Program
We start the main program
@skip int
@until try
by setting the values of the parameters
@skip Real
@until DtNmaxN
and defining the attributes of the boundary, the area of the scatterer and 
the surrounding area respectively.
@skip uint
@until attr0

We specify whether the TE-case or the TM-case is computed.
@skip enum
@until mode

Then we declare the coefficients <tt>alpha0</tt>, <tt>alphaSc</tt>, 
<tt>beta0</tt> and <tt>betaSc</tt> as complex numbers
@skip Cmplx
@until betaSc
and assign their values according to the chosen mode, i.e. in the TE-mode 
we have
@skip if
@until }
and in the TM-mode we write
@skip if
@until }

With these values we define \f$ \alpha \f$ and \f$ \beta \f$ as     
<tt>PiecewiseConstFormula</tt>
@skip concepts
@until beta0

Now we compute the wave number in the exterior domain
@skip const
@until Real
and introduce the incoming field \f$ u^{\rm inc} = e^{{\rm i} k_0 x} \f$ 
and its normal gradient \f$ \nabla u^{\rm inc} \cdot \mathbf{n}\f$ as 
<tt>ParsedFormula</tt>.
@skip std::stringstream
@until concepts::ParsedFormula
@until concepts::ParsedFormula

The mesh is created as a circle surrounded by a ring.
@skip Real
@until endl
Now the mesh is plotted using a scaling factor of 100, a greyscale of 1.0 
and 20 points per edge.
@until drawMeshEPS

The space is built by refining the mesh once and setting the polynomial degree
to 10. Then the elements of the space are built and the space is plotted.
@skip hp2D
@until drawMeshEPS

Now the trace space of the boundary \f$ \partial\Omega \f$ can be built.
@skip hp2D
@until endl

The right hand side
@skip hp1D
@until concepts::Vector
and the system matrix
@skip concepts::SparseMatrix
@until M.addInto
without DtN part are computed.

Now we add the the DtN sum to the system matrix and the right hand side vector.
To this end, we first compute the coefficients 
\f[ -\frac{k_0}{2\pi R} \frac{H_n'(k_0 R)}{H_n(k_0 R)} \f]
using <tt> getHelmholtzDtNCoeff_Circle2D </tt>
@skip concepts::Sequence
@until concepts::getHelmholtzDtNCoeff_Circle2D
where we use
\f[ H_n'(k_0 R) =  H_{n-1}(k_0 R) - \frac{n}{k_0 R}H_n(k_0 R), \f]
c.f. Eq. 9.1.27 in <a href="#references">[1]</a>. Then we add the term 
\f$- \alpha_0\int_{\partial\Omega} B u\,v \;{\rm d}s(\mathbf{x}) \f$ 
to the system matrix <tt> S </tt> and the term
\f$- \alpha_0\int_{\partial\Omega} B u^{\rm inc}\,v \;{\rm d}s(\mathbf{x}) \f$
to the right hand side vector <tt> rhs </tt> using <tt> addExactDtN_Circle2D </tt>
@until concepts::addExactDtN_Circle2D

The system matrix including the DtN part is displayed
@skip std::cout
@until endl
and the system is solved using the direct solver SuperLU.
@skip concepts::SuperLU
@until solver
@until concepts::Vector
@until solver

In order to plot the solution the shape functions are computed on 
equidistant points using the trapezoidal quadrature rule.
@skip hp2D
@until graphics::MatlabBinaryGraphics

Finally, exceptions are caught and a sane return value is given back.
@skip }
@until return
@until return
@until }

@section results Results
Output of the program:
@code 
Mesh: Circle(ncell = 13, cells: Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(0), (Vertex(Key(0)), Vertex(Key(1)), Vertex(Key(2)), Vertex(Key(3))), Attribute(21)), vtx = [<2>(0.5, 0), <2>(0, 0.5), <2>(-0.5, 0), <2>(0, -0.5)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(1), (Vertex(Key(1)), Vertex(Key(0)), Vertex(Key(4)), Vertex(Key(5))), Attribute(21)), vtx = [<2>(0, 0.5), <2>(0.5, 0), <2>(1, 0), <2>(0, 1)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(2), (Vertex(Key(2)), Vertex(Key(1)), Vertex(Key(5)), Vertex(Key(6))), Attribute(21)), vtx = [<2>(-0.5, 0), <2>(0, 0.5), <2>(0, 1), <2>(-1, 0)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(3), (Vertex(Key(3)), Vertex(Key(2)), Vertex(Key(6)), Vertex(Key(7))), Attribute(21)), vtx = [<2>(0, -0.5), <2>(-0.5, 0), <2>(-1, 0), <2>(0, -1)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(4), (Vertex(Key(0)), Vertex(Key(3)), Vertex(Key(7)), Vertex(Key(4))), Attribute(21)), vtx = [<2>(0.5, 0), <2>(0, -0.5), <2>(0, -1), <2>(1, 0)]),

Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(5), (Vertex(Key(5)), Vertex(Key(4)), Vertex(Key(8)), Vertex(Key(9))), Attribute(22)), vtx = [<2>(0, 1), <2>(1, 0), <2>(8, 0), <2>(0, 8)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(6), (Vertex(Key(6)), Vertex(Key(5)), Vertex(Key(9)), Vertex(Key(10))), Attribute(22)), vtx = [<2>(-1, 0), <2>(0, 1), <2>(0, 8), <2>(-8, 0)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(7), (Vertex(Key(7)), Vertex(Key(6)), Vertex(Key(10)), Vertex(Key(11))), Attribute(22)), vtx = [<2>(0, -1), <2>(-1, 0), <2>(-8, 0), <2>(0, -8)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(8), (Vertex(Key(4)), Vertex(Key(7)), Vertex(Key(11)), Vertex(Key(8))), Attribute(22)), vtx = [<2>(1, 0), <2>(0, -1), <2>(0, -8), <2>(8, 0)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(9), (Vertex(Key(9)), Vertex(Key(8)), Vertex(Key(12)), Vertex(Key(13))), Attribute(22)), vtx = [<2>(0, 8), <2>(8, 0), <2>(15, 0), <2>(0, 15)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key( 10), (Vertex(Key(10)), Vertex(Key(9)), Vertex(Key(13)), Vertex(Key(14))), Attribute(22)), vtx = [<2>(-8, 0), <2>(0, 8), <2>(0, 15), <2>(-15, 0)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(11), (Vertex(Key(11)), Vertex(Key(10)), Vertex(Key(14)), Vertex(Key(15))), Attribute(22)), vtx = [<2>(0, -8), <2>(-8, 0), <2>(-15, 0), <2>(0, -15)]), Quad2d(idx = Level<2>(0, 0); 0, 0), cntr = Quad(Key(12), (Vertex(Key(8)), Vertex(Key(11)), Vertex(Key(15)), Vertex(Key(12))), Attribute(22)), vtx = [<2>(8, 0), <2>(0, -8), <2>(0, -15), <2>(15, 0)]))

Space: hpAdaptiveSpaceH1(dim = 5241 (V:57, E:972, I:4212), nelm = 52)

Trace Space: TraceSpace(QuadEdgeFirst(), dim = 5241, nelm = 8)

System Matrix: SparseMatrix(5241x5241, HashedSparseMatrix: 754721 (2.74763%) entries bound.)

Solver: SuperLU(n = 5241)
@endcode

Plot of the refined mesh:
<img src="exactDtN_mesh.png" alt="Mesh output">

Matlab plots of the total field \f$ u \f$:
<img src="exactDtN.png" alt="Matlab output">

@section references References
[1] M. Abramowitz, and I.A. Stegun (ed.), "Handbook of mathematical functions with formulas, graphs, and mathematical tables", National Bureau of Standards, Applied Mathematics Series 55, tenth printing, 1972.

[2] A. Sommerfeld, "Partial Differential Equations in Physics - Lectures on Theoretical Physics Volume VI", Academic Press, New York, 1964.

@section complete Complete Source Code
@author Dirk Klindworth, 2011/2012

#include <iostream>
 
#include "basics.hh"
#include "toolbox.hh"
#include "geometry.hh"
#include "space.hh"
#include "formula.hh"
#include "function.hh"
#include "operator.hh"
#include "graphics.hh"
#include "hp1D.hh"
#include "hp2D.hh"
#include "models/DtNmap2D.hh"
 
using concepts::Real;
using concepts::Cmplx;
 
 
// ******************************************************************** main **
 
int main(int argc, char** argv) {
  try {
 
#ifdef HAS_MPI
    // Initialisation is necessary if parallel version of Mumps is compiled
    MPI_Init(&argc, &argv);
#endif
    
    // ** parameters **
    const Real  R       = 15.0;       // outer radius
    const Real  innerR  =  1.0;       // inner radius (i.e. radius of scatterer)
    const Real  omega   =  1.0;       // frequency
    const Real  eps0    =  1.0;       // permittivity outside the scatterer
    const Real  mu0     =  1.0;       // permeability outside the scatterer
    const Cmplx epsSc     (4.0, 0.0); // permittivity inside the scatterer
    const Real  muSc    =  1.0;       // permeability inside the scatterer
    const uint  DtNmaxN =  5;         // number DtN summands
    
    // ** attributes of the computational domain **
    const uint attrBoundary = 11;     // boundary
    const uint attrSc       = 21;     // scatterer
    const uint attr0        = 22;     // computational domain without scatterer
    
    // ** set mode **
    enum TETM {TE, TM};
    TETM mode = TM;
    
    // ** declare coefficients alpha and beta **
    Cmplx alpha0;
    Cmplx alphaSc;
    Cmplx beta0;
    Cmplx betaSc;
    
    // ** write coefficients alpha and beta (TE-case) **
    if (mode == TE) {
      alphaSc = Cmplx(pow(epsSc,-1.0));
      alpha0  = Cmplx(pow(eps0,-1.0),0.0);
      betaSc  = Cmplx(omega*omega*muSc,0.0);
      beta0   = Cmplx(omega*omega*mu0,0.0);
    }
    
    // ** write coefficients alpha and beta (TM-case) **
    if (mode == TM) {
      alphaSc = Cmplx(1.0, 0.0);
      alpha0  = Cmplx(1.0, 0.0);
      betaSc  = Cmplx(omega*omega*epsSc*muSc);
      beta0   = Cmplx(omega*omega*eps0*mu0,0.0);
    }
    
    // ** piecewise constant formulas for alpha and beta **
    concepts::PiecewiseConstFormula<Cmplx> alpha;
    alpha[concepts::Attribute(attrSc)] = alphaSc;
    alpha[concepts::Attribute(attr0)]  = alpha0;
    concepts::PiecewiseConstFormula<Cmplx> beta;
    beta[concepts::Attribute(attrSc)]  = betaSc;
    beta[concepts::Attribute(attr0)]   = beta0;
    
    // ** wave vector in the exterior **
    const Real k0 = omega*sqrt(mu0*eps0);
    
    // ** incoming wave and its normal gradient **
    std::stringstream uIncReal, uIncImag, uIncGradReal, uIncGradImag;
    uIncReal     << "(cos((" << k0 << ")*x))";
    uIncImag     << "(sin((" << k0 << ")*x))";
    uIncGradReal << "(-(x/(sqrt(x*x+y*y)))*" << k0 << "*sin((" << k0 << ")*x))";
    uIncGradImag << "( (x/(sqrt(x*x+y*y)))*" << k0 << "*cos((" << k0 << ")*x))";
    concepts::ParsedFormula<Cmplx>     uInc(    uIncReal.str(),     uIncImag.str());
    concepts::ParsedFormula<Cmplx> uIncGrad(uIncGradReal.str(), uIncGradImag.str());
    
    // ** create mesh **
    Real ringData[3]     = {innerR,(innerR+R)/2.0,R};
    uint quadAttrData[3] = {attrSc,attr0,attr0};
    uint edgeAttrData[3] = {0,0,attrBoundary};
    concepts::Array<Real> ring(ringData,3);
    concepts::Array<uint> ringQuadAttr(quadAttrData,3);
    concepts::Array<uint> ringEdgeAttr(edgeAttrData,3);
    concepts::Circle msh(ring,ringQuadAttr,ringEdgeAttr);
    std::cout << std::endl << "Mesh: " << msh << std::endl;
    graphics::drawMeshEPS(msh,"exactDtN_mesh.eps",100,1.0,20);
    
    // ** space **
    hp2D::hpAdaptiveSpaceH1 spc(msh, 1, 10);
    spc.rebuild();
    std::cout << std::endl << "Space: " << spc << std::endl;
    graphics::drawMeshEPS(spc,"exactDtN_space.eps",100,1.0,20); 
    
    // ** build trace space at boundary **
    hp2D::TraceSpace tspc(spc,concepts::makeSet<uint>({attrBoundary}),hp2D::TraceSpace::FIRST);
    std::cout << std::endl << "Trace Space: " << tspc << std::endl;
    
    // ** right hand side without DtN part **
    hp1D::Riesz<Cmplx> lform(alpha0*uIncGrad);
    concepts::Vector<Cmplx> rhs(tspc, lform);
        
    // ** system matrix without DtN part **
    concepts::SparseMatrix<Cmplx> S(spc.dim(), spc.dim());
    hp2D::Laplace<Cmplx> la(alpha);
    concepts::SparseMatrix<Cmplx> A(spc, la);
    A.addInto(S, 1.0);
    hp2D::Identity<Cmplx> id(beta);
    concepts::SparseMatrix<Cmplx> M(spc, id);
    M.addInto(S, -1.0);
    
    // ** add DtN parts to system matrix and rhs **
    /*
     * implementation of the DtN map
     * 
     *  \int (Bu)v ds = -k0/(2*\pi*R) * \sum_n H'_n(k0*R)/H_n(k0*R) *
     *                  \int v exp(i*n*\phi(x)) ds * \int u exp(-i*n*\phi(x)) ds
     *
     */ 
    concepts::Sequence<Cmplx> DtNCoeff = 
      concepts::getHelmholtzDtNCoeff_Circle2D(k0,R,DtNmaxN);
    concepts::addExactDtN_Circle2D(S, tspc, DtNCoeff*alpha0, rhs, uInc, DtNCoeff*alpha0);
    
    // ** print system matrix **
    std::cout << std::endl << "System Matrix: " << S << std::endl;
            
    // ** solve **
#ifdef HAS_MUMPS
    concepts::Mumps<Cmplx> solver(S);
#else
    concepts::SuperLU<Cmplx> solver(S);
#endif
 
    std::cout << std::endl << "Solver: " << solver << std::endl;
    concepts::Vector<Cmplx> sol(spc);
    solver(rhs, sol);
    
    // ** plot solution **  
    hp2D::IntegrableQuad::factory().get()->setTensor(concepts::TRAPEZE, true, 10);
    spc.recomputeShapefunctions();
    graphics::MatlabBinaryGraphics(spc, "exactDtN", sol);
    
  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }
 
#ifdef HAS_MPI
  MPI_Finalize();
#endif
 
  return 0;
}