TBIEkernels.m

function [ Skij Dkij] = TBIEkernels(i,j, k, r, dr, drdn, normal )
% calculate the kernels for the Traction Boundary Integral Equation 

global const3 const4 shearMod mu

Skij = shearMod / (2 * pi * (1 - mu) ) * 1 / ( r^2 );
Skij = Skij * ( 2 * drdn * ( const4 * (i==j) * dr(k) + ...
            mu * ( dr(j) * (k==i) + dr(i) * (j==k) ) -...
            4 * dr(i) * dr(j) * dr(k) ) +...
            2 * mu * ( normal(i) * dr(j) * dr(k) +...
            normal(j) * dr(i) * dr(k) ) +...
            const4 * ( 2 * normal(k) * dr(i) * dr(j) + normal(j) * (k==i) +...
            normal(i) * (j==k) )...
            - ( 1 - 4 * mu ) * normal(k) * (i==j) );

Dkij = const3 * (1 / r) * ( const4 * ( -dr(k) * (i==j) + dr(j) * (k==i) + dr(i) * (j==k) )...
					+ 2 * dr(i) * dr(j) * dr(k) );
    
end