diff --git a/chunkie/+chnk/+lap2d/kern.m b/chunkie/+chnk/+lap2d/kern.m index 248e32f..83fcf40 100644 --- a/chunkie/+chnk/+lap2d/kern.m +++ b/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 - diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m index 62a95df..3eedbb3 100644 --- a/chunkie/@kernel/lap2d.m +++ b/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); diff --git a/devtools/test/dint_kernel_integration_ellipse_Test.m b/devtools/test/dint_kernel_integration_ellipse_Test.m new file mode 100644 index 0000000..b1b33b8 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/devtools/test/sint_kernel_integration_ellipse_Test.m b/devtools/test/sint_kernel_integration_ellipse_Test.m new file mode 100644 index 0000000..c336402 --- /dev/null +++ b/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