The Kernel Class

The kernel class provides a convenient way to specify integral operators in classical potential theory. The kernel object has several attributes including parameters required to define the kernel, the type of singularity for determining the quadrature to be used, function handles for evaluating the kernel, and fast multipole method (FMM) acceleration routines if available.

The constructor for the kernel objects take the form

kernel('PDE name', 'kernel name', varargin)

The kernel class documentation lists the information that can be stored in the object

%KERNEL class which describes an integral kernel, usually
% related to the solution of a partial differential equation (PDE).
%
%   K = KERNEL(NAME, TYPE) constructs a kernel of the specified name and
%   type. The currently supported kernels names and types are:
%
%      NAME                                         TYPE
%      ----                                         ----
%      'laplace'    ('lap', 'l')                    's', 'd', 'sp', 'c'
%      'helmholtz'  ('helm', 'h')                   's', 'd', 'sp', 'dp', 'c'
%                                                   'cp'
%      'helmholtz1d'  ('helm1d', 'h1d')             's'
%      'helmholtz difference'                       's', 'd', 'sp', 'dp'
%         ('helmdiff', 'hdiff', 'helm_diff')
%      'elasticity' ('elast', 'e')                  's', 'strac', 'd', 'dalt'
%      'stokes'     ('stok', 's')                   'svel', 'spres', 
%                                                   'strac', 'sgrad'
%                                                   'dvel', 'dpres',
%                                                   'dtrac', 'dgrad'
%                                                   'cvel', 'cpres',
%                                                   'ctrac', 'cgrad'
%      'zeros'       ('zero','z') 
%      'axis sym helmholtz'                         's' 'd' 'sp' 'c'
%         ('axissymh', 'axissymhelm')
%      'axis sym helmholtz difference'              's' 'd' 'sp' 'dp'
%         ('axissymhdiff', 'axissymhelmdiff', 'axissymhelm_diff') 
%      'quasiperiodic helmholtz'                    's', 'd', 'sp', 'dp', 'c'
%         ('helmquas', 'hq')                        'cp'
%   The types may also be written in longer form, e.g. 'single', 'double',
%   'sprime', 'combined', 'svelocity', 'spressure', 'straction',
%   'dvelocity', 'dpressure', 'dtraction'.
%
%   Some kernels may require extra parameters to be specified, e.g. the
%   Helmholtz wavenumber, combined field parameter, or material parameters.
%   Such parameters should be passed as trailing arguments. For instance,
%   the Helmholtz single-layer kernel with wavenumber ZK can be constructed
%   as KERNEL('helmholtz', 's', ZK).
%
%   The kernel class includes an evaluator, K.eval(s,t), which evaluates
%   the kernel K between sources s and targets t specified in the ptinfo
%   format, i.e. a struct with entries s.r, s.d, s.d2, s.n corresponding to
%   the positions, first and second derivatives, and normals (the latter
%   three being only defined for points on a parameterized curve).
%
%   The kernel class can also be used to store other useful methods and
%   data for working with the kernel in a standardized format. For example,
%   for a kernel K:
%
%      - K.opdims: a 2-vector [m n] specifying the dimensions of the
%        kernel values. K.eval(s,t) is an m x n matrix if s and t are a
%        single source and target, respectively. For scalar kernels,
%        opdims = [1 1].
%
%      - K.sing: a string/character array specifying the singularity  
%        strength of the kernel for sources and targets which are on
%        the same curve. Currently recognized singularities are considered
%        as part of a hierarchy
%        - 'smooth' a smooth integral kernel
%        - 'log' sum of above type kernel and phi(s)log(s-t)
%        - 'pv' sum of above type kernels and phi(s)/(s-t)
%        - 'hs' sum of above type kernels and phi(s)/(s-t)^2
%
%      - K.fmm: A function handle which calls the FMM for the corresponding
%        kernel. K.fmm(eps, s, t, sigma) evaluates the kernel with density
%        sigma from sources s to targets t with accuracy eps.
%
%      - K.splitinfo, which provides the quantities needed to evaluate
%        integrals using the Helsing-Ojala kernel-split quadrature
%        technique.

The eval property and the ptinfo format

The eval property in the kernel class is the function which actually evaluates the interaction between sources and targets. It expects two structs, usually denoted s and t, which specify the sources and targets respectively.

Points are specified as structs, with the same naming conventions as the properties of a chunker object. We refer to this specification as the ptinfo format. In particular, if pts is a collection of points in the ptinfo format, then it has the following properties:

• pts.r - 2 x npts array of (x,y) coordinates for each point [required]

• pts.d - 2 x npts array of (x,y) coordinates of derivative for each point [optional]

• pts.d2 - 2 x npts array of (x,y) coordinates of second derivative for each point [optional]

• pts.n - 2 x npts array of (x,y) normal vector for each point [optional]

• pts.data - m x npts array of m-vectors of data for each point [optional, not fully supported]

Aside from the pts.r property, the remaining properties are optional and most are only meaningful for points on a curve.

Suppose that kern is the kernel object corresponding to a scalar kernel function of the form \(K(x,y)\) and that s and t are structs specifying a collection of ns sources and nt targets, respectively. Then, kern.eval(s,t) will return an nt x ns matrix with the entries

\[\begin{split}\begin{bmatrix} K(x^{(1)},y^{(1)}) & K(x^{(1)},y^{(2)}) & \cdots & K(x^{(1)},y^{(\texttt{ns})}) \\ K(x^{(2)},y^{(1)}) & K(x^{(2)},y^{(2)}) & \cdots & K(x^{(2)},y^{(\texttt{ns})}) \\ \vdots & \vdots & & \vdots \\ K(x^{(\texttt{nt})},y^{(1)}) & K(x^{(\texttt{nt})},y^{(2)}) & \cdots & K(x^{(\texttt{nt})},y^{(\texttt{ns})}) \end{bmatrix}\end{split}\]

where \(x^{(i)}\) is the \(i\)th target with coordinates t.r(:,i) and \(y^{(j)}\) is the \(j\)th source with coordinates s.r(:,j).

If \(K\) is a vector- or matrix-valued kernel, then the dimensions of the kernel are stored in the opdims property and the output is a (nt*kern.opdims(1)) x (ns*kern.opdims(2)) matrix. The matrix is organized into nt x ns blocks, each of size kern.opdims(1) x kern.opdims(2). For example, if \(K\) is a \(2 \times 2\) matrix valued kernel, then opdims = [2 2] and kern.eval(s,t) will return a (nt*kern.opdims(1)) x (ns*kern.opdims(2)) matrix with the entries

\[\begin{split}\begin{bmatrix} K_{11}(x^{(1)},y^{(1)}) & K_{12}(x^{(1)},y^{(1)}) & \cdots & K_{11}(x^{(1)},y^{(\texttt{ns})}) & K_{12}(x^{(1)},y^{(\texttt{ns})}) \\ K_{21}(x^{(1)},y^{(1)}) & K_{22}(x^{(1)},y^{(1)}) & \cdots & K_{21}(x^{(1)},y^{(\texttt{ns})}) & K_{22}(x^{(1)},y^{(\texttt{ns})}) \\ \vdots & \vdots & & \vdots & \vdots \\ K_{11}(x^{(\texttt{nt})},y^{(1)}) & K_{12}(x^{(\texttt{nt})},y^{(1)}) & \cdots & K_{11}(x^{(\texttt{nt})},y^{(\texttt{ns})}) & K_{12}(x^{(\texttt{nt})},y^{(\texttt{ns})}) \\ K_{21}(x^{(\texttt{nt})},y^{(1)}) & K_{22}(x^{(\texttt{nt})},y^{(1)}) & \cdots & K_{21}(x^{(\texttt{nt})},y^{(\texttt{ns})}) & K_{22}(x^{(\texttt{nt})},y^{(\texttt{ns})}) \\ \end{bmatrix}\end{split}\]

Defining your own kernel class object

While the integral kernels for many common PDEs are available as built-ins in chunkIE, the utilities will work for a broad range of integral kernels. To define your own kernel class object you can first obtain an empty kernel via kern=kernel() and then update the properties as needed. It is generally simpler to define the eval function in a separate file. The example below implements the kernel evaluator for the kernel \(K(x,y) = \exp(\imath k|x-y|)\) in a file called expkernel.m.

function mat = expkernel(s,t,k)
%EXPKERNEL this implements the kernel evaluator function for the
% kernel exp(ik|x-y|) in the chunkIE kernel format as part of the
% chunkIE guide   
%
% Syntax:
%    mat = expkernel(s,t,k) 
%   

  xs = s.r(1,:);
  ys = s.r(2,:);
  xt = t.r(1,:);
  yt = t.r(2,:);

  mat = exp(1i*k*sqrt((xt.'-xs).^2+(yt.'-ys).^2));

You can then create an empty kernel object and point the eval property to this function:

% start with an empty kernel and add info

kern = kernel();

zk = 1.3; % parameter
kern.eval = @(s,t) expkernel(s,t,zk);
kern.opdims = [1 1]; % for a scalar kernel, dims are [1 1]

Combining kernel objects

Given two kernel objects with compatible dimensions, it is possible to add and scale them.

% start with an empty kernel and add info

kern1 = kernel('lap','d');
kern2 = kernel('lap','s');

% this creates a combined layer kernel
kern = -2*(kern1 + kern2);

For a chunkgraph with ne edges, it is possible to specify interactions between each pair of edges in a ne x ne matrix of kernel objects. Note that most chunkIE routines will require that the opdims(2) property is the same for all kernels within any column of the kernel matrix and the opdims(1) property is the same for all kernels within any row of the kernel matrix.

A simple way to create a matrix of kernel objects is to first build an empty matrix of the desired size and then to assign each entry a kernel.

% suppose edge 1 has a dirichlet condition and edge 2 has 
% a neumann condition. we can use a double layer potential 
% on edge 1 and a single layer potential on edge 2

kern11 = kernel('lap','d'); % influence of edge 1 on itself
kern21 = kernel('lap','dprime'); % influence of edge 1 on edge 2
kern12 = kernel('lap','s'); % influence of edge 2 on edge 1
kern22 = kernel('lap','sprime'); % influence of edge 2 on itself

kernmat(2,2) = kernel;
kernmat(1,1) = kern11; kernmat(2,1) = kern21;
kernmat(1,2) = kern12; kernmat(2,2) = kern22;

A matrix of kernels object can also be used to create matrix-valued kernels. A matrix of kernels can be merged into a single matrix-valued kernel using the kernel constructor.

% for a transmission problem it is common to have the solution defined 
% as the sum of two layer potentials. in the exterior of the obstacle, 
% there is one wave speed, k_2 and in the interior there is another, k_1.
% One version of the boundary condition is that the potential and its
% normal derivative should be continuous across the boundary, thus the
% difference in these quantities is prescribed in terms of the incident 
% field

% parameters 
zk1 = 1.5;
zk2 = 1.1;
zks = [zk1,zk2];

% a transmission boundary condition set of kernels 
% first column -> unknown for double layer
% second column -> unknown for single layer
% first row -> jump in potential
% second row -> jump in normal derivative 

% using the "difference" kernels 
kern11 = kernel.helm2ddiff('d',zks); % D_{k1} - D_{k2}
kern21 = kernel.helm2ddiff('dprime',zks); 
kern12 = kernel.helm2ddiff('s',zks); 
kern22 = kernel.helm2ddiff('sprime',zks); 

kernmat(2,2) = kernel;
kernmat(1,1) = kern11; kernmat(2,1) = kern21;
kernmat(1,2) = kern12; kernmat(2,2) = kern22;
kern = kernel(kernmat); % interleaves into a 2x2 matrix-valued kernel

Built-in Kernels

chunkIE literalincludes built-in kernels that arise in the solution of the following PDEs.

• Laplace kernels

• Helmholtz kernels

• Stokes kernels

• Linear elasticity kernels

Some common notation used below is that \(x=(x_1,x_2)\) will refer to a “target” location and \(y=(y_1,y_2)\) will refer to a “source location. Likewise, \(n(x)\) refers to the curve normal at the target and \(n(y)\) refers to the normal at the source. The notation \(|x-y|\) denotes the distance between the source and target, i.e.

\[|x-y| = \sqrt{ (x_1-y_1)^2 + (x_2-y_2)^2 } \; .\]

Laplace kernels

The free-space Green’s function for Laplace’s equation in two dimensions, denoted by \(G_{0}(x, y)\), is given by

\[G_{0}(x,y) = -\frac{1}{2\pi} \log |x-y| \,,\]

where \(x=(x_{1},x_{2})\), and \(y=(y_{1},y_{2})\).

The kernel \(S(x,y)=G_0(x,y)\) is the integral kernel for the single layer potential. The double layer kernel is then \(D(x,y) = n(y)\cdot \nabla_y G_0(x,y)\). All Laplace integral kernels are derived from these two.

For each of the single layer, double layer, and combined layer (a linear combination of single and double layer), there are also kernels for the normal derivative (with the name “prime”), tangential derivative (with the name “tau”), and the gradient (with the name “grad”) taken at the target location. For a kernel \(K(x,y)\), these are the kernels \(n(x)\cdot \nabla_xK(x,y)\), \(\tau(x)\cdot \nabla_x K(x,y)\), and \(\nabla_x K(x,y)\), respectively.

The PDE name for using Laplace kernels can be ‘Laplace’ or ‘l’. These kernel types can be obtained with either calling sequence below

kern = kernel('l',type)
kern = kernel.lap2d(type)

The documentation of the kernel.lap2d function has details on any expected additional arguments and default values:

function obj = lap2d(type, coefs)
%KERNEL.LAP2D   Construct the Laplace kernel.
%   KERNEL.LAP2D('s') or KERNEL.LAP2D('single') constructs the single-layer
%   Laplace kernel.
%
%   KERNEL.LAP2D('d') or KERNEL.LAP2D('double') constructs the double-layer
%   Laplace kernel.
%
%   KERNEL.LAP2D('sp') or KERNEL.LAP2D('sprime') constructs the
%   normal derivative of the single-layer Laplace kernel.
%
%   KERNEL.LAP2D('sg') or KERNEL.LAP2D('sgrad') constructs the
%   gradient of the single-layer Laplace kernel.
%
%   KERNEL.LAP2D('dp') or KERNEL.LAP2D('dprime') constructs the normal
%   derivative of the double-layer Laplace kernel.
%
%   KERNEL.LAP2D('dg') or KERNEL.LAP2D('dgrad') constructs the gradient
%   of the double-layer Laplace kernel.
%
%   KERNEL.LAP2D('c', coefs) or KERNEL.LAP2D('combined', coefs) constructs
%   the combined-layer Laplace kernel with parameter coefs, i.e.,
%   coefs(1)*KERNEL.LAP2D('d') + coefs(2)*KERNEL.LAP2D('s'). If no
%   value of coefs is specified the default is coefs = [1 1]  
%
%   KERNEL.LAP2D('cp', coefs) or KERNEL.LAP2D('cprime', coefs) constructs
%   the normal derivative of the combined-layer Laplace kernel with
%   parameter coefs
%
%   KERNEL.LAP2D('cg', coefs) or KERNEL.LAP2D('cgrad', coefs) constructs
%   the gradient of the combined-layer Laplace kernel with parameter coefs
%   
%   KERNEL.LAP2D('sint') or KERNEL.LAP2D('s int') constructs the volume
%   integral of the single layer Laplace kernel
%
%   KERNEL.LAP2D('sintt') or KERNEL.LAP2D('s int trans') constructs the transpose of the 
%   volume integral of the single layer Laplace kernel
%   
%   KERNEL.LAP2D('dint') or KERNEL.LAP2D('d int') constructs the volume
%   integral of the double layer Laplace kernel
%
%   KERNEL.LAP2D('cint', coefs) or KERNEL.LAP2D('c int', coefs) constructs the volume
%   integral of the combined-layer Laplace kernel with parameter coefs, 
%   i.e. (coef(1)*  KERNEL.LAP2D('dint') + coef(2)* KERNEL.LAP2D('sint')). If no
%   value of coefs is specified the default is coefs = [1 1] 
%   
%
% See also CHNK.LAP2D.KERN.

Helmholtz kernels

The free-space Green’s function for the Helmholtz equation in two dimensions, denoted by \(G_{k}(x, y)\), is given by

\[G_{k}(x,y) = \frac{\imath}{4} H_0^{(1)}(k |x-y|) \, ,\]

where \(k\) is the wavenumber, \(x=(x_{1},x_{2})\), and \(y=(y_{1},y_{2})\).

The documentation of the kernel.helm2d function has details on any expected additional arguments and default values:

function obj = helm2d(type, zk, coefs)
%KERNEL.HELM2D   Construct the Helmholtz kernel.
%   KERNEL.HELM2D('s', ZK) or KERNEL.HELM2D('single', ZK) constructs the
%   single-layer Helmholtz kernel with wavenumber ZK.
%
%   KERNEL.HELM2D('d', ZK) or KERNEL.HELM2D('double', ZK) constructs the
%   double-layer Helmholtz kernel with wavenumber ZK.
%
%   KERNEL.HELM2D('sp', ZK) or KERNEL.HELM2D('sprime', ZK) constructs the
%   derivative of the single-layer Helmholtz kernel with wavenumber ZK.
%
%   KERNEL.HELM2D('c', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS)
%   constructs the combined-layer Helmholtz kernel with wavenumber ZK and
%   parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('d', ZK) + 
%   COEFS(2)*KERNEL.HELM2D('s', ZK).
%
%   KERNEL.HELM2D('cp', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS)
%   constructs the derivative of the combined-layer Helmholtz kernel with 
%   wavenumber ZK and  parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + 
%   COEFS(2)*KERNEL.HELM2D('sp', ZK).
%
%   KERNEL.HELM2D('all', ZK, COEFS) or KERNEL.HELM2D('tsys', ZK, COEFS) 
%   constructs the (2x2) matrix of kernels [D, S; D' S'] scaled by coefs
%   where coefs is (2x2) matrix, i.e. D part is scaled as
%   COEFS(1,1)*KERNEL.HELM2D('d', zk), and so on.
%
%   KERNEL.HELM2D('trans_rep', ZK, COEFS) or KERNEL.HELM2D('trep', ZK, COEFS) 
%   constructs the transmission repretsentation, i.e. the (1x2) matrix of
%   kernels [D, S] scaled by coefs where coefs is (1x2) matrix, i.e. D part
%   is scaled as COEFS(1)*KERNEL.HELM2D('d', zk), and so on.
%
%   KERNEL.HELM2D('trans_rep_p', ZK, COEFS) or KERNEL.HELM2D('trep_p', ZK, COEFS) 
%   constructs the derivative of the transmission repretsentation, i.e. the
%   (1x2) matrix of kernels [D', S'] scaled by coefs where coefs is (1x2)
%   matrix, i.e. D part is scaled as COEFS(1)*KERNEL.HELM2D('dp', zk), and
%   so on.
%
%   KERNEL.HELM2D('c2trans', ZK, COEFS) or KERNEL.HELM2D('c2t', ZK, COEFS) 
%   evaluates the combined-layer Helmholtz kernel and its derivative, i.e.
%   the (2x1) matrix of kernels [C; C'] scaled by coefs where coefs is
%   (2x2) matrix, i.e. kernel returns 
%   [COEFS(1,1)*KERNEL.HELM2D('d',zk)+COEFS(1,2)*KERNEL.HELM2D('s', zk); 
%   COEFS(2,1)*KERNEL.HELM2D('dp',zk)+COEFS(2,2)*KERNEL.HELM2D('sp', zk)]
% See also CHNK.HELM2D.KERN.

Stokes kernels

The free-space Green’s function for the Stokes equation in two dimensions, known as a Stokeslet and denoted by \({ G}(x, y)\), is a \(2\times 2\) tensor whose \(i,j\) entry is

\[G_{ij}(x,y) = \frac{1}{4\pi \mu} \left ( \frac{(x_i-y_i)(x_j-y_j)}{|x-y|^2} - \delta_{ij} \log |x-y| \right ) \, ,\]

where \(\mu\) is the dynamic viscosity, \(x=(x_{1},x_{2})\), and \(y=(y_{1},y_{2})\). The quantity \({ G} \cdot { f}\) gives the velocity induced by a point charge at \(y\) with strength and orientation determined by \(f = (f_1,f_2)\). The corresponding pressure is \({P}\cdot { f}\) where

\[P_{i}(x,y) = \frac{x_i-y_i}{2\pi |x-y|^2} \, .\]

Let \({ \sigma}({ u},p)\) denote the Cauchy stress tensor for given velocity and pressure fields, i.e.

\[{ \sigma}({ u},p) = 2 \mu { \epsilon} + p I = \mu ( \nabla { u} + \nabla { u}^T) + p I \; ,\]

where \(\epsilon\) is the strain rate.

On a curve in the fluid, the corresponding traction vector is \(t(x) = n(x) \cdot \sigma(u,p)\). Specifying the traction is a common boundary condition for Stokes problems. Kernels corresponding to the traction of some base kernel are denoted by the suffix “trac”, e.g. strac is the name for the traction of the Stokeslet. For reference, the 3-tensor corresponding to the stress of each column of the Stokeslet, denoted by \(T\) and known as a “stresslet”, has the entries

\[T_{ijk}(x,y) = -\frac{(x_i-y_i)(x_j-y_j)(x_k-y_k)}{\pi |x-y|^4} \; ,\]

so that the traction of a Stokeslet is \(n(x) \cdot T(x,y) \cdot f(y)\).

The Stokes double layer kernel is defined in terms of the stresslet as follows

\[D_{ij}(x,y) = -T_{jki}(x,y) n_k(y) \; .\]

function obj = stok2d(type, mu, coefs)
%KERNEL.STOK2D   Construct the Stokes kernel.
%   KERNEL.STOK2D('svel', MU) or KERNEL.ELAST2D('svelocity', MU)
%   constructs the single-layer Stokes kernel for velocity with viscosity
%   MU. KERNEL.STOK2D('s', MU) and KERNEL.STOK2D('single', MU) are
%   equivalent.
%
%   KERNEL.STOK2D('spres', MU) or KERNEL.ELAST2D('spressure', MU)
%   constructs the single-layer Stokes kernel for pressure with viscosity
%   MU.
%
%   KERNEL.STOK2D('strac', MU) or KERNEL.ELAST2D('straction', MU)
%   constructs the single-layer Stokes kernel for traction with viscosity
%   MU.
%
%   KERNEL.STOK2D('sgrad', MU) or KERNEL.ELAST2D('sgradient', MU)
%   constructs the single-layer Stokes kernel for velocity gradient with 
%   viscosity MU.
%
%   KERNEL.STOK2D('dvel', MU) or KERNEL.ELAST2D('dvelocity', MU)
%   constructs the double-layer Stokes kernel for velocity with viscosity
%   MU. KERNEL.STOK2D('d', MU) and KERNEL.STOK2D('double', MU) are
%   equivalent.
%
%   KERNEL.STOK2D('dpres', MU) or KERNEL.ELAST2D('dpressure', MU)
%   constructs the double-layer Stokes kernel for pressure with viscosity
%   MU.
%
%   KERNEL.STOK2D('dtrac', MU) or KERNEL.ELAST2D('dtraction', MU)
%   constructs the double-layer Stokes kernel for traction with viscosity
%   MU.
%
%   KERNEL.STOK2D('dgrad', MU) or KERNEL.ELAST2D('dgradient', MU)
%   constructs the double-layer Stokes kernel for velocity gradient with 
%   viscosity MU.
%  
%   KERNEL.STOK2D('cvel', MU) or KERNEL.ELAST2D('cvelocity', MU)
%   constructs the combined-layer Stokes kernel for velocity with viscosity
%   MU. KERNEL.STOK2D('c', MU) and KERNEL.STOK2D('comb', MU) are
%   equivalent.
%
%   KERNEL.STOK2D('cpres', MU) or KERNEL.ELAST2D('cpressure', MU)
%   constructs the combined-layer Stokes kernel for pressure with viscosity
%   MU.
%
%   KERNEL.STOK2D('ctrac', MU) or KERNEL.ELAST2D('ctraction', MU)
%   constructs the combined-layer Stokes kernel for traction with viscosity
%   MU.
%
%   KERNEL.STOK2D('cgrad', MU) or KERNEL.ELAST2D('cgradient', MU)
%   constructs the combined-layer Stokes kernel for velocity gradient with 
%   viscosity MU.
%
% See also CHNK.STOK2D.KERN.

Linear elasticity kernels

The Lamé parameters \(\lambda\) and \(\mu\) are used to specify the stress-strain relations in linear elasticity, i.e.

\[\sigma = \lambda \operatorname{Trace} (\epsilon) I + 2\mu \epsilon \; .\]

The governing equations for a linear elastic solid, in the absence of body forces, are \(\nabla \cdot \sigma = 0\). The displacement induced by a point force with strength and orientation given by \(f\) is \(G \cdot f\), where

\[G_{ij}(x,y) = \frac{1}{4\pi \mu(\lambda + 2\mu)} \left ( (\lambda + 3\mu) \delta_{ij} \log|x-y| - (\lambda + \mu) \frac{(x_i-y_i)(x_j-y_j)}{|x-y|^2} \right ) \; .\]

The standard traction and double layer operators are analogous to the Stokes case.

function obj = elast2d(type, lam, mu)
%KERNEL.ELAST2D   Construct the elasticity kernel.
%   KERNEL.ELAST2D('s', LAM, MU) or KERNEL.ELAST2D('single', LAM, MU)
%   constructs the single-layer elasticity kernel with Lamé parameters LAM
%   and MU.
%
%   KERNEL.ELAST2D('strac', LAM, MU) or KERNEL.ELAST2D('straction', LAM,
%   MU) constructs the single-layer elasticity kernel for traction with
%   Lamé parameters LAM and MU.
%
%   KERNEL.ELAST2D('sgrad', LAM, MU) or KERNEL.ELAST2D('straction', LAM,
%   MU) constructs the gradient of the single-layer elasticity kernel with
%   Lamé parameters LAM and MU.
%
%   KERNEL.ELAST2D('d', LAM, MU) or KERNEL.ELAST2D('double', LAM, MU)
%   constructs the double-layer elasticity kernel with Lamé parameters LAM
%   and MU.
%
%   KERNEL.ELAST2D('dalt', LAM, MU) constructs the alternative double-layer
%   elasticity kernel with Lamé parameters LAM and MU.
%
%   KERNEL.ELAST2D('dalttrac', LAM, MU) or 
%   KERNEL.ELAST2D('dalttraction', LAM, MU)
%   constructs the traction of the alternative double-layer elasticity 
%   kernel with Lamé parameters LAM and MU.
%
%   KERNEL.ELAST2D('daltgrad', LAM, MU) constructs the gradient of the 
%   alternative double-layer elasticity kernel with Lamé parameters LAM 
%   and MU.
%
% See also CHNK.ELAST2D.KERN