NiHu
2.0
|
This tutorial explains how to use NiHu
to compute a singular double integral of the form
\( \displaystyle I = \int_{S} \int_{S} K({\bf x},{\bf y}) dS_y dS_x, \)
where \( S \) is a square surface
\( S = \left\{{\bf x} = (\xi,\eta,0) : -1 \le \xi \le 1, -1 \le \eta \le 1 \right\} \)
and the integrand kernel \( K \) is the fundamental solution of the Laplace equation in 3D
\( K({\bf x},{\bf y}) = 1 / 4\pi r, \quad r = |{\bf x}-{\bf y}| \).
The integrand contains an \( O(1/r) \) singularity when \( {\bf x} = {\bf y} \). The analytical result is
\( \displaystyle I = \frac{32}{4\pi} \left[ \log \left( 1+\sqrt{2} \right) - \frac{\sqrt{2}-1}{3} \right] \).
We need to include two modules:
We are going to work with dynamically resizeable Eigen matrices. We define two convenient typedefs for double and unsigned matrices.
Our problem domain \( S \) is represented by a mesh consisting of a single square element. For more information on mesh building, refer to the tutorial Mesh building.
We reformulate our original problem in the form
\( \displaystyle I = \sum_{i,j} I_{ij}, \quad \) where \( \quad I_{ij} = \left< w_i, \left(\mathcal{K} w_j\right)_S \right>_S, \)
where the set \( w_i \) denote a piecewise constant or isoparametric function space, and \( \mathcal{K} \) denotes the integral operator based on the Laplace kernel.
Two function spaces are generated from the mesh, in order to test both field generation options:
As the different function space views are of different types, the tester function is written as a template that can receive any function space type as parameter.
The function first allocates the result matrix. The matrix size is determined by the total number of DOF in the funcion space w
, computed by the member function NiHu::function_space_base::get_num_dofs.
After instantiating the integral operator and evaluating the weighted double integral into the matrix, the tester function prints the matrix, its element sum, and compares the sum to the analytical integral.
The complete source of the tutorial is found here: tutorial/laplace_double_integral.cpp