NiHu  2.0
Single integral of the laplace kernel

Introduction

This tutorial explains how to use NiHu to compute a single singular integral of the form

\( \displaystyle I = \int_{S} K({\bf x}_0,{\bf y}) dS_y, \)

where \( S \) is a square surface

\( \displaystyle S = \left\{{\bf x} = (\xi,\eta,0) : 0 \le \xi \le 1, 0 \le \eta \le 1 \right\} \),

\( {\bf x}_0 \) is a point located at the surface's center \( (.5, .5, 0) \), and the integrand kernel \( K \) is the fundamental solution of the Laplace equation in 3D

\( \displaystyle K({\bf x},{\bf y}) = 1 / 4\pi r, \quad r = |{\bf x}-{\bf y}| . \)

The analytical result of the singular integral is

\( \displaystyle I = \log \left(1+\sqrt{2}\right) / \pi . \)

The C++ code

Mesh generation

The mesh of our surface domain will consist of only one 4-noded quadrilateral (NiHu::quad_1_elem) element:

dMatrix nodes(4,3);
nodes << 0., 0., 0., /*|*/ 1., 0., 0., /*|*/ 1., 1., 0., /*|*/ 0., 1., 0.;
uMatrix elements(1, 1+4);
elements << NiHu::quad_1_elem::id, 0, 1, 2, 3;
auto msh = NiHu::create_mesh(nodes, elements, NiHu::quad_1_tag());

Function space and weighted integral

The single integral is generalised as a double weighted integral where one of the weighting functions is a Dirac delta function and the other is constant:

\( \displaystyle I = \int_{S} K({\bf x}_0,{\bf y}) dS_y = \int_{S} \delta({\bf x}_0) \int_{S} K({\bf x},{\bf y}) 1({\bf x}) dS_y dS_x = \left< \delta_{{\bf x}_0}, \left(\mathcal{K}1\right)_S \right>_S \)

The corresponding function spaces are instantiated as

auto const &trial = NiHu::constant_view(msh);
auto const &test = NiHu::dirac(trial);

As our mesh contains one element, and the piecewise constant function space contains one DOF per element, the resulting weighted integral matrix will contain one entry.

dMatrix A(1, 1);
A.setZero();
A << test * K[trial];

Results

The integration result is printed and compared to the analytical solution.

std::cout << "WR matrix: " << A << std::endl;
double anal = std::log(1. + root_two) / pi;
std::cout << "log10 error = " << std::log10(std::abs(A.sum() / anal - 1.)) << std::endl;

The complete source of the tutorial is found here: tutorial/laplace_single_integral.cpp

NiHu::mex::real_matrix
Definition: mex_matrix_separate.hpp:133
NiHu::bessel::K
make_complex< T >::type K(T const &z)
K_nu(z) modified Bessel function.
NiHu::constant_view
const function_space_view< Mesh, field_option::constant, Dimension > & constant_view(Mesh const &msh, Dimension dim=Dimension())
factory function to transform a mesh into a constant function space
Definition: function_space.hpp:263
NiHu::create_mesh
mesh< tmp::vector< typename tag2type< Args >::type... > > create_mesh(nodes_t const &nodes, elements_t const &elements, Args...)
factory function to create a mesh from nodes and elements matrices
Definition: mesh.hpp:287
NiHu::normal_derivative_kernel
Normal derivative of a distance dependent kernel.
Definition: normal_derivative_kernel.hpp:26
NiHu::type2tag
Netafunction assigning a tag to a type.
Definition: type2tag.hpp:17
NiHu::dirac
const dirac_space< FuncSpace > & dirac(function_space_base< FuncSpace > const &space)
factory function to convert a function space into a dirac space
Definition: function_space.hpp:355
NiHu::create_integral_operator
integral_operator< Kernel > create_integral_operator(Kernel &&kernel)
factory function of an integral operator
Definition: integral_operator.hpp:326