Add the integral of lap2d single and double-layer potential

chunkie/+chnk/+lap2d/kern.m

@@ -151,9 +151,53 @@
     submat = coef(1)*reshape(permute(submat,[3,1,2]),2*nt,ns);
     submat = submat+coef(2)*reshape(permute(grad,[3,1,2]),2*nt,ns);
 
+
+
+% integral of the single layer
+case{'sint'}
+    %srcnorm = chnk.normal2d(srcinfo);
+    [s] = chnk.lap2d.green(src,targ);
+    rx = bsxfun(@minus,targ(1,:).',src(1,:));
+    ry = bsxfun(@minus,targ(2,:).',src(2,:));
+    nrt = rx.*targinfo.n(1,:).' + ry.*targinfo.n(2,:).';
+    submat = nrt.*(s/2+1/(8*pi));
+
+% transpose of the previous one.
+case{'sintt'}
+    %srcnorm = chnk.normal2d(srcinfo);
+    [s] = chnk.lap2d.green(src,targ);
+    rx = bsxfun(@minus,targ(1,:).',src(1,:));
+    ry = bsxfun(@minus,targ(2,:).',src(2,:));
+    nrt = rx.*srcinfo.n(1,:) + ry.*srcinfo.n(2,:);
+    submat = -nrt.*(s/2+1/(8*pi));
+
+% integral of the double layer
+case{'dint'}
+    %srcnorm = chnk.normal2d(srcinfo);
+    [s] = chnk.lap2d.green(src,targ);
+    ntsx = bsxfun(@times,targinfo.n(1,:).',srcinfo.n(1,:));
+    ntsy = bsxfun(@times,targinfo.n(2,:).',srcinfo.n(2,:));
+    nts = ntsx+ntsy;
+    submat = -nts.*s;
+% integral of the combined field (coef(1)* integral D + coef(2)*integral S)
+case{'cint'}
+    coef = ones(2,1);
+    if(nargin == 4); coef = varargin{1}; end
+    %srcnorm = chnk.normal2d(srcinfo);
+    [s] = chnk.lap2d.green(src,targ);
+    ntsx = bsxfun(@times,targinfo.n(1,:).',srcinfo.n(1,:));
+    ntsy = bsxfun(@times,targinfo.n(2,:).',srcinfo.n(2,:));
+    nts = ntsx+ntsy;
+    submat = -nts.*s;
+
+    rx = bsxfun(@minus,targ(1,:).',src(1,:));
+    ry = bsxfun(@minus,targ(2,:).',src(2,:));
+    nrt = rx.*targinfo.n(1,:).' + ry.*targinfo.n(2,:).';
+    submat = coef(1)*submat+coef(2)*(nrt.*(s/2+1/(8*pi)));
+
+
 otherwise
     error('Unknown Laplace kernel type ''%s''.', type);
 end
 
 end
-

chunkie/@kernel/lap2d.m

@@ -29,6 +29,21 @@
 %
 %   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.
 
@@ -133,6 +148,35 @@
         obj.sing = 'hs';
         obj.opdims = [2,1];
 
+    case {'sint', 's int'}
+        obj.type = 'sint';
+        obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'sint');
+        obj.sing = 'log';
+        obj.opdims = [1,1];
+
+    case {'sintt', 's int trans'}
+        obj.type = 'sint';
+        obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'sintt');
+        obj.sing = 'log';
+        obj.opdims = [1,1];
+
+    case {'dint', 'd int'}
+        obj.type = 'dint';
+        obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'dint');
+        obj.sing = 'log';
+        obj.opdims = [1,1];
+
+    case {'cint', 'c int'}
+        if ( nargin < 2 )
+            warning('Missing combined layer parameter coefs. Defaulting to [1 1].');
+            coefs = ones(2,1);
+        end
+        obj.type = 'cint';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'cint',coefs);
+        obj.sing = 'log';
+        obj.opdims = [1,1];
+
     otherwise
         error('Unknown Laplace kernel type ''%s''.', type);
 

devtools/test/dint_kernel_integration_ellipse_Test.m

@@ -0,0 +1,60 @@
+dint_kernel_integration_ellipse_Test0();
+% test the volume integral of the double-layer in an ellipse
+%
+
+function dint_kernel_integration_ellipse_Test0()
+
+cparams = [];
+
+cparams.eps = 1e-6;
+cparams.nover = 0;
+cparams.maxchunklen = 0.5;      
+
+
+
+chnkr = chunkerfunc(@(t) ellipse(t,4,1), cparams);
+
+figure(1)                                                   % plot the chunker-object 
+clf
+plot(chnkr, '-x')
+hold on
+quiver(chnkr)
+axis equal
+
+
+fkern = @(s,t) chnk.lap2d.kern(s,t,'d');            % double layer
+
+opts = [];
+opts.sing = 'log';
+
+mat = chunkermat(chnkr, fkern, opts);
+
+lhs = -(1/2).*eye(chnkr.npt) + mat   ;
+
+nx = chnkr.n(1,:); ny = chnkr.n(2,:);
+
+rhs = chnkr.r(1,:).^2 - chnkr.r(2,:).^2;
+
+
+rhs = rhs(:);
+
+rho = lhs\rhs;
+
+
+kern_integral = @(s,t) chnk.lap2d.kern(s, t, 'dint');      % double-layer integration kernel
+
+kern_mat = chunkermat(chnkr, kern_integral, opts);
+integrand = kern_mat*rho;
+
+numeric_result = chunkerintegral(chnkr, integrand);
+analytic_result = 15*pi;
+
+abserr = abs(analytic_result - numeric_result);                       % absolute error
+relerr = abs((numeric_result -analytic_result)/analytic_result);      % rel error
+
+fprintf('absolute error %5.2e\n',abserr);
+fprintf('relative error %5.2e\n',relerr);
+
+assert(abserr<1e-9)
+assert(relerr<1e-9)
+end

devtools/test/sint_kernel_integration_ellipse_Test.m

@@ -0,0 +1,69 @@
+sint_kernel_integration_ellipse_Test0();
+% test the volume integral of the single-layer in an ellipse
+% 
+function sint_kernel_integration_ellipse_Test0()
+cparams = [];
+
+cparams.eps = 1e-6;
+cparams.nover = 0;
+cparams.maxchunklen = 0.5;      
+
+
+
+chnkr = chunkerfunc(@(t) ellipse(t,2,1), cparams);
+
+figure(1)                                                   % plot the chunker-object 
+clf
+plot(chnkr, '-x')
+hold on
+quiver(chnkr)
+axis equal
+
+
+fkern = @(s,t) chnk.lap2d.kern(s,t,'sprime');
+
+opts = [];
+opts.sing = 'log';
+
+mat = chunkermat(chnkr, fkern, opts);
+
+lhs = (1/2).*eye(chnkr.npt) + mat + onesmat(chnkr)  ;
+
+nx = chnkr.n(1,:); ny = chnkr.n(2,:);
+
+rhs = 2.*chnkr.r(1,:).*nx - 2.*chnkr.r(2,:).*ny;
+
+rhs = rhs(:);
+
+rho = lhs\rhs;
+
+eval_kern = @(s,t) chnk.lap2d.kern(s,t, 's');
+
+sol = chunkerkerneval(chnkr, eval_kern, rho,[0;0]);
+
+
+true_sol = 0;                    
+
+uerr = sol - true_sol;                         % figuring out the extra constant in the interior Neumann problem
+
+
+
+kern_integral = @(s,t) chnk.lap2d.kern(s, t, 'sint');              % single layer integration kernel
+
+kern_mat = chunkermat(chnkr, kern_integral, opts);
+integrand = kern_mat*rho;
+
+numeric_result = chunkerintegral(chnkr, integrand);
+analytic_result = 3*pi/2 + uerr*2*pi;
+
+abserr = abs(analytic_result - numeric_result);                            % abs error
+
+relerr = abs((numeric_result -analytic_result)/analytic_result);   % rel error
+
+fprintf('absolute error %5.2e\n',abserr);
+fprintf('relative error %5.2e\n',relerr);
+
+assert(abserr<1e-9)
+assert(relerr<1e-9)
+
+end