diff --git a/chunkie/+chnk/+axissymhelm2d/kern.m b/chunkie/+chnk/+axissymhelm2d/kern.m index 1bd2ea8..e469444 100644 --- a/chunkie/+chnk/+axissymhelm2d/kern.m +++ b/chunkie/+chnk/+axissymhelm2d/kern.m @@ -40,7 +40,7 @@ % type == 'sprime', normal derivative of single % layer S' % type == 'dprime', normal derivative of double layer D' -% type == 'ddiff', D_{k}' - D_{ik}', for this routine k must be real +% type == 'dprimediff', D_{k}' - D_{ik}', for this routine k must be real % type == 'c', combined layer kernel coef(1) D + coef(2) S % % varargin{1} - coef: length 2 array in the combined layer diff --git a/chunkie/+chnk/+flam/kernbyindex.m b/chunkie/+chnk/+flam/kernbyindex.m index b431397..b6b899b 100644 --- a/chunkie/+chnk/+flam/kernbyindex.m +++ b/chunkie/+chnk/+flam/kernbyindex.m @@ -37,8 +37,12 @@ % l2scale - boolean type that determines if we should % rescale the matrix by l2scale. the default value is false. -if isa(kern,'kernel') - kern = kern.eval; +if ~isa(kern,'kernel') + try + kern = kernel(kern); + catch + error('KERNBYINDEX: fourth input kern not of supported type'); + end end if nargin < 7 @@ -145,9 +149,9 @@ if size(kern) == 1 - matuni = kern(srcinfo,targinfo); + matuni = kern.eval(srcinfo,targinfo); else - matuni = kern{itrg,isrc}(srcinfo,targinfo); + matuni = kern(itrg,isrc).eval(srcinfo,targinfo); end % scale matrix by weights diff --git a/chunkie/+chnk/+flam/proxyfun.m b/chunkie/+chnk/+flam/proxyfun.m index de9a9aa..38f6b6c 100644 --- a/chunkie/+chnk/+flam/proxyfun.m +++ b/chunkie/+chnk/+flam/proxyfun.m @@ -41,13 +41,14 @@ % rescale the matrix by l2scale. the default value is false. % -if isa(kern,'kernel') - kern = kern.eval; +if ~isa(kern,'kernel') + try + kern = kernel(kern); + catch + error('PROXYFUN: sixth input kern not of supported type'); + end end -if nargin < 13 - l2scale = false; -end % scaled proxy points and weights (no scaling necessary on tangents) @@ -153,9 +154,9 @@ opdims_trans = opdims_mat(:,j,i); if length(kern) == 1 - matuni = kern(srcinfo,targinfo); + matuni = kern.eval(srcinfo,targinfo); else - matuni = kern{i,j}(srcinfo,targinfo); + matuni = kern(i,j).eval(srcinfo,targinfo); end wsrc = wts(juni); @@ -180,9 +181,9 @@ if ifaddtrans if length(kern) == 1 - matuni2 = kern(targinfo,srcinfo); + matuni2 = kern.eval(targinfo,srcinfo); else - matuni2 = kern{j,i}(targinfo,srcinfo); + matuni2 = kern(j,i).eval(targinfo,srcinfo); end wsrc = pw(:); diff --git a/chunkie/+chnk/chunkerkerneval_smooth.m b/chunkie/+chnk/chunkerkerneval_smooth.m index c3c83e8..4d1fba4 100644 --- a/chunkie/+chnk/chunkerkerneval_smooth.m +++ b/chunkie/+chnk/chunkerkerneval_smooth.m @@ -87,7 +87,7 @@ % nothing to ignore for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); - dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); + dsdtdt = sqrt(sum(chnkr.d(:,:,i).^2,1)); dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); @@ -108,7 +108,7 @@ % ignore interactions in flag array for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); - dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); + dsdtdt = sqrt(sum(chnkr.d(:,:,i).^2,1)); dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); @@ -183,7 +183,7 @@ if ~isempty(flag) for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); - dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); + dsdtdt = sqrt(sum(chnkr.d(:,:,i).^2,1)); dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); diff --git a/chunkie/chunkerflam.m b/chunkie/chunkerflam.m index 9cae0eb..e693e97 100644 --- a/chunkie/chunkerflam.m +++ b/chunkie/chunkerflam.m @@ -143,8 +143,8 @@ error(msg) end -for chnkr = chnkrs - if or(chnkr.nch < 1,chnkr.k < 1) +for i=1:length(chnkrs) + if chnkrs(i).nch < 1 || chnkrs(i).k < 1 warning('empty chunker, doing nothing') return end diff --git a/devtools/test/chunkgrphrcipTransmissionTest.m b/devtools/test/chunkgrphrcipTransmissionTest.m index efb2a80..5bae51c 100644 --- a/devtools/test/chunkgrphrcipTransmissionTest.m +++ b/devtools/test/chunkgrphrcipTransmissionTest.m @@ -48,12 +48,15 @@ ncurve = size(edge2verts,1); fchnks = cell(1,ncurve); +cparams = cell(ncurve,1); for icurve = 1:ncurve fchnks{icurve} = @(t) circulararc(t,cpars{1}); fchnks{icurve} = @(t) sinearc(t,amp,frq); + cparams{icurve}.ta = 0; + cparams{icurve}.tb = 1; end -[cgrph] = chunkgraph(verts,edge2verts,fchnks); +[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams); vstruc = procverts(cgrph); @@ -63,7 +66,7 @@ quiver(cgrph); nregions = 2; -ks = [1.1;2.1]*30; +ks = [1.1;2.1]*10; coefs = [1.0;1.0]; cs(1,1:ncurve) = 1; cs(2,1:ncurve) = 2; diff --git a/devtools/test/complexificationTest.m b/devtools/test/complexificationTest.m new file mode 100644 index 0000000..0daf284 --- /dev/null +++ b/devtools/test/complexificationTest.m @@ -0,0 +1,95 @@ +% test solve_sommerfeld_dens and eval_sommerfeld_dens +% and compare to complexification solution + +%%%% +% +% . . . testing a point charge below +% +%%%% + +eps = 1E-16; +zks = pi*[1,1.3]; +zk1 = zks(1); +zk2 = zks(2); + +close all + +a = 3; +b = 1/a/2; +t0 =-5; +t1 = 5; + +f = @(t) flat_interface(t, a, b, t0, t1); +nch = 50; +xrad = -log(eps)/abs(real(zk1)) + max([abs(t0),abs(t1)]); +xmin = -xrad; +xmax = xrad; +cparams = []; +cparams.ta = xmin; +cparams.tb = xmax; +cparams.ifclosed = 0; +chnkr = chunkerfuncuni(f, nch, cparams); + + +ddiff = kernel('helmdiff', 'd', zks); +sdiff = (-1)*kernel('helmdiff', 's', zks); +dpdiff = kernel('helmdiff', 'dp', zks); +spdiff = (-1)*kernel('helmdiff', 'sp', zks); + +K = kernel([ddiff, sdiff; ... + dpdiff, spdiff]); + +n = 2*chnkr.npt; +sysmat = chunkermat(chnkr, K); +sysmat = sysmat - eye(n); + +%% Analytic solution test + +r = chnkr.r; +r0 = [0.1; -3]; +s = []; +s.r = r0; + +sk1 = kernel('helm', 's', zk1); +dk1 = kernel('helm', 'd', zk1); +skp1 = kernel('helm', 'sp', zk1); +rhs = complex(zeros(n,1)); +rhs(1:2:end) = sk1.eval(s, chnkr); +rhs(2:2:end) = skp1.eval(s, chnkr); + +dens = sysmat \ rhs; + +rt = [0.3; 2]; +targ = []; +targ.r = rt; +Keval = kernel([dk1, (-1)*sk1]); + +opts = []; +opts.forcesmooth = true; +opts.accel = false; +pot = chunkerkerneval(chnkr, Keval, dens, targ, opts); + +uex = sk1.eval(s, targ); +assert(norm(uex - pot) < 1e-10); + + + +function [f,fd,fdd] = flat_interface(t, a, b, t0, t1) + + phi = @(t,u,v,z) u*(t-z).*erfc(u*(t-z))*v - exp(-u^2*(t-z).^2)/sqrt(pi)*v; + phid = @(t,u,v,z) u*erfc(u*(t-z))*v; + phidd = @(t,u,v,z) -u*u*exp(-u^2*(t-z).^2)*2*v/sqrt(pi); + f = zeros([2,size(t)]); + fd = zeros([2,size(t)]); + fdd = zeros([2,size(t)]); + + f(1,:) = t + 1i*(phi(t,a,b,t0) - phi(t,-a,b,t1)); + fd(1,:)= 1 + 1i*(phid(t,a,b,t0) - phid(t,-a,b,t1)); + fdd(1,:) = 1i*(phidd(t,a,b,t0) - phidd(t,-a,b,t1)); + + f(2,:) = 0; + fd(2,:) = 0; + fdd(2,:) = 0; + +end + diff --git a/devtools/test/flamopdimsTest.m b/devtools/test/flamopdimsTest.m new file mode 100644 index 0000000..74f1760 --- /dev/null +++ b/devtools/test/flamopdimsTest.m @@ -0,0 +1,184 @@ + +zk1 = 10; % exterior wave number +zk2 = 15; % coating wave number + + +zkuse = max(real([zk1, zk2])); +cparams = []; +cparams.eps = 1.0e-10; +cparams.nover = 0; +cparams.maxchunklen = 4.0/zkuse; % setting a chunk length helps when the + % frequency is known + +pref = []; +pref.k = 16; +narms = 5; +amp = 0.15; + +chnkr1 = chunkerfunc(@(t) starfish(t,narms,amp, [], [], 1.15),cparams,pref); + +narms2 = 10; +amp2 = 0.03; +chnkr2 = chunkerfunc(@(t) starfish(t, narms2, amp2, [], [], 0.9), cparams, pref); + + +ifflam = true; +opts = []; +opts.ifflam = ifflam; + +%% Get matrix +tic, A = get_linop(chnkr1, chnkr2, [zk1, zk2], opts); toc + +%% Test analytic solution +src1 = zeros(2,2); +src1(1,1) = 0.01; +src1(2,1) = 0.03; + +src1(1,2) = -1; +src1(2,2) = -0.5; + +src2 = zeros(2,2); +src2(:,1) = src1(:,1); +src2(1,2) = 5.3; +src2(2,2) = -3.1; + +err1 = test_analytic_soln([zk1, zk2], chnkr1, chnkr2, src1, src2, A, ifflam); +assert(err1 < 1e-8); + + + +function err1 = test_analytic_soln(zks, chnkr1, chnkr2, src1 ,src2, A, ifflam) + zk1 = zks(1); + zk2 = zks(2); + + dk1 = kernel('helm', 'd', zk1); + sk1 = kernel('helm', 's', zk1); + sk1p = kernel('helm', 'sp', zk1); + + dk2 = kernel('helm', 'd', zk2); + sk2 = kernel('helm', 's', zk2); + sk2p = kernel('helm', 'sp', zk2); + ck2 = kernel('helm', 'c', zk2, [2, -2*1j*zk2]); + + [~, n1] = size(src1); + strengths1 = randn(n1,1) + 1j*randn(n1,1); + + srcinfo1 = []; + srcinfo1.r = src1; + + targinfo1 = []; + targinfo1.r = chnkr1.r(:,:); + targinfo1.n = chnkr1.n(:,:); + + targinfo2 = []; + targinfo2.r = chnkr2.r(:,:); + targinfo2.n = chnkr2.n(:,:); + + u1 = sk1.eval(srcinfo1, targinfo1)*strengths1; + dudn1 = sk1p.eval(srcinfo1, targinfo1)*strengths1; + + + srcinfo2 = []; + srcinfo2.r = src2; + [~, n2] = size(src2); + + strengths2 = randn(n2,1) + 1j*randn(n2,1); + + u2_1 = sk2.eval(srcinfo2, targinfo1)*strengths2; + dudn2_1 = sk2p.eval(srcinfo2, targinfo1)*strengths2; + + u2_2 = sk2.eval(srcinfo2, targinfo2)*strengths2; + + ntot = chnkr1.npt*2 + chnkr2.npt; + n = 2*chnkr1.npt; + rhs = zeros(ntot,1); + rhs(1:2:n) = u1 - u2_1; + rhs(2:2:n) = dudn1 - dudn2_1; + rhs((n+1):end) = u2_2; + + if ~ifflam + sol = A\rhs; + else + sol = rskelf_sv(A, rhs); + end + +% Now test solution at interior and exterior point + + targ1 = src2(:,2); + targ2 = src1(:,2); + targinfo1 = []; + targinfo1.r = targ1; + + targinfo2 = []; + targinfo2.r = targ2; + + uex1 = sk1.eval(srcinfo1, targinfo1)*strengths1; + uex2 = sk2.eval(srcinfo2, targinfo2)*strengths2; + + + K1 = kernel([dk1, -1*sk1]); + u1 = chunkerkerneval(chnkr1, K1, sol(1:n), targinfo1); + + K2 = kernel([dk2, -1*sk2]); + u2 = chunkerkerneval(chnkr1, K2, sol(1:n), targinfo2); + u2 = u2 + chunkerkerneval(chnkr2, ck2, sol((n+1):end), targinfo2); + err1 = norm(u1 - uex1) + norm(u2 - uex2); + + fprintf('error in exterior=%d\n', norm(u1-uex1)); + fprintf('error in coating=%d\n', norm(u2-uex2)); +end + +function [A] = get_linop(chnkr1, chnkr2, zks, opts) + zk1 = zks(1); + zk2 = zks(2); + if nargin < 4 + opts = []; + end + ifflam = false; + if isfield(opts, 'ifflam') + ifflam = opts.ifflam; + end + chnkrs(2,1) = chunker(); + chnkrs(1,1) = chnkr1; + chnkrs(2,1) = chnkr2; + +%% Define kernels + skdiff = kernel('helmdiff', 's', [zk1, zk2]); + skpdiff = kernel('helmdiff', 'sp', [zk1, zk2]); + dkdiff = kernel('helmdiff', 'd', [zk1, zk2]); + dkpdiff = kernel('helmdiff', 'dp', [zk1, zk2]); + + K = kernel([dkdiff, -1*skdiff; ... + dkpdiff, -1*skpdiff]); + + + dk2 = kernel('helm', 'd', zk2); + sk2 = kernel('helm', 's', zk2); + + ck2 = kernel('helm', 'c', zk2, [2, -2*1j*zk2]); + ck2p = kernel('helm', 'cp', zk2, [2, -2*1j*zk2]); + + K2 = kernel([dk2, -1*sk2]); + K3 = -1*kernel([ck2; ck2p]); + + Kmat(2,2) = kernel(); + Kmat(1,1) = K; + Kmat(2,1) = K2; + Kmat(2,2) = ck2; + Kmat(1,2) = K3; + +%% Evaluate matrix + +opts_loc = []; +opts_loc.adaptive_correction = true; + +if ~ifflam + A = chunkermat(chnkrs, Kmat, opts_loc); + ntot = size(A); + A = A + eye(ntot); +else + + A = chunkerflam(chnkrs, Kmat, 1, opts_loc); +end + +end