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/README.md b/README.md index 280853e..d4b11db 100644 --- a/README.md +++ b/README.md @@ -80,7 +80,7 @@ the library will not function without FLAM and its subdirectories included in th chunkie has benefitted from the contributions of several developers: Travis Askham, Manas Rachh, Michael O'Neil, Jeremy Hoskins, Dan Fortunato, Shidong Jiang, -Fredrik Fryklund, Hai Yang Wang, and Hai Zhu. +Fredrik Fryklund, Hai Yang Wang, Hai Zhu, and Tristan Goodwill. James Bremer and Zydrunas Gimbutas provided generalized Gaussian quadrature rules (chunkie/+chnk/+quadggq) 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/+quadadap/buildmat.m b/chunkie/+chnk/+quadadap/buildmat.m index 6946598..ec42364 100644 --- a/chunkie/+chnk/+quadadap/buildmat.m +++ b/chunkie/+chnk/+quadadap/buildmat.m @@ -22,7 +22,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; n = chnkr.n; [t,wts,u] = lege.exps(k); @@ -77,7 +76,7 @@ dt = d(:,:,ibefore); d2t = d2(:,:,ibefore); nt = n(:,:,ibefore); - submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ... + submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ... kern,opdims,t2,w2); imat = 1 + (ibefore-1)*k*opdims(1); @@ -90,7 +89,7 @@ dt = d(:,:,iafter); nt = n(:,:,iafter); d2t = d2(:,:,iafter); - submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ... + submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ... kern,opdims,t2,w2); imat = 1 + (iafter-1)*k*opdims(1); @@ -101,7 +100,7 @@ % self - submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,[],i,kern,opdims,... + submat = chnk.quadggq.diagbuildmat(r,d,n,d2,[],i,kern,opdims,... xs0,wts0,ainterps0kron,ainterps0); imat = 1 + (i-1)*k*opdims(1); @@ -163,7 +162,7 @@ d2t = d2(:,targfix); nt = n(:,targfix); - submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ... + submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ... kern,opdims,t2,w2); imats = bsxfun(@plus,(1:opdims(1)).',opdims(1)*(targfix(:)-1).'); diff --git a/chunkie/+chnk/+quadba/buildmattd.m b/chunkie/+chnk/+quadba/buildmattd.m index dfe942e..07069d5 100644 --- a/chunkie/+chnk/+quadba/buildmattd.m +++ b/chunkie/+chnk/+quadba/buildmattd.m @@ -11,7 +11,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; [~,~,u] = lege.exps(k); diff --git a/chunkie/+chnk/+quadggq/buildmat.m b/chunkie/+chnk/+quadggq/buildmat.m index 9470819..2e24f04 100644 --- a/chunkie/+chnk/+quadggq/buildmat.m +++ b/chunkie/+chnk/+quadggq/buildmat.m @@ -50,7 +50,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; n = chnkr.n; data = []; @@ -101,7 +100,7 @@ if ~isempty(ilist) && ismember(ibefore,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,ibefore,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,ibefore,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); imat = 1 + (ibefore-1)*k*opdims(1); @@ -115,7 +114,7 @@ if ~isempty(ilist) && ismember(iafter,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,iafter,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,iafter,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); imat = 1 + (iafter-1)*k*opdims(1); @@ -129,7 +128,7 @@ if ~isempty(ilist) && ismember(j,ilist) % skip construction if the chunk is in the "bad" chunk list else - submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,data,j,kern,opdims,... + submat = chnk.quadggq.diagbuildmat(r,d,n,d2,data,j,kern,opdims,... xs0,wts0,ainterps0kron,ainterps0); imat = 1 + (j-1)*k*opdims(1); diff --git a/chunkie/+chnk/+quadggq/buildmattd.m b/chunkie/+chnk/+quadggq/buildmattd.m index fd820a6..ee5e01c 100644 --- a/chunkie/+chnk/+quadggq/buildmattd.m +++ b/chunkie/+chnk/+quadggq/buildmattd.m @@ -22,7 +22,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; n = chnkr.n; data = []; @@ -89,7 +88,7 @@ if ~isempty(ilist) && ismember(ibefore,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,ibefore,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,ibefore,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); imat = 1 + (ibefore-1)*k*opdims(1); @@ -106,7 +105,7 @@ if ~isempty(ilist) && ismember(iafter,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,iafter,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,iafter,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); @@ -124,7 +123,7 @@ if ~isempty(ilist) && ismember(j,ilist) % skip construction if the chunk is in the "bad" chunk list else - submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,data,j,kern,opdims,... + submat = chnk.quadggq.diagbuildmat(r,d,n,d2,data,j,kern,opdims,... xs0,wts0,ainterps0kron,ainterps0); imat = 1 + (j-1)*k*opdims(1); diff --git a/chunkie/+chnk/+quadggq/diagbuildmat.m b/chunkie/+chnk/+quadggq/diagbuildmat.m index 8c405ef..b8ceabf 100644 --- a/chunkie/+chnk/+quadggq/diagbuildmat.m +++ b/chunkie/+chnk/+quadggq/diagbuildmat.m @@ -1,10 +1,10 @@ -function submat = diagbuildmat(r,d,n,d2,h,data,i,fkern,opdims,... +function submat = diagbuildmat(r,d,n,d2,data,i,fkern,opdims,... xs0,whts0,ainterps0kron,ainterps0) %CHNK.QUADGGQ.DIAGBUILDMAT % grab specific boundary data -rs = r(:,:,i); ds = d(:,:,i); d2s = d2(:,:,i); hs = h(i); +rs = r(:,:,i); ds = d(:,:,i); d2s = d2(:,:,i); ns = n(:,:,i); if(isempty(data)) dd = []; @@ -29,7 +29,7 @@ d2fine{i} = (ainterps0{i}*(d2s.')).'; dfinenrm = sqrt(sum(dfine{i}.^2,1)); nfine{i} = [dfine{i}(2,:); -dfine{i}(1,:)]./dfinenrm; - dsdt{i} = (dfinenrm(:)).*whts0{i}*hs; + dsdt{i} = (dfinenrm(:)).*whts0{i}; end srcinfo = []; diff --git a/chunkie/+chnk/+quadggq/nearbuildmat.m b/chunkie/+chnk/+quadggq/nearbuildmat.m index 1bf815c..08502d1 100644 --- a/chunkie/+chnk/+quadggq/nearbuildmat.m +++ b/chunkie/+chnk/+quadggq/nearbuildmat.m @@ -1,4 +1,4 @@ -function submat = nearbuildmat(r,d,n,d2,h,data,i,j,fkern,opdims,... +function submat = nearbuildmat(r,d,n,d2,data,i,j,fkern,opdims,... xs1,whts1,ainterp1kron,ainterp1) %CHNKR.QUADGGQ.NEARBUILDMAT @@ -6,7 +6,7 @@ rs = r(:,:,j); ds = d(:,:,j); d2s = d2(:,:,j); ns = n(:,:,j); rt = r(:,:,i); dt = d(:,:,i); d2t = d2(:,:,i); nt = n(:,:,i); -hs = h(j); + if(isempty(data)) dd = []; else @@ -64,7 +64,7 @@ dfinenrm = sqrt(sum(dfine.^2,1)); %dfinenrm = dfine(1,:,:); % for complex contour, by SJ 09/30/21 -dsdt = dfinenrm(:).*whts1(:)*hs; +dsdt = dfinenrm(:).*whts1(:); dsdtndim2 = repmat(dsdt(:).',opdims(2),1); dsdtndim2 = dsdtndim2(:); diff --git a/chunkie/+chnk/+quadnative/buildmat.m b/chunkie/+chnk/+quadnative/buildmat.m index 7dbbd01..f4e9da2 100644 --- a/chunkie/+chnk/+quadnative/buildmat.m +++ b/chunkie/+chnk/+quadnative/buildmat.m @@ -21,7 +21,6 @@ r = chnkr.rstor; d = chnkr.dstor; d2 = chnkr.d2stor; -h = chnkr.hstor; n = chnkr.nstor; [dim,k,~] = size(r); @@ -34,10 +33,10 @@ srcinfo = []; srcinfo.r = rs; srcinfo.d = ds; srcinfo.d2 = d2s; srcinfo.n = ns; targinfo = []; targinfo.r = rt; targinfo.d = dt; targinfo.d2 = d2t; targinfo.n = nt; -hs = h(j); ht = h(i); + dsnrms = sqrt(sum(ds.^2,1)); -ws = kron(hs(:),wts(:)); +ws = repmat(wts(:),length(j),1); dsdt = dsnrms(:).*ws; diff --git a/chunkie/+chnk/+rcip/Rcomp.m b/chunkie/+chnk/+rcip/Rcomp.m deleted file mode 100644 index 4b75f6d..0000000 --- a/chunkie/+chnk/+rcip/Rcomp.m +++ /dev/null @@ -1,68 +0,0 @@ -function [R]=Rcomp(ngl,nedge,ndim,Pbc,PWbc,nsub,... - starL,circL,starS,circS,ilist,... - h0,isstart,fcurve,glxs,fkern) - % carry out the forward recursion for computing the preconditioner R - % in the RCIP method - % - % A more general version of Shidong's rcip version - % - % Function is passed as a handle, number of equations is given by - % ndim - % - % Kernel on input takes in arguments (chnkrlocal,ilistl); - % - % Note that matrix must be scaled to have identity on the diagonal, - % will not work with scaled version of identity - - - pref = []; - pref.k = ngl; - % size of the system matrix - nsys = 3*ngl*nedge*ndim; - - % size of the preconditioner R - nR = 2*ngl*nedge*ndim; - - ts = zeros(4,nedge); - chnkrlocal(1,nedge) = chunker(); - - for level=1:nsub - h = h0/2^(nsub-level); - - for i=1:nedge - if isstart(i) - ts(:,i) = [0, 0.5, 1, 2]*h(i); - else - ts(:,i) = -[2, 1, 0.5, 0]*h(i); - end - end - % construct local chunks around the corner - for i=1:nedge - chnkrlocal(i) = chnk.rcip.chunkerfunclocal(fcurve{i},ts(:,i),pref,glxs); - end -% figure -% clf -% plot(chnkrlocal,'r.') -% axis equal -% pause - % construct the system matrix for local chunks - if level == 1 - ilistl = []; - else - ilistl = ilist; - end - - opts = []; - % test for opdims ~= [1,1] - MAT = chunkermat(chnkrlocal,fkern,opts,ilistl); - - % - MAT = eye(nsys) + MAT; - if level==1 % Dumb and lazy initializer for R, for now - %R=eye(nR); - R = inv(MAT(starL,starL)); - end - R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS); - end - - end diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m index e106bb1..438cd45 100644 --- a/chunkie/+chnk/+rcip/Rcompchunk.m +++ b/chunkie/+chnk/+rcip/Rcompchunk.m @@ -106,7 +106,6 @@ d2 = chnkri.d2(:,:,ie); il = chnkri.adj(1,ie); ir = chnkri.adj(2,ie); - h = chnkri.h(ie); if (il > 0 && ir < 0) nextchunk(i) = il; ileftright(i) = 1; @@ -116,8 +115,8 @@ rcs(:,:,i) = sbcrmat*(r.'); dcs(:,:,i) = u*(d.'); d2cs(:,:,i) = u*(d2.'); - dscal(i) = h*2; - d2scal(i) = h^2*4; + dscal(i) = 2; + d2scal(i) = 4; elseif (il < 0 && ir > 0) nextchunk(i) = ir; ileftright(i) = -1; @@ -126,8 +125,8 @@ rcs(:,:,i) = sbclmat*(r.'); dcs(:,:,i) = u*(d.'); d2cs(:,:,i) = u*(d2.'); - dscal(i) = h*2; - d2scal(i) = h^2*4; + dscal(i) = 2; + d2scal(i) = 4; else error('RCIP: edge chunk not adjacent to one vertex and one neighbor') end @@ -156,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 @@ -197,7 +226,7 @@ chnkrlocal(i).r(:,:,nchi+1) = chnkr(ic).r(:,:,nc)-ctr(:,i); chnkrlocal(i).d(:,:,nchi+1) = chnkr(ic).d(:,:,nc); chnkrlocal(i).d2(:,:,nchi+1) = chnkr(ic).d2(:,:,nc); - chnkrlocal(i).h(nchi+1) = chnkr(ic).h(nc); + if ileftright(i) == -1 chnkrlocal(i).adj(1,nchi+1) = nchi; chnkrlocal(i).adj(2,nchi+1) = -1; @@ -212,6 +241,7 @@ chnkrlocal(i).wts = weights(chnkrlocal(i)); end end + % construct the system matrix for local chunks if level == 1 @@ -221,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/+chnk/+rcip/chunkerfunclocal.m b/chunkie/+chnk/+rcip/chunkerfunclocal.m index 8c21dee..f788463 100644 --- a/chunkie/+chnk/+rcip/chunkerfunclocal.m +++ b/chunkie/+chnk/+rcip/chunkerfunclocal.m @@ -71,10 +71,11 @@ ts = a + (b-a)*(xs+1)/2; [out{:}] = fcurve(ts); + h = (b-a)/2; chnkr.rstor(:,:,i) = reshape(out{1},dim,k); - chnkr.dstor(:,:,i) = reshape(out{2},dim,k); - chnkr.d2stor(:,:,i) = reshape(out{3},dim,k); - chnkr.hstor(i) = (b-a)/2; + chnkr.dstor(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2stor(:,:,i) = reshape(out{3},dim,k)*h*h; + end chnkr.adjstor(:,1:nch) = adjs(:,1:nch); diff --git a/chunkie/+chnk/adapgausskerneval.m b/chunkie/+chnk/adapgausskerneval.m index a1b7235..63f53e3 100644 --- a/chunkie/+chnk/adapgausskerneval.m +++ b/chunkie/+chnk/adapgausskerneval.m @@ -1,4 +1,4 @@ -function [fints,maxrecs,numints,iers] = adapgausskerneval(r,d,n,d2,h,ct,bw,j,... +function [fints,maxrecs,numints,iers] = adapgausskerneval(r,d,n,d2,ct,bw,j,... dens,rt,nt,dt,d2t,kern,opdims,t,w,opts) %CHNK.ADAPGAUSSKERNEVAL adaptive integration for interaction of kernel on chunk % at targets @@ -13,7 +13,6 @@ % r - chnkr nodes % d - chnkr derivatives at nodes % d2 - chnkr 2nd derivatives at nodes -% h - lengths of chunks in parameter space % ct - Legendre nodes at order of chunker % bw - barycentric interpolation weights for Legendre nodes at order of % chunker @@ -66,7 +65,6 @@ ds = d(:,:,j); ns = n(:,:,j); d2s = d2(:,:,j); -hs = h(j); jstart = opdims(2)*k*(j-1)+1; jend = opdims(2)*k*j; densj = reshape(dens(jstart:jend),opdims(2),k); @@ -161,8 +159,6 @@ end -fints = fints*hs; - end function val = oneintp(a,b,rs,ds,ns,d2s,densj,ct,bw, ... diff --git a/chunkie/+chnk/adapgausswts.m b/chunkie/+chnk/adapgausswts.m index 6babbd5..7afc58b 100644 --- a/chunkie/+chnk/adapgausswts.m +++ b/chunkie/+chnk/adapgausswts.m @@ -1,4 +1,4 @@ -function [mat,maxrecs,numints,iers] = adapgausswts(r,d,n,d2,h,ct,bw,j,... +function [mat,maxrecs,numints,iers] = adapgausswts(r,d,n,d2,ct,bw,j,... rt,dt,nt,d2t,kern,opdims,t,w,opts) %CHNK.ADAPGAUSSWTS adaptive integration for interaction of kernel on chunk % at targets @@ -8,7 +8,7 @@ % WARNING: this routine currently assumes that the kernel function is % translation invariant to recenter (improves stability). % -% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,h,ct,bw,j, ... +% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,ct,bw,j, ... % rt,dt,d2t,kern,opdims,t,w,opts) % % Input: @@ -16,7 +16,6 @@ % d - chnkr derivatives at nodes % n - chnkr normals at nodes % d2 - chnkr 2nd derivatives at nodes -% h - lengths of chunks in parameter space % ct - Legendre nodes at order of chunker % bw - barycentric interpolation weights for Legendre nodes at order of % chunker @@ -68,7 +67,6 @@ ds = d(:,:,j); ns = n(:,:,j); d2s = d2(:,:,j); -hs = h(j); stack = zeros(2,maxdepth); vals = zeros(opdims(1)*opdims(2)*k,maxdepth); @@ -167,8 +165,6 @@ end -mat = mat*hs; - end function val = oneintp(a,b,rs,ds,ns,d2s,ct,bw,rt,dt,nt,d2t,kern,opdims,t,w) diff --git a/chunkie/+chnk/chunk_nearparam.m b/chunkie/+chnk/chunk_nearparam.m index efdf2fd..9d670cd 100644 --- a/chunkie/+chnk/chunk_nearparam.m +++ b/chunkie/+chnk/chunk_nearparam.m @@ -24,6 +24,7 @@ maxnewt = 15; thresh0 = 1.0d-14; +iextra = 3; if nargin < 4 || isempty(t) || nargin < 5 || isempty(u) [t,~,u] = lege.exps(chnkr.k); @@ -40,12 +41,6 @@ thresh0 = opts.thresh; end -dist = Inf; - -if nargin < 3 - ich = 1:nch; -end - dim = size(pts,1); npts = numel(pts)/dim; pts = reshape(pts,dim,npts); @@ -92,9 +87,10 @@ rdiff0 = r0(:) - ref; dprime = rdiff0.'*d0(:); - dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dprime2 = (d0(:).'*d0(:) + rdiff0(:).'*d20(:)); newtsuccess = false; + ifend = 0; for iii = 1:maxnewt @@ -102,14 +98,17 @@ deltat = - dprime/dprime2; - t0 = t0+deltat; + t1 = t0 + deltat; - if (t0 > 1.0) - t0 = 1; + if (t1 > 1.0) + t1 = 1; end - if (t0 < -1.0) - t0 = -1; + if (t1 < -1.0) + t1 = -1; end + + deltat = t1-t0; + t0 = t1; all0 = (lege.exev(t0,cfs)); r0 = all0(1:dim); @@ -118,9 +117,13 @@ rdiff0 = r0(:)-ref(:); dprime = rdiff0.'*d0(:); - dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dprime2 = d0(:).'*d0(:) + rdiff0(:).'*d20(:); + + if min(abs(dprime),abs(deltat)) < thresh + ifend = ifend+1; + end - if abs(dprime) < thresh + if (ifend >= iextra) newtsuccess = true; break; end @@ -196,7 +199,7 @@ d20 = d21; rdiff0 = rdiff1; dprime = rdiff0.'*d0(:); - dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dprime2 = d0(:).'*d0(:); dist0 = dist1; lam = lam/3; if abs(dprime) < thresh diff --git a/chunkie/+chnk/chunkerkerneval_smooth.m b/chunkie/+chnk/chunkerkerneval_smooth.m index 884066f..90ea518 100644 --- a/chunkie/+chnk/chunkerkerneval_smooth.m +++ b/chunkie/+chnk/chunkerkerneval_smooth.m @@ -81,7 +81,7 @@ for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:)*chnkr.h(i); + dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); srcinfo = []; srcinfo.r = chnkr.r(:,:,i); @@ -102,7 +102,7 @@ for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:)*chnkr.h(i); + dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); srcinfo = []; srcinfo.r = chnkr.r(:,:,i); @@ -177,7 +177,7 @@ for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:)*chnkr.h(i); + dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); srcinfo = []; srcinfo.r = chnkr.r(:,:,i); @@ -215,4 +215,4 @@ end end % sum(fints) -end \ No newline at end of file +end diff --git a/chunkie/+chnk/funcuni.m b/chunkie/+chnk/funcuni.m index ac93bd3..b11dffb 100644 --- a/chunkie/+chnk/funcuni.m +++ b/chunkie/+chnk/funcuni.m @@ -104,16 +104,18 @@ for i = 1:nch a=ab(1,i); b=ab(2,i); + h = (b-a)/2; ts = a + (b-a)*(xs+1)/2; [out{:}] = fcurve(ts); chnkr.r(:,:,i) = reshape(out{1},dim,k); - chnkr.d(:,:,i) = reshape(out{2},dim,k); - chnkr.d2(:,:,i) = reshape(out{3},dim,k); - chnkr.h(i) = (b-a)/2; + chnkr.d(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2(:,:,i) = reshape(out{3},dim,k)*h*h; end chnkr.adj = adjs(:,1:nch); +chnkr.wtsstor(:,1:nch) = weights(chnkr); +chnkr.nstor(:,:,1:nch) = normals(chnkr); end diff --git a/chunkie/+chnk/pquadwts.m b/chunkie/+chnk/pquadwts.m index 01a95fe..59bc795 100644 --- a/chunkie/+chnk/pquadwts.m +++ b/chunkie/+chnk/pquadwts.m @@ -1,4 +1,4 @@ -function mat = pquadwts(r,d,n,d2,h,ct,bw,j,... +function mat = pquadwts(r,d,n,d2,ct,bw,j,... rt,dt,nt,d2t,kern,opdims,t,w,opts,intp_ab,intp) %CHNK.INTERPQUADWTS product integration for interaction of kernel on chunk % at targets @@ -14,7 +14,6 @@ % d - chnkr derivatives at nodes % n - chnkr normals at nodes % d2 - chnkr 2nd derivatives at nodes -% h - lengths of chunks in parameter space % ct - Legendre nodes at order of chunker % bw - barycentric interpolation weights for Legendre nodes at order of % chunker @@ -34,11 +33,11 @@ % % Helsing-Ojala (interior/exterior?) -h_i = h(j); + xlohi = intp_ab*(r(1,:,j)'+1i*r(2,:,j)'); % panel end points r_i = intp*(r(1,:,j)'+1i*r(2,:,j)'); % upsampled panel -d_i = h_i*(intp*(d(1,:,j)'+1i*d(2,:,j)')); % r' -d2_i = h_i^2*(intp*(d2(1,:,j)'+1i*d2(2,:,j)')); % r'' +d_i = (intp*(d(1,:,j)'+1i*d(2,:,j)')); % r' +d2_i = (intp*(d2(1,:,j)'+1i*d2(2,:,j)')); % r'' sp = abs(d_i); tang = d_i./sp; % speed, tangent n_i = -1i*tang; % normal cur = -real(conj(d2_i).*n_i)./sp.^2; % curvature diff --git a/chunkie/@chunker/arclengthdens.m b/chunkie/@chunker/arclengthdens.m index cecc8d2..470e567 100644 --- a/chunkie/@chunker/arclengthdens.m +++ b/chunkie/@chunker/arclengthdens.m @@ -1,10 +1,9 @@ function ds = arclengthdens(chnkr) %ARCLENGTHDENS arc length density on chunks % -% warning: this takes the length of the ith chunk in parameter space to -% be 2*chnkr.h(i) as opposed to 2. Thus the smooth integration rule +% The smooth integration rule % on the ith chunk is -% ds(:,i).*w*chnkr.h(i) +% ds(:,i).*w % where w are the standard Legendre weights of appropriate order and % ds is the output of this routine % diff --git a/chunkie/@chunker/arclengthfun.m b/chunkie/@chunker/arclengthfun.m index 092da6b..18dd033 100644 --- a/chunkie/@chunker/arclengthfun.m +++ b/chunkie/@chunker/arclengthfun.m @@ -13,10 +13,9 @@ istart = 1; A = lege.intmat(chnkr.k); -[~,w] = lege.exps(chnkr.k); ds = arclengthdens(chnkr); -chunklens = sum((w(:)*(chnkr.h(:).')).*ds,1); -s = (A*ds).*(chnkr.h(:).'); +chunklens = chunklen(chnkr); +s = A*ds; for i = 1:ncomp nch = nchs(i); diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m index ded6cb4..1ad1317 100644 --- a/chunkie/@chunker/chunker.m +++ b/chunkie/@chunker/chunker.m @@ -12,9 +12,6 @@ % npt - returns k*nch, the total number of points on the curve % r - dim x k x nch array, r(:,i,j) gives the coordinates of the ith % node on the jth chunk of the chunker -% h - nch array of scaling factors for chunks. the chunk derivatives are -% scaled as if the coordinates in r(:,:,j) are sampled on an -% interval of length 2*h(j). This affects integration formulas. % d - dim x k x nch array, d(:,i,j) gives the time derivative of the % coordinate at the ith node on the jth chunk of the chunker % d2 - dim x k x nch array, d(:,i,j) gives the 2nd time derivative of the @@ -79,7 +76,6 @@ d d2 adj - h n wts data @@ -89,7 +85,6 @@ dstor d2stor adjstor - hstor nstor wtsstor datastor @@ -130,9 +125,12 @@ p = chunkerpref(p); end k = p.k; + assert(k >= 2,'CHUNKER: order k of panels must be at least 2'); if nargin < 3 [obj.tstor,obj.wstor] = lege.exps(k); else + assert(length(t)==k && length(w)==k,... + 'CHUNKER: precomputed Legendre nodes appear to be wrong order'); obj.tstor = t; obj.wstor = w; end @@ -149,7 +147,6 @@ obj.wtsstor = zeros(k,nchstor); obj.d2stor = zeros(dim,k,nchstor); obj.adjstor = zeros(2,nchstor); - obj.hstor = zeros(nchstor,1); obj.vert = {}; obj.hasdata = false; obj.datastor = []; @@ -168,9 +165,6 @@ function adj = get.adj(obj) adj = obj.adjstor(:,1:obj.nch); end - function h = get.h(obj) - h = obj.hstor(1:obj.nch); - end function n = get.n(obj) n = obj.nstor(:,:,1:obj.nch); end @@ -201,9 +195,6 @@ function obj = set.adj(obj,val) obj.adjstor(:,1:obj.nch) = val; end - function obj = set.h(obj,val) - obj.hstor(1:obj.nch) = val; - end function k = get.k(obj) k = size(obj.rstor,2); end @@ -253,7 +244,6 @@ wtstemp = obj.wts; d2temp = obj.d2; adjtemp = obj.adj; - htemp = obj.h; datatemp = obj.data; obj.rstor = zeros(obj.dim,obj.k,nchstornew); obj.dstor = zeros(obj.dim,obj.k,nchstornew); @@ -261,7 +251,6 @@ obj.wtsstor = zeros(obj.k,nchstornew); obj.d2stor = zeros(obj.dim,obj.k,nchstornew); obj.adjstor = zeros(2,nchstornew); - obj.hstor = zeros(nchstornew,1); obj.datastor = zeros(obj.datadim,obj.k,nchstornew); obj.r = rtemp; obj.d = dtemp; @@ -269,7 +258,6 @@ obj.wts = wtstemp; obj.d2 = d2temp; obj.adj = adjtemp; - obj.h = htemp; obj.data = datatemp; obj.nchstor = nchstornew; end diff --git a/chunkie/@chunker/merge.m b/chunkie/@chunker/merge.m index ddc9587..052e3fb 100644 --- a/chunkie/@chunker/merge.m +++ b/chunkie/@chunker/merge.m @@ -42,7 +42,6 @@ ipos = chnkrtemp.adj > 0; chnkrtemp.adj(ipos) = chnkrtemp.adj(ipos)+nchold; chnkrout.adj(:,istart:iend) = chnkrtemp.adj; - chnkrout.h(istart:iend) = chnkrtemp.h; chnkrout.wts(:,istart:iend) = chnkrtemp.wts; chnkrout.n(:,:,istart:iend) = chnkrtemp.n; end diff --git a/chunkie/@chunker/move.m b/chunkie/@chunker/move.m index 8cbe34a..9484aef 100644 --- a/chunkie/@chunker/move.m +++ b/chunkie/@chunker/move.m @@ -19,4 +19,4 @@ obj.wts = weights(obj); -end \ No newline at end of file +end diff --git a/chunkie/@chunker/nearest.m b/chunkie/@chunker/nearest.m index d585f3e..f87c7d3 100644 --- a/chunkie/@chunker/nearest.m +++ b/chunkie/@chunker/nearest.m @@ -40,6 +40,7 @@ maxnewt = 15; thresh0 = 1.0d-9; iextra = 3; +nch = chnkr.nch; if nargin < 5 || isempty(u) [~,~,u] = lege.exps(chnkr.k); @@ -69,124 +70,15 @@ % grab chunk data ii = ich(i); r = chnkr.rstor(:,:,ii); - d = chnkr.dstor(:,:,ii); - d2 = chnkr.d2stor(:,:,ii); - h = chnkr.hstor(ii); - rc = u*(r.'); - dc = u*(d.'); - d2c = u*(d2.'); - % starting guess on over sampled grid - rover = (ainterpover*(r.')).'; - rdist0 = sqrt(sum((rover-ref(:)).^2,1)); - [disti,ind] = min(rdist0); - tt = xover(ind); - - % check the endpoints - ts = [-1;1]; - rs = (lege.exev(ts,rc)).'; - rsdist = sqrt(sum((rs-ref).^2,1)); - [dists,inds] = min(rsdist); - if dists < disti - % continue if endpoints closer - ri = rs(:,inds); - if dists < dist - dist = dists; - ds = (lege.exev(ts,dc)).'; - d2s = (lege.exev(ts,d2c)).'; - di = ds(:,inds); - d2i = d2s(:,inds); - rn = ri; dn = di; d2n = d2i; - tn = ts(inds); - ichn = ii; - end - continue; - end - - % run levenberg - - t0 = tt; - - ifend = 0; - - r0 = (lege.exev(t0,rc)).'; - - rdiff0 = ref(:) - r0(:); - - drc = lege.derpol(rc); - %d2rc = lege.derpol(drc); - thresh = thresh0; - - for iii = 1:maxnewt - %d0 = (lege.exev(t0,dc)).'*h; - d0 = (lege.exev(t0,drc)).'; % actual derivative of chunk -% d20 = (lege.exev(t0,d2rc)).'; % actual 2nd derivative of chunk - - % newton info -% -% dprime = rdiff0.'*d0(:); -% dprime2 = (d0(:).'*d0(:) + d0(:).'*d20(:)); - - % levenberg instead - - deltat = d0(:)\(rdiff0(:)); - t1 = t0+deltat; - - if (t1 > 1.0) - t1 = (1.0+t0)/2; - end - if (t1 < -1.0) - t1 = (-1.0+t0)/2; - end - - r1 = (lege.exev(t1,rc)).'; - rdiff1 = ref(:)-r1(:); - - err = min(abs(deltat),abs(sum(d0(:).*rdiff0(:)))); - if err < thresh - ifend = ifend+1; - end - - %sum(rdiff0(:).*d0(:)) - %abs(deltat) - - t0 = t1; - rdiff0 = rdiff1; - -% if (err < thresh) -% ifend = ifend + 1; -% end -% - %t0 = t1; - - if (ifend >= iextra) - break; - end - end - - % done levenberg - - ri = (lege.exev(t0,rc)).'; - disti = sqrt(sum((ri(:)-ref(:)).^2,1)); - if (ifend < iextra) - warning('newton failed in nearest'); -% deltat -% disti -% jac -% t0 -% rhs -% ifend -% sum(jac.*rhs) - else -% fprintf('success\n') - end - + [ti,ri,di,d2i,disti] = chnk.chunk_nearparam(r,ref,opts,chnkr.tstor,u); + disti = sqrt(disti); if disti < dist dist = disti; rn = ri; - dn = (lege.exev(t1,dc)).'*h; - d2n = (lege.exev(t1,d2c)).'*h*h; - tn = t1; + dn = di; + d2n = d2i; + tn = ti; ichn = ii; end end diff --git a/chunkie/@chunker/refine.m b/chunkie/@chunker/refine.m index f1ce045..6d96e6f 100644 --- a/chunkie/@chunker/refine.m +++ b/chunkie/@chunker/refine.m @@ -18,10 +18,7 @@ % (10*chnkr.nch) % opts.lvlr = level restriction flag ('a'), if 'a', enforce that no % two adjacent chunks should differ in length by more -% than a factor of approx 2 (see lvlrfac). if 't', -% enforce that the -% length in parameter space of two adjacent chunks -% differs by no more than a factor of approx 2. if 'n' +% than a factor of approx 2 (see lvlrfac). if 'n' % don't enforce (not recommended) % opts.lvlrfac = (2.0) factor for enforcing level restriction % opts.maxchunklen = maximum chunk length (Inf). enforce that no @@ -85,9 +82,6 @@ % splitting type stype = 'a'; -if strcmpi(lvlr,'t') - stype = 't'; -end % chunks told to split @@ -165,7 +159,7 @@ % level restriction -if (strcmpi(lvlr,'a') || strcmpi(lvlr,'t')) +if (strcmpi(lvlr,'a')) for ijk = 1:maxiter_lvlr nchold=chnkr.nch; @@ -175,45 +169,24 @@ i1=chnkr.adj(1,i); i2=chnkr.adj(2,i); - if (strcmpi(lvlr,'t')) - rlself = chnkr.h(i); - rl1 = rlself; - rl2 = rlself; - if (i1 > 0) - rl1 = chnkr.h(i1); - end - if (i2 > 0) - rl2 = chnkr.h(i2); - end - if (i1 < 0) - % meets at vertex (parameter space ref not recommended) - rl1 = min(chnkr.h(vert{-i1})); - end - if (i2 < 0) - rl2 = min(chnkr.h(vert{-i2})); - end - - else - - rlself = chunklens(i); + rlself = chunklens(i); - rl1=rlself; - rl2=rlself; + rl1=rlself; + rl2=rlself; - if (i1 > 0) - rl1 = chunklens(i1); - end - if (i2 > 0) - rl2 = chunklens(i2); - end - if (numel(vert) ~= 0) - if (i1 < 0) - rl1 = min(chunklens(vert{-i1})); - end - if (i2 < 0) - rl2 = min(chunklens(vert{-i2})); - end - end + if (i1 > 0) + rl1 = chunklens(i1); + end + if (i2 > 0) + rl2 = chunklens(i2); + end + if (numel(vert) ~= 0) + if (i1 < 0) + rl1 = min(chunklens(vert{-i1})); + end + if (i2 < 0) + rl2 = min(chunklens(vert{-i2})); + end end % only check if self is larger than either of adjacent blocks, @@ -265,7 +238,7 @@ % split chunk i now, and recalculate nodes, d, etc if (chnkr.nch + 1 > nchmax) - error('too many chunks') + error('CHUNKER.REFINE nchmax=%d exceeded during oversample',nchmax) end chnkr = split(chnkr,i,[],x,w,u,stype); diff --git a/chunkie/@chunker/sort.m b/chunkie/@chunker/sort.m index 6bad4c4..8e9887f 100644 --- a/chunkie/@chunker/sort.m +++ b/chunkie/@chunker/sort.m @@ -32,7 +32,6 @@ chnkr.r = chnkr.r(:,:,inds); chnkr.d = chnkr.d(:,:,inds); chnkr.d2 = chnkr.d2(:,:,inds); -chnkr.h = chnkr.h(inds); chnkr.n = chnkr.n(:,:,inds); chnkr.wts = chnkr.wts(:,inds); diff --git a/chunkie/@chunker/split.m b/chunkie/@chunker/split.m index 3f36283..995d309 100644 --- a/chunkie/@chunker/split.m +++ b/chunkie/@chunker/split.m @@ -66,7 +66,7 @@ if strcmpi(stype1,'a') % first construct dsdt -dsdt = sqrt(sum(d.^2,1))*chnkr.hstor(ich); +dsdt = sqrt(sum(d.^2,1)); dsdt = dsdt(:); rltot = dot(dsdt,w); @@ -124,8 +124,6 @@ d2_2 = lege.exev(ts2,cd2); d2_2 = d2_2.'; -hold=chnkr.hstor(ich); - % update chnkr %i1=chnkr.adjstor(1,ich); @@ -133,19 +131,17 @@ chnkr = chnkr.addchunk(); -chnkr.hstor(ich) = hold*(t1+1)/2; -chnkr.hstor(nch+1) = hold*(1-t1)/2; +h1 = (t1+1)/2; +h2 = (1-t1)/2; chnkr.rstor(:,:,ich) = r_1; chnkr.rstor(:,:,nch+1) = r_2; -chnkr.dstor(:,:,ich) = d_1; -chnkr.wtsstor(:,ich) = (sqrt(sum(d_1.^2,1)).') .* chnkr.wstor * ... - chnkr.hstor(ich); -chnkr.dstor(:,:,nch+1) = d_2; -chnkr.wtsstor(:,nch+1) = (sqrt(sum(d_2.^2,1)).') .* chnkr.wstor * ... - chnkr.hstor(nch+1); -chnkr.d2stor(:,:,ich) = d2_1; -chnkr.d2stor(:,:,nch+1) = d2_2; +chnkr.dstor(:,:,ich) = d_1*h1; +chnkr.wtsstor(:,ich) = (sqrt(sum(d_1.^2,1)).') .* chnkr.wstor*h1; +chnkr.dstor(:,:,nch+1) = d_2*h2; +chnkr.wtsstor(:,nch+1) = (sqrt(sum(d_2.^2,1)).') .* chnkr.wstor*h2; +chnkr.d2stor(:,:,ich) = d2_1*h1*h1; +chnkr.d2stor(:,:,nch+1) = d2_2*h2*h2; chnkr.adjstor(2,ich)=nch+1; chnkr.adjstor(1,nch+1)=ich; diff --git a/chunkie/@chunker/upsample.m b/chunkie/@chunker/upsample.m index 5aeb177..25d6f77 100644 --- a/chunkie/@chunker/upsample.m +++ b/chunkie/@chunker/upsample.m @@ -36,7 +36,6 @@ chnkrup.d2= permute(reshape(upmat*reshape( ... permute(chnkr.d2,[2,1,3]),k,nn),kup,dim,nch),[2,1,3]); chnkrup.adj= chnkr.adj; -chnkrup.h= chnkr.h; chnkrup.n = normals(chnkrup); chnkrup.wts = weights(chnkrup); diff --git a/chunkie/@chunker/weights.m b/chunkie/@chunker/weights.m index af6afe9..fcce774 100644 --- a/chunkie/@chunker/weights.m +++ b/chunkie/@chunker/weights.m @@ -24,6 +24,6 @@ nch = chnkr.nch; w = chnkr.wstor; wts = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch); - wts = wts.*bsxfun(@times,w(:),(chnkr.h(:)).'); + wts = wts.*w(:); end diff --git a/chunkie/@chunker/whts.m b/chunkie/@chunker/whts.m index 1da1df0..43bd30c 100644 --- a/chunkie/@chunker/whts.m +++ b/chunkie/@chunker/whts.m @@ -1,12 +1,8 @@ -function wchnk = whts(chnkr) +function wts = whts(chnkr) warning('whts is deprecated and will be removed, use weights instead'); - k = chnkr.k; - nch = chnkr.nch; - [~,w] = lege.exps(k); - wchnk = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch); - wchnk = wchnk.*bsxfun(@times,w(:),(chnkr.h(:)).'); + wts = weights(chnkr); end diff --git a/chunkie/@chunkgraph/balance.m b/chunkie/@chunkgraph/balance.m index dfa9bca..111e7f2 100644 --- a/chunkie/@chunkgraph/balance.m +++ b/chunkie/@chunkgraph/balance.m @@ -39,7 +39,6 @@ parcl = zeros([numel(vedge),1]); pinds = zeros([numel(vedge),1]); - [~,wleg16,~,~] = lege.exps(16); for ii=1:numel(vedge) if (sign(vsign(ii)) == -1) ds = obj.echnks(vedge(ii)).d(:,:,1); @@ -49,15 +48,11 @@ pinds(ii) = size(obj.echnks(vedge(ii)).r,3); end - h = obj.echnks(vedge(ii)).h(pinds(ii)); k = obj.echnks(vedge(ii)).k; - if (k ~=16) - [~,wleg,~,~] = lege.exps(k); - else - wleg = wleg16; - end - arc = sum(sqrt(ds(1,:).^2+ds(2,:).^2).*wleg'*h); + wleg = echnks(vedge(ii)).wstor; + + arc = sum(sqrt(ds(1,:).^2+ds(2,:).^2).*wleg'); parcl(ii) = arc; end 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/elast2d.m b/chunkie/@kernel/elast2d.m index ea6991b..3940faf 100644 --- a/chunkie/@kernel/elast2d.m +++ b/chunkie/@kernel/elast2d.m @@ -90,4 +90,9 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + end diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m index 73cc8ce..057b5d3 100644 --- a/chunkie/@kernel/helm2d.m +++ b/chunkie/@kernel/helm2d.m @@ -86,4 +86,10 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + + end diff --git a/chunkie/@kernel/helm2ddiff.m b/chunkie/@kernel/helm2ddiff.m index a1392b9..9dae103 100644 --- a/chunkie/@kernel/helm2ddiff.m +++ b/chunkie/@kernel/helm2ddiff.m @@ -92,4 +92,10 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + + end diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index 30803a6..08a72e9 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,23 +83,32 @@ 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 elseif ( isa(kern, 'function_handle') ) obj.eval = kern; elseif ( isa(kern, 'kernel') ) - if ( numel(kern) == 1 ) - obj = kern; - else - % The input is a matrix of kernels. - % Create a single kernel object by interleaving the - % outputs of each sub-kernel's eval() and fmm() routines. - % TODO: Check that opdims are consistent - obj = interleave(kern); + if (numel(kern) == 1) + obj = kern; + else + obj = interleave(kern); end + elseif isa(kern,'cell') + sz = size(kern); assert(length(sz)==2,'KERNEL: first input not of a supported type'); + obj(sz(1),sz(2)) = kernel(); + for j = 1:sz(2) + for i = 1:sz(1) + obj(i,j) = kernel(kern{i,j}); + end + end else - error('First input must be a string or function handle.'); + error('KERNEL: first input not of a supported type'); end end @@ -109,6 +122,8 @@ obj = helm2ddiff(varargin); obj = stok2d(varargin); obj = elast2d(varargin); + obj = axissymhelm2d(varargin); + obj = axissymhelm2ddiff(varargin); obj = zeros(varargin); end @@ -154,6 +169,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 +269,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/lap2d.m b/chunkie/@kernel/lap2d.m index 8db13bb..1eb548a 100644 --- a/chunkie/@kernel/lap2d.m +++ b/chunkie/@kernel/lap2d.m @@ -85,4 +85,9 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + end diff --git a/chunkie/@kernel/stok2d.m b/chunkie/@kernel/stok2d.m index c8f2327..e017192 100644 --- a/chunkie/@kernel/stok2d.m +++ b/chunkie/@kernel/stok2d.m @@ -84,4 +84,10 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + + end 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/chunkerfit.m b/chunkie/chunkerfit.m new file mode 100644 index 0000000..c3b060d --- /dev/null +++ b/chunkie/chunkerfit.m @@ -0,0 +1,166 @@ +function chnkr = chunkerfit(xy, opts) +%CHUNKERFIT Create a chunker by fitting a curve to a set of points. +% +% Syntax: chnkr = CHUNKERFIT(xy, opts) +% +% Input: +% xy - (2,:) array of points used to fit a curve +% +% Optional input: +% opts - options structure (defaults) +% opts.method = string ('spline') +% Method to use for curve fitting. +% opts.ifclosed = boolean (true) +% Flag for fitting a closed (true) or open (false) curve. +% opts.splitatpoints = boolean (false) +% If true, start the adaptive splitting process from the set of +% splits implied by the given points. This ensures that the given +% points will always coincide with panel endpoints in the resulting +% chunker object. If false, the adaptive splitting process will +% start from scratch. +% opts.cparams = struct +% Curve parameters structure to be passed to chunkerfunc after +% fitting is done. See chunkerfunc for parameters. +% opts.pref = chunkerpref object or struct +% Chunker preferences to be passed to chunkerfunc. See chunkerpref +% for parameters. +% +% Output: +% chnkr - chunker object discretizing the fitted curve +% +% Examples: +% r = chnk.curves.bymode(sort(2*pi*rand(20,1)), [2 0.5 0.2 0.7]); +% chnkr = chunkerfit(r); % Closed curve +% opts = []; opts.ifclosed = false; +% chnkr = chunkerfit(r(:,1:10), opts); % Open curve +% +% See also CHUNKERFUNC, CHUNKERPREF, CHUNKER. + +if ( size(xy, 1) ~= 2 ) + error('CHUNKIE:CHUNKERFIT:size', 'Points must be specified as a 2xN matrix.'); +end + +x = xy(1,:).'; +y = xy(2,:).'; + +% Set default options +defaults = []; +defaults.method = 'spline'; +defaults.ifclosed = true; +defaults.splitatpoints = false; +defaults.cparams = []; +defaults.pref = []; + +if ( nargin < 2 ) + opts = defaults; +else + for field = fieldnames(defaults).' + if ( ~isfield(opts, field{1}) ) + opts.(field{1}) = defaults.(field{1}); + end + end +end + +if ( ~strcmpi(opts.method, 'spline') ) + error('CHUNKIE:CHUNKERFIT:method', 'Unsupported method ''%s''.', opts.method); +end + +if ( opts.ifclosed ) + % Wrap the nodes around + x = [x; x(1)]; + y = [y; y(1)]; +end + +% Choose the parametric grid to be polygonal arc length +t = [0; cumsum(sqrt(diff(x).^2 + diff(y).^2))]; + +ppx = myspline(t, x, opts.ifclosed); +ppy = myspline(t, y, opts.ifclosed); +ppx_dx = ppdiff(ppx, 1); +ppy_dy = ppdiff(ppy, 1); +ppx_dxx = ppdiff(ppx, 2); +ppy_dyy = ppdiff(ppy, 2); + +% Pack polynomial coefficients into a tensor to vectorize ppval +cfs = zeros([6 size(ppx.coefs)]); +cfs(1,:,:) = ppx.coefs; +cfs(2,:,:) = ppy.coefs; +cfs(3,:,:) = ppx_dx.coefs; +cfs(4,:,:) = ppy_dy.coefs; +cfs(5,:,:) = ppx_dxx.coefs; +cfs(6,:,:) = ppy_dyy.coefs; +breaks = ppx.breaks.'; + + function [r, d, d2] = splinefunc(tt) + vals = myppval(breaks, cfs, tt); + xs = vals(1,:); + ys = vals(2,:); + dxs = vals(3,:); + dys = vals(4,:); + d2xs = vals(5,:); + d2ys = vals(6,:); + r = [xs(:).' ; ys(:).' ]; + d = [dxs(:).' ; dys(:).' ]; + d2 = [d2xs(:).'; d2ys(:).']; + end + +cparams = opts.cparams; +cparams.ifclosed = opts.ifclosed; +cparams.ta = t(1); +cparams.tb = t(end); +if ( opts.splitatpoints ) + cparams.tsplits = t; +end +chnkr = chunkerfunc(@splinefunc, cparams, opts.pref); + +end + +function y = myppval(breaks, cfs, x) + [~, idx] = histc(x, [-inf; breaks(2:end-1); inf]); + p = reshape((x-breaks(idx)).^(3:-1:0), [1 size(idx,1) size(cfs,3)]); + y = sum(p .* cfs(:,idx,:), 3); +end + +function pp = ppdiff(pp, n) + deg = pp.order-1; + D = diag(deg:-1:1, 1); + for k = 1:n + pp.coefs = pp.coefs * D^n; + end +end + +function pp = myspline(x, y, ifclosed) + n = length(y); + dx = diff(x); + dy = diff(y); + dydx = dy ./ dx; + b = zeros(n, 1); + A = spdiags([ [dx(2:n-1) ; 0 ; 0], ... + 2*[0 ; dx(2:n-1)+dx(1:n-2) ; 0], ... + [0 ; 0 ; dx(1:n-2)] ], ... + [-1 0 1], n, n); + + if ( ifclosed ) + % Periodic conditions: match first and second derivatives at the first + % and last points + A(1,[1 n]) = [1 -1]; + A(n,1:2) = dx(n-1)*[2 1]; + A(n,n-1:n) = A(n,n-1:n) + dx(1)*[1 2]; + b(2:n-1,:) = 3*(dx(2:n-1).*dydx(1:n-2) + dx(1:n-2).*dydx(2:n-1)); + b(n,:) = 3*(dx(n-1)*dydx(1) + dx(1)*dydx(n-1)); + else + % Not-a-knot conditions: match third derivatives at the first and last + % interior breaks + A(1,1:2) = [dx(2), x(3)-x(1)]; + A(n,n-1:n) = [x(n)-x(n-2), dx(n-2)]; + b(1,:) = ((dx(1)+2*(x(3)-x(1)))*dx(2)*dydx(1) + dx(1)^2*dydx(2)) / (x(3)-x(1)); + b(2:n-1,:) = 3*(dx(2:n-1).*dydx(1:n-2) + dx(1:n-2).*dydx(2:n-1)); + b(n,:) = (dx(n-1)^2*dydx(n-2) + ((2*(x(n)-x(n-2))+dx(n-1))*dx(n-2))*dydx(n-1)) / (x(n)-x(n-2)); + end + + % Solve linear system for the coefficients + c = A \ b; + c4 = (c(1:n-1)+c(2:n)-2*dydx(1:n-1)) ./ dx; + c3 = (dydx(1:n-1)-c(1:n-1))./dx - c4; + pp = mkpp(x, [c4./dx c3 c(1:n-1) y(1:n-1)]); +end diff --git a/chunkie/chunkerfunc.m b/chunkie/chunkerfunc.m index 4367785..1925403 100644 --- a/chunkie/chunkerfunc.m +++ b/chunkie/chunkerfunc.m @@ -67,11 +67,12 @@ pref = chunkerpref(pref); end +chnkr = chunker(pref); % empty chunker ta = 0.0; tb = 2*pi; ifclosed=true; chsmall = Inf; nover = 0; eps = 1.0e-6; -lvlr = 'a'; maxchunklen = Inf; lvlrfac = 2.0; +lvlr = 'a'; maxchunklen = Inf; lvlrfac = 2; nout = 1; if isfield(cparams,'ta') @@ -152,11 +153,11 @@ tsplits = sort(unique(tsplits),'ascend'); lab = length(tsplits); if (lab-1 > nchmax) - error(['nchmax exceeded in chunkerfunc on initial splits.\n ',... + error(['CHUNKERFUNC: nchmax exceeded in chunkerfunc on initial splits.\n ',... 'try increasing nchmax in preference struct']); end if (any(tsplits > tb) || any(tsplits < ta)) - error(['tsplits outside of interval of definition.\n', ... + error(['CHUNKERFUNC: tsplits outside of interval of definition.\n', ... 'check definition of splits, ta and tb']); end @@ -258,7 +259,7 @@ % . . . if here, not resolved % divide - first update the adjacency list if (nch +1 > nchmax) - error('too many chunks') + error('CHUNKERFUNC: nchmax=%d exceeded. Unable to resolve curve.',nchmax) end ifprocess(ich)=0; @@ -366,7 +367,7 @@ % split chunk i now, and recalculate nodes, ders, etc if (nch + 1 > nchmax) - error('too many chunks') + error('CHUNKERFUNC: nchmax=%d exceeded during level restriction (curve resolved).',nchmax); end @@ -436,7 +437,7 @@ fbkm1 = fbk; bkm1 = bk; - for iter = 1:200 + for iter = 1:52 m = (ak+bk)/2; s = m; if fbkm1 ~= fbk @@ -486,7 +487,7 @@ end if (nch + 1 > nchmax) - error('too many chunks') + error('CHUNKERFUNC: nchmax=%d exceeded while oversampling by nover=%d',nchmax,nover); end adjs(1,nch+1)=i; @@ -507,7 +508,6 @@ % . . . finally evaluate the k nodes on each chunk, along with % derivatives and chunk lengths -chnkr = chunker(pref); % empty chunker chnkr = chnkr.addchunk(nch); @@ -520,10 +520,10 @@ for j = nout+1:3 out{j} = out{j-1}*dermat*(2/(b-a)); end + h = (b-a)/2; chnkr.rstor(:,:,i) = reshape(out{1},dim,k); - chnkr.dstor(:,:,i) = reshape(out{2},dim,k); - chnkr.d2stor(:,:,i) = reshape(out{3},dim,k); - chnkr.hstor(i) = (b-a)/2; + chnkr.dstor(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2stor(:,:,i) = reshape(out{3},dim,k)*h*h; end chnkr.adjstor(:,1:nch) = adjs(:,1:nch); diff --git a/chunkie/chunkerfuncuni.m b/chunkie/chunkerfuncuni.m index b58765d..edfd459 100644 --- a/chunkie/chunkerfuncuni.m +++ b/chunkie/chunkerfuncuni.m @@ -123,10 +123,10 @@ for j = nout+1:3 out{j} = out{j-1}*dermat*(2/(b-a)); end + h = (b-a)/2; chnkr.rstor(:,:,i) = reshape(out{1},dim,k); - chnkr.dstor(:,:,i) = reshape(out{2},dim,k); - chnkr.d2stor(:,:,i) = reshape(out{3},dim,k); - chnkr.hstor(i) = (b-a)/2; + chnkr.dstor(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2stor(:,:,i) = reshape(out{3},dim,k)*h*h; end chnkr.adjstor(:,1:nch) = adjs(:,1:nch); diff --git a/chunkie/chunkerintegral.m b/chunkie/chunkerintegral.m index 1a239fe..ed723ef 100644 --- a/chunkie/chunkerintegral.m +++ b/chunkie/chunkerintegral.m @@ -66,7 +66,7 @@ for i = 1:nch rci = rc(:,:,i); dci = dc(:,:,i); - fint = fint + chnk.intchunk.fhandle(f,rci,dci)*chnkr.h(i); + fint = fint + chnk.intchunk.fhandle(f,rci,dci); end else fint = 0.0; @@ -74,7 +74,7 @@ for i = 1:nch dci = dc(:,:,i); fci = fc(:,i); - fint = fint + chnk.intchunk.fcoefs(fci,dci)*chnkr.h(i); + fint = fint + chnk.intchunk.fcoefs(fci,dci); end end diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m index cde9605..f0c29fe 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,34 @@ 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.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 +141,8 @@ [chnkr2] = upsample(chnkr,nlegnew); + + icont = false; if usefmm_final try diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index 0246a75..4906252 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -179,7 +179,6 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; % interpolation matrix intp = lege.matrin(k,t); % interpolation from k to 2*k @@ -207,7 +206,7 @@ % Helsing-Ojala (interior/exterior?) - mat1 = chnk.pquadwts(r,d,n,d2,h,ct,bw,j,targs(:,ji), ... + mat1 = chnk.pquadwts(r,d,n,d2,ct,bw,j,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts,intp_ab,intp); % depends on kern, different mat1? fints(ji) = fints(ji) + mat1*dens(idxjmat); @@ -268,14 +267,12 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; - if isempty(flag) % do all to all adaptive for i = 1:nch for j = 1:nt - fints1 = chnk.adapgausskerneval(r,d,n,d2,h,ct,bw,i,dens,targs(:,j), ... + fints1 = chnk.adapgausskerneval(r,d,n,d2,ct,bw,i,dens,targs(:,j), ... targd(:,j),targn(:,j),targd2(:,j),kerneval,opdims,t,w,opts); indj = (j-1)*opdims(1); @@ -286,7 +283,7 @@ else % do only those flagged for i = 1:nch [ji] = find(flag(:,i)); - [fints1,maxrec,numint,iers] = chnk.adapgausskerneval(r,d,n,d2,h,ct,bw,i,dens,targs(:,ji), ... + [fints1,maxrec,numint,iers] = chnk.adapgausskerneval(r,d,n,d2,ct,bw,i,dens,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kerneval,opdims,t,w,opts); indji = (ji-1)*opdims(1); diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m index a697ac2..7a5f593 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/chunkie/chunkerkernevalmat.m @@ -1,4 +1,4 @@ -function mat = chunkerkernevalmat(chnkr,kern,targobj,opts) + function mat = chunkerkernevalmat(chnkr,kern,targobj,opts) %CHUNKERKERNEVALMAT compute the matrix which maps density values on % the chunk geometry to the value of the convolution of the given % integral kernel with the density at the specified target points @@ -273,28 +273,13 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; - h = chnkr.h; for i = 1:nch jmat = 1 + (i-1)*k*opdims(2); jmatend = i*k*opdims(2); - mat(:,jmat:jmatend) = chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs, ... - targd,targd2,kern,opdims,t,w,opts); - - js1 = jmat:jmatend; - js1 = repmat( (js1(:)).',1,opdims(1)*numel(ji)); - - indji = (ji-1)*opdims(1); - indji = repmat( (indji(:)).', opdims(1),1) + ( (1:opdims(1)).'); - indji = indji(:); - indji = repmat(indji,1,opdims(2)*k); - - iend = istart+numel(mat1)-1; - is(istart:iend) = indji(:); - js(istart:iend) = js1(:); - vs(istart:iend) = mat1(:); - istart = iend+1; + mat(:,jmat:jmatend) = chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs, ... + targd,targn,targd2,kern,opdims,t,w,opts); end else @@ -310,13 +295,12 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; - h = chnkr.h; for i = 1:nch jmat = 1 + (i-1)*k*opdims(2); jmatend = i*k*opdims(2); [ji] = find(flag(:,i)); - mat1 = chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs(:,ji), ... + mat1 = chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts); js1 = jmat:jmatend; @@ -388,7 +372,6 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; - h = chnkr.h; % interpolation matrix intp = lege.matrin(k,t); % interpolation from k to 2*k @@ -401,7 +384,7 @@ [ji] = find(flag(:,i)); % Helsing-Ojala (interior/exterior?) - mat1 = chnk.pquadwts(r,d,n,d2,h,ct,bw,i,targs(:,ji), ... + mat1 = chnk.pquadwts(r,d,n,d2,ct,bw,i,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts,intp_ab,intp); % depends on kern, different mat1? js1 = jmat:jmatend; diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index d190ea5..79f2fa5 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -95,33 +95,30 @@ % opts.sing provides a default value for singularities if not % defined for kernels -if isa(kern,'function_handle') - kern2 = kernel(kern); - kern = kern2; -elseif isa(kern,'cell') - sz = size(kern); - kern2(sz(1),sz(2)) = kernel(); - for j = 1:sz(2) - for i = 1:sz(1) - if isa(kern{i,j},'function_handle') - kern2(i,j) = kernel(kern{i,j}); - elseif isa(kern{i,j},'kernel') - kern2(i,j) = kern{i,j}; - else - msg = "Second input is not a kernel object, function handle, " ... - + "or cell array"; - error(msg); - end - end + +% Flag for determining whether input object is a chunkergraph +icgrph = 0; + +if (class(chnkobj) == "chunker") + chnkrs = chnkobj; + npttot = chnkobj.npt; +elseif(class(chnkobj) == "chunkgraph") + icgrph = 1; + chnkrs = chnkobj.echnks; + npttot = chnkobj.npt; +else + msg = "CHUNKERMAT: first input is not a chunker or chunkgraph object"; + error(msg) +end + +if ~isa(kern,'kernel') + try + kern = kernel(kern); + catch + error('CHUNKERMAT: second input kern not of supported type'); end - kern = kern2; - -elseif ~isa(kern,'kernel') - msg = "Second input is not a kernel object, function handle, " ... - + "or cell array"; - error(msg); end - + if nargin < 3 opts = []; end @@ -177,21 +174,6 @@ adaptive_correction = opts.adaptive_correction; end -% Flag for determining whether input object is a chunkergraph -icgrph = 0; - -if (class(chnkobj) == "chunker") - chnkrs = chnkobj; - npttot = chnkobj.npt; -elseif(class(chnkobj) == "chunkgraph") - icgrph = 1; - chnkrs = chnkobj.echnks; - npttot = chnkobj.npt; -else - msg = "First input is not a chunker or chunkgraph object"; - error(msg) -end - nchunkers = length(chnkrs); opdims_mat = zeros(2,nchunkers,nchunkers); @@ -728,14 +710,13 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; for i = 1:nch jmat = 1 + (i-1)*k*opdims(2); jmatend = i*k*opdims(2); [ji] = find(flag(:,i)); - mat1 = chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs(:,ji), ... + mat1 = chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kernev,opdims,t,w,opts); js1 = jmat:jmatend; diff --git a/chunkie/bieapply.m b/chunkie/chunkermatapply.m similarity index 97% rename from chunkie/bieapply.m rename to chunkie/chunkermatapply.m index 58b3cdf..8b91580 100644 --- a/chunkie/bieapply.m +++ b/chunkie/chunkermatapply.m @@ -1,7 +1,7 @@ -function u = bieapply(chnkr,kern,dval,dens,cormat,opts) -%BIEAPPLY - apply BIE system on chunker defined by kern +function u = chunkermatapply(chnkr,kern,dval,dens,cormat,opts) +%CHUNKERMATAPPLY - apply chunkermat system on chunker defined by kern % -% Syntax: u = bieapply(chnkr,kern,dval,dens,sigma,opts) +% Syntax: u = chunkermatapply(chnkr,kern,dval,dens,sigma,opts) % % Input: % chnkobj - chunker object describing boundary @@ -240,4 +240,4 @@ end end -end \ No newline at end of file +end diff --git a/chunkie/chunkerpoly.m b/chunkie/chunkerpoly.m index 69f8199..016d760 100644 --- a/chunkie/chunkerpoly.m +++ b/chunkie/chunkerpoly.m @@ -194,15 +194,15 @@ val1 = edgevals(:,i); chnkr.data(:,:,nch) = bsxfun(@times,val1,onek(:).'); end + h = (l-w2-w1)/2.0; chnkr.r(:,:,nch) = r1 + bsxfun(@times,v(:),(ts(:)).'); - chnkr.d(:,:,nch) = repmat(v(:),1,k); - chnkr.d2(:,:,nch) = zeros(dim,k); + chnkr.d(:,:,nch) = repmat(v(:),1,k)*h; + chnkr.d2(:,:,nch) = zeros(dim,k)*h*h; chnkr.adj(1,nch) = 0; chnkr.adj(2,nch) = 0; if nch > 1 chnkr.adj(1,nch) = nch-1; chnkr.adj(2,nch-1) = nch; end - chnkr.h(nch) = (l-w2-w1)/2.0; if or(i < nv-1,ifclosed) if rounded @@ -245,11 +245,9 @@ chnkr.adj(:,nch+1:nch+ncht) = chnkrt.adj+nch; chnkr.adj(2,nch) = nch+1; chnkr.adj(1,nch+1) = nch; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; if nvals > 0 - ds = bsxfun(@times,reshape(sum((rotmat*chnkrt.d(:,:)).^2,1),k,ncht), ... - (chnkrt.h(:)).'); + ds = reshape(sum((rotmat*chnkrt.d(:,:)).^2,1),k,ncht); dsw = bsxfun(@times,w(:),ds); dssums = sum(dsw,1); dssums2 = cumsum([0,dssums(1:end-1)]); @@ -286,10 +284,11 @@ ncht = depth + 1; chnkrt = chunker(pref); chnkrt = chnkrt.addchunk(ncht); + hall = hadap*w2/2; hall = repmat(hall(:).',k,1); + hall=hall(:).'; chnkrt.r = bsxfun(@times,v,tadap)*w2 + r2; - chnkrt.d = repmat(v,1,k,ncht); - chnkrt.d2 = zeros(dim,k,ncht); - chnkrt.h = hadap*w2/2; + chnkrt.d(:,:) = repmat(v,1,k*ncht).*hall; + chnkrt.d2(:,:) = zeros(dim,k*ncht).*(hall.^2); chnkrt.adj = [0, 1:(ncht-1); 2:ncht, 0]; chnkrt = reverse(chnkrt); chnkrt = sort(chnkrt); @@ -303,7 +302,6 @@ chnkr.adj(2,nch) = nch+1; chnkr.adj(1,nch+1) = nch; chnkr.adj(2,nch+ncht) = 0; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; if nvals > 0 chnkr.data(:,:,nch+1:nch+ncht) = repmat(val1,1,k,ncht); @@ -317,9 +315,9 @@ chnkrt = chunker(pref); chnkrt = chnkrt.addchunk(ncht); chnkrt.r = bsxfun(@times,v2,tadap)*w2 + r2; - chnkrt.d = repmat(v2,1,k,ncht); - chnkrt.d2 = zeros(dim,k,ncht); - chnkrt.h = hadap*w2/2; + chnkrt.d(:,:) = repmat(v2,1,k*ncht).*hall; + chnkrt.d2(:,:) = zeros(dim,k*ncht).*(hall.^2); + chnkrt.adj = [0, 1:(ncht-1); 2:ncht, 0]; chnkrt = sort(chnkrt); @@ -332,7 +330,6 @@ chnkr.adj(:,nch+1:nch+ncht) = chnkrt.adj+nch; chnkr.adj(2,nch) = nch+1; chnkr.adj(1,nch+1) = nch; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; chnkr = chnkr.addvert([nch,nch+1]); 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/chunkie/demo/demo_barbell.m b/chunkie/demo/demo_barbell.m index 96cc7ce..65965c8 100644 --- a/chunkie/demo/demo_barbell.m +++ b/chunkie/demo/demo_barbell.m @@ -21,12 +21,13 @@ % parameters for curve rounding/chunking routine to oversample boundary cparams = []; cparams.widths = 0.1*ones(size(verts,2),1);% width to cut about each corner -cparams.eps = 1e-6; % tolerance at which to resolve curve +cparams.eps = 1e-12; % tolerance at which to resolve curve cparams.rounded = true; % call smoothed-polygon chunking routine % a smoothed version of edgevals is returned in % chnkr.data + chnkr = chunkerpoly(verts,cparams,[],edgevals); chnkr = chnkr.refine(); diff --git a/chunkie/roundedpolygon.m b/chunkie/roundedpolygon.m index 867e3ca..5ba0c08 100644 --- a/chunkie/roundedpolygon.m +++ b/chunkie/roundedpolygon.m @@ -53,210 +53,4 @@ % % See also CHUNKERFUNC, CHUNKER, CHUNKERPREF, REFINE -[dimv,nv] = size(verts); - -assert(dimv > 1) - -if nargin < 2 - cparams = []; -end -if nargin < 3 - p = []; - p.dim = dimv; - pref = chunkerpref(p); -else - pref = chunkerpref(pref); - if pref.dim ~= dimv - warning('dimensions dont match, overwriting with vertex dim'); - pref.dim = dimv; - end -end -if nargin < 4 - edgevals = []; -end - -autowidths = false; -autowidthsfac = 0.1; -ifclosed = true; -eps = 1e-6; -nover = 0; - -if isfield(cparams,'ifclosed') - ifclosed = cparams.ifclosed; -end -if isfield(cparams,'eps') - eps = cparams.eps; -end -if isfield(cparams,'nover') - nover = cparams.nover; -end - -if (ifclosed) - verts2 = [verts(:,end), verts, verts(:,1)]; - edges2 = sqrt(sum(diff(verts2,1,2).^2,1)); -else - edges2 = [0, sqrt(sum(diff(verts,1,2).^2,1)), 0]; -end - -widths_not_set = true; - -if isfield(cparams,'widths') - widths = cparams.widths; - widths_not_set = false; -end -if isfield(cparams,'autowidths') - autowidths = cparams.autowidths; -end -if isfield(cparams,'autowidthsfac') - autowidthsfac = cparams.autowidthsfac; -end - -if (autowidths || widths_not_set) - widths = autowidthsfac*... - min(edges2(1:end-1),edges2(2:end)); -end - -if ifclosed - widths = [widths(:); widths(1)]; - verts = [verts, verts(:,1)]; - nv = size(verts,2); -end - -nedge = nv - 1; -nvals = 0; -if numel(edgevals) ~= 0 - assert(rem(numel(edgevals),nedge) == 0, ... - 'number of edge values should be multiple of number of edges'); - nvals = numel(edgevals)/nedge; - edgevals = reshape(edgevals,nvals,nedge); -end - -chnkr = chunker(pref); -chnkr = chnkr.makedatarows(nvals); -k = chnkr.k; dim = chnkr.dim; -[t,w] = lege.exps(k); -onek = ones(k,1); - - -if nvals > 0 - aint = lege.intmat(chnkr.k); -end - -for i = 1:nv-1 - % grab vertices - r1 = verts(:,i); r2 = verts(:,i+1); - w1 = widths(i); w2 = widths(i+1); - l = sqrt(sum((r1-r2).^2)); - assert(l > w1+w2+2*eps(1)*l,'widths too large for side'); - - % make chunk in middle - v = (r2-r1)/l; - ts = w1+(l-w2-w1)*(t+1)/2.0; - chnkr = chnkr.addchunk(); - nch = chnkr.nch; - if nvals > 0 - val1 = edgevals(:,i); - chnkr.data(:,:,nch) = bsxfun(@times,val1,onek(:).'); - end - chnkr.r(:,:,nch) = r1 + bsxfun(@times,v(:),(ts(:)).'); - chnkr.d(:,:,nch) = repmat(v(:),1,k); - chnkr.d2(:,:,nch) = zeros(dim,k); - chnkr.adj(1,nch) = 0; chnkr.adj(2,nch) = 0; - if nch > 1 - chnkr.adj(1,nch) = nch-1; - chnkr.adj(2,nch-1) = nch; - end - chnkr.h(nch) = (l-w2-w1)/2.0; - - if or(i < nv-1,ifclosed) - % chunk up smoothed corner made by three verts - if (i==nv-1) - if nvals > 0 - val2 = edgevals(:,1); - end - r3 = verts(:,2); - else - if nvals > 0 - val2 = edgevals(:,i+1); - end - r3 = verts(:,i+2); - end - l2 = sqrt(sum((r2-r3).^2)); - v = -v; - v2 = (r3-r2)/l2; - cosphi = dot(v2,v); %cosine of angle between edges - sinphi = sqrt(1-cosphi*cosphi); - trange = w2*sqrt((1-cosphi)/2); %parameter space range of gaussian - hbell = abs(trange)/8.0; % width parameter of gaussian - m = sinphi/(1-cosphi); % slope of abs approx - cpt.ta = -trange; - cpt.tb = trange; - cpt.eps = eps; cpt.levrestr = 0; cpt.ifclosed = 0; - chnkrt = sort(chunkerfunc(@(t)fround(t,m,hbell,dim),cpt,pref)); - - % do optimal procrustes match of left and right ends - rl = fround(-trange,m,hbell,dim); rr = fround(trange,m,hbell,dim); - [um,~,vm] = svd([w2*v w2*v2]*([rl rr].')); - rotmat = um*vm.'; - - % copy in rotated and translated chunks - ncht = chnkrt.nch; - chnkr = chnkr.addchunk(ncht); - chnkr.r(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.r(:,:),dim,k,ncht)+r2; - chnkr.d(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.d(:,:),dim,k,ncht); - chnkr.d2(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.d2(:,:),dim,k,ncht); - chnkr.adj(:,nch+1:nch+ncht) = chnkrt.adj+nch; - chnkr.adj(2,nch) = nch+1; - chnkr.adj(1,nch+1) = nch; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; - - if nvals > 0 - ds = bsxfun(@times,reshape(sum((rotmat*chnkrt.d(:,:)).^2,1),k,ncht), ... - (chnkrt.h(:)).'); - dsw = bsxfun(@times,w(:),ds); - dssums = sum(dsw,1); - dssums2 = cumsum([0,dssums(1:end-1)]); - dsint = aint*ds; - dsint = bsxfun(@plus,dsint,dssums2); - lencorner = sum(dsw(:)); - ss = -lencorner/2.0 + dsint; - ss = ss/lencorner*16; - erfss = erf(ss); - datass = reshape((val2(:)-val1(:))/2*((erfss(:)).'+1) ... - +val1(:),nvals,k,ncht); - chnkr.data(:,:,nch+1:nch+ncht) = datass; - end - end - -end - -if ifclosed - nch = chnkr.nch; - chnkr.adj(1,1) = nch; - chnkr.adj(2,nch) = 1; -end - -chnkr = chnkr.refine(cparams); - -% compute normals - -chnkr.n = normals(chnkr); - -end - -function [r,d,d2] = fround(t,m,h,dim) - -[y,dy,d2y] = chnk.spcl.absconvgauss(t,m,0.0,h); - -r = zeros(dim,length(t)); -d = zeros(dim,length(t)); -d2 = zeros(dim,length(t)); - -r(1,:) = t; -r(2,:) = y; -d(1,:) = 1.0; -d(2,:) = dy; -d2(1,:) = 0.0; -d2(2,:) = d2y; - -end +chnkr = chunkerpoly(verts,cparams,pref,edgevals); \ No newline at end of file diff --git a/devtools/test/adapgausswtsTest.m b/devtools/test/adapgausswtsTest.m index 54471c9..02d565b 100644 --- a/devtools/test/adapgausswtsTest.m +++ b/devtools/test/adapgausswtsTest.m @@ -78,7 +78,7 @@ opdims = zeros(2,1); opdims(1) = 1; opdims(2) = 1; -r = chnkr.r; d = chnkr.d; h = chnkr.h; d2 = chnkr.d2; n =chnkr.n; +r = chnkr.r; d = chnkr.d; d2 = chnkr.d2; n =chnkr.n; k = chnkr.k; [t,w] = lege.exps(k); bw = lege.barywts(k); @@ -98,7 +98,7 @@ ntimes = 100; start = tic; for i = 1:ntimes - [mat,maxrecs,numints,iers] = chnk.adapgausswts(r,d,n,d2,h,t,bw,j, ... + [mat,maxrecs,numints,iers] = chnk.adapgausswts(r,d,n,d2,t,bw,j, ... rt,dt,nt,d2t,fkern,opdims,t2,w2,opts); end t1 = toc(start); diff --git a/devtools/test/arclengthfunTest.m b/devtools/test/arclengthfunTest.m new file mode 100644 index 0000000..1d11d70 --- /dev/null +++ b/devtools/test/arclengthfunTest.m @@ -0,0 +1,36 @@ +%ARCLENGTHFUNTEST tests the arclengthfun routine for computing distance +% along length of curve +% + +clearvars; close all; +seed = 8675309; +rng(seed); +addpaths_loc(); + +% geometry parameters and construction + + +cparams = []; +cparams.eps = 1.0e-9; +pref = []; +pref.k = 16; +narms = 0; +amp = 0.0; +start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); + +[s, ~, chnkr] = arclengthfun(chnkr); +% Compare to analytic arclength for circle +ts = squeeze(atan2(chnkr.r(2,:,:), chnkr.r(1,:,:))); +ts(ts<0) = ts(ts<0) + 2*pi; +assert(norm(s-ts) < 1e-12); + +% Now test two circles +chnkrs(2) = chunker(); +chnkrs(1) = chnkr; +rfac = 1.1; +chnkrs(2) = move(chnkr, [0;0], [3;0], 0, rfac); +chnkrtotal = merge(chnkrs); +[s, nchs, ~] = arclengthfun(chnkrtotal); + +assert(norm(s(:,1:nchs(1)) - ts) < 1e-12); +assert(norm(s(:,nchs(1)+1:end) - rfac*ts) < 1e-12); \ No newline at end of file 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/chunker_nearestTest.m b/devtools/test/chunker_nearestTest.m new file mode 100644 index 0000000..54ab6c5 --- /dev/null +++ b/devtools/test/chunker_nearestTest.m @@ -0,0 +1,24 @@ +%chunker.nearest test +% + +iseed = 1234; +rng(iseed); + +% compare closest computed point on circle geometry to exact answer + +circ = chunkerfunc(@(t) chnk.curves.bymode(t,1)); +nt = 1000; +thetas = -pi+2*pi*rand(1,nt); + +% the 0.1 here keeps targets away from coordinate singularity of circle +scal = 0.1+2*rand(1,nt); +targs = [cos(thetas); sin(thetas)].*scal; + +for j = 1:nt + [rn,dn,dist,tn,ichn] = nearest(circ,targs(:,j)); + th1 = atan2(targs(2,j),targs(1,j)); + th2 = atan2(rn(2),rn(1)); + err = min(abs(th1-th2),abs(th1-th2+2*pi)); + err = min(err,abs(th1-th2-2*pi)); + assert(err < 1e-12); +end \ No newline at end of file diff --git a/devtools/test/chunkerfitTest.m b/devtools/test/chunkerfitTest.m new file mode 100644 index 0000000..a0cef29 --- /dev/null +++ b/devtools/test/chunkerfitTest.m @@ -0,0 +1,24 @@ + +%CHUNKERFITTEST + +clearvars; close all; +addpaths_loc(); + +% Sample a smooth curve at random points +rng(0) +n = 20; +tt = sort(2*pi*rand(n,1)); +r = chnk.curves.bymode(tt, [2 0.5 0.2 0.7]); + +opts = []; +opts.ifclosed = true; +opts.cparams = []; +opts.cparams.eps = 1e-6; +opts.pref = []; +opts.pref.k = 16; +chnkr = chunkerfit(r, opts); +assert(checkadjinfo(chnkr) == 0); + +opts.ifclosed = false; +chnkr = chunkerfit(r(:,1:10), opts); +assert(checkadjinfo(chnkr) == 0); diff --git a/devtools/test/chunkerfuncTest.m b/devtools/test/chunkerfuncTest.m index 67de6ca..1820d76 100644 --- a/devtools/test/chunkerfuncTest.m +++ b/devtools/test/chunkerfuncTest.m @@ -83,6 +83,9 @@ a = area(chnkr); assert(abs(a - pi*r^2) < 1e-12,'area wrong for circle domain') +chnkr = refine(chnkr,struct('nover',1)); +a = area(chnkr); +assert(abs(a - pi*r^2) < 1e-12,'area wrong for circle domain') % subplot(1,3,3) diff --git a/devtools/test/chunkerinteriorTest.m b/devtools/test/chunkerinteriorTest.m index 2538be6..281d205 100644 --- a/devtools/test/chunkerinteriorTest.m +++ b/devtools/test/chunkerinteriorTest.m @@ -82,5 +82,11 @@ fprintf('%5.2e s : time for chunkerinterior (chunker, with fmm)\n',t5); assert(all(in5(:) == 1)); - - +% test axissym option +chnkr = chunkerfunc(@(t) starfish(t),struct('ta',-pi/2,'tb',pi/2,'ifclosed',0)); +nt = 1000; +ttarg = -pi/2+pi*rand(nt,1); scal = 2*rand(1,nt); +targs = starfish(ttarg).*scal; + +in = chunkerinterior(chnkr,targs,struct('axissym',true)); +assert(all(in(:) == (scal(:) <= 1))); \ No newline at end of file diff --git a/devtools/test/chunkerkernevalmat_greenlapTest.m b/devtools/test/chunkerkernevalmat_greenlapTest.m index ed4913f..4ffa782 100644 --- a/devtools/test/chunkerkernevalmat_greenlapTest.m +++ b/devtools/test/chunkerkernevalmat_greenlapTest.m @@ -6,7 +6,7 @@ clearvars; close all; seed = 8675309; rng(seed); -%addpaths_loc(); +addpaths_loc(); doadap = false; @@ -81,17 +81,26 @@ start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); Du = Dmat*densu; toc(start) +opts=[]; opts.forceadap=true; +start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); +Du2 = Dmat*densu; +toc(start) +opts=[]; start=tic; Smat = chunkerkernevalmat(chnkr,kerns,targets,opts); Sun = Smat*densun; toc(start) utarg2 = Sun-Du; +utarg3 = Sun-Du2; % relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro'); +assert(relerr < 1e-11); fprintf('relative frobenius error %5.2e\n',relerr); +relerr = norm(utarg-utarg3,'fro')/norm(utarg,'fro'); assert(relerr < 1e-11); + diff --git a/devtools/test/chunkermat_helm2dfmmTest.m b/devtools/test/chunkermat_axissymhelm2dTest.m similarity index 50% rename from devtools/test/chunkermat_helm2dfmmTest.m rename to devtools/test/chunkermat_axissymhelm2dTest.m index dcea14f..3b4985b 100644 --- a/devtools/test/chunkermat_helm2dfmmTest.m +++ b/devtools/test/chunkermat_axissymhelm2dTest.m @@ -1,7 +1,7 @@ -%CHUNKERMAT_HELM2DFMMTEST +%CHUNKERMAT_AXISSYMHELM2DTEST % -% test the matrix builder and the fmm for post processing and do a basic solve +% test the matrix builder and do a basic solve clearvars; close all; iseed = 8675309; @@ -9,34 +9,55 @@ 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 = 32; -narms = 3; +pref.k = 16; +narms = 0; amp = 0.25; -start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); + + +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 = 1.1; +zk = 80.1; % sources ns = 10; -ts = 0.0+2*pi*rand(ns,1); +ts = -pi/2+pi*rand(ns,1); sources = starfish(ts,narms,amp); -sources = 3.0*sources; +sources = 0.5*sources; +sources(1,:) = sources(1,:) + ctr(1); strengths = randn(ns,1); % targets -nt = 10000; -ts = 0.0+2*pi*rand(nt,1); +nt = 100; +ts = -pi/2+pi*rand(nt,1); targets = starfish(ts,narms,amp); -targets = targets.*repmat(rand(1,nt),2,1)*0.5; +targets = targets .* (1 + 0.5*repmat(rand(1, nt), 2, 1)); +targets(1,:) = targets(1,:) + ctr(1); % plot geo and sources @@ -52,9 +73,10 @@ scatter(targets(1,:),targets(2,:),'x') axis equal + % -kerns = @(s,t) chnk.helm2d.kern(zk,s,t,'s'); +kerns = @(s,t) chnk.axissymhelm2d.kern(zk,s,t, [0,0], 's'); % eval u on bdry @@ -75,18 +97,20 @@ % -% build helmholtz dirichlet matrix +% build laplace dirichlet matrix -fkern = @(s,t) chnk.helm2d.kern(zk,s,t,'D'); +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; +sys = 0.5*eye(chnkr.k*chnkr.nch) + D; rhs = ubdry; rhs = rhs(:); -start = tic; sol = gmres(sys,rhs,[],1e-14,100); t1 = toc(start); + +start = tic; sol = gmres(sys,rhs,[],1e-14,200); t1 = toc(start); fprintf('%5.2e s : time for dense gmres\n',t1) @@ -99,22 +123,21 @@ fprintf('difference between direct and iterative %5.2e\n',err) % evaluate at targets and compare -wchnkr = chnkr.wts; -% evaluate at targets using FMM and compare -iffmm = 0; +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) -if(iffmm) - sol_use = sol2.*wchnkr(:); - eps = 1e-6; - pgt = 1; - pot = chnk.helm2d.fmm(eps,zk,chnkr,targets,'D',sol_use); +% - relerr = norm(utarg-pot,'fro')/(sqrt(chnkr.nch)*norm(utarg,'fro')); - relerr2 = norm(utarg-pot,'inf')/dot(abs(sol(:)),wchnkr(:)); - fprintf('relative frobenius error %5.2e\n',relerr); - fprintf('relative l_inf/l_1 error %5.2e\n',relerr2); +wchnkr = chnkr.wts; - assert(relerr < eps); -end +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/chunkermat_truepolygonTest.m b/devtools/test/chunkermat_truepolygonTest.m index e023ff7..e12be30 100644 --- a/devtools/test/chunkermat_truepolygonTest.m +++ b/devtools/test/chunkermat_truepolygonTest.m @@ -119,9 +119,7 @@ % evaluate at targets and compare -opts.usesmooth=false; opts.verb=false; -opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; start=tic; Ssol = chunkerkerneval(chnkr,kerns,sol,targets,opts); t1 = toc(start); fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) diff --git a/devtools/test/bieapplyTest.m b/devtools/test/chunkermatapplyTest.m similarity index 94% rename from devtools/test/bieapplyTest.m rename to devtools/test/chunkermatapplyTest.m index 9ff6b7c..c488754 100644 --- a/devtools/test/bieapplyTest.m +++ b/devtools/test/chunkermatapplyTest.m @@ -1,4 +1,4 @@ -%BIEAPPLYTEST test the routines for a matrix free apply of the system +%CHUNKERMATAPPLYTEST test the routines for a matrix free apply of the system %matrix % % @@ -39,7 +39,7 @@ udense = sys*dens; cormat = chunkermat(chnkr,fkern,struct("corrections",true)); -sysapply = @(sigma) bieapply(chnkr,fkern,-0.5,sigma,cormat); +sysapply = @(sigma) chunkermatapply(chnkr,fkern,-0.5,sigma,cormat); start = tic; u = sysapply(dens); t1 = toc(start); fprintf('%5.2e s : time for matrix free apply\n',t1) @@ -87,7 +87,7 @@ udense = sys*bdry_data; cormat = chunkermat(chnkr,kerns,struct("corrections",true)); -sysapply = @(sigma) bieapply(chnkr,kerns,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(chnkr,kerns,1,sigma,cormat); start = tic; u = sysapply(bdry_data); t1 = toc(start); fprintf('%5.2e s : time for matrix free apply\n',t1) @@ -151,7 +151,7 @@ udense = sys*dens; cormat = chunkermat(cgrph,fkern,struct("corrections",true)); -sysapply = @(sigma) bieapply(cgrph,fkern,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(cgrph,fkern,1,sigma,cormat); start = tic; u = sysapply(dens); t1 = toc(start); fprintf('%5.2e s : time for matrix free apply\n',t1) @@ -199,7 +199,7 @@ udense = sys*bdry_data; cormat = chunkermat(cgrph,kerns,struct("corrections",true)); -sysapply = @(sigma) bieapply(cgrph,kerns,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat); start = tic; u = sysapply(bdry_data); t1 = toc(start); fprintf('%5.2e s : time for matrix free apply\n',t1) diff --git a/devtools/test/chunkgrphrcipTransmissionTest.m b/devtools/test/chunkgrphrcipTransmissionTest.m index 7f1c4cf..4b7f6d0 100644 --- a/devtools/test/chunkgrphrcipTransmissionTest.m +++ b/devtools/test/chunkgrphrcipTransmissionTest.m @@ -84,13 +84,14 @@ opts.charges = cell(1,2); opts.charges{1} = charges{1}; opts.charges{2} = charges{2}; -[kerns,bdry_data] = chnk.helm2d.transmission_helper(cgrph,ks,cs,coefs,opts); + +[kerns, bdry_data] = chnk.helm2d.transmission_helper(cgrph,ks,cs,coefs,opts); opts = []; opts.nonsmoothonly = false; opts.rcip = true; -[sysmat] = chunkermat(cgrph,kerns,opts); +[sysmat] = chunkermat(cgrph, kerns, opts); sysmat = sysmat + eye(size(sysmat,2)); 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); diff --git a/devtools/test/mixedbcTest.m b/devtools/test/mixedbcTest.m index c1bc5d1..89a7b3f 100644 --- a/devtools/test/mixedbcTest.m +++ b/devtools/test/mixedbcTest.m @@ -70,7 +70,7 @@ sol1 = sys\bdrydata; cormat = chunkermat(cgrph,kerns,struct("corrections",true)); -sysapply = @(sigma) bieapply(cgrph,kerns,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat); xapply1 = sys*bdrydata; xapply2 = sysapply(bdrydata); @@ -194,7 +194,7 @@ sol1 = sys\bdrydata; cormat = chunkermat(cgrph,kerns,struct("corrections",true)); -sysapply = @(sigma) bieapply(cgrph,kerns,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat); xapply1 = sys*bdrydata; xapply2 = sysapply(bdrydata); @@ -259,4 +259,4 @@ relerr2 = max(abs(uscat2)); relerr = max(relerr1, relerr2); fprintf('relative field error %5.2e\n',relerr); -assert(relerr < 1e-4) \ No newline at end of file +assert(relerr < 1e-4) diff --git a/devtools/test/rcip_analyticTest.m b/devtools/test/rcip_analyticTest.m deleted file mode 100644 index 4ee6b19..0000000 --- a/devtools/test/rcip_analyticTest.m +++ /dev/null @@ -1,340 +0,0 @@ - -%RCIPTEST -% -% This file tests the rcip routines for solving the exterior dirichlet -% problem on a domain defined by two arcs of circles meeting at two vertices - -clearvars; close all; -addpaths_loc(); - -ncurve = 2; -chnkr(1,ncurve) = chunker(); - -% set wave number -zk = 1.1; - - - - -nch = 5*ones(1,ncurve); - -a = -1.0; -b = 1.0; - -% set angles for top and bottom curve -theta = [pi/2.2 pi/2.4]; - -% set parametrization start and end points for chunkie objects -tab = [-theta(1)/2 -theta(2)/2; theta(1)/2 theta(2)/2]; - -% set parameters for curve parametrization -vert = [a b; 0 0]; -cpars = cell(1,ncurve); -cpars{1}.v0 = vert(:,1); cpars{1}.v1 = vert(:,2); -cpars{1}.theta = theta(1); cpars{1}.ifconvex = 0; cpars{1}.islocal = -1; - -cpars{2}.v0 = vert(:,1); cpars{2}.v1 = vert(:,2); -cpars{2}.theta = theta(2); cpars{2}.ifconvex = 2; cpars{2}.islocal = -1; - -% number of gauss legendre nodes -ngl = 16; -pref = []; -pref.k = ngl; - -cparams = cell(1,ncurve); -for i=1:ncurve - cparams{i}.ta = tab(1,i); - cparams{i}.tb = tab(2,i); - cparams{i}.ifclosed = false; -end - -inds = [0, cumsum(nch)]; - -[glxs,glwts,glu,glv] = lege.exps(ngl); - -fcurve = cell(1,ncurve); -figure(1) -clf; hold on; -for icurve = 1:ncurve - fcurve{icurve} = @(t) circulararc(t,cpars{icurve}); - chnkr(icurve) = chunkerfuncuni(fcurve{icurve},nch(icurve),cparams{icurve},pref); - plot(chnkr(icurve)); quiver(chnkr(icurve)); -end - -iedgechunks = [1 2; 1 chnkr(2).nch]; -isstart = [1 0]; -nedge = 2; -ndim=1; -[Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... - nedge,isstart); - -fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'D'); - -%% - -% sources - -ns = 10; -ts = 0.0+2*pi*rand(1,ns); -sources = [cos(ts); sin(ts)]; -sources = 3.0*sources; -strengths = randn(ns,1); - -% targets - -nt = 20; -ts = 0.0+2*pi*rand(1,nt); -targets = [cos(ts);sin(ts)]; -targets = 0.2*targets; - -scatter(sources(1,:),sources(2,:),'o'); -scatter(targets(1,:),targets(2,:),'x'); -axis equal - - - - -chnkrtotal = merge(chnkr); -fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'D'); -np = chnkrtotal.k*chnkrtotal.nch; -start = tic; D = chunkermat(chnkrtotal,fkern); -t1 = toc(start); - - - -kerns = @(s,t) chnk.helm2d.kern(zk,s,t,'s'); - -% eval u on bdry - -targs = chnkrtotal.r; targs = reshape(targs,2,chnkrtotal.k*chnkrtotal.nch); -targstau = tangents(chnkrtotal); -targstau = reshape(targstau,2,chnkrtotal.k*chnkrtotal.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; - - - -fprintf('%5.2e s : time to assemble matrix\n',t1) - - -sysmat = D; - -RG = speye(np); -ncorner = 2; -corners= cell(1,ncorner); -R = cell(1,ncorner); - - -corners{1}.clist = [1,2]; -corners{1}.isstart = [1,0]; -corners{1}.nedge = 2; - - -corners{2}.clist = [1,2]; -corners{2}.isstart = [0 1]; -corners{2}.nedge = 2; - - -ndim = 1; - -nsub = 100; - -opts = []; - -for icorner=1:ncorner - clist = corners{icorner}.clist; - isstart = corners{icorner}.isstart; - nedge = corners{icorner}.nedge; - cparslocal = cell(1,nedge); - fcurvelocal = cell(1,nedge); - h0 = zeros(1,nedge); - starind = []; - for i=1:nedge - cparslocal{i} = cpars{clist(i)}; - cparslocal{i}.islocal = isstart(i); - fcurvelocal{i} = @(t) circulararc(t,cparslocal{i}); - if(isstart(i)) - h0(i) = chnkr(clist(i)).h(1); - starind = [starind inds(clist(i))*ngl*ndim+(1:2*ngl*ndim)]; - else - h0(i) = chnkr(clist(i)).h(end); - starind = [starind inds(clist(i)+1)*ngl*ndim-fliplr(0:2*ngl*ndim-1)]; - end - end - h0 = h0*2; - - [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... - nedge,isstart); - opts_use = []; - tic; R{icorner} = chnk.rcip.Rcomp(ngl,nedge,ndim,Pbc, ... - PWbc,nsub,starL,circL, ... - starS,circS,ilist,h0,isstart,fcurvelocal,glxs,fkern); - toc - - sysmat(starind,starind) = inv(R{icorner}) - eye(2*ngl*nedge*ndim); -end -sysmat = sysmat + eye(np); - -sol = gmres(sysmat,ubdry,np,eps*20,np); - -opts.usesmooth=true; -opts.verb=false; -opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; -start=tic; Dsol = chunkerkerneval(chnkrtotal,fkern,sol,targets,opts); -t1 = toc(start); -fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) - -% - -wchnkr = chnkrtotal.wts; - -relerr = norm(utarg-Dsol,'fro')/(sqrt(chnkrtotal.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); - - -%% - - -%% -%---------------------------------------- -% -% -% Auxiliary routines for generating boundary -% - -function [r,d,d2] = circulararc(t,cpars) -%%circulararc -% return position, first and second derivatives of a circular arc -% that passes through points (x,y)=(a,0) and (x,y)=(b,0) with opening -% angle theta0. -% -% Inputs: -% t - paramter values (-theta0/2,theta0/2) to evaluate these quantities -% -% Outputs: -% r - coordinates -% d - first derivatives w.r.t. t -% d2 - second derivatives w.r.t. t - -v0 = cpars.v0; -v1 = cpars.v1; -theta0 = cpars.theta; -ifconvex = cpars.ifconvex; - -a = v0(1); -b = v1(1); - -islocal = -1; -if isfield(cpars,'islocal') - islocal = cpars.islocal; -end - -if islocal == -1 - cx = (a+b)/2; - r0 = (b-a)/2/sin(theta0/2); - - if ifconvex - cy = v0(2)+(b-a)/2/tan(theta0/2); - theta = 3*pi/2+t; - else - cy = v0(2)-(b-a)/2/tan(theta0/2); - theta = pi/2+t; - end - - - xs = r0*cos(theta); - ys = r0*sin(theta); - - xp = -ys; - yp = xs; - - xpp = -xs; - ypp = -ys; - - xs = cx + xs; - ys = cy + ys; - -else - - r0 = (b-a)/2/sin(theta0/2); - if ~ifconvex - if islocal == 0 - cx = r0*cos(pi/2+theta0/2); - sx = r0*sin(pi/2+theta0/2); - elseif islocal == 1 - cx = r0*cos(pi/2-theta0/2); - sx = r0*sin(pi/2-theta0/2); - end - elseif ifconvex - if islocal == 0 - cx = r0*cos(3*pi/2+theta0/2); - sx = r0*sin(3*pi/2+theta0/2); - elseif islocal == 1 - cx = r0*cos(3*pi/2-theta0/2); - sx = r0*sin(3*pi/2-theta0/2); - end - end - -% t2 = t.*t; -% t3 = t2.*t; -% t4 = t3.*t; -% t5 = t4.*t; -% t6 = t5.*t; -% t7 = t6.*t; -% t8 = t7.*t; -% t9 = t8.*t; -% t10 = t9.*t; -% -% -% ct = -t2/2 + t4/24 - t6/720 + t8/40320 - t10/3628800; -% st = t - t3/6 + t5/120 - t7/5040 + t9/362880; - -% better to use sin(t) directly instead of its series expansion because -% 1. the evaluation of sin(t) is accurate at t=0; -% 2. one needs to determine how many terms are needed in series expansion, -% which looks simple but it could be either inefficient or inaccurate -% depending on the value of t; 3. if we write "if" statements, it's -% difficult to vectorize the evaluation. - -% same reason that we should use -2*sin(t/2).^2 instead of its Taylor -% expansion, but of course one should NOT use cos(t)-1. - ct = -2*sin(t/2).^2; - st = sin(t); - - ctp = -st; - stp = cos(t); - - ctpp = -stp; - stpp = -st; - - xs = -sx*st + cx*ct; - ys = cx*st + sx*ct; - - xp = -sx*stp + cx*ctp; - yp = cx*stp + sx*ctp; - - xpp = -sx*stpp + cx*ctpp; - ypp = cx*stpp + sx*ctpp; -end - - -r = [(xs(:)).'; (ys(:)).']; -d = [(xp(:)).'; (yp(:)).']; -d2 = [(xpp(:)).'; (ypp(:)).']; - -end - - - diff --git a/devtools/test/refine_trefinementTest.m b/devtools/test/refine_trefinementTest.m deleted file mode 100644 index fde86aa..0000000 --- a/devtools/test/refine_trefinementTest.m +++ /dev/null @@ -1,83 +0,0 @@ -%REFINE_TREFINEMENTTEST -% -% get chunker description of a starfish domain. check if that domain -% satisfies a level restriction in the underlying parameter space. -% call a refinement routine to fix. - -clearvars; close all; -iseed = 8675309; -rng(iseed); - -addpaths_loc(); - -% define curve - -cparams = []; -cparams.eps = 1.0e-10; -cparams.nover = 0; -pref = []; -pref.k = 16; -narms = 5; -amp = 0.25; - -% form curve without strict enforcement of level restriction -% in h - -start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); -t1 = toc(start); -fprintf('%5.2e s : time to build geo\n',t1) -fprintf('originally, nch = %d\n',chnkr.nch) - -% check for violations of level restriction - -happy = true; -for i = 1:chnkr.nch - hself = chnkr.h(i); - i1 = chnkr.adj(1,i); - i2 = chnkr.adj(2,i); - h1 = hself; h2 = hself; - if (i1 > 0) - h1 = chnkr.h(i1); - end - if (i2 > 0) - h2 = chnkr.h(i2); - end - if (hself > 2*h1) - happy = false; - end - if (hself > 2*h2) - happy = false; - end -end - -% enforce level restriction using refine code with -% specific options - -opts.lvlr = 't'; % refine in t rather than arclength -% standard is 2.05 which is too loose for our purposes -opts.lvlrfac = 1.99; - -start = tic; chnkr = refine(chnkr,opts); t1 = toc(start); -fprintf('%5.2e s : time to refine geo\n',t1) -fprintf('now, nch = %d\n',chnkr.nch) - -happy = true; -for i = 1:chnkr.nch - hself = chnkr.h(i); - i1 = chnkr.adj(1,i); - i2 = chnkr.adj(2,i); - h1 = hself; h2 = hself; - if (i1 > 0) - h1 = chnkr.h(i1); - end - if (i2 > 0) - h2 = chnkr.h(i2); - end - if (hself > 2*h1) - happy = false; - end - if (hself > 2*h2) - happy = false; - end -end -