diff --git a/chunkie/+chnk/+helm2d/transmission_helper.m b/chunkie/+chnk/+helm2d/transmission_helper.m index 1cc20e5b..cbc6a4e0 100644 --- a/chunkie/+chnk/+helm2d/transmission_helper.m +++ b/chunkie/+chnk/+helm2d/transmission_helper.m @@ -124,24 +124,31 @@ c1 = coefs(cs(1,i)); c2 = coefs(cs(2,i)); - cc1(1,1) = -alpha1(i)*c1; - cc2(1,1) = -alpha1(i)*c2; - - cc1(1,2) = -alpha1(i); - cc2(1,2) = -alpha1(i); - - cc1(2,1) = alpha2(i); - cc2(2,1) = alpha2(i); + cc1(1,1) = alpha1(i)*c1; + cc1(1,2) = alpha1(i); + cc1(2,1) = -alpha2(i); + cc1(2,2) = -alpha2(i)/c1; + + cc2(1,1) = alpha1(i)*c2; + cc2(1,2) = alpha1(i); + cc2(2,1) = -alpha2(i); + cc2(2,2) = -alpha2(i)/c2; - cc1(2,2) = alpha2(i)/c1; - cc2(2,2) = alpha2(i)/c2; + % Build the correct system kernel cc1*K_{k1} - cc2*K_{k2}. + % Using helmdiff naively with cat(3,cc1,cc2) gives + % cc1*(K_{k1}-K_Lap) - cc2*(K_{k2}-K_Lap) + % = cc1*K_{k1} - cc2*K_{k2} - (cc1-cc2)*K_Lap, + % missing the Laplace correction when cc1~=cc2 (i.e. c1~=c2). + % Instead use the full Helmholtz 'all' kernel directly. + k1_cap = k1(i); k2_cap = k2(i); + cc1_cap = cc1; cc2_cap = cc2; for j=1:ncurve - zks = [k1(i), k2(i)]; - cc_use = cat(3, -cc1, -cc2); - kerns(i,j) = kernel('helmdiff', 'all', zks, cc_use); -% kerns{i,j} = @(s,t) -(chnk.helm2d.kern(k1(i),s,t,'all',cc1)- ... -% chnk.helm2d.kern(k2(i),s,t,'all',cc2)); - end + kerns(i,j) = kernel(); + kerns(i,j).eval = @(s,t) chnk.helm2d.kern(k1_cap, s, t, 'all', cc1_cap) - ... + chnk.helm2d.kern(k2_cap, s, t, 'all', cc2_cap); + kerns(i,j).sing = 'log'; + kerns(i,j).opdims = [2,2]; + end end diff --git a/chunkie/@chunker/refine.m b/chunkie/@chunker/refine.m index 8d9f501d..b9e1c651 100644 --- a/chunkie/@chunker/refine.m +++ b/chunkie/@chunker/refine.m @@ -107,7 +107,13 @@ error('too many chunks') end - chnkr = split(chnkr,ii,[],x,w,u,stype); + + if max(abs(imag(chnkr.r(:,:,ii))), [], 'all') > 1e-16 + chnkr = split(chnkr,ii,[],x,w,u,'t'); + else + chnkr = split(chnkr,ii,[],x,w,u,stype); + end + % update chunklens @@ -145,7 +151,13 @@ error('too many chunks') end - chnkr = split(chnkr,i,[],x,w,u,stype); + + if max(abs(imag(chnkr.r(:,:,i))), [], 'all') > 1e-16 + chnkr = split(chnkr,i,[],x,w,u,'t'); + else + chnkr = split(chnkr,i,[],x,w,u,stype); + end + % update chunklens @@ -215,7 +227,13 @@ error('too many chunks') end - chnkr = split(chnkr,i,[],x,w,u,stype); + + if max(abs(imag(chnkr.r(:,:,i))), [], 'all') > 1e-16 + chnkr = split(chnkr,i,[],x,w,u,'t'); + else + chnkr = split(chnkr,i,[],x,w,u,stype); + end + % update chunklens @@ -254,7 +272,12 @@ error('CHUNKER.REFINE nchmax=%d exceeded during oversample',nchmax) end - chnkr = split(chnkr,i,[],x,w,u,stype); + + if max(abs(imag(chnkr.r(:,:,i))), [], 'all') > 1e-16 + chnkr = split(chnkr,i,[],x,w,u,'t'); + else + chnkr = split(chnkr,i,[],x,w,u,stype); + end % update chunklens diff --git a/devtools/test/chunkgrphrcipTransmissionTest.m b/devtools/test/chunkgrphrcipTransmissionTest.m index 04da2442..f9557816 100644 --- a/devtools/test/chunkgrphrcipTransmissionTest.m +++ b/devtools/test/chunkgrphrcipTransmissionTest.m @@ -59,7 +59,7 @@ function chunkgrphrcipTransmissionTest0() nregions = 2; ks = [1.1;2.1]*10; -coefs = [1.0;1.0]; +coefs = [1.1;1.0]; cs(1,1:ncurve) = 1; cs(2,1:ncurve) = 2; opts = [];