diff --git a/chunkie/+chnk/pquadwts.m b/chunkie/+chnk/pquadwts.m index 59bc795..242c8fa 100644 --- a/chunkie/+chnk/pquadwts.m +++ b/chunkie/+chnk/pquadwts.m @@ -1,56 +1,78 @@ -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 +function [varargout] = pquadwts(r,d,n,d2,wts,j,... + rt,t,w,opts,intp_ab,intp,types) +%CHNK.pquadwts product integration for interaction of kernel on chunk % at targets % % WARNING: this routine is not designed to be user-callable and assumes % a lot of precomputed values as input % -% Syntax: mat = interpquadwts(r,d,d2,h,ct,bw,j, ... -% rt,dt,d2t,kern,opdims,t,w,opts) +% Syntax: [varargout] = pquadwts(r,d,n,d2,wts,j,... +% rt,t,w,opts,intp_ab,intp,types) % % Input: % r - chnkr nodes % d - chnkr derivatives at nodes % n - chnkr normals at nodes % d2 - chnkr 2nd derivatives at nodes -% ct - Legendre nodes at order of chunker -% bw - barycentric interpolation weights for Legendre nodes at order of -% chunker +% wts - chnkr integration weights for smooth functions % j - chunk of interest -% rt,dt,nt,d2t - position, derivative, normal,second derivative of select -% target points. if any are not used by kernel (or not well -% defined, e.g. when not on curve), a dummy array -% of the appropriate size should be supplied -% kern - kernel function of form kern(srcinfo,targinfo) -% opdims - dimensions of kernel -% t - (Legendre) integration nodes for adaptive integration -% w - integration nodes for adaptive integrator (t and w not necessarily -% same order as chunker order) +% rt - position of target points. if any are not used by kernel +% (or not well defined, e.g. when not on curve), a +% dummy array of the appropriate size should be supplied +% t - (Legendre) integration nodes +% w - (Legendre) integration weights +% opts - opts.side +% intp_ab - panel endpoints interpolation matrix +% types - specified singularity types +% [0 0 0 0] - smooth quadr +% [1 0 0 0] - log singularity +% [0 0 -1 0] - cauchy singularity +% each [a, b, c, d] +% corresponds to (log(z-w))^a (log(zc-wc))^b (z-w)^c (zc-wc)^d % % Output -% mat - integration matrix +% varargout - integration matrices for specified singularity types % % Helsing-Ojala (interior/exterior?) - 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 = (intp*(d(1,:,j)'+1i*d(2,:,j)')); % r' -d2_i = (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'' +wts_i = wts(:,j)'; % 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 -wxp_i = w.*d_i; % complex speed weights (Helsing's wzp) +wxp_i = w.*d_i; % complex speed weights (Helsing's wzp) + +varargout = cell(size(types)); + +% [As, Ad, A1, A2, A3, A4] = SDspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),... +% struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side); +[As, Ad] = SDspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),... + struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side); +As = As*intp; +Ad = Ad*intp; -mat_ho_slp = Sspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),... - struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),'e'); -mat_ho_dlp = Dspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),... - struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),'e'); +for j = 1:length(types) + type0 = types{j}; -mat = (mat_ho_slp+real(mat_ho_dlp))*intp; % depends on kern, different mat1? + if (all(type0 == [0 0 0 0])) + varargout{j} = ones(size(rt,2),numel(wts_i)).*wts_i; + elseif (all(type0 == [1 0 0 0])) + % varargout{j} = Sspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),... + % struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side)*intp; + varargout{j} = As; + elseif (all(type0 == [0 0 -1 0])) + % varargout{j} = Dspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),... + % struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side)*intp; + varargout{j} = Ad; + else + error("split panel quad type " + join(string([1 2 3]),",") + " not available"); + end end +end function [A, A1, A2] = Sspecialquad(t,s,a,b,side) % SSPECIALQUAD - SLP val+grad close-eval Helsing "special quadrature" matrix @@ -66,6 +88,7 @@ % Efficient only if multiple targs, since O(p^3). % See Helsing-Ojala 2008 (special quadr Sec 5.1-2), Helsing 2009 mixed (p=16), % and Helsing's tutorial demo11b.m LGIcompRecFAS() +% See https://github.com/ahbarnett/BIE2D/blob/master/panels/LapSLP_closepanel.m if nargin<5, side = 'i'; end % interior or exterior zsc = (b-a)/2; zmid = (b+a)/2; % rescaling factor and midpoint of src segment y = (s.x-zmid)/zsc; x = (t.x-zmid)/zsc; % transformed src nodes, targ pts @@ -146,6 +169,7 @@ % Efficient only if multiple targs, since O(p^3). % See Helsing-Ojala 2008 (special quadr Sec 5.1-2), Helsing 2009 mixed (p=16), % and Helsing's tutorial demo11b.m M1IcompRecFS() +% See https://github.com/ahbarnett/BIE2D/blob/master/panels/LapDLP_closepanel.m if nargin<5, side = 'i'; end % interior or exterior zsc = (b-a)/2; zmid = (b+a)/2; % rescaling factor and midpoint of src segment y = (s.x-zmid)/zsc; x = (t.x-zmid)/zsc; % transformed src nodes, targ pts @@ -203,3 +227,76 @@ end end end + +function [As, A, A1, A2, A3, A4] = SDspecialquad(t,s,a,b,side) +% S+D together... (with Sn) should be the same as requesting S or D +% more expensive if request Dn... +% https://github.com/ahbarnett/BIE2D/blob/master/panels/LapDLP_closepanel.m +%%% d = 100; +zsc = (b-a)/2; zmid = (b+a)/2; % rescaling factor and midpoint of src segment +y = (s.x-zmid)/zsc; x = (t.x-zmid)/zsc; % transformed src nodes, targ pts +%figure; plot(x,'.'); hold on; plot(y,'+-'); plot([-1 1],[0 0],'ro'); % debug +N = numel(x); % # of targets +p = numel(s.x); % assume panel order is # nodes +if N*p==0 + As = 0; A = 0; A1=0; A2=0; + return +end +c = (1-(-1).^(1:p))./(1:p); % Helsing c_k, k = 1..p. +V = ones(p,p); for k=2:p, V(:,k) = V(:,k-1).*y; end % Vandermonde mat @ nodes +P = zeros(p+1,N); % Build P, Helsing's p_k vectorized on all targs... +d = 1.1; inr = abs(x)<=d; ifr = abs(x)>d; % near & far treat separately +%gam = 1i; +gam = exp(1i*pi/4); % smaller makes cut closer to panel. barnett 4/17/18 +if side == 'e', gam = conj(gam); end % note gam is a phase, rots branch cut +P(1,:) = log(gam) + log((1-x)./(gam*(-1-x))); % init p_1 for all targs int + +% upwards recurrence for near targets, faster + more acc than quadr... +% (note rotation of cut in log to -Im; so cut in x space is lower unit circle) +if N>1 || (N==1 && inr==1) % Criterion added by Bowei Wu 03/05/15 to ensure inr not empty + for k=1:p + P(k+1,inr) = x(inr).'.*P(k,inr) + c(k); + end % recursion for p_k +end +% fine quadr (no recurrence) for far targets (too inaccurate cf downwards)... +Nf = numel(find(ifr)); wxp = s.wxp/zsc; % rescaled complex speed weights +if Nf>0 % Criterion added by Bowei Wu 03/05/15 to ensure ifr not empty + P(end,ifr) = sum(((wxp.*(V(:,end).*y(:)))*ones(1,Nf))./bsxfun(@minus,y,x(ifr).')); % int y^p/(y-x) + for ii = p:-1:2 + P( ii,ifr) = (P(ii+1,ifr)-c(ii))./x(ifr).'; + end +end + +Q = zeros(p,N); % compute q_k from p_k via Helsing 2009 eqn (18)... (p even!) +% Note a rot ang appears here too... 4/17/18 +%gam = exp(1i*pi/4); % 1i; % moves a branch arc as in p_1 +%if side == 'e', gam = conj(gam); end % note gam is a phase, rots branch cut +Q(1:2:end,:) = P(2:2:end,:) - repmat(log((1-x.').*(-1-x.')),[ceil(p/2) 1]); % guessed! +% (-1)^k, k odd, note each log has branch cut in semicircle from -1 to 1 +% Q(2:2:end,:) = P(3:2:end,:) - repmat(log((1-x.')./((-1-x.'))),[floor(p/2) 1]); % same cut as for p_1 +Q(2:2:end,:) = P(3:2:end,:) - repmat(log(gam) + log((1-x.')./(gam*(-1-x.'))),[floor(p/2) 1]); % same cut as for p_1 +Q = Q.*repmat(1./(1:p)',[1 N]); % k even +As = real((V.'\Q).'.*repmat((1i*s.nx)',[N 1])*zsc)/(2*pi*abs(zsc)); +As = As*abs(zsc) - log(abs(zsc))/(2*pi)*repmat(abs(s.wxp)',[N 1]); % unscale, yuk +warning('off','MATLAB:nearlySingularMatrix'); +% A = real((V.'\P).'*(1i/(2*pi))); % solve for special quadr weights +A = ((V.'\P(1:p,:)).'*(1i/(2*pi))); % do not take real for the eval of Stokes DLP non-laplace term. Bowei 10/19/14 +%A = (P.'*inv(V))*(1i/(2*pi)); % equiv in exact arith, but not bkw stable. +if nargout>2 + R = -(kron(ones(p,1),1./(1-x.')) + kron((-1).^(0:p-1).',1./(1+x.'))) +... + repmat((0:p-1)',[1 N]).*[zeros(1,N); P(1:p-1,:)]; % hypersingular kernel weights of Helsing 2009 eqn (14) + Az = (V.'\R).'*(1i/(2*pi*zsc)); % solve for targ complex-deriv mat & rescale + A1 = Az; + if nargout > 3 + S = -(kron(ones(p,1),1./(1-x.').^2) - kron((-1).^(0:p-1).',1./(1+x.').^2))/2 +... + repmat((0:p-1)',[1 N]).*[zeros(1,N); R(1:p-1,:)]/2; % supersingular kernel weights + Azz = (V.'\S).'*(1i/(2*pi*zsc^2)); + if nargout > 4 + A1 = real(Az); A2 = -imag(Az); % note sign for y-deriv from C-deriv + A3 = real(Azz); A4 = -imag(Azz); + else + A1 = Az; A2 = Azz; + end + end +end +end diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m index 057b5d3..80e3c2a 100644 --- a/chunkie/@kernel/helm2d.m +++ b/chunkie/@kernel/helm2d.m @@ -40,12 +40,20 @@ obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 's'); obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 's', sigma); obj.sing = 'log'; + obj.splitinfo = []; + obj.splitinfo.type = {[0 0 0 0],[1 0 0 0]}; + obj.splitinfo.action = {'r','r'}; + obj.splitinfo.functions = @(s,t) helm2d_s_split(zk,s,t); case {'d', 'double'} obj.type = 'd'; obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'd'); obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'd', sigma); obj.sing = 'log'; + obj.splitinfo = []; + obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]}; + obj.splitinfo.action = {'r','r','r'}; + obj.splitinfo.functions = @(s,t) helm2d_d_split(zk,s,t); case {'sp', 'sprime'} obj.type = 'sp'; @@ -69,6 +77,10 @@ obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'c', coefs); obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs); obj.sing = 'log'; + obj.splitinfo = []; + obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]}; + obj.splitinfo.action = {'r','r','r'}; + obj.splitinfo.functions = @(s,t) helm2d_c_split(zk,s,t,coefs); case {'cp', 'cprime'} if ( nargin < 3 ) @@ -93,3 +105,37 @@ end + +function f = helm2d_s_split(zk,s,t) +seval = chnk.helm2d.kern(zk, s, t, 's'); +dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'); +logeval = log(abs(dist)); +f = cell(2,1); +f{1} = seval + 1/(2*pi)*4*logeval.*imag(seval); +f{2} = 4*imag(seval); +end + +function f = helm2d_d_split(zk,s,t) +deval = chnk.helm2d.kern(zk, s, t, 'd'); +dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'); +logeval = log(abs(dist)); +cauchyeval = (s.n(1,:)+1i*s.n(2,:))./(dist); +f = cell(3, 1); +f{1} = deval + 1/(2*pi)*4*logeval.*imag(deval) + 1/(2*pi)*real(cauchyeval); +f{2} = 4*imag(deval); +f{3} = ones(size(t.r,2),size(s.r,2)); +end + +function f = helm2d_c_split(zk,s,t,coefs) +seval = chnk.helm2d.kern(zk, s, t, 's'); +deval = chnk.helm2d.kern(zk, s, t, 'd'); +dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'); +logeval = log(abs(dist)); +cauchyeval = (s.n(1,:)+1i*s.n(2,:))./(dist); +f = cell(3, 1); +f{1} = coefs(1) * (deval + 2/pi*logeval.*imag(deval) + 1/(2*pi)*real(cauchyeval)) + ... + coefs(2) * (seval + 2/pi*logeval.*imag(seval)); +f{2} = coefs(1) * (4*imag(deval)) + coefs(2) * (4*imag(seval)); +f{3} = coefs(1)*ones(size(t.r,2),size(s.r,2)); +end + diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m index e7317bd..df7a5dd 100644 --- a/chunkie/@kernel/lap2d.m +++ b/chunkie/@kernel/lap2d.m @@ -42,12 +42,20 @@ obj.eval = @(s,t) chnk.lap2d.kern(s, t, 's'); obj.fmm = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 's', sigma); obj.sing = 'log'; + obj.splitinfo = []; + obj.splitinfo.type = {[1 0 0 0]}; + obj.splitinfo.action = {'r'}; + obj.splitinfo.functions = @(s,t) lap2d_s_split(s,t); case {'d', 'double'} obj.type = 'd'; obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'd'); obj.fmm = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 'd', sigma); obj.sing = 'smooth'; + obj.splitinfo = []; + obj.splitinfo.type = {[0 0 -1 0]}; + obj.splitinfo.action = {'r'}; + obj.splitinfo.functions = @(s,t) lap2d_d_split(s,t); case {'sp', 'sprime'} obj.type = 'sp'; @@ -91,6 +99,10 @@ obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'c', coefs); obj.fmm = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 'c', sigma, coefs); obj.sing = 'log'; + obj.splitinfo = []; + obj.splitinfo.type = {[1 0 0 0],[0 0 -1 0]}; + obj.splitinfo.action = {'r','r'}; + obj.splitinfo.functions = @(s,t) lap2d_c_split(s,t,coefs); otherwise error('Unknown Laplace kernel type ''%s''.', type); @@ -103,3 +115,20 @@ end end + +function f = lap2d_s_split(s,t) +f = cell(1,1); +f{1} = ones(size(t.r,2),size(s.r,2)); +end + +function f = lap2d_d_split(s,t) +f = cell(1,1); +f{1} = ones(size(t.r,2),size(s.r,2)); +end + +function f = lap2d_c_split(s,t,coefs) +f = cell(2,1); +f{1} = coefs(2)*ones(size(t.r,2),size(s.r,2)); +f{2} = coefs(1)*ones(size(t.r,2),size(s.r,2)); +end + diff --git a/chunkie/@kernel/stok2d.m b/chunkie/@kernel/stok2d.m index 236cbdb..b5bb102 100644 --- a/chunkie/@kernel/stok2d.m +++ b/chunkie/@kernel/stok2d.m @@ -78,6 +78,10 @@ obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'svel', sigma); obj.opdims = [2, 2]; obj.sing = 'log'; + obj.splitinfo = []; + obj.splitinfo.type = {[1 0 0 0],[0 0 -1 0],[0 0 -1 0]}; + obj.splitinfo.action = {'r','r','i'}; + obj.splitinfo.functions = @(s,t) stok2d_s_split(mu, s, t); case {'spres', 'spressure'} obj.type = 'spres'; @@ -188,3 +192,27 @@ end + +function f = stok2d_s_split(mu,s,t) +dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'); +distr = real(dist); +disti = imag(dist); +f = cell(3, 1); +ntarg = numel(t.r(1,:)); +nsrc = numel(s.r(1,:)); +sn1mat = repmat(s.n(1,:),[ntarg 1]); +sn2mat = repmat(s.n(2,:),[ntarg 1]); +f{1} = zeros(2*ntarg,2*nsrc); +f{2} = zeros(2*ntarg,2*nsrc); +f{3} = zeros(2*ntarg,2*nsrc); +f{1}(1:2:end,1:2:end) = 1/(2*mu); +f{1}(2:2:end,2:2:end) = 1/(2*mu); +f{2}(1:2:end,1:2:end) = -sn1mat.*distr/(2*mu); +f{2}(1:2:end,2:2:end) = -sn2mat.*distr/(2*mu); +f{2}(2:2:end,1:2:end) = -sn1mat.*disti/(2*mu); +f{2}(2:2:end,2:2:end) = -sn2mat.*disti/(2*mu); +f{3}(1:2:end,1:2:end) = -sn2mat.*distr/(2*mu); +f{3}(1:2:end,2:2:end) = sn1mat.*distr/(2*mu); +f{3}(2:2:end,1:2:end) = -sn2mat.*disti/(2*mu); +f{3}(2:2:end,2:2:end) = sn1mat.*disti/(2*mu); +end diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index d12ae14..46d6aea 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -116,6 +116,7 @@ if isfield(opts,'forcesmooth'); opts_use.forcesmooth = opts.forcesmooth; end if isfield(opts,'forceadap'); opts_use.forceadap = opts.forceadap; end if isfield(opts,'forcepquad'); opts_use.forcepquad = opts.forcepquad; end +if isfield(opts,'side'); opts_use.side = opts.side; end if isfield(opts,'flam') opts_use.flam = opts.flam; end @@ -266,7 +267,12 @@ fints(irow0) = fints(irow0) + chnk.chunkerkerneval_smooth(chnkr0,kern0,opdims0,dens0,targinfo0, ... flag,opts_use); - fints(irow0) = fints(irow0) + chunkerkerneval_ho(chnkr0,kern0,opdims0,dens0, ... + if ~isfield(opts_use,'side') + msg = "Error: for pquad, set opts.side to 'i' or 'e' "; + error(msg) + end + + fints(irow0) = fints(irow0) + chunkerkerneval_pquad(chnkr0,kern0,opdims0,dens0, ... targinfo0,flag,opts_use); return @@ -298,13 +304,18 @@ end -function fints = chunkerkerneval_ho(chnkr,kern,opdims,dens, ... +function fints = chunkerkerneval_pquad(chnkr,kern,opdims,dens, ... targinfo,flag,opts) +if ~isa(kern,'kernel') || isempty(kern.splitinfo) + error('Helsing-Ojala quad only available for kernel class objects with splitinfo defined'); +end + % target [~,nt] = size(targinfo.r); fints = zeros(opdims(1)*nt,1); -k = size(chnkr.r,2); +k = chnkr.k; +nch = chnkr.nch; [t,w] = lege.exps(2*k); ct = lege.exps(k); bw = lege.barywts(k); @@ -312,6 +323,10 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; +wts = chnkr.wts; + +assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size') +dens = reshape(dens,opdims(2),k,nch); % interpolation matrix intp = lege.matrin(k,t); % interpolation from k to 2*k @@ -319,30 +334,53 @@ targs = targinfo.r; -targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); -targn = zeros(chnkr.dim,nt); -if isfield(targinfo, 'd') - targd = targinfo.d; -end - -if isfield(targinfo, 'd2') - targd2 = targinfo.d2; -end - -if isfield(targinfo, 'n') - targn = targinfo.n; -end for j=1:size(chnkr.r,3) [ji] = find(flag(:,j)); if(~isempty(ji)) idxjmat = (j-1)*k+(1:k); + targinfoji = []; + targinfoji.r = targinfo.r(:,ji); + if isfield(targinfo, 'd') + targinfoji.d = targinfo.d(:,ji); + end - % Helsing-Ojala (interior/exterior?) - 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? + if isfield(targinfo, 'd2') + targinfoji.d2 = targinfo.d2(:,ji); + end + + if isfield(targinfo, 'n') + targinfoji.n = targinfo.n(:,ji); + end + + srcinfo = []; + srcinfo.r = r(:,:,j); + srcinfo.d = d(:,:,j); + srcinfo.d2 = d2(:,:,j); + srcinfo.n = n(:,:,j); + densj = dens(:,:,j); + indji = (ji-1)*opdims(1); + ind = (indji(:)).' + (1:opdims(1)).'; - fints(ji) = fints(ji) + mat1*dens(idxjmat); + % Helsing-Ojala (interior/exterior?) + allmats = cell(size(kern.splitinfo.type)); + [allmats{:}] = chnk.pquadwts(r,d,n,d2,wts,j,targs(:,ji), ... + t,w,opts,intp_ab,intp,kern.splitinfo.type); + + funs = kern.splitinfo.functions(srcinfo,targinfoji); + for l = 1:length(allmats) + switch kern.splitinfo.action{l} + case 'r' % real part + mat0 = real(allmats{l}); + case 'i' % imaginary part + mat0 = imag(allmats{l}); + case 'c' % complex + mat0 = allmats{l}; + end + mat0opdim = kron(mat0,ones(opdims')); + mat0xsplitfun = mat0opdim.*funs{l}; + fints(ind(:)) = fints(ind(:)) + mat0xsplitfun*densj(:); + end end end end diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m index 4e8da6f..a47ba62 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/chunkie/chunkerkernevalmat.m @@ -174,8 +174,8 @@ if forcepquad optsflag = []; optsflag.fac = fac; flag = flagnear(chnkr,targinfo.r,optsflag); - spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ... - targinfo,flag,optsadap); + spmat = chunkerkernevalmat_pquad(chnkr,kern,opdims, ... + targinfo,flag,opts); mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... flag,opts); mat = mat + spmat; @@ -343,7 +343,7 @@ end -function mat = chunkerkernevalmat_ho(chnkr,kern,opdims, ... +function mat = chunkerkernevalmat_pquad(chnkr,kern,opdims, ... targinfo,flag,opts) k = chnkr.k; @@ -391,6 +391,7 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; + wts = chnkr.wts; % interpolation matrix intp = lege.matrin(k,t); % interpolation from k to 2*k @@ -402,10 +403,36 @@ [ji] = find(flag(:,i)); + targinfoji = []; + targinfoji.r = targinfo.r(:,ji); + + srcinfo = []; + srcinfo.r = r(:,:,i); + srcinfo.d = d(:,:,i); + srcinfo.d2 = d2(:,:,i); + srcinfo.n = n(:,:,i); + % Helsing-Ojala (interior/exterior?) - 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? - + allmats = cell(size(kern.splitinfo.type)); + [allmats{:}] = chnk.pquadwts(r,d,n,d2,wts,i,targs(:,ji), ... + t,w,opts,intp_ab,intp,kern.splitinfo.type); + + mat1 = zeros(size(allmats{1})); + funs = kern.splitinfo.functions(srcinfo,targinfoji); + for l = 1:length(allmats) + switch kern.splitinfo.action{l} + case 'r' + mat0 = real(allmats{l}); + case 'i' + mat0 = imag(allmats{l}); + case 'c' + mat0 = allmats{l}; + end + mat0opdim = kron(mat0,ones(opdims)); + mat0xsplitfun = mat0opdim.*funs{l}; + mat1 = mat1 + mat0xsplitfun; + end + js1 = jmat:jmatend; js1 = repmat( (js1(:)).',opdims(1)*numel(ji),1); diff --git a/devtools/test/chunkermat_quadadap_closetotouchingTest.m b/devtools/test/chunkermat_quadadap_closetotouchingTest.m index 450b309..3f1bafa 100644 --- a/devtools/test/chunkermat_quadadap_closetotouchingTest.m +++ b/devtools/test/chunkermat_quadadap_closetotouchingTest.m @@ -72,20 +72,23 @@ % use adaptive routine to build matrix (self done by ggq, nbor by adaptive) eta = 1; -fkern = @(s,t) chnk.lap2d.kern(s,t,'C',[1,eta]); +fkern = kernel('laplace','c',[1,eta]); type = 'log'; opts = []; opts.robust = true; +opts.adaptive_correction = true; opdims = [1 1]; start = tic; -mata = chnk.quadadap.buildmat(chnkr,fkern,opdims,type,opts); +mata = chunkermat(chnkr,fkern,opts); t1 = toc(start); fprintf('%5.2e s : time to assemble matrix (adaptive)\n',t1) +opts.adaptive_correction = false; start = tic; mato = chunkermat(chnkr,fkern); t1 = toc(start); fprintf('%5.2e s : time to assemble matrix (original GGQ)\n',t1) start = tic; opts.forcepquad = true; targs_bd = [chnkr.r(1,:);chnkr.r(2,:)]; +opts.side = 'e'; matho = chunkerkernevalmat(chnkr,fkern,targs_bd,opts); t1 = toc(start); fprintf('%5.2e s : time to assemble matrix (Helsing-Ojala)\n',t1) diff --git a/devtools/test/chunkermat_stok2dTest.m b/devtools/test/chunkermat_stok2dTest.m index cc18225..0d1b2ed 100644 --- a/devtools/test/chunkermat_stok2dTest.m +++ b/devtools/test/chunkermat_stok2dTest.m @@ -38,6 +38,7 @@ ts = 0.0+2*pi*rand(nt,1); targets = starfish(ts,narms,amp); targets = targets.*repmat(rand(1,nt),2,1)*0.75; +% targets = targets.*repmat(rand(1,nt),2,1)*0.99; plot(chnkr, 'r.'); hold on; plot(targets(1,:), targets(2,:), 'kx') @@ -93,6 +94,19 @@ t1 = toc(start); fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) +% just test stokes slp velocity pquad... +% targets = targets.*repmat(rand(1,nt),2,1)*0.99 seems to break below stokes kernel test... +fkerns = kernel('stok', 'svel', mu); +Ssol = chunkerkerneval(chnkr,fkerns,sol,targets,opts); +opts.forcepquad=true; +opts.side = 'i'; +Ssol_pquad = chunkerkerneval(chnkr,fkerns,sol,targets,opts); +err = abs(Ssol - Ssol_pquad); +assert(norm(err) < 1e-10); +% figure(1),clf, +% plot(chnkr), hold on +% scatter(targets(1,:),targets(2,:),[],log10(err(1:2:end))); colorbar +opts.forcepquad=false; % wchnkr = chnkr.wts; diff --git a/devtools/test/pquadTest.m b/devtools/test/pquadTest.m new file mode 100644 index 0000000..aa6d9a1 --- /dev/null +++ b/devtools/test/pquadTest.m @@ -0,0 +1,119 @@ +%DEMO_BARBELL_2 +% +% Define a polygonal "barbell" shaped domain with +% prescribed constant boundary data on each edge. +% Solves the corresponding Helmholtz Dirichlet problem +% using corner rounding and smoothing the boundary data +% across the rounded corners +% + +% profile clear +% profile on + +clearvars; close all; +iseed = 8675309; +rng(iseed,'twister'); +addpaths_loc(); + +% planewave vec + +kvec = 10*[1;-1.5]; + +% + +zk = norm(kvec); + +% define geometry and boundary conditions +% (vertices defined by another function) +verts = chnk.demo.barbell(2.0,2.0,1.0,1.0); +nv = size(verts,2); +edgevals = randn(1,nv); + +% 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-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 = refine(chnkr,struct('nover',1)); +% +% plot smoothed geometry and data + +figure(1) +clf +plot(chnkr,'-x') +hold on +quiver(chnkr) +axis equal + +figure(2) +clf +nchplot = 1:chnkr.nch; +plot3(chnkr,1) + +% solve and visualize the solution + +% build laplace dirichlet matrix + +% fkern = kernel('laplace','d'); +fkern = kernel('helmholtz','c',zk,[1.5,-2i]); +opts = []; +start = tic; C = chunkermat(chnkr,fkern,opts); +t1 = toc(start); + +fprintf('%5.2e s : time to assemble matrix\n',t1) + +% + +sys = -0.5*eye(chnkr.npt) + C; + +% rhs = chnkr.data(1,:); rhs = rhs(:); +rhs = besselh(0,zk*abs((1+1.25*1i)-(chnkr.r(1,:)+1i*chnkr.r(2,:)))); rhs = rhs(:); +start = tic; sol = gmres(sys,rhs,[],1e-14,100); t1 = toc(start); + +fprintf('%5.2e s : time for dense gmres\n',t1) + +% evaluate at targets and plot + +rmin = min(chnkr); rmax = max(chnkr); +nplot = 300; +hx = (rmax(1)-rmin(1))/nplot; +hy = (rmax(2)-rmin(2))/nplot; +xtarg = linspace(rmin(1)+hx/2,rmax(1)-hx/2,nplot); +ytarg = linspace(rmin(2)+hy/2,rmax(2)-hy/2,nplot); +[xxtarg,yytarg] = meshgrid(xtarg,ytarg); +targets = zeros(2,length(xxtarg(:))); +targets(1,:) = xxtarg(:); targets(2,:) = yytarg(:); + +start = tic; in = chunkerinterior(chnkr,targets); t1 = toc(start); +fprintf('%5.2e s : time to find points in domain\n',t1) + +% compute layer potential at interior points +opts.forcepquad = true; +opts.side = 'i'; % 'i' for interior, 'e' for exterior, for positively oriented curve. +start = tic; +Csolpquad = chunkerkerneval(chnkr,fkern,sol,targets(:,in),opts); +% below check slp and dlp +% fkernd = kernel('helmholtz','d',zk); +% fkerns = kernel('helmholtz','s',zk); +% Ssolpquad = chunkerkerneval(chnkr,fkerns,sol,targets(:,in),opts); +% Dsolpquad = chunkerkerneval(chnkr,fkernd,sol,targets(:,in),opts); +% Csolpquad = Dsolpquad + (-zk*1i)*Ssolpquad; +t1 = toc(start); +fprintf('%5.2e s : time for kerneval (Helsing-Ojala for near)\n',t1); + +start = tic; +Csol = chunkerkerneval(chnkr,fkern,sol,targets(:,in)); t1 = toc(start); +fprintf('%5.2e s : time for kerneval (adaptive for near)\n',t1); + +% Compare with reference solution Dsol +error = max(abs(Csol-Csolpquad))/max(abs(Csol)); +fprintf('%5.2e : Relative max error\n',error); +% + +assert(error < 1e-10) +% profile viewer \ No newline at end of file