Introduction

The purpose of this tutorial is to demonstrate how a new family of partial differential equations (PDE) can be introduced into the NiHu framework. In a boundary element context, the introduction of a new PDE is equivalent to the implementation of its fundamental solutions. The demonstrative example is linear isotropic elastostatics.

Theory

PDE

3D linear isotropic elasticity is governed by the PDE

\( \displaystyle \sigma_{ij,j} = \delta_{ij} u_{k,ki} + \mu \left( u_{i,jj} + u_{j,ij} \right) \\ \displaystyle \mu u_{i,jj}(x) + \left(\mu + \lambda\right) u_{j,ij}(x) = 0, \quad x \in \Omega \subset \mathbb{R}^{3} \)

with the boundary conditions

\( \displaystyle t_i(x) = \bar{t}_i(x), x \in \Gamma_t, \\ u_i(x) = \bar{u}_i(x), x \in \Gamma_u \)

where \( u_i(x) \) denotes the displacement vector field, \( t_i(x) \) denotes the traction vector field, \( \mu \) and \( \lambda \) are the Lamé-coefficients, \( \Omega \) denotes the solution domain and \( \Gamma_{\alpha} \) stands for its boundary with prescribed traction or displacement boundary conditions.

BIE

The equivalent boundary integral representation of the PDE is

\( \displaystyle \int_{\Gamma} t^*_{ij}(x,y) u_j(y) \mathrm{d} y - \int_{\Gamma} u^*_{ij}(x,y) t_{j}(y) \mathrm{d} y = u_i(x) \)

where the displacement and traction fundamental solutions are

\( \displaystyle u^*_{ij} = \frac{(3-4\nu) \delta_{ij} + r,_i r,_j}{16 \pi \mu (1-\nu) r},\\ \displaystyle t^*_{ij} = \frac{-r,_n ((1-2\nu)\delta_{ij} + 3 r,_i r,_j) + (1-2\nu) (r,_i n_j - r,_j n_i)}{8 \pi (1-\nu) r^2} \)

where \( \nu \) and \( \mu \) denote the Poisson's number and the shear modulus, respectively. The displacement kernel contains a weak \( O(1/r) \) singularity, and can be integrated with blind singular quadratures. The traction kernel is strongly singular, and its \( O(1/r^2) \)-type singularity needs to be handled in a CPV sense.

Implementing the fundamental solution

In order to be able to define expressions for the fundamental solutions, we need to define three classes:

The kernel's test input \( x \)
The kernel's trial input \( y \)
The kernel's parameters

For the case of the displacement kernel, both inputs are simple locations. For the case of the traction kernel, the test input is location, while the trial input is location with normal vector. The kernel parameter class needs to encapsulate the two material properties:

The Poisson's ratio \( \nu \)
The Shear modulus \( \mu \)

As the shear modulus serves only as a scaling factor, it can be omitted in the implementation without the loss of generality.

The parameter class

First the kernel parameter class is defined that encapsulates a simple double as the Poisson's ratio:

The kernel functor

As the next step, a function object (functor) named Ukernel is defined that returns the value of the displacement fundamental solution.

As it is clear from the definition, the functor not only returns the kernel's expression but defines the kernel's return type (3d double matrix) as well.

The kernel class

The next step is the definition of the final kernel class that is going to be used in numerical integrations. This is done in three steps.

The kernel class is declared
Compile time properties of the kernel class are defined in its traits class
The kernel class is defined, derived from NiHu::kernel_base using the CRTP pattern.

Step (1): The kernel class is forward declared as

Step(2): The kernel traits are defined by specialising the traits class NiHu::kernel_traits as follows

The following traits are defined:

members test_input_t and trial_input_t select the kernel's test and trial inputs, as mentioned before.
member data_t denotes the kernel's parameter class, and member output_t denotes the kernel's output class. NiHu is capable of generating combined kernels optimised for parallel evaluations, if the kernel's inputs, parameters and outputs are defined in terms of building blocks (bricks). This feature is not exploited in the present example, meaning that the output type is a single brick wall, defined by the functor class Ukernel.
member result_dimensions indicates that the kernel is tensor valued (vector PDE)
member far_field_behaviour indicates that the kernel's far field behaviour is \( O(1/r) \).
member quadrature_family_t indicates that the kernel is integrated using Gaussian quadratures in the far field.
members is_symmetric and is_singular indicate that the kernel is symmetric and singular.

As the kernel is defined as singular, additional kernel traits need to be defined by specialising template NiHu::singular_kernel_traits

The singular traits define the singularity type as \( O(1/r) \). In 3D, such a singularity is weak, and can be cancelled out by blind quadrature methods. Member singular_quadrature_order defines that the singular quadratures are going to be evaluated with a 7-th order blind singular quadrature rule.

Step (3): Finally, the kernel class can be defined:

The kernel is derived from NiHu::kernel_base, and a constructor is provided that constructs the kernel object from a single Poisson's ratio value.

The traction kernel

The traction kernel is defined similarly. Without detailed explanation, we simply print the corresponding code snippets as follows:

Compared to the displacement kernel, the only significant difference is the singular behaviour that is defined as \( O(1/r^2) \). As this singularity is strong, it can not be handled by blind quadratures. Instead, we specialise Guiggiani's method to integrate the strong singularity.

Table of Contents