Bembel is a header-only template library. On this page we explain fundamental design considerations. This is intended to improve a fundamental understanding of the code base and provide starting points for extensions.

Curiously Recurring Template Pattern

The three essential components LinearOperator, LinearForm and Potential of the BEM are implemented in Bembel according to the Curiously Recurring Template Pattern.

Linear Operator

If you want to add an implementation of an LinearOperator, then LinearOperatorTraits must be defined and an implementation for the integral evaluation must be given.

Traits

The LinearOperatorTraits are explained using the example of the LaplaceSingleLayerOperator. These traits are passed to the implementation via template parameters and fundamentally define the properties of the software.

template <>
struct LinearOperatorTraits<LaplaceSingleLayerOperator> {
  typedef Eigen::VectorXd EigenType;
  typedef Eigen::VectorXd::Scalar Scalar;
  enum {
    OperatorOrder = -1,
    Form = DifferentialForm::Discontinuous,
    NumberOfFMMComponents = 1
  };
};

The explanation of the traits is

EigenType defines the type of the vectors.
Scalar defines if the operator is real or complex valued.
OperatorOrder defines the order of the integral operator which is necessary for the quadrature.
Form defines which B-spline basis is chosen.
NumberOfFMMComponents is used in the H2Matrix.

There are three options for the parameter Form: DifferentialForm::Continuous, DifferentialForm::DivConforming and DifferentialForm::Discontinuous. This selects a scalar continuous, vector-valued and div conforming or scalar discontinuous B-spline basis for the discretization.

Evaluate Integrand

Furthermore, there must be an implementation of the function evaluateIntegrand_impl().

void evaluateIntegrand_impl(
    const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2,
    Eigen::Matrix<
        typename LinearOperatorTraits<LaplaceSingleLayerOperator>::Scalar,
        Eigen::Dynamic, Eigen::Dynamic> *intval);

The SuperSpace manages the basis functions and interactions thereof. Since a BEM always evaluates double integrals, two SurfacePoints are provided, which also contain the information about derivatives. Finally, the result must be added to intval.

Potential

Just as above, as with the LinearOperator, we define and traits and the implementation of the integral evaluation.

Traits

template <typename LinOp>
struct PotentialTraits<LaplaceSingleLayerPotential<LinOp>> {
  typedef Eigen::VectorXd::Scalar Scalar;
  static constexpr int OutputSpaceDimension = 1;
};

Scalar defines if the potential is real or complex valued.
OutputSpaceDimension defines if the potential is scalar or vector valued.

Evaluate Integrand

Eigen::Matrix<typename PotentialReturnScalar<
                  typename LinearOperatorTraits<LinOp>::Scalar, double>::Scalar,
              1, 1>
evaluateIntegrand_impl(const FunctionEvaluator<LinOp> &fun_ev,
                       const ElementTreeNode &element,
                       const Eigen::Vector3d &point, const SurfacePoint &p);

For the return type the Scalar type must be derived from the LinearOperator and Potential. Furthermore, the number of rows in the returned matrix corresponds to the dimension of the output space. The FunctionEvaluator fun_ev takes care over the evaluation of the surface quantity, which needs an ElementTreeNode element as input. For the SurfacePoint p the contribution of the potential in the given point is returned.

Linear Form

We also define the scalar type and the evaluation of the integral for the right-hand side integrated using the LinearForm.

Traits

template <typename ScalarT>
struct LinearFormTraits<DirichletTrace<ScalarT>> {
  typedef ScalarT Scalar;
};

Scalar defines if the linear form is real or complex valued.

Evaluate Integrand

void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p,

Eigen::Matrix<Scalar, Eigen::Dynamic, 1> *intval);

The SuperSpace manages the basis functions. Integral is evaluated at the SurfacePoint and the value is returned in intval.