# Introduction

In this tutorial the implementation of the so-called BGT-0 boundary condition (Bayliss-Gunzburger-Turkel boundary condition of order 0) is shown. The BGT-0 boundary condition deals as an approximation to the Sommerfeld radiation condition of scattering problems on a disc.

The equation solved is the 2D Helmholtz equation

on the disc of radius , where the coefficient functions and depend on whether the TE or TM mode is considered. In the TE-case, i.e. corresponds to the third component of the magnetic field, and are given by

with the permittivity , the permeability and the frequency . On the other hand, in the TM-case, i.e. corresponds to the third component of the electric field, and are given by

In the exterior domain the total field can be split into a known incoming field and an unknown scattered field , i.e.

The coefficient functions and are assumed to be constant in . We will denote the constant values of the coefficient functions by and and they are defined according to the definitions above and where denotes the the vacuum permittivity and denotes the vacuum permeability.

Now we introduce the BGT-0 boundary condition [1]

where denotes the unit outward normal and the wave number in the exterior domain.

Finally we assume the total field and its co-normal derivative to be continuous over .

# Variational Formulation

Now we derive the corresponding variational formulation of the introduced problem: find such that

Considering its continuity over the co-normal of the total field can be written in the form

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 BGT-0 boundary condition gives

Now we use the continuity of the total field over to obtain

Incorporating this into the variational formulation yields

# 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

we do not need to prepend concepts:: to Real and Cmplx everytime.

## 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

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 and as PiecewiseConstFormula and as ConstFormula.

We write the values of , and 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);
}

The wave number

const Real k0 = omega*sqrt(mu0*eps0);

in the exterior domain and the BGT-0 coefficient

const Cmplx BGTcoeff(0.0,k0);

are introduced and the incoming field , its normal gradient and its product with the BGT-0 coefficient are set as ParsedFormula.

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))";
uIncBGTReal << "(" << real(BGTcoeff) << "*" << uIncReal.str() << "-" << imag(BGTcoeff) << "*" << uIncImag.str() << ")";
uIncBGTImag << "(" << real(BGTcoeff) << "*" << uIncImag.str() << "+" << imag(BGTcoeff) << "*" << uIncReal.str() << ")";
concepts::ParsedFormula<Cmplx> uInc( uIncReal.str(), uIncImag.str());
concepts::ParsedFormula<Cmplx> uIncBGT( uIncBGTReal.str(), uIncBGTImag.str());

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

Real ringData[3] = {innerR,(innerR+R)/2.0,R};
uint edgeAttrData[3] = {0,0,attrBoundary};
concepts::Array<Real> ring(ringData,3);
concepts::Array<uint> ringEdgeAttr(edgeAttrData,3);
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,"BGT_0_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::Space spc(msh, 1, 10);
spc.rebuild();
std::cout << std::endl << "Space: " << spc << std::endl;
graphics::drawMeshEPS(spc,"BGT_0_space.eps",100,1.0,20);

Now the trace space of the boundary can be built.

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

The right hand side

concepts::Vector<Cmplx> rhs(tspc, lform);

and the system matrix

concepts::SparseMatrix<Cmplx> S(spc.dim(), spc.dim());
hp1D::Identity<Cmplx> id1D(alpha0);
std::cout << std::endl << "System Matrix: " << S << std::endl;

are computed.

We solve the equation using the direct solver SuperLU.

#endif
std::cout << std::endl << "Solver: " << solver << std::endl;
solver(rhs, sol);

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

spc.recomputeShapefunctions();
graphics::MatlabBinaryGraphics(spc, "BGT_0", sol);

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

}
std::cout << e << std::endl;
return 1;
}
#ifdef HAS_MPI
MPI_Finalize();
#endif
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)]),
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: Space(dim = 2485, nelm = 52)
Trace Space: TraceSpace(QuadEdgeFirst(), dim = 2485, nelm = 8)
System Matrix: SparseMatrix(2485x2485, HashedSparseMatrix: 228397 (3.6986%) entries bound.)
Solver: SuperLU(n = 2485)

Plot of the refined mesh:

Matlab plots of the total field :

# References

[1] A. Bayliss, M. Gunzburger, and E. Turkel, "Boundary conditions for the numerical solution of elliptic equations in exterior domains", SIAM J. Appl. Math., vol. 42, no. 2, pp. 430-451, 1982.

# Complete Source Code

#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"
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
// ** 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 **
// ** 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);
// ** BGT-0 coefficient **
const Cmplx BGTcoeff(0.0,k0);
// ** incoming wave, its normal gradient and its product with the BGT-0 coefficient **
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))";
uIncBGTReal << "(" << real(BGTcoeff) << "*" << uIncReal.str() << "-" << imag(BGTcoeff) << "*" << uIncImag.str() << ")";
uIncBGTImag << "(" << real(BGTcoeff) << "*" << uIncImag.str() << "+" << imag(BGTcoeff) << "*" << uIncReal.str() << ")";
concepts::ParsedFormula<Cmplx> uInc( uIncReal.str(), uIncImag.str());
concepts::ParsedFormula<Cmplx> uIncBGT( uIncBGTReal.str(), uIncBGTImag.str());
// ** create mesh **
Real ringData[3] = {innerR,(innerR+R)/2.0,R};
uint edgeAttrData[3] = {0,0,attrBoundary};
concepts::Array<Real> ring(ringData,3);
concepts::Array<uint> ringEdgeAttr(edgeAttrData,3);
std::cout << std::endl << "Mesh: " << msh << std::endl;
graphics::drawMeshEPS(msh,"BGT_0_mesh.eps",100,1.0,20);
// ** space **
hp2D::Space spc(msh, 1, 10);
spc.rebuild();
std::cout << std::endl << "Space: " << spc << std::endl;
graphics::drawMeshEPS(spc,"BGT_0_space.eps",100,1.0,20);
// ** build trace space of boundary **
hp2D::TraceSpace tspc(spc,concepts::makeSet<uint>({attrBoundary}),hp2D::TraceSpace::FIRST);
std::cout << std::endl << "Trace Space: " << tspc << std::endl;
// ** right hand side **
concepts::Vector<Cmplx> rhs(tspc, lform);
// ** system matrix **
hp1D::Identity<Cmplx> id1D(alpha0);
std::cout << std::endl << "System Matrix: " << S << std::endl;
// ** solve **
#ifdef HAS_MUMPS
#else
#endif
std::cout << std::endl << "Solver: " << solver << std::endl;
solver(rhs, sol);
std::cout << " ... solved " << std::endl;
// ** plot solution **
graphics::MatlabBinaryGraphics(spc, "BGT_0", sol);
}
std::cout << e << std::endl;
return 1;
}
#ifdef HAS_MPI
MPI_Finalize();
#endif
return 0;
}
A function class to calculate element matrices for the mass matrix.
Definition: bilinearForm.hh:64
Builds the trace space of an FE space.
Definition: traces.hh:52
Piecewise constant function defined by the attribute of a cell.
Definition: formula.hh:84
static std::unique_ptr< concepts::QuadRuleFactoryTensor2d > & factory()
@ TRAPEZE
Definition: defines.hh:13
Direct sparse solver for unsymmetric matrices.
Definition: superLU.hh:70
Class for a constant formula.
Definition: constFormula.hh:26
Base class for exceptions.
Definition: exceptions.hh:86
MeshEPS< Real > drawMeshEPS(concepts::Mesh &msh, std::string filename, const Real scale=100, const Real greyscale=1.0, const unsigned int nPoints=2)
Trampoline function to create a MeshEPS.
void addInto(Matrix< H > &dest, const I fact, const uint rowoffset=0, const uint coloffset=0) const
This matrix is added as block to the given matrix dest.
A function class to calculate element matrices for the mass matrix.
Definition: bf_identity.hh:58
Linear form on edges in nD.
Definition: linearForm.hh:67
Space rebuild()
Rebuilds the space.
Class that allows to store graphical infomations in .mat files to use them in Matlab.
std::complex< Real > Cmplx
Type for a complex number. It also depends on the setting of Real.
Definition: typedefs.hh:39
virtual uint dim() const
Definition: space.hh:393
void recomputeShapefunctions()
Recompute shape functions, e.g.
double beta(const double a, const double b)
A function class to calculate element matrices for the Laplacian.
Definition: bf_laplace.hh:91
A 2D hp FEM space with continuous, piecewise polynomial basis functions.
Definition: space.hh:88
MUMPS : MUltifrontal Massively Parallel sparse direct Solver.
Definition: mumps.hh:72
Mesh for a circle.
Definition: circle.hh:47
Attributes for elements of the topology.
Definition: connector.hh:22
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich