Contents
- Introduction
- Variational Formulation
- Commented Program
- Results
- References
- Complete Source Code
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
and concepts files
#include "formula.hh"
#include "function.hh"
#include "graphics.hh"
#include "operator.hh"
#include "space.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 Cmplx epsSc (4.0, 0.0);
and defining the attributes of the boundary, the area of the scatterer and the surrounding area respectively.
const uint attrBoundary = 11;
const uint attrSc = 21;
const uint attr0 = 22;
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) {
alpha0 =
Cmplx(pow(eps0,-1.0),0.0);
}
and in the TM-case we write
if (mode == TM) {
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
.
std::stringstream uIncReal, uIncImag, uIncGradReal, uIncGradImag, uIncBGTReal, uIncBGTImag;
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() << ")";
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};
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.
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.
spc.rebuild();
std::cout << std::endl << "Space: " << spc << std::endl;
Now the trace space of the boundary
can be built.
std::cout << std::endl << "Trace Space: " << tspc << std::endl;
The right hand side
and the system matrix
A.addInto(S, 1.0);
M2D.addInto(S, -1.0);
M1D.addInto(S, -BGTcoeff);
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();
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)]),
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: 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: ![Mesh output](BGT_0_mesh.png)
Matlab plots of the total field
: ![Matlab output](BGT_0.png)
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
- Author
- Dirk Klindworth, 2011
#include <iostream>
#include "formula.hh"
#include "function.hh"
#include "graphics.hh"
#include "operator.hh"
#include "space.hh"
int main(int argc, char** argv) {
try {
#ifdef HAS_MPI
MPI_Init(&argc, &argv);
#endif
const Cmplx epsSc (4.0, 0.0);
const uint attrBoundary = 11;
const uint attrSc = 21;
const uint attr0 = 22;
enum TETM {TE, TM};
TETM mode = TM;
if (mode == TE) {
alpha0 =
Cmplx(pow(eps0,-1.0),0.0);
}
if (mode == TM) {
alpha0 =
Cmplx(1.0, 0.0);
}
const Real k0 = omega*sqrt(mu0*eps0);
const Cmplx BGTcoeff(0.0,k0);
std::stringstream uIncReal, uIncImag, uIncGradReal, uIncGradImag, uIncBGTReal, uIncBGTImag;
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() << ")";
Real ringData[3] = {innerR,(innerR+R)/2.0,R};
uint quadAttrData[3] = {attrSc,attr0,attr0};
uint edgeAttrData[3] = {0,0,attrBoundary};
std::cout << std::endl << "Mesh: " << msh << std::endl;
std::cout << std::endl << "Space: " << spc << std::endl;
std::cout << std::endl << "Trace Space: " << tspc << std::endl;
std::cout << std::endl << "System Matrix: " << S << std::endl;
#ifdef HAS_MUMPS
#else
#endif
std::cout << std::endl << "Solver: " << solver << std::endl;
solver(rhs, sol);
std::cout << " ... solved " << std::endl;
}
std::cout << e << std::endl;
return 1;
}
#ifdef HAS_MPI
MPI_Finalize();
#endif
return 0;
}