Home | Doxygen Documentation | Tutorials | Developer Tools (restricted)
exactDtN.cc
  1. Introduction
  2. Variational Formulation
  3. Dirichlet-to-Neumann Mapping
  4. Commented Program
  5. Results
  6. References
  7. Complete Source Code

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

\[ - \nabla \cdot \alpha \nabla u - \beta u = 0 \]

on the disc $ \Omega = \{\mathbf{x}=(x,y)\in\mathbb{R}^2 \mid r=\sqrt{x^2+y^2} < R\} $ of radius $ R $, where the coefficient functions $ \alpha(x) $ and $ \beta(x) $ depend on whether the TE or TM mode is considered. In the TE-case, i.e. $ u $ corresponds to the third component of the magnetic field, $ \alpha $ and $ \beta $ are given by

\[ \alpha = \frac{1}{\epsilon}, \qquad \beta = \omega^2 \mu, \]

with the permittivity $ \epsilon $, the permeability $ \mu $ and the frequency $ \omega $. On the other hand, in the TM-case, i.e. $ u $ corresponds to the third component of the electric field, $ \alpha $ and $ \beta $ are given by

\[ \alpha \equiv 1, \qquad \beta = \omega^2 \epsilon \mu. \]

In the exterior domain $ \Omega^{\rm ext} = \mathbb{R}^2\setminus\Omega $ the total field $ u^{\rm ext} $ can be split into a known incoming field $ u^{\rm inc} $ and an unknown scattered field $ u^{\rm sc} $, i.e.

\[ u^{\rm ext} = u^{\rm inc} + u^{\rm sc}. \]

The coefficient functions $ \alpha $ and $ \beta $ are assumed to be constant in $ \Omega^{\rm ext} $. We will denote the constant values of the coefficient functions by $ \alpha_0 $ and $ \beta_0 $ and they are defined according to the definitions above $ \epsilon \equiv \epsilon_0 $ and $ \mu \equiv \mu_0 $ where $ \epsilon_0 $ denotes the the vacuum permittivity and $ \mu_0 $ denotes the vacuum permeability.

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

We will suppose for now that there exists an appropriate Dirichlet-to-Neumann (DtN) mapping $ B $ for the scattered field $ u^{sc} $ on the boundary $ \partial\Omega $, i.e.

\[ \nabla u^{\rm sc} \cdot \mathbf{n} = B u^{\rm sc} \qquad \text{on }\partial\Omega, \]

In the tutorial BGT_0.cc this mapping was given by the BGT-0 boundary condition $ B = {\rm i} k_0 $. In the section Dirichlet-to-Neumann Mapping we will see how to construct such a linear DtN mapping that can be regarded as "exact".

Variational Formulation

Now we derive the corresponding variational formulation of the introduced problem: find $ u\in H^1(\Omega) $ such that

\[ \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). \]

Considering its continuity over $ \partial\Omega $ the co-normal of the total field $ u $ can be written in the form

\[ \alpha\nabla u\cdot \mathbf{n} = \alpha_0\nabla u^{\rm ext} \cdot \mathbf{n} \qquad \text{on }\partial\Omega. \]

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 $ B $ gives

\[ \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. \]

Now we use the linearity of the DtN mapping $ B $ and the continuity of the total field $ u $ over $ \partial\Omega $ to obtain

\[ \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. \]

Incorporating this into the variational formulation yields

\[ \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 \;{\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). \]

Dirichlet-to-Neumann Mapping

In this section we want derive the exact DtN mapping for scattering problems on a disc. Since the total field $ u $ and the incoming field $ u^{\rm inc} $ are supposed to satisfy the Helmholtz equation in the exterior domain $ \Omega^{\rm ext} $, the scattered field $ u^{\rm sc} $ also satisfies the Helmholtz equation in $ \Omega^{\rm ext} $, i.e. we have

\[ - \Delta u^{\rm sc} - k_0^2 u^{\rm sc} = 0 \qquad \text{in }\Omega^{\rm ext} \]

with the wave number $ k_0 = \omega\sqrt{\epsilon_0\mu_0} $ for both, the TE-case and the TM-case.

Rewriting this equation in polar coordinates $ (r,\phi) $ gives

\[ -\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. \]

Using the approach $ u^{\rm sc}(r,\phi) = \sum_{n\in\mathbb{Z}}u^{\rm sc}_n(r)e^{{\rm i}n\phi} $, that is based on the Fourier expansion in $ \phi $, we get

\[ \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, \]

i.e.

\[ \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, \]

Since $ e^{{\rm i}n\phi} $, $ n\in\mathbb{Z} $, are mutually linear independent, we have

\[ -\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} \]

and thus

\[ 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}. \]

By substituting $ \rho = k_0 r $ we arrive at the Bessel differential equation

\[ \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}. \]

The solutions to this equation define the Bessel functions $ J_n $ of the first kind and the Bessel functions $ Y_n $ of the second kind.

As the Bessel differential equation is linear, any linear combination of the Bessel functions $ J_n $ and $ Y_n $ solve the equation. In particular, we can choose the Hankel functions $ H^{(1)}_n \equiv J_n + {\rm i} Y_n $. In the limit $ r\rightarrow\infty $ the Hankel functions behave like $ \sqrt{\frac{2}{\pi k_0 r}}e^{{\rm i} k_0 r - n\pi/2-\pi/4} $ [1] and thus, they satisfy the Sommerfeld radiation condition

\[ \lim_{r\rightarrow\infty}\sqrt{r}\left(\frac{\partial u^{\rm sc}}{\partial r} - {\rm i}k_0 u^{\rm sc}\right)=0, \]

c.f. [2].

Taking the Hankel functions $ H^{(1)}_n $ as basis for $ u_n^{\rm sc} $ we can write the scattered field in the form

\[ u^{\rm sc}(r,\phi) = \sum_{n\in\mathbb{Z}} a_n H_n(k_0 r) e^{{\rm i} n\phi} \]

and its derivative

\[ \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}. \]

In particular, we have

\[ u^{\rm sc}(R,\phi) = \sum_{n\in\mathbb{Z}} a_n H_n(k_0 R) e^{{\rm i} n\phi} \]

and by multiplying with $ e^{-{\rm i}m\phi} $, $ m\in\mathbb{Z} $ and integrating over $ [0,2\pi] $ we arrive at

\[ \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) \]

where we used the orthogonality of $ e^{{\rm i}n\phi} $, $ n\in\mathbb{Z} $. Thus, we obtain

\[ 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, \]

and $ \frac{\partial}{\partial r} u^{\rm sc}(R,\phi) $ can written as

\[ \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}. \]

This is our sought DtN mapping which we can incorporate into the boundary integral

\[ \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}). \]

Commented Program

First, system files 

#include <iostream>
and concepts files 
#include "basics.hh"
#include "formula.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "hp1D.hh"
#include "hp2D.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"
are included. With the following using directives 
using concepts::Real;
using concepts::Cmplx;
we do not need to prepend concepts:: to Real and Cmplx everytime.

Class integrandF

Before starting with the main program we define a class integrandF derived from concepts::Formula in order to evaluate $ e^{{\rm i}n\phi(\mathbf{x})} $ at $ \mathbf{x}\in\partial\Omega $ for different values of $ n\in\mathbb{Z} $. Since $ |n| $ is usually very small (in the implementation below we have $ -5 \leq n\leq 5 $) we evaluate the value of $ e^{{\rm i}n\phi(\mathbf{x})} $ by simply computing $ 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. $ 

class integrandF : public concepts::Formula<Cmplx> {
public:
  // ** Constructor **
  integrandF(const int n = 0) 
    : n_(n) {}
  // ** () opterator **
  virtual Cmplx operator() (const Real p,    const Real t = 0.0) const;
  virtual Cmplx operator() (const concepts::Real2d& p, const Real t = 0.0) const;
  virtual Cmplx operator() (const concepts::Real3d& p, const Real t = 0.0) const;
  // ** clone **
  virtual integrandF* clone() const {
    return new integrandF(n_);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const;
private:
  // ** exponent of the function **
  const int n_;
};
  
Cmplx integrandF::operator() (const Real p, const Real t) const {
  std::stringstream error;
  error << "integrandF::operator() not defined for input parameter of type Real";
  throw concepts::MissingFeature(error.str());
  return 0.0;
}
  
Cmplx integrandF::operator() (const concepts::Real2d& p, const Real t) const {
  Real x = p[0];
  Real y = p[1];
  Real r = sqrt(x*x+y*y);
  Cmplx c(x/r,y/r);
  return pow(c,n_);
}
  
Cmplx integrandF::operator() (const concepts::Real3d& p, const Real t) const {
  const concepts::Real2d p2(p[0], p[1]);
  return (*this)(p2,t);
}

std::ostream& integrandF::info(std::ostream& os) const {
  return os << "integrandF(" << n_ << ")";
}

Main Program

We start the main program 

int main(int argc, char** argv) {
  try {
by setting the values of the 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 int   DtNmaxN =  5;         // number DtN summands
and defining the attributes of the boundary, the area of the scatterer and the surrounding area respectively. 
    const uint attrBoundary = 11;     // boundary
    const uint attrSc       = 21;     // scatterer
    const uint attr0        = 22;     // computational domain without scatterer

We specify whether the TE-case or the TM-case is computed. 

    enum TETM {TE, TM};
    TETM mode = TM;

Then we declare $ \alpha $ and $ \beta $ as PiecewiseConstFormula and $ \alpha_0 $ as ConstFormula. 

    concepts::PiecewiseConstFormula<Cmplx> alpha;
    concepts::PiecewiseConstFormula<Cmplx> beta;
    concepts::ConstFormula<Cmplx>          alpha0;

We write the values of $ \alpha $, $ \beta $ and $ \alpha_0 $ according to the defined mode. If the TE-case was chosen we write 

    if (mode == TE) {
      alpha[concepts::Attribute(attrSc)] = Cmplx(pow(epsSc,-1.0));
      alpha[concepts::Attribute(attr0)]  = Cmplx(pow(eps0,-1.0),0.0);
      beta[concepts::Attribute(attrSc)] = Cmplx(omega*omega*muSc,0.0);
      beta[concepts::Attribute(attr0)]  = Cmplx(omega*omega*mu0,0.0);
      alpha0 = Cmplx(pow(eps0,-1.0),0.0);
    }
and in the TM-case we write 
    if (mode == TM) {
      alpha[concepts::Attribute(attrSc)] = Cmplx(1.0, 0.0);
      alpha[concepts::Attribute(attr0)]  = Cmplx(1.0, 0.0);
      beta[concepts::Attribute(attrSc)] = Cmplx(omega*omega*epsSc*muSc);
      beta[concepts::Attribute(attr0)]  = Cmplx(omega*omega*eps0*mu0,0.0);
      alpha0 = Cmplx(1.0, 0.0);
    }

Now we compute the wave number 

    const Real k0 = omega*sqrt(mu0*eps0);
and introduce the incoming field $ u^{\rm inc} = e^{{\rm i} k_0 x} $ and its normal gradient $ \nabla u^{\rm inc} \cdot \mathbf{n}$ as ParsedFormula. 
    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());

The mesh is created as a circle surrounded by a ring. 

    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;
Now the mesh is plotted using a scaling factor of 100, a greyscale of 1.0 and 20 points per edge. 
    graphics::drawMeshEPS(msh,"exactDTN_mesh.eps",100,1.0,20);

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. 

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

Now the trace space of the boundary $ \partial\Omega $ can be built. 

    hp2D::TraceSpace tspc(spc,concepts::makeSet<uint>(1,attrBoundary),hp2D::TraceSpace::FIRST);
    std::cout << std::endl << "Trace Space: " << tspc << std::endl;

The right hand side 

    hp1D::Riesz<Cmplx> lform(alpha0*uIncGrad);
    concepts::Vector<Cmplx> rhs(tspc, lform);
and the system matrix 
    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);
without DtN part are computed.

Now we perform the DtN sum. First we define the DtN factor $ \frac{\omega}{2\pi R} $ 

    const Real DtNfactor = omega/(2*M_PI*R);
and declare the complex Bessel coefficient. 
    Cmplx besselCoeff;
Then we sum from -DtNmaxN to DtNmaxN. 
    for (int n = -abs(DtNmaxN); n <= abs(DtNmaxN); ++n) {
      // v-integral
      integrandF exp_p_inphi(n);      
      hp1D::Riesz<Cmplx> vLform(alpha0*exp_p_inphi);
      concepts::Vector<Cmplx> vVector(tspc, vLform);
      // u-integral
      integrandF exp_m_inphi(-n);
      hp1D::Riesz<Cmplx> uLform(exp_m_inphi);
      concepts::Vector<Cmplx> uVector(tspc, uLform);
      // rank-1-product of v- and u-integral
      concepts::SparseMatrix<Cmplx> DtNmatrix(tspc,vVector,uVector);
      // ratio of derivative of Hankel function and Hankel function
      besselCoeff = Cmplx(concepts::besselJn(omega*R,n-1),concepts::besselYn(omega*R,n-1))
                    / Cmplx(concepts::besselJn(omega*R,n),concepts::besselYn(omega*R,n))   
                    - n/(omega*R);
      DtNmatrix.addInto(S, -DtNfactor*besselCoeff);
      // integral of uInc * exp_m_inphi
      Cmplx integral = concepts::integrate(tspc,uInc*exp_m_inphi);
      rhs = rhs - vVector*(DtNfactor*besselCoeff*integral);
    }

The system matrix including the DtN part is displayed 

    std::cout << std::endl << "System Matrix: " << S << std::endl;
and the system is solved using the direct solver SuperLU. 
    concepts::SuperLU<Cmplx> solver(S);
    std::cout << std::endl << "Solver: " << solver << std::endl;
    concepts::Vector<Cmplx> sol(spc);
    solver(rhs, sol);

In order to plot the solution the shape functions are computed on equidistant points using the trapezoidal quadrature rule. 

    hp2D::IntegrableQuad::rule().set(concepts::TRAPEZE, true, 10);
    spc.recomputeShapefunctions();
    graphics::MatlabGraphics(spc, "exactDtN", sol);

Finally, exceptions are caught and a sane return value is given back. 

  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }
  return 0;
}

Results

Output of the program: 

    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)

Plot of the refined mesh:

exactDtN_mesh.png

Matlab plots of the total field $ u $:

exactDtN.png

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.

Complete Source Code

Author:
Dirk Klindworth, 2011
#include <iostream>

#include "basics.hh"
#include "formula.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "hp1D.hh"
#include "hp2D.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"

using concepts::Real;
using concepts::Cmplx;

// ************************************************************** integrandF **
/*
 * Class for evaluating the integrand function of the exact DtN map
   
   The order \c n has to be given to evaluate
     \f[(x+iy)^n\f],
   where \f$x\f$ and \f$y\f$ are cartesian coordinates.
 */
class integrandF : public concepts::Formula<Cmplx> {
public:
  // ** Constructor **
  integrandF(const int n = 0) 
    : n_(n) {}
  // ** () opterator **
  virtual Cmplx operator() (const Real p,    const Real t = 0.0) const;
  virtual Cmplx operator() (const concepts::Real2d& p, const Real t = 0.0) const;
  virtual Cmplx operator() (const concepts::Real3d& p, const Real t = 0.0) const;
  // ** clone **
  virtual integrandF* clone() const {
    return new integrandF(n_);
  }
protected:
  virtual std::ostream& info(std::ostream& os) const;
private:
  // ** exponent of the function **
  const int n_;
};
  
Cmplx integrandF::operator() (const Real p, const Real t) const {
  std::stringstream error;
  error << "integrandF::operator() not defined for input parameter of type Real";
  throw concepts::MissingFeature(error.str());
  return 0.0;
}
  
Cmplx integrandF::operator() (const concepts::Real2d& p, const Real t) const {
  Real x = p[0];
  Real y = p[1];
  Real r = sqrt(x*x+y*y);
  Cmplx c(x/r,y/r);
  return pow(c,n_);
}
  
Cmplx integrandF::operator() (const concepts::Real3d& p, const Real t) const {
  const concepts::Real2d p2(p[0], p[1]);
  return (*this)(p2,t);
}

std::ostream& integrandF::info(std::ostream& os) const {
  return os << "integrandF(" << n_ << ")";
}

// ******************************************************************** main **

int main(int argc, char** argv) {
  try {
    
    // ** 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 int   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 piecewise constant formulas for alpha and beta **
    concepts::PiecewiseConstFormula<Cmplx> alpha;
    concepts::PiecewiseConstFormula<Cmplx> beta;
    concepts::ConstFormula<Cmplx>          alpha0;
    
    // ** Piecewise constant alpha and beta (TE-case) **
    if (mode == TE) {
      alpha[concepts::Attribute(attrSc)] = Cmplx(pow(epsSc,-1.0));
      alpha[concepts::Attribute(attr0)]  = Cmplx(pow(eps0,-1.0),0.0);
      beta[concepts::Attribute(attrSc)] = Cmplx(omega*omega*muSc,0.0);
      beta[concepts::Attribute(attr0)]  = Cmplx(omega*omega*mu0,0.0);
      alpha0 = Cmplx(pow(eps0,-1.0),0.0);
    }
    
    // ** Piecewise constant alpha and beta (TM-case) **
    if (mode == TM) {
      alpha[concepts::Attribute(attrSc)] = Cmplx(1.0, 0.0);
      alpha[concepts::Attribute(attr0)]  = Cmplx(1.0, 0.0);
      beta[concepts::Attribute(attrSc)] = Cmplx(omega*omega*epsSc*muSc);
      beta[concepts::Attribute(attr0)]  = Cmplx(omega*omega*eps0*mu0,0.0);
      alpha0 = Cmplx(1.0, 0.0);
    }
    
    // ** 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 of boundary **
    hp2D::TraceSpace tspc(spc,concepts::makeSet<uint>(1,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);
    
    // ** DtN part of system matrix and rhs
    /*
     * implementation of the DtN map
     * 
     *  \int (Bu)v ds = \omega/(2*\pi*R) * \sum H'_n(\omega*R)/H_n(\omega*R) *
     *                  \int v exp(i*n*\phi(x)) ds * \int u exp(-i*n*\phi(x)) ds
     *
     */ 
    const Real DtNfactor = omega/(2*M_PI*R);
    Cmplx besselCoeff;
    for (int n = -abs(DtNmaxN); n <= abs(DtNmaxN); ++n) {
      // v-integral
      integrandF exp_p_inphi(n);      
      hp1D::Riesz<Cmplx> vLform(alpha0*exp_p_inphi);
      concepts::Vector<Cmplx> vVector(tspc, vLform);
      // u-integral
      integrandF exp_m_inphi(-n);
      hp1D::Riesz<Cmplx> uLform(exp_m_inphi);
      concepts::Vector<Cmplx> uVector(tspc, uLform);
      // rank-1-product of v- and u-integral
      concepts::SparseMatrix<Cmplx> DtNmatrix(tspc,vVector,uVector);
      // ratio of derivative of Hankel function and Hankel function
      besselCoeff = Cmplx(concepts::besselJn(omega*R,n-1),concepts::besselYn(omega*R,n-1))
                    / Cmplx(concepts::besselJn(omega*R,n),concepts::besselYn(omega*R,n))   
                    - n/(omega*R);
      DtNmatrix.addInto(S, -DtNfactor*besselCoeff);
      // integral of uInc * exp_m_inphi
      Cmplx integral = concepts::integrate(tspc,uInc*exp_m_inphi);
      rhs = rhs - vVector*(DtNfactor*besselCoeff*integral);
    }
    
    // ** print system matrix and rhs vector **
    std::cout << std::endl << "System Matrix: " << S << std::endl;
        
    // ** solve **
    concepts::SuperLU<Cmplx> solver(S);
    std::cout << std::endl << "Solver: " << solver << std::endl;
    concepts::Vector<Cmplx> sol(spc);
    solver(rhs, sol);
    
    // ** plot solution **  
    hp2D::IntegrableQuad::rule().set(concepts::TRAPEZE, true, 10);
    spc.recomputeShapefunctions();
    graphics::MatlabGraphics(spc, "exactDtN", sol);
    
  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }
  return 0;
}
Home | Doxygen Documentation | Tutorials | Developer Tools (restricted)