NiHu  2.0
Integral operators and weighted residuals

Introduction

This tutorial explains how NiHu implements weighted double integral matrices of the form

\( \displaystyle W_{ij} = \left<t,\left(\mathcal{K}d\right)_S\right>_F = \int_{F} t_i({\bf x}) \int_{S} K({\bf x}, {\bf y}) d_j({\bf y}) \mathrm{d}S_y \mathrm{d}S_x \).

Integral operators

Integral operators are defined by their kernel functions. NiHu creates integral operators using the function NiHu::create_integral_operator, called with a kernel function instance. The lines below, for example

#include "library/laplace_kernel.hpp"
auto K = NiHu::create_integral_operator(NiHu::laplace_3d_DLP_kernel());

instantiate the double layer potential kernel of the Laplace equation in 3D from the library, and transform the kernel to an integral operator.

A special integral operator is the identity operator with the kernel function \( I({\bf x}, {\bf y}) = \delta({\bf y}-{\bf x})\). This is instantiated without argument, as

auto I = identity_integral_operator();

Integral operations

Integral transform

The most important operation with an integral operator is letting it act on a function (or function space) \( d({\bf y}) \) as:

\( \displaystyle p({\bf x}) = (\mathcal{K}d)_S({\bf x}) \)

NiHu implements this operation by indexing the integral operator with a function space. The integration domain \( S \) (the mesh) is included in the function space's definition:

auto const &d = NiHu::constant_view(my_mesh);
auto K = NiHu::create_integral_operator(my_kernel_instance);
auto p = K[d];

Inner product

The inner product

\( \displaystyle W = \left< t, p \right> _F= \int_{F} t({\bf x}) p({\bf x}) \mathrm{d} F_x \)

is implemented as a multiplication between a function space and an integral transform:

auto const &t = NiHu::constant_view(other_mesh);
auto W = t * K[d]; // or W = t * p;

where the integration domain \( F \) is contained by the function space t;

Finally, the result of the weighted double integral can be evaluated into a matrix using the << operator:

myMatrix << ( t * K[d] ); // or myMatrix << W;
Note
Although we deeply believe in the do-not-memorise-precedence-use-parentheses-instead theorem, we mention here that the form
myMatrix << t * K[d];

is also valid, as << is weaker than *.

The inner product multiplication only works between a function space and an integral transform. If we need to compose the inner product of two function spaces (defined on the same mesh, of course), then the identity operator needs to be applied:

auto const &v = constant_view(my_mesh);
auto const &w = isoparametric_view(my_mesh);
auto I = identity_integral_operator();
matrix << v * I[w];