BEM++  2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
Classes | Functions
Maxwell equations in 3D

Classes

class  Bempp::Maxwell3dDoubleLayerPotentialOperator< BasisFunctionType_ >
 Double-layer potential operator for the Maxwell equations in 3D. More...
 
class  Bempp::Maxwell3dFarFieldDoubleLayerPotentialOperator< BasisFunctionType_ >
 Far-field limit of the double-layer potential operator for the Maxwell equations in 3D. More...
 
class  Bempp::Maxwell3dFarFieldSingleLayerPotentialOperator< BasisFunctionType_ >
 Far-field limit of the single-layer potential operator for the Maxwell equations in 3D. More...
 
class  Bempp::Maxwell3dIdentityOperator< BasisFunctionType_, ResultType_ >
 "Identity operator" for Maxwell equations in 3D. More...
 
class  Bempp::Maxwell3dSingleLayerPotentialOperator< BasisFunctionType_ >
 Single-layer potential operator for the Maxwell equations in 3D. More...
 

Functions

template<typename BasisFunctionType >
BoundaryOperator
< BasisFunctionType, typename
ScalarTraits
< BasisFunctionType >
::ComplexType > 
Bempp::maxwell3dDoubleLayerBoundaryOperator (const shared_ptr< const Context< BasisFunctionType, typename ScalarTraits< BasisFunctionType >::ComplexType > > &context, const shared_ptr< const Space< BasisFunctionType > > &domain, const shared_ptr< const Space< BasisFunctionType > > &range, const shared_ptr< const Space< BasisFunctionType > > &dualToRange, typename ScalarTraits< BasisFunctionType >::ComplexType waveNumber, const std::string &label="", int symmetry=NO_SYMMETRY, bool useInterpolation=false, int interpPtsPerWavelength=DEFAULT_HELMHOLTZ_INTERPOLATION_DENSITY)
 Construct a double-layer boundary operator for Maxwell equations in 3D. More...
 
template<typename BasisFunctionType >
BoundaryOperator
< BasisFunctionType, typename
ScalarTraits
< BasisFunctionType >
::ComplexType > 
Bempp::maxwell3dSingleLayerBoundaryOperator (const shared_ptr< const Context< BasisFunctionType, typename ScalarTraits< BasisFunctionType >::ComplexType > > &context, const shared_ptr< const Space< BasisFunctionType > > &domain, const shared_ptr< const Space< BasisFunctionType > > &range, const shared_ptr< const Space< BasisFunctionType > > &dualToRange, typename ScalarTraits< BasisFunctionType >::ComplexType waveNumber, const std::string &label="", int symmetry=NO_SYMMETRY, bool useInterpolation=false, int interpPtsPerWavelength=DEFAULT_HELMHOLTZ_INTERPOLATION_DENSITY)
 Construct a single-layer boundary operator associated with the Maxwell equations in 3D. More...
 

Detailed Description

This submodule contains classes relevant to the implementation of boundary operators and potential operators related to the time-harmonic Maxwell equations in 3D. Inside a region with a uniform scalar permittivity and permeability, these equations can be reduced to the form

\[ \boldsymbol\nabla \times \boldsymbol\nabla \times \boldsymbol u - k^2 \boldsymbol u = 0, \]

where $\boldsymbol u$ stands either for the electric or the magnetic field and the number $k \neq 0$ is referred to as the wave number.

The treatment of Maxwell equations in BEM++ closely follows that of A. Buffa and R. Hiptmair, "Galerkin Boundary Element Methods for Electromagnetic Scattering", in "Computational Methods in Wave Propagation", M. Ainsworth, ed., Springer, New York 2003, pp. 85-126. In the following we cite the principal definitions relevant to the implementation of BEM in the Maxwell case.

The starting point is the Stratton-Chu theorem (Theorem 6 in Buffa and Hiptmair), which plays a part analogous to the Green's representation theorem for the Laplace equation. It states that any solution $u$ of the Maxwell equations in a bounded domain $\Omega$ with boundary $\Gamma$ has the representation

\[ \boldsymbol u = \boldsymbol \Psi_{\mathrm{SL}}(\gamma_{\mathrm{N,int}}\boldsymbol u) + \boldsymbol \Psi_{\mathrm{DL}}(\gamma_{\mathrm{D,int}}\boldsymbol u). \]

(all the symbols will be defined shortly). Furthermore, any solution $u$ that satisfies the Maxwell equations in a domain $\Omega^{\mathrm c}$ extending to infinity and meets the Silver-Mueller radiation conditions at infinity has the representation

\[ \boldsymbol u = -\boldsymbol \Psi_{\mathrm{SL}}(\gamma_{\mathrm{N,ext}}\boldsymbol u) - \boldsymbol \Psi_{\mathrm{DL}}(\gamma_{\mathrm{D,ext}}\boldsymbol u). \]

The action of the single-layer potential operator is defined by

\[ (\boldsymbol \Psi_{\mathrm{SL}} \boldsymbol v)(\boldsymbol x) \equiv \mathrm i k \int_{\Gamma} G(\boldsymbol x, \boldsymbol y) \boldsymbol v(\boldsymbol y) \mathrm\Gamma(\boldsymbol y) - \frac{1}{\mathrm i k} \boldsymbol\nabla_{\boldsymbol x} \int_{\Gamma} G(\boldsymbol x, \boldsymbol y) (\boldsymbol \nabla_\Gamma \cdot \boldsymbol v)(\boldsymbol y) \mathrm\Gamma(\boldsymbol y), \]

where

\[ G(\boldsymbol x, \boldsymbol y) \equiv \frac{\exp(\mathrm i k \lvert \boldsymbol x - \boldsymbol y \rvert)} {4\pi \lvert \boldsymbol x - \boldsymbol y \rvert} \]

is the Green's function of the Helmholtz equation. The action of the double-layer potential operator, in turn, is defined by

\[ (\Psi_{\mathrm{DL}} \boldsymbol v)(\boldsymbol x) \equiv \boldsymbol\nabla_{\boldsymbol x} \times \int_{\Gamma} G(\boldsymbol x, \boldsymbol y) \boldsymbol v(\boldsymbol y) \mathrm\Gamma(\boldsymbol y). \]

The interior Dirichlet trace $\gamma_{\mathrm{D,int}}\boldsymbol u$ at a point $\boldsymbol x \in \Gamma$ is defined as

\[(\gamma_{\mathrm{D,int}}\boldsymbol u)(\boldsymbol x) \equiv \boldsymbol u\rvert_{\Gamma,\mathrm{int}}(\boldsymbol x) \times \boldsymbol n(\boldsymbol x),\]

where $\boldsymbol n$ is the outward unit vector normal to $\Gamma$ at $\boldsymbol x$ and $\boldsymbol u\rvert_{\Gamma,\mathrm{int}}(\boldsymbol x)$ is the limit of $\boldsymbol u(\boldsymbol y)$ as $\boldsymbol y$ approaches $\boldsymbol x$ from within $\Omega$. The interior Neumann trace $\gamma_{\mathrm{N,int}}\boldsymbol u$ at $\boldsymbol x \in \Gamma$ is defined as

\[(\gamma_{\mathrm{N,int}}\boldsymbol u)(\boldsymbol x) \equiv \frac{1}{\mathrm i k} (\boldsymbol \nabla \times \boldsymbol u)\rvert_{\Gamma,\mathrm{int}} (\boldsymbol x) \times \boldsymbol n(\boldsymbol x).\]

Note that compared to Buffa and Hiptmair, we include an additional $\mathrm i$ factor in the denominator of the Neumann trace and in the numerator of the single-layer potential operator, as this makes both quantities real-valued in the special case of $k$ purely imaginary (and $\boldsymbol u$ purely real).

The exterior Dirichlet and Neumann traces are defined analogously, with the relevant quantities assumed to approach the point $\boldsymbol x$ from within the complement of $\Omega$.

The single-layer boundary operator $\boldsymbol S$ is defined as the average of the interior and exterior Dirichlet traces of the single-layer potential operator:

\[ \boldsymbol S = \frac12 (\gamma_{\mathrm{D,int}} \boldsymbol \Psi_{\mathrm{SL}} + \gamma_{\mathrm{D,ext}} \boldsymbol \Psi_{\mathrm{SL}}). \]

The double-layer boundary operator $\boldsymbol C$ is defined, analogously, as the average of the interior and exterior Dirichlet traces of the double-layer potential operator.

To exploit the symmetry between the electric and magnetic field in Maxwell equations, it is advantageous to use the following antisymmetric pseudo-inner product in the Galerkin discretisation of the above boundary operators:

\[ \langle \boldsymbol u, \boldsymbol v \rangle_{\boldsymbol\tau,\Gamma} \equiv \int_\Gamma \boldsymbol u^* \cdot (\boldsymbol v \times \boldsymbol n) \mathrm d\Gamma, \]

where $*$ denotes complex conjugation. The weak forms of the single- and double-layer boundary operators with respect to this pseudo-inner product are given by

\[ \langle \boldsymbol u, \boldsymbol S \boldsymbol v\rangle_{\boldsymbol\tau,\Gamma} = \int_\Gamma \int_\Gamma G(\boldsymbol x, \boldsymbol y) \biggl[-\mathrm i k \boldsymbol u^*(\boldsymbol x) \cdot \boldsymbol v(\boldsymbol y) -\frac{1}{\mathrm i k} (\boldsymbol \nabla_\Gamma \cdot \boldsymbol u^*)(\boldsymbol x) (\boldsymbol \nabla_\Gamma \cdot \boldsymbol v)(\boldsymbol y) \biggr] \mathrm{d}\Gamma(\boldsymbol x)\,\mathrm{d}\Gamma(\boldsymbol y) \]

and

\[ \langle \boldsymbol u, \boldsymbol C \boldsymbol v\rangle_{\boldsymbol\tau,\Gamma} = \int_\Gamma \int_\Gamma \boldsymbol\nabla_{\boldsymbol x}G(\boldsymbol x, \boldsymbol y) \cdot [\boldsymbol u^*(\boldsymbol x) \times \boldsymbol v(\boldsymbol y)] \mathrm{d}\Gamma(\boldsymbol x)\,\mathrm{d}\Gamma(\boldsymbol y). \]

Solution of particular Maxwell problems relies on the following boundary integral equations imposed on $\Gamma$:

\begin{align*} \biggl(-\frac12 \boldsymbol I + \boldsymbol C_{\mathrm{int}}\biggr) \gamma_{\mathrm{D,int}}\boldsymbol u + \boldsymbol S_{\mathrm{int}} \gamma_{\mathrm{N,int}}\boldsymbol u &= 0 \\ -\boldsymbol S_{\mathrm{int}} \gamma_{\mathrm{D,int}}\boldsymbol u + \biggl(-\frac12 \boldsymbol I + \boldsymbol C_{\mathrm{int}}\biggr) \gamma_{\mathrm{N,int}}\boldsymbol u &= 0 \\ \biggl(\frac12 \boldsymbol I + \boldsymbol C_{\mathrm{ext}}\biggr) \gamma_{\mathrm{D,ext}}\boldsymbol u + \boldsymbol S_{\mathrm{ext}} \gamma_{\mathrm{N,ext}}\boldsymbol u &= 0 \\ -\boldsymbol S_{\mathrm{ext}} \gamma_{\mathrm{D,ext}}\boldsymbol u + \biggl(\frac12 \boldsymbol I + \boldsymbol C_{\mathrm{ext}}\biggr) \gamma_{\mathrm{N,ext}}\boldsymbol u &= 0. \end{align*}

Here, the symbol $\boldsymbol I$ denotes the identity operator. It is crucial to remember that the weak forms of $\boldsymbol S$ and $\boldsymbol C$ are defined with respect to the pseudo-inner product rather than the standard sesquilinear inner product. Thus, in practical implementation an operator whose weak form represents the identity operator under the pseudo-inner product is indispensable. Such an operator is provided by the Maxwell3dIdentityOperator class.

Function Documentation

template<typename BasisFunctionType >
BoundaryOperator< BasisFunctionType, typename ScalarTraits< BasisFunctionType >::ComplexType > Bempp::maxwell3dDoubleLayerBoundaryOperator ( const shared_ptr< const Context< BasisFunctionType, typename ScalarTraits< BasisFunctionType >::ComplexType > > &  context,
const shared_ptr< const Space< BasisFunctionType > > &  domain,
const shared_ptr< const Space< BasisFunctionType > > &  range,
const shared_ptr< const Space< BasisFunctionType > > &  dualToRange,
typename ScalarTraits< BasisFunctionType >::ComplexType  waveNumber,
const std::string &  label = "",
int  symmetry = NO_SYMMETRY,
bool  useInterpolation = false,
int  interpPtsPerWavelength = DEFAULT_HELMHOLTZ_INTERPOLATION_DENSITY 
)

Construct a double-layer boundary operator for Maxwell equations in 3D.

This function constructs a BoundaryOperator object representing the double-layer boundary operator $\boldsymbol C$ for the Maxwell equations in 3D, as defined in Maxwell equations in 3D.

Parameters
[in]contextA Context object that will be used to build the weak form of the boundary operator when necessary.
[in]domainFunction space being the domain of the boundary operator.
[in]rangeFunction space being the range of the boundary operator.
[in]dualToRangeFunction space dual to the the range of the boundary operator.
[in]waveNumberWave number. See Maxwell equations in 3D for its definition.
[in]labelTextual label of the operator. If empty, a unique label is generated automatically.
[in]symmetrySymmetry of the weak form of the operator. Can be any combination of the flags defined in the enumeration type Symmetry.
[in]useInterpolationIf set to false (default), the standard exp() function will be used to evaluate the exponential factor occurring in the kernel. If set to true, the exponential factor will be evaluated by piecewise-cubic interpolation of values calculated in advance on a regular grid. This normally speeds up calculations, but might result in a loss of accuracy. This is an experimental feature: use it at your own risk.
[in]interpPtsPerWavelengthIf useInterpolation is set to true, this parameter determines the number of points per "effective wavelength" (defined as $2\pi/|k|$, where $k$ = waveNumber) used to construct the interpolation grid. The default value (5000) is normally enough to reduce the relative or absolute error, whichever is smaller, below 100 * machine precision. If useInterpolation is set to false, this parameter is ignored.

None of the shared pointers may be null and the spaces range and dualToRange must be defined on the same grid, otherwise an exception is thrown.

See Also
Maxwell3dDoubleLayerPotentialOperator()
template<typename BasisFunctionType >
BoundaryOperator< BasisFunctionType, typename ScalarTraits< BasisFunctionType >::ComplexType > Bempp::maxwell3dSingleLayerBoundaryOperator ( const shared_ptr< const Context< BasisFunctionType, typename ScalarTraits< BasisFunctionType >::ComplexType > > &  context,
const shared_ptr< const Space< BasisFunctionType > > &  domain,
const shared_ptr< const Space< BasisFunctionType > > &  range,
const shared_ptr< const Space< BasisFunctionType > > &  dualToRange,
typename ScalarTraits< BasisFunctionType >::ComplexType  waveNumber,
const std::string &  label = "",
int  symmetry = NO_SYMMETRY,
bool  useInterpolation = false,
int  interpPtsPerWavelength = DEFAULT_HELMHOLTZ_INTERPOLATION_DENSITY 
)

Construct a single-layer boundary operator associated with the Maxwell equations in 3D.

This function constructs a BoundaryOperator object representing the single-layer boundary operator $\boldsymbol S$ for the Maxwell equations in 3D, as defined in Maxwell equations in 3D.

Parameters
[in]contextA Context object that will be used to build the weak form of the boundary operator when necessary.
[in]domainFunction space being the domain of the boundary operator.
[in]rangeFunction space being the range of the boundary operator.
[in]dualToRangeFunction space dual to the the range of the boundary operator.
[in]waveNumberWave number. See Maxwell equations in 3D for its definition.
[in]labelTextual label of the operator. If empty, a unique label is generated automatically.
[in]symmetrySymmetry of the weak form of the operator. Can be any combination of the flags defined in the enumeration type Symmetry.
[in]useInterpolationIf set to false (default), the standard exp() function will be used to evaluate the exponential factor occurring in the kernel. If set to true, the exponential factor will be evaluated by piecewise-cubic interpolation of values calculated in advance on a regular grid. This normally speeds up calculations, but might result in a loss of accuracy. This is an experimental feature: use it at your own risk.
[in]interpPtsPerWavelengthIf useInterpolation is set to true, this parameter determines the number of points per "effective wavelength" (defined as $2\pi/|k|$, where $k$ = waveNumber) used to construct the interpolation grid. The default value (5000) is normally enough to reduce the relative or absolute error, whichever is smaller, below 100 * machine precision. If useInterpolation is set to false, this parameter is ignored.

None of the shared pointers may be null and the spaces range and dualToRange must be defined on the same grid, otherwise an exception is thrown.

This operator supports local-mode ACA assembly. It does not support hybrid-mode assembly, though.

If local-mode ACA assembly is requested (see AcaOptions::mode), after discretization, the weak form of this operator is stored as

\[ A = -ik \sum_{i=1}^3 P_i A_{\textrm{d}} Q_i - \frac{1}{ik} R_i A_{\textrm{d}} S_i, \]

where:

  • $k$ is the wave number
  • $A_{\textrm{d}}$ is the weak form of the single-layer modified Helmholtz boundary operator discretised with test and trial functions being the restrictions of the basis functions of domain and range to individual elements;
  • $P_i$ is the sparse matrix whose transpose represents the expansion of the $i$th component of the basis functions of dualToRange in the single-element test functions mentioned above;
  • $Q_i$ is the sparse matrix representing the expansion of the $i$th component of the basis functions of domain in the single-element trial functions;
  • $R_i$ is the sparse matrix whose transpose represents the expansion of the surface divergence of the basis functions of dualToRange in the single-element test functions;
  • $S_i$ is the sparse matrix representing the expansion of the surface divergence of the basis functions of domain in the single-element trial functions.
Template Parameters
BasisFunctionTypeType of the values of the basis functions into which functions acted upon by the operator are expanded. It can take the following values: float, double, std::complex<float> and std::complex<double>.
See Also
maxwell3dSingleLayerPotentialOperator()

References Bempp::AssemblyOptions::ACA, Bempp::AssemblyOptions::acaOptions(), Bempp::AssemblyOptions::assemblyMode(), and Bempp::AcaOptions::mode.