Double integral of the Laplace kernel

Table of Contents

• Introduction
• The C++ code
  • Libraries and typedefs
  • The main function
  • The testing function

Introduction

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] \).

The C++ code

Libraries and typedefs

We need to include two modules:

• core/weighted_residual.hpp is the main module of the NiHu C++ core
• library/laplace_kernel.hpp defines our kernel.

#include "library/lib_element.hpp"
#include "core/weighted_residual.hpp"
#include "library/laplace_kernel.hpp"

We are going to work with dynamically resizeable Eigen matrices. We define two convenient typedefs for double and unsigned matrices.

// dynamically resizeable double matrix
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dMatrix;
// dynamically resizeable unsigned matrix
typedef Eigen::Matrix<unsigned, Eigen::Dynamic, Eigen::Dynamic> uMatrix;

The main function

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.

    // nodal coordinates in 3D, 3 nodes
    dMatrix nodes(4,3);
    nodes <<
        -1.0, -1.0, 0.0,
         1.0, -1.0, 0.0,
        -1.0,  1.0, 0.0,
         1.0,  1.0, 0.0;
 
    // element nodal indices
    uMatrix elements(1, 1+4);
    elements << NiHu::quad_1_elem::id, 0, 1, 3, 2;
 
    // create the mesh
    auto msh = NiHu::create_mesh(nodes, elements, NiHu::tria_1_tag(), NiHu::quad_1_tag());

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:

    // create a piecewise constant function space and call the tester
    std::cout << "Testing with constant field" << std::endl;
    tester(NiHu::constant_view(msh));
 
    // create a piecewise linear function space and call the tester
    std::cout << "Testing with isoparametric field" << std::endl;
    tester(NiHu::isoparametric_view(msh));

The testing function

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.

template <class func_space>
void tester(func_space const &w)
{
    using namespace boost::math::double_constants;
    
    // compute number of DOF and allocate result matrix
    size_t nDOF = w.get_num_dofs();
    dMatrix I(nDOF, nDOF);
    I.setZero();
 
    // create integral operator from kernel and perform weighted double integral
    auto K = NiHu::create_integral_operator(NiHu::laplace_3d_SLP_kernel());
    I <<  w * K[w];
 
    // Display matrix elements and their sum
    std::cout << "WR matrix:\n" << I << std::endl;
    std::cout << "sum of elements: " << I.sum() << std::endl;
 
    // Compare to analytical solution
    double anal = 32.0 * (std::log(1.0+root_two)-(root_two-1.0)/3.0) / 4.0 / pi;
    std::cout << "log10 error = " << std::log10(std::abs(I.sum() / anal - 1.0)) << std::endl;
}

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