diff --git a/.gitignore b/.gitignore index d572bd6..3138d41 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,6 @@ chunkie/FLAM chunkie/chunkergraph/test_chnkgraph_local.m chunkie/chunkergraph/local_tests/test_chnkgraph_local.m chunkie/chunkergraph/local_tests/test_chnkgraph_local2.m + +*k2test* +*ktest* diff --git a/chunkie/+chnk/+axissymhelm2d/asym_helm_data.mat b/chunkie/+chnk/+axissymhelm2d/asym_helm_data.mat new file mode 100644 index 0000000..cebb0a8 Binary files /dev/null and b/chunkie/+chnk/+axissymhelm2d/asym_helm_data.mat differ diff --git a/chunkie/+chnk/+axissymhelm2d/dalph.m b/chunkie/+chnk/+axissymhelm2d/dalph.m new file mode 100644 index 0000000..c268ed7 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/dalph.m @@ -0,0 +1,48 @@ +function [dout] = dalph(r,dr,z) + + dout = {}; + + den = 2*r.^2+2*r.*dr+dr.^2+z.^2; + num = 2*r.*dr+dr.^2+z.^2; + val = num./den.^2; + ar = -2*(r+dr).*val; + num = -2*r.*dr-dr.^2+z.^2; + val = num./den.^2; + adr= -2*r.*val; + az = 4*(r+dr).*r.*z./den.^2; + + dout.ar = ar; + dout.adr= adr; + dout.az = az; + + ardr = dr.^4+4*dr.^3.*r-8*dr.*r.^3-4*r.^4-z.^4; + ardr = 2*ardr./den.^3; + + arz = 4*(r+dr).*z.*(dr.^2+2*dr.*r-2*r.^2+z.^2)./den.^3; + adrz= -4*r.*z.*(3*dr.^2+6*dr.*r+2*r.^2-z.^2)./den.^3; + azz = 4*r.*(r+dr).*(dr.^2+2*dr.*r+2*r.^2-3*z.^2)./den.^3; + + dout.ardr = ardr; + dout.arz = arz; + dout.adrz = adrz; + dout.azz = azz; + + r0 = sqrt(r.^2+(r+dr).^2+z.^2); + kr = r./r0; + kdr = (r+dr)./r0; + kz = z./r0; + krdr = -r.*(r+dr)./r0.^3; + krz = -r.*z./r0.^3; + kdrz = -(r+dr).*z./r0.^3; + kzz = (2*r.^2+2*r.*dr+dr.^2)./r0.^3; + + dout.kr = kr; + dout.kdr = kdr; + dout.kz = kz; + dout.krdr= krdr; + dout.krz = krz; + dout.kdrz= kdrz; + dout.kzz = kzz; + +end + diff --git a/chunkie/+chnk/+axissymhelm2d/der_ak_to_grad.m b/chunkie/+chnk/+axissymhelm2d/der_ak_to_grad.m new file mode 100644 index 0000000..15eb3f1 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/der_ak_to_grad.m @@ -0,0 +1,39 @@ +function [dout] = der_ak_to_grad(r,q,z,i,ida,idk,... + idaa,idak,idkk) + + dout = {}; + dout.int = i; + + ds = chnk.axissymhelm2d.dalph(r,q,z); + + intdr = ds.ar.*ida +ds.kr.*idk; + intdq = ds.adr.*ida+ds.kdr.*idk; + intdz = ds.az.*ida +ds.kz.*idk; + + intdrq= ds.ardr.*ida + ds.krdr.*idk + ... + ds.ar.*ds.adr.*idaa + ds.kr.*ds.kdr.*idkk + ... + ds.ar.*ds.kdr.*idak+ds.kr.*ds.adr.*idak; + + intdrz= ds.arz.*ida + ds.krz.*idk + ... + ds.ar.*ds.az.*idaa + ds.kr.*ds.kz.*idkk + ... + ds.ar.*ds.kz.*idak+ds.kr.*ds.az.*idak; + + intdqz= ds.adrz.*ida + ds.kdrz.*idk + ... + ds.adr.*ds.az.*idaa + ds.kdr.*ds.kz.*idkk + ... + ds.adr.*ds.kz.*idak+ds.kdr.*ds.az.*idak; + + intdzz= ds.azz.*ida + ds.kzz.*idk + ... + ds.az.*ds.az.*idaa + ds.kz.*ds.kz.*idkk + ... + ds.az.*ds.kz.*idak+ds.kz.*ds.az.*idak; + + dout.intdr = intdr; + dout.intdq = intdq; + dout.intdz = intdz; + + dout.intdrq = intdrq; + dout.intdrz = intdrz; + dout.intdqz = intdqz; + dout.intdzz = intdzz; + +end + diff --git a/chunkie/+chnk/+axissymhelm2d/div_by_kap.m b/chunkie/+chnk/+axissymhelm2d/div_by_kap.m new file mode 100644 index 0000000..05f48b3 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/div_by_kap.m @@ -0,0 +1,36 @@ +function [dout] = div_by_kap(r,q,z,dfunc) + + r0 = sqrt(r.^2+(r+q).^2+z.^2); + i0 = 1./r0; + i0da = 0; + i0dk = -1./r0.^2; + i0daa= 0; + i0dak= 0; + i0dkk= 2./r0.^3; + [dout] = chnk.axissymhelm2d.der_ak_to_grad(r,q,z,i0,i0da,i0dk,... + i0daa,i0dak,i0dkk); + + i = dfunc.int.*i0; + idr = dfunc.int.*dout.intdr + dfunc.intdr.*dout.int; + idq = dfunc.int.*dout.intdq + dfunc.intdq.*dout.int; + idz = dfunc.int.*dout.intdz + dfunc.intdz.*dout.int; + idqr = dfunc.intdrq.*dout.int + dfunc.int.*dout.intdrq + ... + dfunc.intdr.*dout.intdq + dfunc.intdq.*dout.intdr; + idrz = dfunc.intdrz.*dout.int + dfunc.int.*dout.intdrz + ... + dfunc.intdr.*dout.intdz + dfunc.intdz.*dout.intdr; + idqz = dfunc.intdqz.*dout.int + dfunc.int.*dout.intdqz + ... + dfunc.intdq.*dout.intdz + dfunc.intdz.*dout.intdq; + idzz = dfunc.intdzz.*dout.int + dfunc.int.*dout.intdzz + ... + 2*dfunc.intdz.*dout.intdz; + + dout.int = i; + dout.intdr = idr; + dout.intdq = idq; + dout.intdz = idz; + dout.intdrq = idqr; + dout.intdrz = idrz; + dout.intdqz = idqz; + dout.intdzz = idzz; + +end + diff --git a/chunkie/+chnk/+axissymhelm2d/green.m b/chunkie/+chnk/+axissymhelm2d/green.m new file mode 100644 index 0000000..fe57615 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/green.m @@ -0,0 +1,107 @@ +function [val, grad, hess] = green(k, src, targ, origin, ifdiff) +%CHNK.AXISSYMHELM2D.GREEN evaluate the helmholtz green's function +% for the given sources and targets. +% +% Note: that the first coordinate is r, and the second z. +% The code relies on precomputed tables and hence loops are required for +% computing various pairwise interactions. +% Finally, the code is not efficient in the sense that val, grad, hess +% are always internally computed independent of nargout +% +% Since the kernels are not translationally invariant in r, the size +% of the gradient is 3, for d/dr, d/dr', and d/dz +% +% Similarly the hessian is of size 6 and ordered as +% d_{rr}, d_{r'r'}, d_{zz}, d_{rr'}, d_{rz}, d_{r'z} + + +zr = real(k); zi = imag(k); +if abs(zr*zi) > eps + error('AXISSYMHELM2D.green: Axissymmetric greens function not supported for arbitrary complex', ... + 'wavenumbers, please input purely real or imaginary wave numbers\n'); +end + +if nargin == 4 + ifdiff = 0; +end + + +% Load precomputed tables +persistent asym_tables +if isempty(asym_tables) + asym_tables = chnk.axissymhelm2d.load_asym_tables(); +end + +[~, ns] = size(src); +[~, nt] = size(targ); + +vtmp = zeros(nt, ns); +gtmp = zeros(nt, ns, 3); +htmp = zeros(nt, ns, 6); + +% Determine whether to use +ifun = 1; +if abs(zr) < eps + ifun = 2; +end + +if ifdiff + if abs(zi) > eps + error('AXISSYMHELM2D.green: Difference of greens function only supported for real wave numbers\n'); + end + ifun = 3; +end + +over2pi = 1/2/pi; +kabs = abs(k); + +rt = repmat(targ(1,:).',1,ns); +rs = repmat(src(1,:),nt,1); +dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns); +r = (rt + origin(1))*kabs; +dz = dz*kabs; +dr = (rs-rt)*kabs; +dout = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, asym_tables,ifun); +int = dout.int; +intdz = dout.intdz; +intdq = dout.intdq; +intdr = dout.intdr; +intdzz = dout.intdzz; +intdrq = dout.intdrq; +intdrz = dout.intdrz; +intdqz = dout.intdqz; +for j = 1:ns + darea = src(1,j) + origin(1); + for i = 1:nt + if (i >1 && size(int,1)<=1) + disp("catastrophic error") + end + + vtmp(i,j) = int(i,j) * over2pi * kabs * darea; + + gtmp(i,j,1) = intdr(i,j) * over2pi * kabs.^2 * darea; + gtmp(i,j,2) = intdq(i,j) * over2pi * kabs.^2 * darea; + gtmp(i,j,3) = -intdz(i,j) * over2pi * kabs.^2 * darea; + +% Fix hessian scalings +% No drr, dr'r', or dr r' currently available + %htmp(i,j,1) = 0; + %htmp(i,j,2) = 0; + htmp(i,j,3) = intdzz(i,j) * over2pi * kabs.^3 * darea; + htmp(i,j,4) = intdrq(i,j) * over2pi * kabs.^3 * darea; + htmp(i,j,5) = intdrz(i,j) * over2pi * kabs.^3 * darea; + htmp(i,j,6) = intdqz(i,j) * over2pi * kabs.^3 * darea; + end +end + +if nargout > 0 + val = vtmp; +end + +if nargout > 1 + grad = gtmp; +end + +if nargout > 2 + hess = htmp; +end diff --git a/chunkie/+chnk/+axissymhelm2d/green_neu_all.m b/chunkie/+chnk/+axissymhelm2d/green_neu_all.m new file mode 100644 index 0000000..7916326 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/green_neu_all.m @@ -0,0 +1,87 @@ +function [valk, gradk, hessk, valik, gradik, ... + hessik, valdiff, graddiff, hessdiff] = green_neu_all(k, src, targ, origin) +%CHNK.AXISSYMHELM2D.GREEN evaluate the helmholtz green's function +% for the given sources and targets. +% +% Note: that the first coordinate is r, and the second z. +% The code relies on precomputed tables and hence loops are required for +% computing various pairwise interactions. +% Finally, the code is not efficient in the sense that val, grad, hess +% are always internally computed independent of nargout +% +% Since the kernels are not translationally invariant in r, the size +% of the gradient is 3, for d/dr, d/dr', and d/dz +% +% Similarly the hessian is of size 6 and ordered as +% d_{rr}, d_{r'r'}, d_{zz}, d_{rr'}, d_{rz}, d_{r'z} + +% Load precomputed tables +persistent asym_tables +if isempty(asym_tables) + asym_tables = chnk.axissymhelm2d.load_asym_tables(); +end + +[~, ns] = size(src); +[~, nt] = size(targ); + +over2pi = 1/2/pi; +kabs = abs(k); + +rt = repmat(targ(1,:).',1,ns); +rs = repmat(src(1,:),nt,1); +dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns); +r = (rt + origin(1))*kabs; +dz = dz*kabs; +dr = (rs-rt)*kabs; +ifun = 4; +[doutk, doutik, doutdiff] = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, asym_tables,ifun); + +pfac = over2pi * kabs; +gfac = pfac * kabs; +hfac = gfac * kabs; +[valk, gradk, hessk] = dout_to_vgh(doutk, ns, nt, src, origin, pfac, gfac, hfac); +[valik, gradik, hessik] = dout_to_vgh(doutik, ns, nt, src, origin, pfac, gfac, hfac); +[valdiff, graddiff, hessdiff] = dout_to_vgh(doutdiff, ns, nt, src, origin, pfac, gfac, hfac); + +end + +function [vtmp, gtmp, htmp] = dout_to_vgh(dout, ns, nt, src, origin, pfac, gfac, hfac) + + vtmp = zeros(nt, ns); + gtmp = zeros(nt, ns, 3); + htmp = zeros(nt, ns, 6); + + int = dout.int; + intdz = dout.intdz; + intdq = dout.intdq; + intdr = dout.intdr; + intdzz = dout.intdzz; + intdrq = dout.intdrq; + intdrz = dout.intdrz; + intdqz = dout.intdqz; + for j = 1:ns + darea = src(1,j) + origin(1); + pfacsc = pfac * darea; + gfacsc = gfac * darea; + hfacsc = hfac * darea; + for i = 1:nt + if (i >1 && size(int,1)<=1) + disp("catastrophic error") + end + + vtmp(i,j) = int(i,j) * pfacsc; + + gtmp(i,j,1) = intdr(i,j) * gfacsc; + gtmp(i,j,2) = intdq(i,j) * gfacsc; + gtmp(i,j,3) = -intdz(i,j) * gfacsc; + +% Fix hessian scalings +% No drr, dr'r', or dr r' currently available + + htmp(i,j,3) = intdzz(i,j) * hfacsc; + htmp(i,j,4) = intdrq(i,j) * hfacsc; + htmp(i,j,5) = intdrz(i,j) * hfacsc; + htmp(i,j,6) = intdqz(i,j) * hfacsc; + end + end +end diff --git a/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m b/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m new file mode 100644 index 0000000..a5dd026 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/helm_axi_all.m @@ -0,0 +1,198 @@ +function [varargout] = helm_axi_all(rs,drs,dzs,htables,ifun) + + r0s = sqrt(rs.^2+(rs+drs).^2+dzs.^2); + alphs = (drs.^2+dzs.^2)./r0s.^2; + + int = zeros(size(alphs)); + + iflag_rk = (ifun == 1 || ifun == 4); + iflag_ik = (ifun == 2 || ifun == 4); + iflag_dk = (ifun == 3 || ifun == 4); + + if (iflag_rk) + rk = cell(6); + for ii=1:6 + rk{ii} = int; + end + end + + if (iflag_ik) + ik = cell(6); + for ii=1:6 + ik{ii} = int; + end + end + + if (iflag_dk) + dk = cell(6); + for ii=1:6 + dk{ii} = int; + end + end + + iclose = find(alphs<0.2); + ifar = 1:numel(alphs); + ifar(iclose) = []; + + alph_far = alphs(ifar); + r0s_far = r0s(ifar); + [i_mid] = find((1-alph_far).*r0s_far < 150); + imid = ifar(i_mid); + ifar(i_mid) = []; + + alph_mid = alphs(imid); + r0s_mid = r0s(imid); + [i_midnear] = find((1-alph_mid).*r0s_mid < 40); + imidnear = imid(i_midnear); + imid(i_midnear) = []; + + alph_midnear = alphs(imidnear); + r0s_midnear = r0s(imidnear); + [i_midnearnear] = find((1-alph_midnear).*r0s_midnear < 40); + imidnearnear = imidnear(i_midnearnear); + imidnear(i_midnearnear) = []; + + + if (iflag_rk) + s ... + = chnk.axissymhelm2d.helm_axi_close_table(r0s(iclose),alphs(iclose),1,htables); + for ii=1:6 + rk{ii}(iclose) = s{ii}; + end + end + + if (iflag_ik) + s ... + = chnk.axissymhelm2d.helm_axi_close_table(r0s(iclose),alphs(iclose),2,htables); + for ii=1:6 + ik{ii}(iclose) = s{ii}; + end + end + + if (iflag_dk) + s ... + = chnk.axissymhelm2d.helm_axi_close_table(r0s(iclose),alphs(iclose),3,htables); + for ii=1:6 + dk{ii}(iclose) = s{ii}; + end + end + + + s ... + = chnk.axissymhelm2d.helm_axi_smooth(r0s(imidnearnear),alphs(imidnearnear),... + ifun,htables.xlege_midnear,htables.wlege_midnear); + + if (iflag_rk) + for ii=1:6 + rk{ii}(imidnearnear) = s.rk{ii}; + end + end + + if(iflag_ik) + for ii=1:6 + ik{ii}(imidnearnear) = s.ik{ii}; + end + end + + if(iflag_dk) + for ii=1:6 + dk{ii}(imidnearnear) = s.dk{ii}; + end + end + + s ... + = chnk.axissymhelm2d.helm_axi_smooth(r0s(imidnear),alphs(imidnear),... + ifun,htables.xlege_midnear,htables.wlege_midnear); + + if (iflag_rk) + for ii=1:6 + rk{ii}(imidnear) = s.rk{ii}; + end + end + + if (iflag_ik) + for ii=1:6 + ik{ii}(imidnear) = s.ik{ii}; + end + end + + if (iflag_dk) + for ii=1:6 + dk{ii}(imidnear) = s.dk{ii}; + end + end + + s ... + = chnk.axissymhelm2d.helm_axi_smooth(r0s(imid),alphs(imid),... + ifun,htables.xlege_mid,htables.wlege_mid); + + if(iflag_rk) + for ii=1:6 + rk{ii}(imid) = s.rk{ii}; + end + end + + if (iflag_ik) + for ii=1:6 + ik{ii}(imid) = s.ik{ii}; + end + end + + if (iflag_dk) + for ii=1:6 + dk{ii}(imid) = s.dk{ii}; + end + end + + s ... + = chnk.axissymhelm2d.helm_axi_smooth(r0s(ifar),alphs(ifar),... + ifun,htables.xlege,htables.wlege); + + if (iflag_rk) + for ii=1:6 + rk{ii}(ifar) = s.rk{ii}; + end + end + + if (iflag_ik) + for ii=1:6 + ik{ii}(ifar) = s.ik{ii}; + end + end + + if (iflag_dk) + for ii=1:6 + dk{ii}(ifar) = s.dk{ii}; + end + end + + + + if(iflag_rk) + [dsk] = proc_kerns(rs,drs,dzs,rk); + varargout{1} = dsk; + if (ifun == 4) + [dsik] = proc_kerns(rs,drs,dzs,ik); + [dsdiff] = proc_kerns(rs,drs,dzs,dk); + varargout{2} = dsik; + varargout{3} = dsdiff; + end + end + if(ifun == 2) + [dsik] = proc_kerns(rs,drs,dzs,ik); + varargout{1} = dsik; + end + if(ifun==3) + [dsdiff] = proc_kerns(rs,drs,dzs,dk); + varargout{1} = dsdiff; + end + +end + +function [ds] = proc_kerns(rs,drs,dzs,s) + [dout] = chnk.axissymhelm2d.der_ak_to_grad(rs,drs,dzs,s{1},s{2},... + s{3},s{6},s{5},s{4}); + [ds] = chnk.axissymhelm2d.div_by_kap(rs,drs,dzs,dout); + ds.intdrz = -ds.intdrz; + ds.intdzz = -ds.intdzz; +end \ No newline at end of file diff --git a/chunkie/+chnk/+axissymhelm2d/helm_axi_close_table.m b/chunkie/+chnk/+axissymhelm2d/helm_axi_close_table.m new file mode 100644 index 0000000..feefa3c --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/helm_axi_close_table.m @@ -0,0 +1,83 @@ +function [sout] ... + = helm_axi_close_table(r0s,alphs,ifun,htables) + int = zeros(size(alphs)); + kerns = int; + kernsdk = int; + kernsda = int; + kernsdkk= int; + kernsdak= int; + kernsdaa= int; + + for ii=1:numel(alphs) + + alph = alphs(ii); + r0 = r0s(ii); + + + ik = ceil(r0/(pi)); + ia = ceil(-log(alph*5)/log(2)); + + if (ia >111 || ik >152) + kern = 0; + kernda = 0; + kerndk = 0; + kerndak = 0; + kerndkk = 0; + kerndaa = 0; + else + krel = (r0-pi*(ik-1))*2/pi-1; + arel = ((alph - 0.2*2^(-ia))/(0.2*2^(-ia))-0.5)*2; + + tas = cos((0:(htables.ncheb-1))*acos(arel)); + tks = cos((0:(htables.ncheb-1))*acos(krel)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + vv = squeeze(htables.allvs(:,:,ifun,ia,ik)); + kern = tas*vv*tks.'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + vdk = squeeze(htables.alldk(:,:,ifun,ia,ik)); + kerndk = tas*vdk*tks.'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + vda = squeeze(htables.allda(:,:,ifun,ia,ik)); + kernda = tas*vda*tks.'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + vdak = squeeze(htables.alldak(:,:,ifun,ia,ik)); + kerndak = tas*vdak*tks.'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + vdkk = squeeze(htables.alldkk(:,:,ifun,ia,ik)); + kerndkk = tas*vdkk*tks.'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + vdaa = squeeze(htables.alldaa(:,:,ifun,ia,ik)); + kerndaa = tas*vdaa*tks.'; + + end + kerns(ii) = kern; + kernsdk(ii) = kerndk; + kernsda(ii) = kernda; + kernsdaa(ii)= kerndaa; + kernsdak(ii)= kerndak; + kernsdkk(ii)= kerndkk; + + end + + sout = {}; + sout{1} = kerns; + sout{3} = kernsdk; + sout{2} = kernsda; + sout{4} = kernsdkk; + sout{5} = kernsdak; + sout{6} = kernsdaa; + +end + diff --git a/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m b/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m new file mode 100644 index 0000000..740631f --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m @@ -0,0 +1,146 @@ +function [sout] ... + = helm_axi_smooth(r0s,alphs,ifun,xlege,wlege) + + +%%%% +% +% On output, a cell array with entries +% *kerns +% *kernsda +% *kernsdk +% *kernsdkk +% *kernsdak +% *kernsdaa +% +%%%% + + + [alp,xle] = meshgrid(alphs,xlege); + xle = cos(xle); + [r0t,wle] = meshgrid(r0s,wlege); + rts = sqrt(1-(1-alp).*xle); + + if (ifun == 2) + r0t = r0t*1i; + end + + efac = exp(1i*r0t.*rts)./rts.*wle; + r2 = rts.^2; + k2 = 1i*r0t.*rts; + + kerns = sum(asymint_v(rts,r0t,efac),1); + kernsdk = sum(asymintdk_v(rts,r0t,efac),1); + kernsda = sum(asymintda_v(xle,rts,r0t,efac),1); + kernsdaa = sum(asymintdaa_v(xle,rts,r0t,efac,k2,r2),1); + kernsdak = sum(asymintdak_v(xle,rts,r0t,efac),1); + kernsdkk = sum(asymintdkk_v(xle,rts,r0t,efac),1); + + if (ifun == 2) + kernsdk = kernsdk*1i; + kernsdak = kernsdak*1i; + kernsdkk = -kernsdkk; + end + + if (ifun == 3) + r0t = r0t*1i; + efac = exp(1i*r0t.*rts)./rts.*wle; + r2 = rts.^2; + k2 = 1i*r0t.*rts; + kern2_v = sum(asymint_v(rts,r0t,efac),1); + kern2dk_v = sum(asymintdk_v(rts,r0t,efac),1); + kern2da_v = sum(asymintda_v(xle,rts,r0t,efac),1); + kern2daa_v = sum(asymintdaa_v(xle,rts,r0t,efac,k2,r2),1); + kern2dak_v = sum(asymintdak_v(xle,rts,r0t,efac),1); + kern2dkk_v = sum(asymintdkk_v(xle,rts,r0t,efac),1); + kerns = kerns - kern2_v; + kernsdk = kernsdk - 1i*kern2dk_v; + kernsda = kernsda - kern2da_v; + kernsdak = kernsdak - 1i*kern2dak_v; + kernsdaa = kernsdaa - kern2daa_v; + kernsdkk = kernsdkk + kern2dkk_v; + end + + s ={}; + s{1} = kerns; + s{3} = kernsdk; + s{2} = kernsda; + s{5} = kernsdak; + s{6} = kernsdaa; + s{4} = kernsdkk; + + if (ifun == 1 || ifun==4) + sout.rk = s; + end + if (ifun ==2) + sout.ik = s; + end + if (ifun ==3) + sout.dk = s; + end + + if (ifun == 4) + r0t = r0t*1i; + efac = exp(1i*r0t.*rts)./rts.*wle; + r2 = rts.^2; + k2 = 1i*r0t.*rts; + kern2_v = sum(asymint_v(rts,r0t,efac),1); + kern2dk_v = sum(asymintdk_v(rts,r0t,efac),1); + kern2da_v = sum(asymintda_v(xle,rts,r0t,efac),1); + kern2daa_v = sum(asymintdaa_v(xle,rts,r0t,efac,k2,r2),1); + kern2dak_v = sum(asymintdak_v(xle,rts,r0t,efac),1); + kern2dkk_v = sum(asymintdkk_v(xle,rts,r0t,efac),1); + + ik = {}; + ik{1} = kern2_v; + ik{2} = kern2da_v; + ik{3} = 1i*kern2dk_v; + ik{6} = kern2daa_v; + ik{5} = 1i*kern2dak_v; + ik{4} = -kern2dkk_v; + sout.ik = ik; + + dk = {}; + for ii=1:6 + dk{ii} = sout.rk{ii}-sout.ik{ii}; + end + sout.dk = dk; + end + +end + +function [val] = asymint_v(r,k,efac) + val = efac; +end + +function [val] = asymintda_v(x,r,k,efac) + val = x.*efac./(r.^2).*(1i*k.*r-1)/2; + %val = exp(1i*k*sqrt(1-(1-a)*cos(x)))./sqrt(1-(1-a)*cos(x)); +end + +function [val] = asymintdaa_v(x,r,k,efac,k2,r2) + rdiv = r2.^2; + prefac = 0.25*(x.^2).*efac./rdiv; + fac = k2.^2-3*k2+3; + val = fac.*prefac; + %val = -exp(1i*k*sqrt(1-(1-a)*cos(x))).*sqrt(1-(1-a)*cos(x)); +end + +function [val] = asymintdak_v(x,r,k,efac) + val = -(k/2).*x.*efac; +end + + + +function [val] = asymintdk_v(r,k,efac) + val = 1i*efac.*r; +end + + +function [val] = asymintdkk_v(x,r,k,efac) + val = -efac.*(r.^2); +end + + + + + diff --git a/chunkie/+chnk/+axissymhelm2d/kern.m b/chunkie/+chnk/+axissymhelm2d/kern.m new file mode 100644 index 0000000..1bd2ea8 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/kern.m @@ -0,0 +1,421 @@ +function submat = kern(zk, srcinfo, targinfo, origin, type, varargin) +%CHNK.AXISSYMHELM2D.KERN axissymmetric Helmholtz layer potential kernels in 2D +% +% Syntax: submat = chnk.axissymhelm2d.kern(zk,srcinfo,targingo,type,htables) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). Note that the normal information is obtained +% by taking the perpendicular to the provided tangential deriviative +% info and normalizing +% +% Here the first and second components correspond to the r and z +% coordinates respectively. +% +% Kernels based on G(x,y) = \int_{0}^{\pi} e^{i k d(t)}/(d(t)) \, dt \, +% where d(t) = \sqrt(r^2 + r'^2 - 2rr' \cos(t) + (z-z')^2) with +% x = (r,z), and y = (r',z') +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% S'(x,y) = \nabla_{n_x} G(x,y) +% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y) +% +% Input: +% zk - complex number, Helmholtz wave number, must be purely real, or purely +% imaginary, doesn't support other complex wavenumbers yet +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% targinfo - description of targets in ptinfo struct format, +% if info not relevant (d/d2) it doesn't need to +% be provided. sprime requires tangent info in +% targinfo.d +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% 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 == 'c', combined layer kernel coef(1) D + coef(2) S +% +% varargin{1} - coef: length 2 array in the combined layer +% formula, 2x2 matrix for all kernels +% otherwise does nothing +% +% Output: +% submat - the evaluation of the selected kernel for the +% provided sources and targets. the number of +% rows equals the number of targets and the +% number of columns equals the number of sources +% +% + +src = srcinfo.r; +targ = targinfo.r; + +[~, ns] = size(src); +[~, nt] = size(targ); + +if strcmpi(type, 'd') + srcnorm = srcinfo.n; + [~, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin); + nx = repmat(srcnorm(1,:), nt, 1); + ny = repmat(srcnorm(2,:), nt, 1); + % Due to lack of translation invariance in r, no sign flip needed, + % as gradient is computed with repsect to r' + submat = (grad(:,:,2).*nx - grad(:,:,3).*ny); + fker = @(x, s, t, rns) fdlp(x, zk, s, t, rns, origin); + for j=1:ns + for i=1:nt + rt = targ(1,i) + origin(1); + dr = (src(1,j) - targ(1,i)); + dz = (src(2,j) - targ(2,i)); + r0 = sqrt(rt^2+(rt+dr)^2+dz^2); + alph = (dr^2+dz^2)/r0^2; + if alph > 2e-4 && alph < 0.2 + [x0, w0] = get_grid(zk, rt, dr, dz); + fvals = fker(x0, src(:, j), targ(:,i), srcnorm(:,j)); + submat(i,j) = 2*w0.'*fvals; + end + end + end +end + +if strcmpi(type, 'sprime') + targnorm = targinfo.n; + [~, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin); + + nx = repmat((targnorm(1,:)).',1,ns); + ny = repmat((targnorm(2,:)).',1,ns); + submat = (grad(:,:,1).*nx + grad(:,:,3).*ny); + + + fker = @(x, s, t, rnt) fsprime(x, zk, s, t, rnt, origin); + for j=1:ns + for i=1:nt + rt = targ(1,i) + origin(1); + dr = (src(1,j) - targ(1,i)); + dz = (src(2,j) - targ(2,i)); + r0 = sqrt(rt^2+(rt+dr)^2+dz^2); + alph = (dr^2+dz^2)/r0^2; + if alph > 2e-4 && alph < 0.2 + [x0, w0] = get_grid(zk, rt, dr, dz); + fvals = fker(x0, src(:, j), targ(:,i), targnorm(:,i)); + submat(i,j) = 2*w0.'*fvals; + end + end + end + +end + +if strcmpi(type, 's') + submat = chnk.axissymhelm2d.green(zk, src, targ, origin); + fker = @(x, s, t) fslp(x, zk, s, t, origin); + for j=1:ns + for i=1:nt + rt = targ(1,i) + origin(1); + dr = (src(1,j) - targ(1,i)); + dz = (src(2,j) - targ(2,i)); + r0 = sqrt(rt^2+(rt+dr)^2+dz^2); + + alph = (dr^2 + dz^2)/r0^2; + + if alph > 2e-4 && alph < 0.2 + [x0, w0] = get_grid(zk, rt, dr, dz); + fvals = fker(x0, src(:, j), targ(:,i)); + submat(i,j) = 2*w0.'*fvals; + end + end + end + +end + + +if strcmpi(type, 'sdiff') + ifdiff = 1; + submat = chnk.axissymhelm2d.green(zk, src, targ, origin, ifdiff); +end + +if strcmpi(type, 'c') + srcnorm = srcinfo.n; + coef = ones(2,1); + if (nargin == 6); coef = varargin{1}; end + nx = repmat(srcnorm(1,:), nt, 1); + ny = repmat(srcnorm(2,:), nt, 1); + [submats, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin); + % Due to lack of translation invariance in r, no sign flip needed, + % as gradient is computed with respect to r' + submat = coef(1)*(grad(:,:,2).*nx - grad(:,:,3).*ny) + coef(2)*submats; + + fker = @(x, s, t, rns) coef(1)*fdlp(x, zk, s, t, rns, origin) + ... + coef(2)*fslp(x, zk, s, t, origin); + for j=1:ns + for i=1:nt + rt = targ(1,i) + origin(1); + dr = (src(1,j) - targ(1,i)); + dz = (src(2,j) - targ(2,i)); + r0 = sqrt(rt^2+(rt+dr)^2+dz^2); + alph = (dr^2+dz^2)/r0^2; + if alph > 2e-4 && alph < 0.2 + [x0, w0] = get_grid(zk, rt, dr, dz); + fvals = fker(x0, src(:, j), targ(:,i), srcnorm(:,j)); + submat(i,j) = 2*w0.'*fvals; + end + end + end +end + + +if strcmpi(type, 'dprime') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin); + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + submat = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... + - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg; +end + + + +if strcmpi(type, 'dprimediff') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + ifdiff = 1; + [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin, ifdiff); + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + submat = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... + - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg; + + fker = @(x, s, t, rns, rnt) fdprimediff(x, zk, s, t, rns, rnt, origin); + for j=1:ns + for i=1:nt + rt = targ(1,i) + origin(1); + dr = (src(1,j) - targ(1,i)); + dz = (src(2,j) - targ(2,i)); + r0 = sqrt(rt^2+(rt+dr)^2+dz^2); + alph = (dr^2+dz^2)/r0^2; + if alph > 2e-4 && alph < 0.2 + [x0, w0] = get_grid(zk, rt, dr, dz); + fvals = fker(x0, src(:, j), targ(:,i), srcnorm(:,j), ... + targnorm(:,i)); + submat(i,j) = 2*w0.'*fvals; + end + end + end + +end + +if strcmpi(type, 'neu_rpcomb') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + [~,gk,~,sikmat,gik,~,~,~,hessdiff] = ... + chnk.axissymhelm2d.green_neu_all(zk, src, targ, origin); + alpha = 1; + if (nargin == 6); alpha = varargin{1}; end + if (size(alpha) > 1) + warning('Incorrect dimensions for coefs, using first component'); + alpha = alpha(1); + end + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + spmat = (gk(:,:,1).*nxtarg + gk(:,:,3).*nytarg); + spikmat = (gik(:,:,1).*nxtarg + gik(:,:,3).*nytarg); + dkdiffmat = hessdiff(:,:,4).*nxsrc.*nxtarg ... + - hessdiff(:,:,5).*nysrc.*nxtarg ... + - hessdiff(:,:,6).*nxsrc.*nytarg + hessdiff(:,:,3).*nysrc.*nytarg; + + rt = repmat(targ(1,:).',1,ns); rt = rt(:); + rs = repmat(src(1,:),nt,1); rs = rs(:); + dr = rs - rt; + rt = rt + origin(1); + dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns); dz = dz(:); + r0 = rt.^2 + (rt + dr).^2 + dz.^2; + alphs = (dr.^2 + dz.^2)./r0; + iflag = (alphs > 2e-4).*(alphs < 0.2); + if sum(iflag) + + % Assumes r,z are specified in meters, and, k is appropriately scaled + rtmax = max(targ(1,:)); + rsmax = max(src(1,:)); + rmax = max(rtmax,rsmax); + dr0 = 2e-4; + dz0 = 2e-4; + ppw = 10; + [x0, w0] = get_grid(zk, rmax, dr0, dz0, ppw); + iflag = reshape(iflag, [nt, ns]); + sxhalf = sin(x0/2); + sxhalf2 = sxhalf.*sxhalf; + cx = 1- 2*sxhalf2; + for j=1:ns + for i=1:nt + if iflag(i,j) + [fkp, fik, fikp, fkdiff] = get_neu_kers(zk, cx, sxhalf2, ... + src(:,j), targ(:,i), srcnorm(:,j), targnorm(:,i), origin); + spmat(i,j) = 2.*w0'*fkp; + sikmat(i,j) = 2.*w0'*fik; + spikmat(i,j) = 2.*w0'*fikp; + dkdiffmat(i,j) = 2.*w0'*fkdiff; + end + end + end + end + submat = zeros(3*nt, 3*ns); + c1 = -1.0/(0.5 + 0.25*1i*alpha); + c2 = 1i*alpha*c1; + submat(1:3:end, 1:3:end) = c1*spmat; + submat(1:3:end, 2:3:end) = c2*dkdiffmat; + submat(1:3:end, 3:3:end) = c2*spikmat; + submat(2:3:end, 1:3:end) = -sikmat; + submat(3:3:end, 1:3:end) = -spikmat; + + + + + + + + +end + +end + + + +function f = fslp (x, zk, s, t, o) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sin(x/2).^2); + f = exp(1j*zk*r)/4/pi./r.*(rs + o(1)); +end + + + +function f = fdlp (x, zk, s, t, rns, o) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + rnd = ((rt +o(1)).*cos(x) - (rs + o(1))).*rns(1) + (zt - zs).*rns(2); + + r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sin(x/2).^2); + f = rnd.*(1 - 1j*zk*r).*exp(1j*zk*r)/4/pi./r.^3.*(rs + o(1)); +end + + + +function f = fsprime (x, zk, s, t, rnt, o) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + rnd = ((rt + o(1)) - (rs + o(1)).*cos(x)).*rnt(1) + (zt - zs).*rnt(2); + r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs + o(1)).*(rt + o(1)).*sin(x/2).^2); + f = -rnd.*(1-1j*zk*r).*exp(1j*zk*r)/4/pi./r.^3.*(rs + o(1)); +end + + + +function f = fdprime (x, zk, s, t, rns, rnt, o) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + sxhalf = sin(x/2); + sxhalf2 = sxhalf.*sxhalf; + cx = 1-2*sxhalf2; + + rndt = ((rt + o(1)) - (rs + o(1)).*cx).*rnt(1) + (zt - zs).*rnt(2); + rnds = ((rt + o(1)).*cx - (rs + o(1))).*rns(1) + (zt - zs).*rns(2); + rnsnt = rns(1)*rnt(1).*cx + rns(2)*rnt(2); + + r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sxhalf2); + f = -(rnsnt.*(1j*zk.*r-1)./r.^3 + ... + rndt.*rnds.*(-zk^2.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs + o(1)); +end + +function [fkp, fik, fikp, fkdiff] = get_neu_kers(zk, cx, sxhalf2, s, t, rns, rnt, o) + + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + + rndt = ((rt + o(1)) - (rs + o(1)).*cx).*rnt(1) + (zt - zs).*rnt(2); + rnds = ((rt + o(1)).*cx - (rs + o(1))).*rns(1) + (zt - zs).*rns(2); + rnsnt = rns(1)*rnt(1).*cx + rns(2)*rnt(2); + + r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sxhalf2); + rinv = 1.0./r; + rinv2 = rinv.*rinv; + rinv3 = rinv.*rinv2; + rinv5 = rinv3.*rinv2; + + afac = 1/4/pi.*(rs + o(1)); + efac = exp(1j*zk*r).*afac; + efac_i = exp(-zk*r).*afac; + + fkp = -rndt.*(1-1j*zk*r).*rinv3.*efac; + + fik = efac_i.*rinv; + fikp = -rndt.*(1 + zk*r).*rinv3.*efac_i; + + + fkdiff = -(rnsnt.*(1j*zk.*r-1).*rinv3 + ... + rndt.*rnds.*(-zk^2.*r.^2 - 3*1j*zk.*r + 3).*rinv5).*efac; + fkdiff = fkdiff + (rnsnt.*(-zk.*r-1).*rinv3 + ... + rndt.*rnds.*(zk^2.*r.^2 + 3*zk.*r + 3).*rinv5).*efac_i; +end + + +function f = fsdiff (x, zk, s, t, o) + f = fslp(x, zk, s, t, o) - fslp(x, 1j*zk, s, t, o); +end + + + +function f = fdprimediff (x, zk, s, t, rns, rnt, o) + f1 = fdprime(x, zk, s, t, rns, rnt, o); + f2 = fdprime(x, 1j*zk, s, t, rns, rnt, o); + f = f1 - f2; +end + + +function [xlegs, wlegs] = get_grid(zk, rt, dr, dz, ppw) + + if (nargin <= 4); ppw = 10; end + persistent xlegloc wlegloc + k = 16; + if isempty(xlegloc) && isempty(wlegloc) + [xlegloc,wlegloc] = lege.exps(k); + end + + rs = rt + dr; + drdiff = sqrt((rt+rs)^2 + dz^2) - sqrt(dr^2 + dz^2); + npan = max(ceil(abs(zk)*drdiff*ppw/2/pi/k),3); + + h = pi/npan; + tends = h:h:pi; + + dr0 = sqrt(dr^2 + dz^2); + nref = max(ceil(-log(dr0/h)/log(2)),2); + tends = [2.^(-nref:-1)*h tends]; + tstarts = [0 tends(1:end-1)]; + + xlegs = tstarts + (xlegloc+1)/2*(tends-tstarts); + wlegs = wlegloc/2*(tends-tstarts); + xlegs = xlegs(:); + wlegs = wlegs(:); +end + diff --git a/chunkie/+chnk/+axissymhelm2d/load_asym_tables.m b/chunkie/+chnk/+axissymhelm2d/load_asym_tables.m new file mode 100644 index 0000000..fb5c523 --- /dev/null +++ b/chunkie/+chnk/+axissymhelm2d/load_asym_tables.m @@ -0,0 +1,95 @@ +function asym_tables = load_asym_tables() +%CHNK.AXISSYMHELM2D.load_asym_tables loads the precomputed tables +% for axissymmetric Helmholtz kernels + p = mfilename('fullpath'); + dirname = dir([p ,'.m']).folder; + fname = [dirname '/asym_helm_data.mat']; + load(fname,'allvs'); + msizes = size(allvs); + ncheb = msizes(1); + nks = msizes(5); + nas = msizes(4); + nkers = msizes(3); + xcheb = (0:(ncheb-1))/(ncheb-1)*pi; + + [N,X]=meshgrid(0:(ncheb-1),xcheb); + Tc2v = cos(N.*X); + Tv2c = inv(Tc2v); + Tfull = kron(Tv2c,Tv2c); + Tc2vd = N.*sin(N.*X)./sqrt(1-cos(X).^2); + Tc2vd(1,:) = (0:(ncheb-1)).^2; + Tc2vd(end,:) = (-1).^(1:ncheb).*(0:(ncheb-1)).^2; + Tc2cd = Tv2c*Tc2vd; + + xcheb = cos(xcheb)'; + + + allda = allvs; + alldk = allvs; + alldaa = allvs; + alldak = allvs; + alldkk = allvs; + + for kk=1:nks + for jj=1:nas + for ii=1:nkers +% vmat = squeeze(allvs(:,:,ii,jj,kk)); + if (ii ==1) + vmat = squeeze(allvs(:,:,ii,jj,kk)); + else + vmat = squeeze(allvs(:,:,ii,jj,kk)); + vmat = reshape(Tfull*vmat(:),[12,12]); + allvs(:,:,ii,jj,kk) = vmat; + end + ia = jj; + vda = (2)^(ia-1)*20*Tc2cd*vmat; + vdaa= (2)^(ia-1)*20*Tc2cd*vda; + vdk = 2/(pi)*Tc2cd*vmat.'; + vdkk= (2)/pi*Tc2cd*vdk; + vdak = 2/(pi)*Tc2cd*vda.'; + allda(:,:,ii,jj,kk) = vda; + alldaa(:,:,ii,jj,kk)= vdaa; + alldk(:,:,ii,jj,kk) = vdk.'; + alldak(:,:,ii,jj,kk)= vdak.'; + alldkk(:,:,ii,jj,kk)= vdkk.'; + end + end + end + + asym_tables =[]; + asym_tables.allvs = allvs; + asym_tables.allda = allda; + asym_tables.alldaa= alldaa; + asym_tables.alldk = alldk; + asym_tables.alldak= alldak; + asym_tables.alldkk= alldkk; + asym_tables.ncheb = ncheb; + + nlege = 500; + [xlege,wlege,~,~] = lege.exps(nlege); + xlege = (pi*(xlege+1)/2); + wlege = wlege*pi/2; + asym_tables.xlege = xlege; + asym_tables.wlege = wlege; + + nlege = 100; + [xlege,wlege,~,~] = lege.exps(nlege); + xlege = (pi*(xlege+1)/2); + wlege = wlege*pi/2; + asym_tables.xlege_mid = xlege; + asym_tables.wlege_mid = wlege; + + nlege = 50; + [xlege,wlege,~,~] = lege.exps(nlege); + xlege = (pi*(xlege+1)/2); + wlege = wlege*pi/2; + asym_tables.xlege_midnear = xlege; + asym_tables.wlege_midnear = wlege; + + nlege = 20; + [xlege,wlege,~,~] = lege.exps(nlege); + xlege = (pi*(xlege+1)/2); + wlege = wlege*pi/2; + asym_tables.xlege_midnearnear = xlege; + asym_tables.wlege_midnearnear = wlege; +end diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m index 477967c..438cd45 100644 --- a/chunkie/+chnk/+rcip/Rcompchunk.m +++ b/chunkie/+chnk/+rcip/Rcompchunk.m @@ -155,7 +155,37 @@ nR = 2*k*nedge*ndim; ts = cell(nedge,1); -chnkrlocal(1,nedge) = chunker(); +pref = []; +pref.k = k; +chnkrlocal(1,nedge) = chunker(pref); + + + +if(size(fkern)==1) + fkernlocal = fkern; + if isa(fkern, 'kernel') + if isa(fkern.shifted_eval, 'function_handle') + fkernlocal.eval = @(s,t) fkern.shifted_eval(s, t, ctr(:,1)); + end + end +else + fkernlocal(nedge,nedge) = kernel(); + for i=1:nedge + ici = iedgechunks(1,i); + for j=1:nedge + icj = iedgechunks(1,j); + fkernlocal(i,j) = fkern(ici,icj); + if isa(fkern(ici,icj), 'kernel') + if isa(fkern(ici,icj).shifted_eval, 'function_handle') + fkernlocal(i,j).eval = ... + @(s,t) fkern(ici,icj).shifted_eval(s,t,ctr(:,i)); + end + end + end + end +end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % begin recursion proper @@ -211,6 +241,7 @@ chnkrlocal(i).wts = weights(chnkrlocal(i)); end end + % construct the system matrix for local chunks if level == 1 @@ -220,7 +251,7 @@ end % test for opdims ~= [1,1] - [MAT,opts] = chunkermat(chnkrlocal,fkernlocal,opts,ilistl); + [MAT,opts] = chunkermat(chnkrlocal, fkernlocal, opts, ilistl); % diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m index 7552ebe..1d6ca2d 100644 --- a/chunkie/@chunkgraph/chunkgraph.m +++ b/chunkie/@chunkgraph/chunkgraph.m @@ -34,7 +34,7 @@ end methods - function obj = chunkgraph(verts,edgesendverts,fchnks,cparams) + function obj = chunkgraph(verts, edgesendverts, fchnks, cparams, pref) if (nargin == 0) return end @@ -90,10 +90,16 @@ end end + if nargin <= 4 + pref = []; + pref.nchmax = 10000; + pref.k = 16; + end + + if ~isfield(pref, 'nchmax') + pref.nchmax = 10000; + end - pref = []; - pref.nchmax = 10000; - pref.k = 16; echnks = chunker.empty(); for i=1:size(edgesendverts,2) diff --git a/chunkie/@kernel/axissymhelm2d.m b/chunkie/@kernel/axissymhelm2d.m new file mode 100644 index 0000000..21b589e --- /dev/null +++ b/chunkie/@kernel/axissymhelm2d.m @@ -0,0 +1,99 @@ +function obj = axissymhelm2d(type, zk, coefs) +%KERNEL.AXISSYMHELM2D Construct the axissymmetric Helmholtz kernel. +% KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) +% constructs the axissymmetric single-layer Helmholtz kernel with +% wavenumber ZK. +% +% KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) +% constructs the axissymmetric double-layer Helmholtz kernel with +% wavenumber ZK. +% +% KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) +% constructs the normal derivative of the axissymmetric single-layer +% Helmholtz kernel with wavenumber ZK. +% +% KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or +% KERNEL.AXISSYMHELM2D('combined', ZK, COEFS) +% constructs the combined-layer axissymmetric Helmholtz kernel with +% wavenumber ZK and parameter COEFS, +% i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + +% COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK). +% +% NOTES: The axissymetric kernels are currently supported only for purely +% real or purely imaginary wave numbers +% +% See also CHNK.AXISSYMHELM2D.KERN. + +if ( nargin < 1 ) + error('Missing Helmholtz kernel type.'); +end + +if ( nargin < 2 ) + error('Missing Helmholtz wavenumber.'); +end + + +zr = real(zk); zi = imag(zk); +if abs(zr)*abs(zi) > eps + error('Only purely real or purely imaginary wavenumbers supported'); +end + +obj = kernel(); +obj.name = 'axissymhelmholtz'; +obj.params.zk = zk; +obj.opdims = [1 1]; + +switch lower(type) + + case {'s', 'single'} + obj.type = 's'; + obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 's'); + obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 's'); + obj.fmm = []; + obj.sing = 'log'; + + case {'d', 'double'} + obj.type = 'd'; + obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'd'); + obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'd'); + obj.fmm = []; + obj.sing = 'log'; + + case {'sp', 'sprime'} + obj.type = 'sp'; + obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'sprime'); + obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'sprime'); + obj.fmm = []; + obj.sing = 'log'; + + case {'c', 'combined'} + if ( nargin < 3 ) + warning('Missing combined layer coefficients. Defaulting to [1,1i].'); + coefs = [1,1i]; + end + obj.type = 'c'; + obj.params.coefs = coefs; + obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'c', coefs); + obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'c', coefs); + obj.fmm = []; + obj.sing = 'log'; + case {'neu_rpcomb'} + obj.type = 'neu_rpcomb'; + if ( nargin < 3 ) + warning('Missing coefficient of 1i*D. Defaulting to 1.'); + coefs = 1; + end + obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'neu_rpcomb', coefs); + obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'neu_rpcomb', coefs); + obj.fmm = []; + obj.sing = 'log'; + obj.opdims = [3 3]; + obj.params.c1 = -1.0/(0.5 + 0.25*1i*coefs); + obj.params.c2 = -1i*coefs/(0.5 + 0.25*1i*coefs); + + otherwise + error('Unknown axissym Helmholtz kernel type ''%s''.', type); + +end + +end diff --git a/chunkie/@kernel/axissymhelm2ddiff.m b/chunkie/@kernel/axissymhelm2ddiff.m new file mode 100644 index 0000000..3c7b1a7 --- /dev/null +++ b/chunkie/@kernel/axissymhelm2ddiff.m @@ -0,0 +1,108 @@ +function obj = axissymhelm2ddiff(type, zks, coefs) +%KERNEL.AXISSYMHELM2DDIFF Construct the difference of +% axissymmetric Helmholtz kernels. +% KERNEL.AXISSYMHELM2DDIFF('dp', ZKS) or +% KERNEL.AXISSYMHELM2DDIFF('dprime', ZKS) +% constructs the difference of Neumann data for double-layer +% axissymmetric Helmholtz kernels with wavenumbers ZKS(1) and ZKS(2), +% scaled by COEFS, i.e. +% COEFS(1)*KERNEL.AXISSYMHELM2D('dp', ZKS(1)) - +% COEFS(2)*KERNEL.AXISSYMHELM2D('dp', ZKS(2)) +% +% COEFS is optional. If not given then COEFS is set to [1 1]. +% +% NOTES: currently only supports coefs = [c c], and zks(1), real, and +% zks(2) = I*zks(1) +% +% See also CHNK.AXISSYMHELM2D.KERN. + +if ( nargin < 1 ) + error('Missing Helmholtz kernel type.'); +end + +if ( nargin < 2 ) + error('Missing Helmholtz wavenumber.'); +end + +zr1 = real(zks(1)); zi1 = imag(zks(1)); +zr2 = real(zks(2)); zi2 = imag(zks(2)); + +if zi1 ~=0 || zr2 ~=0 || zr1 ~= zi2 + error('Wave numbers must be of the form zk, 1i*zk with zk real'); +end + +obj = kernel(); +obj.name = 'axissymhelmholtzdiff'; +obj.params.zks = zks; +obj.opdims = [1 1]; +if(nargin < 3) + coefs = [1 1]; +end + + +obj.params.coefs = coefs; + +switch lower(type) + + case {'s', 'single'} + obj.type = 's'; + obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(zks(1), ... + s, t, [0,0], 's') - ... + coefs(2)*chnk.axissymhelm2d.kern(zks(2), ... + s, t, [0,0], 's'); + obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(zks(1), ... + s, t, o, 's') - ... + coefs(2)*chnk.axissymhelm2d.kern(zks(2), ... + s, t, o, 's'); + obj.fmm = []; + obj.sing = 'log'; + + case {'d', 'double'} + obj.type = 'd'; + obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(zks(1), ... + s, t, [0,0], 'd') - ... + coefs(2)*chnk.axissymhelm2d.kern(zks(2), ... + s, t, [0,0], 'd'); + obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(zks(1), ... + s, t, o, 'd') - ... + coefs(2)*chnk.axissymhelm2d.kern(zks(2), ... + s, t, o, 'd'); + obj.fmm = []; + obj.sing = 'log'; + + case {'sp', 'sprime'} + obj.type = 'sp'; + obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(zks(1), ... + s, t, [0,0], 'sprime') - ... + coefs(2)*chnk.axissymhelm2d.kern(zks(2), ... + s, t, [0,0], 'sprime'); + obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(zks(1), ... + s, t, o, 'sprime') - ... + coefs(2)*chnk.axissymhelm2d.kern(zks(2), ... + s, t, o, 'sprime'); + obj.fmm = []; + obj.sing = 'log'; + + case {'dp', 'dprime'} + + if coefs(1) ~= coefs(2) + error('Coefs must be of form [c c]'); + end + obj.type = 'dp'; + obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(real(zks(1)), ... + s, t, [0,0], 'dprimediff'); + obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(real(zks(1)), ... + s, t, o, 'dprimediff'); + obj.fmm = []; + if ( abs(coefs(1)-coefs(2)) < eps ) + obj.sing = 'log'; + else + obj.sing = 'hs'; + end + + otherwise + error('Unknown axissymmetric Helmholtz difference kernel type ''%s''.', type); + +end + +end diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index 30803a6..a25b09f 100644 --- a/chunkie/@kernel/kernel.m +++ b/chunkie/@kernel/kernel.m @@ -15,7 +15,10 @@ % 'stokes' ('stok', 's') 'svel', 'spres', 'strac', % 'dvel', 'dpres', 'dtrac' % 'zeros' ('zero','z') -% +% 'axis sym helmholtz' 's' 'd' 'sp' 'c' +% ('axissymh', 'axissymhelm') +% 'axis sym helmholtz difference' 's' 'd' 'sp' 'dp' +% ('axissymhdiff', 'axissymhelmdiff') % The types may also be written in longer form, e.g. 'single', 'double', % 'sprime', 'combined', 'svelocity', 'spressure', 'straction', % 'dvelocity', 'dpressure', 'dtraction'. @@ -50,6 +53,7 @@ type % Type of the kernel params % Structure of kernel parameters eval % Function handle for kernel evaluation + shifted_eval % Function handle for evaluating translated kernel evaluation fmm % Function handle for kernel FMM sing % Singularity type splitinfo % Kernel-split information @@ -79,6 +83,11 @@ obj = kernel.elast2d(varargin{:}); case {'zeros', 'zero', 'z'} obj = kernel.zeros(varargin{:}); + case {'axis sym helmholtz', 'axissymh', 'axissymhelm'} + obj = kernel.axissymhelm2d(varargin{:}); + case {'axis sym helmholtz difference', 'axissymhdiff' ... + 'axissymhelmdiff'} + obj = kernel.axissymhelm2ddiff(varargin{:}); otherwise error('Kernel ''%s'' not found.', kern); end @@ -109,6 +118,8 @@ obj = helm2ddiff(varargin); obj = stok2d(varargin); obj = elast2d(varargin); + obj = axissymhelm2d(varargin); + obj = axissymhelm2ddiff(varargin); obj = zeros(varargin); end @@ -154,6 +165,37 @@ end + + function out = shifted_eval_(s, t, o) + + [~, ns] = size(s.r); + [~, nt] = size(t.r); + + % Compute interleaved indices + ridx = cell(m, 1); + cidx = cell(n, 1); + for k = 1:m + ridx{k} = (rowstarts(k)+1:rowstarts(k+1)).' + (0:nt-1)*opdims(1); + ridx{k} = ridx{k}(:).'; + end + for l = 1:n + cidx{l} = ((colstarts(l)+1):colstarts(l+1)).' + (0:ns-1)*opdims(2); + cidx{l} = cidx{l}(:).'; + end + + % Evaluate each sub-kernel and assign the resulting block to the + % output matrix using interleaved indices + out = zeros(opdims(1)*nt, opdims(2)*ns); + for k = 1:m + for l = 1:n + out(ridx{k},cidx{l}) = kerns(k,l).shifted_eval(s,t,o); + end + end + + end + + + function varargout = fmm_(eps, s, t, sigma) [~, ns] = size(s.r); @@ -223,11 +265,24 @@ if ( any(strcmpi(sings, 'pv')) ), K.sing = 'pv'; end if ( any(strcmpi(sings, 'hs')) ), K.sing = 'hs'; end +% Set params +K.params = cell(m, n); +for kk=1:m + for ll=1:n + K.params{kk,ll} = kerns(kk,ll).params; + end +end + % The new kernel has eval() only if all sub-kernels have eval() if ( all(cellfun('isclass', {kerns.eval}, 'function_handle')) ) K.eval = @eval_; end +% The new kernel has shifted_eval() only if all sub-kernels have eval() +if ( all(cellfun('isclass', {kerns.shifted_eval}, 'function_handle')) ) + K.shifted_eval = @shifted_eval_; +end + % The new kernel has fmm() only if all sub-kernels have fmm() if ( all(cellfun('isclass', {kerns.fmm}, 'function_handle')) ) K.fmm = @fmm_; diff --git a/chunkie/@kernel/times.m b/chunkie/@kernel/times.m index 8f10ed7..db7ace3 100644 --- a/chunkie/@kernel/times.m +++ b/chunkie/@kernel/times.m @@ -16,6 +16,12 @@ f.eval = []; end + if(isa(f.shifted_eval, 'function_handle')) + f.shifted_eval = @(varargin) g*f.shifted_eval(varargin{:}); + else + f.shifted_eval = []; + end + if(isa(f.fmm, 'function_handle')) f.fmm = @(varargin) g*f.fmm(varargin{:}); else diff --git a/chunkie/@kernel/zeros.m b/chunkie/@kernel/zeros.m index 1982952..4c108c3 100644 --- a/chunkie/@kernel/zeros.m +++ b/chunkie/@kernel/zeros.m @@ -23,6 +23,12 @@ out = zeros(m*nt, n*ns); end + function out = shifted_eval_(s, t, o) + [~, ns] = size(s.r); + [~, nt] = size(t.r); + out = zeros(m*nt, n*ns); + end + function varargout = fmm_(eps, s, t, sigma) if ( isstruct(t) ) @@ -45,6 +51,7 @@ obj.opdims = [m n]; obj.sing = 'smooth'; obj.eval = @eval_; +obj.shifted_eval = @shifted_eval_; obj.fmm = @fmm_; end diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m index cde9605..4a23291 100644 --- a/chunkie/chunkerinterior.m +++ b/chunkie/chunkerinterior.m @@ -18,6 +18,7 @@ % opts - options structure with entries: % opts.fmm = boolean, use FMM % opts.flam = boolean, use FLAM routines +% opts.axissym = boolean, chunker is axissymmetric % Note on the default behavior: % by default it tries to use the fmm if it exists, if it doesn't % then unless explicitly set to false, it tries to use flam @@ -36,7 +37,6 @@ grid = false; - if nargin < 3 opts = []; end @@ -80,6 +80,35 @@ useflam = opts.flam; end +axissym = false; +if isfield(opts,'axissym') + axissym = opts.axissym; +end + +if axissym + nch = chnkr.nch; + istart = nch+1; + iend = 2*nch; + chnkr = sort(chnkr); + chnkr = chnkr.addchunk(nch); + chnkr.r(:,:,istart:iend) = fliplr(chnkr.r(:,:,1:nch)); + chnkr.r(1,:,istart:iend) = -chnkr.r(1,:,istart:iend); + chnkr.d(:,:,istart:iend) = fliplr(chnkr.d(:,:,1:nch)); + chnkr.d(2,:,istart:iend) = -chnkr.d(2,:,istart:iend); + + chnkr.d2(:,:,istart:iend) = fliplr(chnkr.d2(:,:,1:nch)); + chnkr.d2(1,:,istart:iend) = chnkr.d2(1,:,istart:iend); + + chnkr.h(istart:iend) = chnkr.h(1:nch); + chnkr.wts = weights(chnkr); + chnkr.n = normals(chnkr); + chnkr.adj(1,:) = 0:chnkr.nch-1; + chnkr.adj(2,:) = 2:chnkr.nch+1; + chnkr.adj(1,1) = chnkr.nch; + chnkr.adj(2,chnkr.nch) = 1; +end + + usefmm_final = false; useflam_final = false; @@ -113,6 +142,8 @@ [chnkr2] = upsample(chnkr,nlegnew); + + icont = false; if usefmm_final try diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index b10b495..46a2ed1 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -281,7 +281,6 @@ warning(sprintf(fstr)) end end - if strcmpi(imethod,'direct') % do dense version diff --git a/chunkie/demo/demo_axissym_neumann.m b/chunkie/demo/demo_axissym_neumann.m new file mode 100644 index 0000000..db289f4 --- /dev/null +++ b/chunkie/demo/demo_axissym_neumann.m @@ -0,0 +1,453 @@ +addpaths_loc(); +clear(); + +rng(pi) + +zk = 6*pi; + +type = 'chnkr-star'; +% type = 'chnkr-torus'; +type = 'cgrph'; +% type = 'cgrph-sphere'; + +irep = 'rpcomb'; + +pref = []; +pref.k = 16; +ns = 10; +nt = 100; +ppw = 10; % points per wavelength; +maxchunklen = pref.k/ppw/real(zk)*2*pi; + +[chnkr, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen); +wts = chnkr.wts; wts = wts(:); + +fprintf('Done building geometry\n'); + + +ngrid = 100; +grid = linspace(-1,1,ngrid); +[xt,yt] = meshgrid(grid, grid); +targs = [xt(:).'; yt(:).']; + +opts_in = []; +opts_in.axissym = 1; +opts_in.fmm = true; +[in] = chunkerinterior(chnkr, targs, opts_in); + + +plot(chnkr,'kx'); hold on; +plot(targs(1,in), targs(2,in), 'b.'); + + +% source strengths +strengths = randn(ns, 1); + + +% Plot everything + +figure(1) +clf +hold off +plot(chnkr, 'k.'); +hold on +quiver(chnkr); +scatter(sources(1,:), sources(2,:), 'o') +scatter(targets(1,:), targets(2,:), 'x') +axis equal + + +% For solving exterior the Neumann boundary value problem, we use the +% following integral equation +% +% u = \beta(S_{k} + i \alpha D_{k} S_{ik}) \sigma +% +% with \beta = -1.0/(0.5 + 1i*0.25*alpha) +% +% On imposing the boundary conditions, we get the following integral +% equation +% +% du/dn = I + \beta S_{k}'[\sigma] + +% i\beta \alpha(D'_{k} - D'_{ik})(S_{ik}[\sigma]) + +% i\beta \alpha(S_{ik}')^2 [\sigma]; +% +% Setting -S_{ik}[\sigma] + \tau = 0, and - S_{ik}'[\sigma] + \mu=0, +% we have the following system of integral equations + + +% Set up boundary data +Skp = kernel('axissymhelm', 'sprime', zk); +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:); +kernmats = Skp.eval(srcinfo, targinfo); +ubdry = kernmats*strengths; + + +% [Skpmat, Sikmat, Sikpmat, Dkdiffmat] = get_neumann_matrices(zk, chnkr); +[Skpmat, Sikmat, Sikpmat, Dkdiffmat] = get_neumann_matrices_fast(zk, chnkr); + + +alpha = 1; +c1 = -1/(0.5 + 1i*alpha*0.25); +c2 = c1*1i*alpha; +c3 = -1; + +n = chnkr.npt; +start = tic; +M = c1*(Skpmat + 1i*alpha*(Dkdiffmat*Sikmat + Sikpmat*Sikpmat)) + eye(n); +t1 = toc(start); +fprintf('Time taken for matrix matrix product:%d \n',t1); + +rhs = ubdry; +start = tic; sol = M\rhs; t1 = toc(start); +fprintf('Time taken in solve: %d\n',t1); + +start = tic; Minv = inv(M); t1 = toc(start); +fprintf('Time taken in computing inverse: %d\n',t1); + +start = tic; sol2 = Minv*rhs; t1 = toc(start); +fprintf('Time taken in computing solution with precomp inverse: %d\n',t1); + +siksol = Sikmat*sol; + +%% Begin postprocessing + + +Sk = kernel('axissymhelm', 's', zk); +Dk = kernel('axissymhelm', 'd', zk); + +% Compute exact solution +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = targets; +kernmatstarg = Sk.eval(srcinfo, targinfo); +utarg = kernmatstarg*strengths; + + +Keval = c1*kernel([Sk 1i*alpha*Dk]); + + +soluse = zeros(2*n,1); +soluse(1:2:end) = sol; +soluse(2:2:end) = siksol; + +if isa(chnkr, 'chunkgraph') + % Collapse cgrph into chnkrtotal + chnkrs = chnkr.echnks; + chnkrtotal = merge(chnkrs); +else + chnkrtotal = chnkr; +end + + +start = tic; +Dsol = chunkerkerneval(chnkrtotal, Keval, soluse, targets); +t2 = toc(start); +fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n', t2); + +wchnkr = chnkrtotal.wts; wchnkr = wchnkr(:); +relerr = norm(utarg-Dsol) / (sqrt(chnkrtotal.nch)*norm(utarg)); +relerr2 = norm(utarg-Dsol, 'inf') / dot(abs(sol(:)), wchnkr(:)); +fprintf('relative frobenius error %5.2e\n', relerr); +fprintf('relative l_inf/l_1 error %5.2e\n', relerr2); + + + + +function [Skpmat, Sikmat, Sikpmat, Dkdiffmat] = get_neumann_matrices(zk, chnkr) + + nn = chnkr.npt; + Skp = kernel('axissymhelm', 'sprime', zk); + Sik = kernel('axissymhelm', 's', 1i*zk); + Sikp = kernel('axissymhelm', 'sprime', 1i*zk); + Dkdiff = kernel('axissymhelmdiff', 'dprime', [zk 1i*zk]); + + opts = []; + opts.nonsmoothonly = true; + opts.rcip = false; + opts.nsub_or_tol = 20; + start = tic; + + kernels = cell(1,4); + kernels{1} = Skp; + kernels{2} = Sik; + kernels{3} = Sikp; + kernels{4} = Dkdiff; + + spmats = cell(1,4); + + parfor imat=1:4 + spmats{imat} = chunkermat(chnkr, kernels{imat}, opts); + end + Skp_spmat = spmats{1}; + Sik_spmat = spmats{2}; + Sikp_spmat = spmats{3}; + Dkdiff_spmat = spmats{4}; + + + + t1 = toc(start); + fprintf('Time taken in sparse matrix generation: %d\n',t1); + + nthd = maxNumCompThreads*2; + nbsize = ceil(chnkr.npt/nthd); + Skp_cellmat = cell(nthd,1); + Sik_cellmat = cell(nthd,1); + Sikp_cellmat = cell(nthd,1); + Dkdiff_cellmat = cell(nthd,1); + l2scale = false; + start = tic; + nn = chnkr.npt; + parfor i=1:nthd + opdims = [1 1]; + iind = ((i-1)*nbsize+1):min(nn,i*nbsize); + Skp_cellmat{i} = chnk.flam.kernbyindex(iind, 1:nn, chnkr, ... + Skp, opdims, Skp_spmat); + + + Sik_cellmat{i} = chnk.flam.kernbyindex(iind, 1:nn, chnkr, ... + Sik, opdims, Sik_spmat); + + + Sikp_cellmat{i} = chnk.flam.kernbyindex(iind, 1:nn, chnkr, ... + Sikp, opdims, Sikp_spmat); + + + Dkdiff_cellmat{i} = chnk.flam.kernbyindex(iind, 1:nn, chnkr, ... + Dkdiff, opdims, Dkdiff_spmat); + end + + Skpmat = vertcat(Skp_cellmat{:}); + Sikmat = vertcat(Sik_cellmat{:}); + Sikpmat = vertcat(Sikp_cellmat{:}); + Dkdiffmat = vertcat(Dkdiff_cellmat{:}); + t1 = toc(start); + fprintf('Time taken in matrix entry generation: %d\n',t1); +end + + +function [Skpmat, Sikmat, Sikpmat, Dkdiffmat] = get_neumann_matrices_fast(zk, chnkr) + + + Nkerns = kernel('axissymhelm', 'neu_rpcomb', zk); + + opts = []; + opts.nonsmoothonly = true; + opts.rcip = false; + opts.nsub_or_tol = 20; + start = tic; + + spmat = chunkermat(chnkr, Nkerns, opts); + + + t1 = toc(start); + fprintf('Time taken in sparse matrix generation (new): %d\n',t1); + + nsys = 3*chnkr.npt; + nthd = maxNumCompThreads*2; + nbsize = ceil(nsys/nthd); + N_cellmat = cell(nthd,1); + start = tic; + opdims = Nkerns.opdims; + parfor i=1:nthd + iind = ((i-1)*nbsize+1):min(nsys,i*nbsize); + N_cellmat{i} = chnk.flam.kernbyindex(iind, 1:nsys, chnkr, ... + Nkerns, opdims, spmat); + end + + Nmat = vertcat(N_cellmat{:}); + Skpmat = Nmat(1:3:end, 1:3:end)/Nkerns.params.c1; + Dkdiffmat = Nmat(1:3:end, 2:3:end)/Nkerns.params.c2; + Sikpmat = Nmat(1:3:end, 3:3:end)/Nkerns.params.c2; + Sikmat = -Nmat(2:3:end, 1:3:end); + t1 = toc(start); + fprintf('Time taken in matrix entry generation (new): %d\n',t1); +end + + + + + +function [chnkobj, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen) + +if nargin == 0 + type = 'chnkr'; +end + +if nargin <= 1 + pref = []; + pref.k = 16; +end + +if nargin <= 2 + ns = 10; +end + +if nargin <= 3 + nt = 10; +end + + +if nargin <= 4 + maxchunklen = 1.0; +end + + +if strcmpi(type, 'cgrph') + + + nverts = 3; + verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1)); + verts = [real(verts);imag(verts)]; + + + iind = 1:(nverts-1); + jind = 1:(nverts-1); + + iind = [iind iind]; + jind = [jind jind + 1]; + jind(jind>nverts) = 1; + svals = [-ones(1,nverts-1) ones(1,nverts-1)]; + edge2verts = sparse(iind, jind, svals, nverts-1, nverts); + + amp = 0.1; + frq = 2; + fchnks = cell(1,size(edge2verts,1)); + for icurve = 1:size(edge2verts,1) + fchnks{icurve} = @(t) sinearc(t, amp, frq); + end + cparams = []; + cparams.nover = 2; + cparams.maxchunklen = maxchunklen; + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + ts = -pi/2 + pi*rand(ns,1); + sources = 0.2*[cos(ts)';sin(ts)']; + + ts = -pi/2 + pi*rand(nt,1); + targets = 3.0*[cos(ts)'; sin(ts)']; + + + +elseif strcmpi(type,'cgrph-sphere') + + + nverts = 2; + verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1)); + verts = [real(verts);imag(verts)]; + + + iind = 1:(nverts-1); + jind = 1:(nverts-1); + + iind = [iind iind]; + jind = [jind jind + 1]; + jind(jind>nverts) = 1; + svals = [-ones(1,nverts-1) ones(1,nverts-1)]; + edge2verts = sparse(iind, jind, svals, nverts-1, nverts); + + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = false; + cparams.ta = -pi/2; + cparams.tb = pi/2; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.0; + + fchnks{1} = @(t) circle(t); + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.1*sources; + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + + + +elseif strcmpi(type,'chnkr-star') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = false; + cparams.ta = -pi/2; + cparams.tb = pi/2; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.5*sources; + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + +elseif strcmpi(type,'chnkr-torus') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = true; + cparams.ta = 0; + cparams.tb = 2*pi; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + ctr = [3 0]; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.5*sources; + sources(1,:) = sources(1,:) + ctr(1); + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + targets(1,:) = targets(1,:) + ctr(1); +end + + + +end + + +function [r,d,d2] = sinearc(t,amp,frq) +xs = t; +ys = amp*sin(frq*t); +xp = ones(size(t)); +yp = amp*frq*cos(frq*t); +xpp = zeros(size(t)); +ypp = -frq*frq*amp*sin(t); + +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end + + +function [r,d,d2] = circle(t) +xs = cos(-pi/4 + pi/2*t); +ys = sin(-pi/4 + pi/2*t); +xp = -pi*ys/2; +yp = pi*xs/2; +xpp = -pi*pi*xs/4; +ypp = -pi*pi*ys/4; +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end + + + diff --git a/devtools/test/axissymkernTest.m b/devtools/test/axissymkernTest.m new file mode 100644 index 0000000..4658df4 --- /dev/null +++ b/devtools/test/axissymkernTest.m @@ -0,0 +1,200 @@ + +clearvars; close all; +addpaths_loc(); +cparams = []; +cparams.eps = 1.0e-4; +pref = []; +pref.k = 16; +narms = 0; +amp = 0.25; +cparams.maxchunklen = 0.5; +cparams.ta = -pi/2; +cparams.tb = pi/2; + +start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); +t1 = toc(start); + +[~,~,info] = sortinfo(chnkr); + +np = 100; +isstart = ceil(chnkr.npt/3); +isend = isstart + np; +isind = isstart:isend; + +isind = 1:chnkr.npt; +% isind = 81:95; + +itstart = ceil(2*chnkr.npt/3); +itstart = isstart+np+1; +itend = itstart + np; +itind = itstart:itend; +itind = 1:chnkr.npt; +% itind = 477:495; + +srcinfo = []; +srcinfo.r = chnkr.r(:,isind); +srcinfo.d = chnkr.d(:,isind); +srcinfo.d2 = chnkr.d2(:,isind); +srcinfo.n = chnkr.n(:,isind); + +targinfo = []; +targinfo.r = chnkr.r(:,itind); +targinfo.d = chnkr.d(:,itind); +targinfo.d2 = chnkr.d2(:,itind); +targinfo.n = chnkr.n(:,itind); + +% Real k tests +zk = 40.1; + +type = 'sprime'; +origin = [0 0]; +start = tic; submat = chnk.axissymhelm2d.kern(zk, srcinfo, targinfo, origin, type); +t1 = toc(start); + + +fprintf('First call to axissymhelm2d.kern time: %d\n',t1); + +v = get_exact_kernels(zk, srcinfo, targinfo, type); + + +err1 = norm(v(:) - submat(:))/norm(submat(:)); +fprintf('Error in kernel %s = %d \n',type,err1); +fprintf('ratios = %d\n', max(abs(v(1:end)./submat(1:end)-1))); + + + +% Now test the shifted kernel bit +origin_new = [3.1 2.1]; +chnkr_shift = chnkr; +chnkr_shift.r(1,:) = chnkr.r(1,:) + origin_new(1); +chnkr_shift.r(2,:) = chnkr.r(2,:) + origin_new(2); + +K = kernel('axissymhelm', type, zk); + +start = tic; submat = chnk.axissymhelm2d.kern(zk, srcinfo, targinfo, origin_new, type); +submat = K.shifted_eval(srcinfo, targinfo, origin_new); +t1 = toc(start); + +fprintf('First call to axissymhelm2d.kern time%d\n',t1); + + +srcinfo.r(1,:) = srcinfo.r(1,:) + origin_new(1); +srcinfo.r(2,:) = srcinfo.r(2,:) + origin_new(2); + +targinfo.r(1,:) = targinfo.r(1,:) + origin_new(1); +targinfo.r(2,:) = targinfo.r(2,:) + origin_new(2); + +v = get_exact_kernels(zk, srcinfo, targinfo, type); + +err1 = norm(v(:) - submat(:)); +fprintf('Error in shifted kernel %s = %d \n',type,err1); +fprintf('ratios shifted kernel= %d\n', max(abs(v(1:end)./submat(1:end)-1))); + + + + + + + + +function [v] = get_exact_kernels(zk, srcinfo, targinfo, type) + +src = srcinfo.r; +targ = targinfo.r; + +srcnorm = srcinfo.n; +targnorm = targinfo.n; + +[~, ns] = size(src); +[~, nt] = size(targ); + +v = zeros(nt,ns); + +if strcmpi(type,'s') + fker = @(x, s, t, rns, rnt) fslp(x, zk, s, t, rns, rnt); +elseif strcmpi(type, 'd') + fker = @(x, s, t, rns, rnt) fdlp(x, zk, s, t, rns, rnt); +elseif strcmpi(type, 'sprime') + fker = @(x, s, t, rns, rnt) fsprime(x, zk, s, t, rns, rnt); +elseif strcmpi(type, 'dprime') + fker = @(x, s, t, rns, rnt) fdprime(x, zk, s, t, rns, rnt); +elseif strcmpi(type, 'sdiff') + fker = @(x, s, t, rns, rnt) fsdiff(x, zk, s, t, rns, rnt); +elseif strcmpi(type, 'dprimediff') + fker = @(x, s, t, rns, rnt) fdprimediff(x, zk, s, t, rns, rnt); +end + +for j=1:ns + for i=1:nt + r = sqrt((src(1,j) - targ(1,i))^2 + (src(2,j)-targ(2,i))^2); + if r > eps + v(i,j) = integral(@(x) fker(x, src(:,j), targ(:,i), srcnorm(:,j), targnorm(:,i)), 0, 2*pi,'AbsTol',1E-12); + else + v(i,j) = 0; + end + end +end + + +end + + +function f = fslp (x, zk, s, t, rns, rnt) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + r = sqrt(rs.^2 + rt.^2 - 2*rs.*rt.*cos(x) + (zs-zt).^2); + f = exp(1j*zk*r)/4/pi./r.*rs; +end + + + +function f = fdlp (x, zk, s, t, rns, rnt) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + rnd = (rt.*cos(x) - rs).*rns(1) + (zt - zs).*rns(2); + + r = sqrt(rs.^2 + rt.^2 - 2*rs.*rt.*cos(x) + (zs-zt).^2); + f = rnd.*(1 - 1j*zk*r).*exp(1j*zk*r)/4/pi./r.^3.*rs; +end + + + +function f = fsprime (x, zk, s, t, rns, rnt) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + rnd = (rt - rs.*cos(x)).*rnt(1) + (zt - zs).*rnt(2); + + r = sqrt(rs.^2 + rt.^2 - 2*rs.*rt.*cos(x) + (zs-zt).^2); + f = -rnd.*(1-1j*zk*r).*exp(1j*zk*r)/4/pi./r.^3.*rs; +end + + + +function f = fdprime (x, zk, s, t, rns, rnt) + rs = s(1); zs = s(2); + rt = t(1); zt = t(2); + + rndt = (rt - rs.*cos(x)).*rnt(1) + (zt - zs).*rnt(2); + rnds = (rt.*cos(x) - rs).*rns(1) + (zt - zs).*rns(2); + rnsnt = rns(1)*rnt(1).*cos(x) + rns(2)*rnt(2); + + r = sqrt(rs.^2 + rt.^2 - 2*rs.*rt.*cos(x) + (zs-zt).^2); + f = -(rnsnt.*(1j*zk.*r-1).*exp(1j*zk*r)/4/pi./r.^3 + ... + rndt.*rnds.*(-zk^2.*r.^2 - 3*1j*zk.*r + 3).*exp(1j*zk*r)/4/pi./r.^5).*rs; +end + + +function f = fsdiff (x, zk, s, t, rns, rnt) + f = fslp(x, zk, s, t, rns, rnt) - fslp(x, 1j*zk, s, t, rns, rnt); +end + + + +function f = fdprimediff (x, zk, s, t, rns, rnt) + f1 = fdprime(x, zk, s, t, rns, rnt); + f2 = fdprime(x, 1j*zk, s, t, rns, rnt); + f = f1 - f2; +end diff --git a/devtools/test/axissymkernel_sphereTest.m b/devtools/test/axissymkernel_sphereTest.m new file mode 100644 index 0000000..105d43c --- /dev/null +++ b/devtools/test/axissymkernel_sphereTest.m @@ -0,0 +1,217 @@ +clearvars; close all; +iseed = 8675309; +rng(iseed); + +addpaths_loc(); + +zk = 40.1; + +type = 'chnkr-star'; + +pref = []; +pref.k = 16; +ns = 10; +nt = 100; +ppw = 80; % points per wavelength; +maxchunklen = pref.k/ppw/real(zk)*2*pi; +maxchunklen = 0.5; + +[chnkr, ~, ~] = get_geometry(type, pref, ns, nt, maxchunklen); +wts = chnkr.wts; wts = wts(:); + +fprintf('Done building geometry\n'); + +% Plot everything + +figure(1) +clf +hold off +plot(chnkr) +axis equal + +type = 'dprimediff'; +K = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], type); + + + +% mode number +nmode = 20; + +if strcmpi(type, 's') + zfac = 1j*zk*sphericalbesselj(nmode,zk)*sphericalbesselh(nmode,zk); +elseif strcmpi(type, 'd') + zfac = 1j*zk^2*sphericalbesselj_der(nmode,zk)*sphericalbesselh(nmode,zk) - 0.5; +elseif strcmpi(type, 'sprime') + zfac = 1j*zk^2*sphericalbesselj(nmode,zk)*sphericalbesselh_der(nmode,zk) + 0.5; +elseif strcmpi(type, 'dprimediff') + zfac = 1j*zk^3*sphericalbesselj_der(nmode,zk)*sphericalbesselh_der(nmode,zk); + zfac = zfac - 1j*(1j*zk)^3*sphericalbesselj_der(nmode,1j*zk)*sphericalbesselh_der(nmode,1j*zk); +end + + + +cos_t = chnkr.r(2,:)./sqrt(chnkr.r(2,:).^2 + chnkr.r(1,:).^2); +ubdry = legendre(nmode,chnkr.r(2,:)); +ubdry = ubdry(1,:); +ubdry = ubdry(:); + +% Form matrix +opts = []; +tic, A = chunkermat(chnkr, K, opts); toc + +pot = A*ubdry; +potex = zfac*ubdry; + +err1 = sqrt(sum(abs(pot - potex).^2.*wts(:))); +fprintf('error in layer pot on sphere=%d\n', err1); + + + +function [chnkobj, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen) + +if nargin == 0 + type = 'chnkr'; +end + +if nargin <= 1 + pref = []; + pref.k = 16; +end + +if nargin <= 2 + ns = 10; +end + +if nargin <= 3 + nt = 10; +end + + +if nargin <= 4 + maxchunklen = 1.0; +end + + +if strcmpi(type, 'cgrph') + + + nverts = 3; + verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1)); + verts = [real(verts);imag(verts)]; + + + iind = 1:(nverts-1); + jind = 1:(nverts-1); + + iind = [iind iind]; + jind = [jind jind + 1]; + jind(jind>nverts) = 1; + svals = [-ones(1,nverts-1) ones(1,nverts-1)]; + edge2verts = sparse(iind, jind, svals, nverts-1, nverts); + + amp = 0.1; + frq = 2; + fchnks = cell(1,size(edge2verts,1)); + for icurve = 1:size(edge2verts,1) + fchnks{icurve} = @(t) sinearc(t, amp, frq); + end + cparams = []; + cparams.nover = 2; + cparams.maxchunklen = maxchunklen; + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + ts = 0.0+2*pi*rand(ns,1); + sources = 3.0*[cos(ts)';sin(ts)']; + + ts = 0.0+2*pi*rand(nt,1); + targets = 0.2*[cos(ts)'; sin(ts)']; + + +elseif strcmpi(type,'chnkr-star') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = false; + cparams.ta = -pi/2; + cparams.tb = pi/2; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.0; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.5*sources; + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + +elseif strcmpi(type,'chnkr-torus') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = true; + cparams.ta = 0; + cparams.tb = 2*pi; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + ctr = [3 0]; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.5*sources; + sources(1,:) = sources(1,:) + ctr(1); + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + targets(1,:) = targets(1,:) + ctr(1); +end + + + +end + + +function [r,d,d2] = sinearc(t,amp,frq) +xs = t; +ys = amp*sin(frq*t); +xp = ones(size(t)); +yp = amp*frq*cos(frq*t); +xpp = zeros(size(t)); +ypp = -frq*frq*amp*sin(t); + +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end + +function jn = sphericalbesselj(n,z) + jn = sqrt(pi/2./z).*besselj(n+0.5,z); +end + + +function jnp = sphericalbesselj_der(n,z) + jnp = 0.5*(sphericalbesselj(n-1,z) - (sphericalbesselj(n,z) + ... + z.*sphericalbesselj(n+1,z))./z); +end + + +function hn = sphericalbesselh(n,z) + hn = sqrt(pi/2./z).*besselh(n+0.5,1,z); +end + + +function hnp = sphericalbesselh_der(n,z) + hnp = 0.5*(sphericalbesselh(n-1,z) - (sphericalbesselh(n,z) + ... + z.*sphericalbesselh(n+1,z))./z); +end + + diff --git a/devtools/test/chunkermat_axissymhelm2dTest.m b/devtools/test/chunkermat_axissymhelm2dTest.m new file mode 100644 index 0000000..3b4985b --- /dev/null +++ b/devtools/test/chunkermat_axissymhelm2dTest.m @@ -0,0 +1,143 @@ + +%CHUNKERMAT_AXISSYMHELM2DTEST +% +% test the matrix builder and do a basic solve + +clearvars; close all; +iseed = 8675309; +rng(iseed); + +addpaths_loc(); + +iftorus = 0; + +cparams = []; +cparams.eps = 1.0e-10; +cparams.nover = 1; + +if ~iftorus + cparams.ta = -pi/2; + cparams.tb = pi/2; + ctr = [0 0]; + cparams.ifclosed = false; +else + cparams.ta = 0; + cparams.tb = 2*pi; + ctr = [3 0]; + cparams.ifclosed = true; +end + + +cparams.maxchunklen = 0.5; +pref = []; +pref.k = 16; +narms = 0; +amp = 0.25; + + +start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp,ctr),cparams,pref); +t1 = toc(start); + +fprintf('%5.2e s : time to build geo\n',t1) + +zk = 80.1; + +% sources + +ns = 10; +ts = -pi/2+pi*rand(ns,1); +sources = starfish(ts,narms,amp); +sources = 0.5*sources; +sources(1,:) = sources(1,:) + ctr(1); +strengths = randn(ns,1); + +% targets + +nt = 100; +ts = -pi/2+pi*rand(nt,1); +targets = starfish(ts,narms,amp); +targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); +targets(1,:) = targets(1,:) + ctr(1); + +% plot geo and sources + +xs = chnkr.r(1,:,:); xmin = min(xs(:)); xmax = max(xs(:)); +ys = chnkr.r(2,:,:); ymin = min(ys(:)); ymax = max(ys(:)); + +figure(1) +clf +hold off +plot(chnkr) +hold on +scatter(sources(1,:),sources(2,:),'o') +scatter(targets(1,:),targets(2,:),'x') +axis equal + + +% + +kerns = @(s,t) chnk.axissymhelm2d.kern(zk,s,t, [0,0], 's'); + +% eval u on bdry + +targs = chnkr.r; targs = reshape(targs,2,chnkr.k*chnkr.nch); +targstau = tangents(chnkr); +targstau = reshape(targstau,2,chnkr.k*chnkr.nch); + +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = targs; +kernmats = kerns(srcinfo,targinfo); +ubdry = kernmats*strengths; + +% eval u at targets + +targinfo = []; targinfo.r = targets; +kernmatstarg = kerns(srcinfo,targinfo); +utarg = kernmatstarg*strengths; + +% + +% build laplace dirichlet matrix + +fkern = kernel('axissymh', 'c', zk); +fkern = @(s,t) chnk.axissymhelm2d.kern(zk,s,t,[0 0],'c', [1 -1j]); +start = tic; D = chunkermat(chnkr,fkern); +t1 = toc(start); + +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = 0.5*eye(chnkr.k*chnkr.nch) + D; + +rhs = ubdry; rhs = rhs(:); + +start = tic; sol = gmres(sys,rhs,[],1e-14,200); t1 = toc(start); + +fprintf('%5.2e s : time for dense gmres\n',t1) + +start = tic; sol2 = sys\rhs; t1 = toc(start); + +fprintf('%5.2e s : time for dense backslash solve\n',t1) + +err = norm(sol-sol2,'fro')/norm(sol2,'fro'); + +fprintf('difference between direct and iterative %5.2e\n',err) + +% evaluate at targets and compare + +opts.usesmooth=false; +opts.verb=false; +opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; +start=tic; Dsol = chunkerkerneval(chnkr,fkern,sol2,targets,opts); +t1 = toc(start); +fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) + +% + +wchnkr = chnkr.wts; + +relerr = norm(utarg-Dsol,'fro')/(sqrt(chnkr.nch)*norm(utarg,'fro')); +relerr2 = norm(utarg-Dsol,'inf')/dot(abs(sol(:)),wchnkr(:)); +fprintf('relative frobenius error %5.2e\n',relerr); +fprintf('relative l_inf/l_1 error %5.2e\n',relerr2); + +assert(relerr < 1e-7); diff --git a/devtools/test/chunkermat_axissymhelm_neumannTest.m b/devtools/test/chunkermat_axissymhelm_neumannTest.m new file mode 100644 index 0000000..335b21e --- /dev/null +++ b/devtools/test/chunkermat_axissymhelm_neumannTest.m @@ -0,0 +1,399 @@ +clearvars; close all; +iseed = 8675309; +rng(iseed); + +addpaths_loc(); + +zk = 1.1; + +type = 'chnkr-star'; +% type = 'chnkr-torus'; +% type = 'cgrph'; +% type = 'cgrph-sphere'; + +irep = 'rpcomb'; +% irep = 'sk'; + +pref = []; +pref.k = 16; +ns = 10; +nt = 100; +ppw = 10; % points per wavelength; +maxchunklen = pref.k/ppw/real(zk)*2*pi; + +[chnkr, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen); +wts = chnkr.wts; wts = wts(:); + +l2scale = false; +fprintf('Done building geometry\n'); + +% source strengths +strengths = randn(ns, 1); + +% targets + + +% Plot everything + +figure(1) +clf +hold off +plot(chnkr, 'k.'); +hold on +quiver(chnkr); +scatter(sources(1,:), sources(2,:), 'o') +scatter(targets(1,:), targets(2,:), 'x') +axis equal + + +% For solving exterior the Neumann boundary value problem, we use the +% following integral equation +% +% u = \beta(S_{k} + i \alpha D_{k} S_{ik}) \sigma +% +% with \beta = -1.0/(0.5 + 1i*0.25*alpha) +% +% On imposing the boundary conditions, we get the following integral +% equation +% +% du/dn = I + \beta S_{k}'[\sigma] + +% i\beta \alpha(D'_{k} - D'_{ik})(S_{ik}[\sigma]) + +% i\beta \alpha(S_{ik}')^2 [\sigma]; +% +% Setting -S_{ik}[\sigma] + \tau = 0, and - S_{ik}'[\sigma] + \mu=0, +% we have the following system of integral equations + +% Set up kernels +Skp = kernel('axissymhelm', 'sprime', zk); +Sk = kernel('axissymhelm', 's', zk); +Dk = kernel('axissymhelm', 'd', zk); + +Z = kernel.zeros(); + + +if strcmpi(irep,'rpcomb') + Sik = kernel('axissymhelm', 's', 1i*zk); + Sikp = kernel('axissymhelm', 'sprime', 1i*zk); + Dkdiff = kernel('axissymhelmdiff', 'dprime', [zk 1i*zk]); + alpha = 1; + c1 = -1/(0.5 + 1i*alpha*0.25); + c2 = -1i*alpha/(0.5 + 1i*alpha*0.25); + c3 = -1; + K = [ c1*Skp c2*Dkdiff c2*Sikp ; + c3*Sik Z Z ; + c3*Sikp Z Z ]; + K = kernel(K); + Keval = c1*kernel([Sk 1i*alpha*Dk Z]); +else + K = -2*Skp; + Keval = -2*Sk; +end + +% Set up boundary data + +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:); +kernmats = Skp.eval(srcinfo, targinfo); +ubdry = kernmats*strengths; + +% rvals = chnkr.r(1,:); +% ubdry = ubdry.*rvals(:); + +npts = chnkr.npt; +nsys = K.opdims(1)*npts; +rhs = zeros(nsys, 1); + + +if(l2scale) + rhs(1:K.opdims(1):end) = ubdry.*sqrt(wts); +else + rhs(1:K.opdims(1):end) = ubdry; +end + +% Form matrix +opts = []; +opts.l2scale = l2scale; +opts.rcip = true; +opts.nsub_or_tol = 20; +tic, A = chunkermat(chnkr, K, opts) + eye(nsys); toc +start = tic; +sol = gmres(A, rhs, [], 1e-14, 200); +t1 = toc(start); + +% Compute exact solution +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = targets; +kernmatstarg = Sk.eval(srcinfo, targinfo); +utarg = kernmatstarg*strengths; + +% Compute solution using chunkerkerneval +% evaluate at targets and compare + +opts.forcesmooth = false; +opts.verb = false; +opts.quadkgparams = {'RelTol', 1e-8, 'AbsTol', 1.0e-8}; + +if(l2scale) + wts_rep = repmat(wts(:).', K.opdims(1),1); + wts_rep = wts_rep(:); + sol = sol./sqrt(wts_rep); +end + +if isa(chnkr, 'chunkgraph') + % Collapse cgrph into chnkrtotal + chnkrs = chnkr.echnks; + chnkrtotal = merge(chnkrs); +else + chnkrtotal = chnkr; +end + + + +start = tic; +Dsol = chunkerkerneval(chnkrtotal, Keval, sol, targets, opts); +t2 = toc(start); +fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n', t2) + + +wchnkr = chnkrtotal.wts; +wchnkr = repmat(wchnkr(:).', K.opdims(1), 1); +relerr = norm(utarg-Dsol) / (sqrt(chnkrtotal.nch)*norm(utarg)); +relerr2 = norm(utarg-Dsol, 'inf') / dot(abs(sol(:)), wchnkr(:)); +fprintf('relative frobenius error %5.2e\n', relerr); +fprintf('relative l_inf/l_1 error %5.2e\n', relerr2); + +return + + + + +% Test fast direct solver interfaces + +% build sparse tridiag part + +opts.nonsmoothonly = true; +opts.rcip = true; +start = tic; spmat = chunkermat(chnkr, K, opts); t1 = toc(start); +fprintf('%5.2e s : time to build tridiag\n',t1) + + +spmat = spmat + speye(nsys); + +% test matrix entry evaluator +start = tic; +opdims = K.opdims; +sys2 = chnk.flam.kernbyindex(1:nsys, 1:nsys, chnkr, K, opdims, ... + spmat, l2scale); + + +t1 = toc(start); + +fprintf('%5.2e s : time for mat entry eval on whole mat\n',t1) + +err2 = norm(sys2-A,'fro')/norm(A,'fro'); +fprintf('%5.2e : fro error of build \n',err2); + +% test fast direct solver +opts.ifproxy = false; +F = chunkerflam(chnkr,K,1.0,opts); + +start = tic; sol2 = rskelf_sv(F,rhs); t1 = toc(start); + +if(l2scale) + wts_rep = repmat(wts(:).', K.opdims(1),1); + wts_rep = wts_rep(:); + sol2 = sol2./sqrt(wts_rep); +end + + +fprintf('%5.2e s : time for rskelf_sv \n',t1) + +err = norm(sol-sol2,'fro')/norm(sol,'fro'); + +fprintf('difference between fast-direct and iterative %5.2e\n',err) + + +function [chnkobj, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen) + +if nargin == 0 + type = 'chnkr'; +end + +if nargin <= 1 + pref = []; + pref.k = 16; +end + +if nargin <= 2 + ns = 10; +end + +if nargin <= 3 + nt = 10; +end + + +if nargin <= 4 + maxchunklen = 1.0; +end + + +if strcmpi(type, 'cgrph') + + + nverts = 3; + verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1)); + verts = [real(verts);imag(verts)]; + + + iind = 1:(nverts-1); + jind = 1:(nverts-1); + + iind = [iind iind]; + jind = [jind jind + 1]; + jind(jind>nverts) = 1; + svals = [-ones(1,nverts-1) ones(1,nverts-1)]; + edge2verts = sparse(iind, jind, svals, nverts-1, nverts); + + amp = 0.1; + frq = 2; + fchnks = cell(1,size(edge2verts,1)); + for icurve = 1:size(edge2verts,1) + fchnks{icurve} = @(t) sinearc(t, amp, frq); + end + cparams = []; + cparams.nover = 2; + cparams.maxchunklen = maxchunklen; + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + ts = -pi/2 + pi*rand(ns,1); + sources = 0.2*[cos(ts)';sin(ts)']; + + ts = -pi/2 + pi*rand(nt,1); + targets = 3.0*[cos(ts)'; sin(ts)']; + + + +elseif strcmpi(type,'cgrph-sphere') + + + nverts = 2; + verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1)); + verts = [real(verts);imag(verts)]; + + + iind = 1:(nverts-1); + jind = 1:(nverts-1); + + iind = [iind iind]; + jind = [jind jind + 1]; + jind(jind>nverts) = 1; + svals = [-ones(1,nverts-1) ones(1,nverts-1)]; + edge2verts = sparse(iind, jind, svals, nverts-1, nverts); + + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = false; + cparams.ta = -pi/2; + cparams.tb = pi/2; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.0; + + fchnks{1} = @(t) circle(t); + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.1*sources; + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + + + +elseif strcmpi(type,'chnkr-star') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = false; + cparams.ta = -pi/2; + cparams.tb = pi/2; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.5*sources; + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + +elseif strcmpi(type,'chnkr-torus') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = true; + cparams.ta = 0; + cparams.tb = 2*pi; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + ctr = [3 0]; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = 0.5*sources; + sources(1,:) = sources(1,:) + ctr(1); + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); + targets(1,:) = targets(1,:) + ctr(1); +end + + + +end + + +function [r,d,d2] = sinearc(t,amp,frq) +xs = t; +ys = amp*sin(frq*t); +xp = ones(size(t)); +yp = amp*frq*cos(frq*t); +xpp = zeros(size(t)); +ypp = -frq*frq*amp*sin(t); + +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end + + +function [r,d,d2] = circle(t) +xs = cos(-pi/4 + pi/2*t); +ys = sin(-pi/4 + pi/2*t); +xp = -pi*ys/2; +yp = pi*xs/2; +xpp = -pi*pi*xs/4; +ypp = -pi*pi*ys/4; +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end + + + diff --git a/devtools/test/chunkermat_axissymhelm_transmissionTest.m b/devtools/test/chunkermat_axissymhelm_transmissionTest.m new file mode 100644 index 0000000..23d3a3b --- /dev/null +++ b/devtools/test/chunkermat_axissymhelm_transmissionTest.m @@ -0,0 +1,349 @@ +clearvars; close all; +iseed = 8675309; +rng(iseed); + +addpaths_loc(); + +zk = 10.1; + +type = 'chnkr-star'; +% type = 'chnkr-torus'; + +pref = []; +pref.k = 16; +ns = 10; +nt = 10; +ppw = 80; % points per wavelength; +maxchunklen = pref.k/ppw/real(zk)*2*pi; +maxchunklen = 0.5; + +[chnkr, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen); +wts = chnkr.wts; wts = wts(:); + +l2scale = false; +fprintf('Done building geometry\n'); + +% source strengths +strengths_in = randn(ns, 1); +strengths_out = randn(nt, 1); + + +% targets + + +% Plot everything + +figure(1) +clf +hold off +plot(chnkr) +hold on +scatter(sources(1,:), sources(2,:), 'o') +scatter(targets(1,:), targets(2,:), 'x') +axis equal + + +% For solving the transmission boundary value problem, we +% use the repesentation +% u_{i} = D_{k_{i}} [\sigma] - S_{k_{i}}[\tau] +% +% and impose jumps in u and du/dn as exterior-interior respectively +% + +% Set up kernels +Skdiff = kernel('axissymhelmdiff', 's', [zk 1i*zk]); +Dkdiff = kernel('axissymhelmdiff', 'd', [zk 1i*zk]); +Skpdiff = kernel('axissymhelmdiff', 'sprime', [zk 1i*zk]); +Dkpdiff = kernel('axissymhelmdiff', 'dprime', [zk 1i*zk]); + +K = [Dkdiff -1*Skdiff; + Dkpdiff -1*Skpdiff]; +K = kernel(K); +Dk = kernel('axissymhelm', 'd', zk); +Sk = kernel('axissymhelm', 's', zk); +Skp = kernel('axissymhelm', 'sprime', zk); + +Dik = kernel('axissymhelm', 'd', 1i*zk); +Sik = kernel('axissymhelm', 's', 1i*zk); +Sikp = kernel('axissymhelm', 'sprime', 1i*zk); + +Kouteval = kernel([Dk, -1*Sk]); +Kineval = kernel([Dik, -1*Sik]); + + +% Set up boundary data + +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:); +kernmats = Sk.eval(srcinfo, targinfo); +kernmatsp = Skp.eval(srcinfo, targinfo); +ubdry = kernmats*strengths_in; +dudnbdry = kernmatsp*strengths_in; + + + +srcinfo = []; srcinfo.r = targets; +targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.n = chnkr.n(:,:); +kernmats = Sik.eval(srcinfo, targinfo); +kernmatsp = Sikp.eval(srcinfo, targinfo); +ubdry = ubdry - kernmats*strengths_out; +dudnbdry = dudnbdry - kernmatsp*strengths_out; + + +npts = chnkr.npt; +nsys = K.opdims(1)*npts; +rhs = zeros(nsys, 1); + + +if(l2scale) + rhs(1:2:end) = ubdry.*sqrt(wts); + rhs(2:2:end) = dudnbdry.*sqrt(wts); +else + rhs(1:2:end) = ubdry; + rhs(2:2:end) = dudnbdry; +end + +% Form matrix +opts = []; +opts.l2scale = l2scale; +tic, A = chunkermat(chnkr, K, opts) + eye(nsys); toc +start = tic; +sol = gmres(A, rhs, [], 1e-14, 200); +t1 = toc(start); + +% Compute exact solution +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = targets; +kernmatstarg = Sk.eval(srcinfo, targinfo); +utarg_out = kernmatstarg*strengths_in; + +% Compute exact solution +srcinfo = []; srcinfo.r = targets; +targinfo = []; targinfo.r = sources; +kernmatstarg = Sik.eval(srcinfo, targinfo); +utarg_in = kernmatstarg*strengths_out; + + + +% Compute solution using chunkerkerneval +% evaluate at targets and compare + +opts.forceadap = true; +opts.verb = false; +opts.quadkgparams = {'RelTol', 1e-8, 'AbsTol', 1.0e-8}; + + +if(l2scale) + wts_rep = repmat(wts(:).', K.opdims(1),1); + wts_rep = wts_rep(:); + sol = sol./sqrt(wts_rep); +end + +if isa(chnkr, 'chunkgraph') + % Collapse cgrph into chnkrtotal + chnkrs = chnkr.echnks; + chnkrtotal = merge(chnkrs); +else + chnkrtotal = chnkr; +end + + + +start = tic; +Dsol_out = chunkerkerneval(chnkrtotal, Kouteval, sol, targets, opts); +Dsol_in = chunkerkerneval(chnkrtotal, Kineval, sol, sources, opts); +t2 = toc(start); +fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n', t2) + + +wchnkr = chnkrtotal.wts; +wchnkr = repmat(wchnkr(:).', K.opdims(1), 1); +relerr = norm(utarg_out-Dsol_out) / (sqrt(chnkrtotal.nch)*norm(utarg_out)); +relerr2 = norm(utarg_out-Dsol_out, 'inf') / dot(abs(sol(:)), wchnkr(:)); +fprintf('relative frobenius error in exterior %5.2e\n', relerr); +fprintf('relative l_inf/l_1 error in exterior %5.2e\n', relerr2); + +relerr = norm(utarg_in-Dsol_in) / (sqrt(chnkrtotal.nch)*norm(utarg_in)); +relerr2 = norm(utarg_in-Dsol_in, 'inf') / dot(abs(sol(:)), wchnkr(:)); +fprintf('relative frobenius error in interior %5.2e\n', relerr); +fprintf('relative l_inf/l_1 error in interior %5.2e\n', relerr2); + + +return + + + + +% Test fast direct solver interfaces + +% build sparse tridiag part + + +opts.nonsmoothonly = true; +opts.rcip = true; +start = tic; spmat = chunkermat(chnkr, K, opts); t1 = toc(start); +fprintf('%5.2e s : time to build tridiag\n',t1) + + +spmat = spmat + speye(nsys); + +% test matrix entry evaluator +start = tic; +opdims = K.opdims; +sys2 = chnk.flam.kernbyindex(1:nsys, 1:nsys, chnkr, K, opdims, ... + spmat, l2scale); + + +t1 = toc(start); + +fprintf('%5.2e s : time for mat entry eval on whole mat\n',t1) + +err2 = norm(sys2-A,'fro')/norm(A,'fro'); +fprintf('%5.2e : fro error of build \n',err2); + +% test fast direct solver +opts.ifproxy = false; +F = chunkerflam(chnkr,K,1.0,opts); + +start = tic; sol2 = rskelf_sv(F,rhs); t1 = toc(start); + +if(l2scale) + wts_rep = repmat(wts(:).', K.opdims(1),1); + wts_rep = wts_rep(:); + sol2 = sol2./sqrt(wts_rep); +end + + +fprintf('%5.2e s : time for rskelf_sv \n',t1) + +err = norm(sol-sol2,'fro')/norm(sol,'fro'); + +fprintf('difference between fast-direct and iterative %5.2e\n',err) + + +function [chnkobj, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen) + +if nargin == 0 + type = 'chnkr'; +end + +if nargin <= 1 + pref = []; + pref.k = 16; +end + +if nargin <= 2 + ns = 10; +end + +if nargin <= 3 + nt = 10; +end + + +if nargin <= 4 + maxchunklen = 1.0; +end + + +if strcmpi(type, 'cgrph') + + + nverts = 3; + verts = exp(-1i*pi/2 + 1i*pi*(0:(nverts-1))/(nverts-1)); + verts = [real(verts);imag(verts)]; + + + iind = 1:(nverts-1); + jind = 1:(nverts-1); + + iind = [iind iind]; + jind = [jind jind + 1]; + jind(jind>nverts) = 1; + svals = [-ones(1,nverts-1) ones(1,nverts-1)]; + edge2verts = sparse(iind, jind, svals, nverts-1, nverts); + + amp = 0.1; + frq = 2; + fchnks = cell(1,size(edge2verts,1)); + for icurve = 1:size(edge2verts,1) + fchnks{icurve} = @(t) sinearc(t, amp, frq); + end + cparams = []; + cparams.nover = 2; + cparams.maxchunklen = maxchunklen; + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + ts = 0.0+2*pi*rand(ns,1); + sources = 3.0*[cos(ts)';sin(ts)']; + + ts = 0.0+2*pi*rand(nt,1); + targets = 0.2*[cos(ts)'; sin(ts)']; + + +elseif strcmpi(type,'chnkr-star') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = false; + cparams.ta = -pi/2; + cparams.tb = pi/2; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = sources .* (0.7 - 0.4.*repmat(rand(1,ns), 2, 1)); + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1.2 + 3.0*repmat(rand(1, nt), 2, 1)); + +elseif strcmpi(type,'chnkr-torus') + cparams = []; + cparams.eps = 1.0e-10; + cparams.nover = 1; + cparams.ifclosed = true; + cparams.ta = 0; + cparams.tb = 2*pi; + cparams.maxchunklen = maxchunklen; + narms = 0; + amp = 0.25; + ctr = [3 0]; + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); + chnkobj = sort(chnkobj); + + ts = -pi/2 + pi*rand(ns, 1); + sources = starfish(ts, narms, amp); + sources = sources .* (0.7 - 0.4.*repmat(rand(1,ns), 2, 1)); + sources(1,:) = sources(1,:) + ctr(1); + + ts = -pi/2 + pi*rand(nt, 1); + targets = starfish(ts, narms, amp); + targets = targets .* (1.2 + 3.0*repmat(rand(1, nt), 2, 1)); + targets(1,:) = targets(1,:) + ctr(1); +end + + + +end + + +function [r,d,d2] = sinearc(t,amp,frq) +xs = t; +ys = amp*sin(frq*t); +xp = ones(size(t)); +yp = amp*frq*cos(frq*t); +xpp = zeros(size(t)); +ypp = -frq*frq*amp*sin(t); + +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end + diff --git a/devtools/test/flamutilitiesTest.m b/devtools/test/flamutilitiesTest.m index 6290065..6e13270 100644 --- a/devtools/test/flamutilitiesTest.m +++ b/devtools/test/flamutilitiesTest.m @@ -33,8 +33,6 @@ cparams.nover = 2; [cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams); -vstruc = procverts(cgrph); -rgns = findregions(cgrph); cgrph = balance(cgrph);