diff --git a/chunkie/+chnk/+flam/proxyfunr.m b/chunkie/+chnk/+flam/proxyfunr.m index 482c89f..6d59a1a 100644 --- a/chunkie/+chnk/+flam/proxyfunr.m +++ b/chunkie/+chnk/+flam/proxyfunr.m @@ -102,8 +102,8 @@ srcinfo = []; srcinfo.r = pxy; srcinfo.d = ptau; if (size(rslf,1) == 2) srcinfo.n = chnk.perp(ptau); + targinfo.n = chnk.perp(dslf); end - Kpxy = kern(srcinfo,targinfo); Kpxy = Kpxy(islfuni2,:); diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m index d7eceb0..4a50233 100644 --- a/chunkie/+chnk/+helm2d/fmm.m +++ b/chunkie/+chnk/+helm2d/fmm.m @@ -78,6 +78,10 @@ pg = 0; pgt = min(nargout, 2); +switch lower(type) + case {'sprime', 'dprime', 'cprime'} + pgt = max(pgt, 2); +end U = hfmm2d(eps, zk, srcuse, pg, targuse, pgt); % Assign potentials diff --git a/chunkie/+chnk/chunkerkerneval_smooth.m b/chunkie/+chnk/chunkerkerneval_smooth.m new file mode 100644 index 0000000..90ea518 --- /dev/null +++ b/chunkie/+chnk/chunkerkerneval_smooth.m @@ -0,0 +1,218 @@ +function fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens, ... + targinfo,flag,opts) + +if isa(kern,'kernel') + kerneval = kern.eval; +else + kerneval = kern; +end + +flam = false; +accel = true; +forcefmm = false; + +if nargin < 6 + flag = []; +end +if nargin < 7 + opts = []; +end +if isfield(opts,'flam'); flam = opts.flam; end +if isfield(opts,'accel'); accel = opts.accel; end +if isfield(opts,'forcefmm'); forcefmm = opts.forcefmm; end + +k = chnkr.k; +nch = chnkr.nch; + +assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size') +dens = reshape(dens,opdims(2),k,nch); + +[~,w] = lege.exps(k); +[~,nt] = size(targinfo.r); + +fints = zeros(opdims(1)*nt,1); + +% assume smooth weights are good enough + +% Sequence of checks, first see ifflam is set as it supercedes +% everything, if not flam, then check to see if the fmm +% exists and whether it should be used +% The number of sources set to 200 is currently a hack, +% must be set based on opdims, accuracy, and kernel type +% considerations + +imethod = 'direct'; +if flam + imethod = 'flam'; +elseif isa(kern,'kernel') && ~isempty(kern.fmm) + if forcefmm + imethod = 'fmm'; + elseif accel + if nt > 200 || chnkr.npt > 200 + imethod = 'fmm'; + end + end +end + +diamsrc = max(abs(chnkr.r(:))); +diamtarg = max(abs(targinfo.r(:))); +diam = max(diamsrc, diamtarg); + +% flag targets that are within 1e-14 of sources +flagslf = chnk.flagself(targinfo.r, chnkr.r, 1e-14*diam); +if isempty(flagslf) + selfzero = sparse(opdims(1)*size(targinfo.r(:,:),2), ... + opdims(2)*chnkr.npt); +else + tmp = repmat((1:opdims(1))',opdims(2),1) + opdims(1)*(flagslf(1,:)-1); + flagslftarg = tmp(:); + tmp = repmat((1:opdims(2)),opdims(1),1); + tmp = tmp(:) + opdims(2)*(flagslf(2,:)-1); + flagslfsrc = tmp(:); + + selfzero = sparse(flagslftarg,flagslfsrc, 1e-300, ... + opdims(1)*size(targinfo.r(:,:),2), opdims(2)*chnkr.npt); +end + +if strcmpi(imethod,'direct') + % do dense version + if isempty(flag) + % nothing to ignore + for i = 1:nch + densvals = dens(:,:,i); densvals = densvals(:); + dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); + dsdtdt = dsdtdt(:).*w(:); + dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); + densvals = densvals.*(dsdtdt(:)); + srcinfo = []; srcinfo.r = chnkr.r(:,:,i); + srcinfo.n = chnkr.n(:,:,i); + srcinfo.d = chnkr.d(:,:,i); srcinfo.d2 = chnkr.d2(:,:,i); + kernmat = kerneval(srcinfo,targinfo); + + selfzeroch = selfzero(:, opdims(2)*k*(i-1) + (1:opdims(2)*k)); + [isp,jsp,~] = find(selfzeroch); + linsp = isp + (jsp-1)*size(selfzeroch,1); + kernmat(linsp) = 0; + + fints = fints + kernmat*densvals; + % sum(fints) + end + else + % ignore interactions in flag array + for i = 1:nch + densvals = dens(:,:,i); densvals = densvals(:); + dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); + dsdtdt = dsdtdt(:).*w(:); + dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); + densvals = densvals.*(dsdtdt(:)); + srcinfo = []; srcinfo.r = chnkr.r(:,:,i); + srcinfo.n = chnkr.n(:,:,i); + srcinfo.d = chnkr.d(:,:,i); srcinfo.d2 = chnkr.d2(:,:,i); + kernmat = kerneval(srcinfo,targinfo); + + rowkill = find(flag(:,i)); + rowkill = (opdims(1)*(rowkill(:)-1)).' + (1:opdims(1)).'; + kernmat(rowkill,:) = 0; + + selfzeroch = selfzero(:, opdims(2)*k*(i-1) + (1:opdims(2)*k)); + [isp,jsp,~] = find(selfzeroch); + linsp = isp + (jsp-1)*size(selfzeroch,1); + kernmat(linsp) = 0; + + fints = fints + kernmat*densvals; + end + end +else + + wts = chnkr.wts; + wts = wts(:); + + if strcmpi(imethod,'flam') + xflam1 = chnkr.r(:,:); + xflam1 = repmat(xflam1,opdims(2),1); + xflam1 = reshape(xflam1,chnkr.dim,numel(xflam1)/chnkr.dim); + + targinfo_flam = []; + targinfo_flam.r = repelem(targinfo.r(:,:),1,opdims(1)); + if isfield(targinfo, 'd') + targinfo_flam.d = repelem(targinfo.d(:,:),1,opdims(1)); + end + + if isfield(targinfo, 'd2') + targinfo_flam.d2 = repelem(targinfo.d2(:,:),1,opdims(1)); + end + + if isfield(targinfo, 'n') + targinfo_flam.n = repelem(targinfo.n(:,:),1,opdims(1)); + end + +% TODO: Pull through data? + + matfun = @(i,j) chnk.flam.kernbyindexr(i, j, targinfo_flam, ..., + chnkr, kerneval, opdims, selfzero); + + + width = max(abs(max(chnkr)-min(chnkr)))/3; + tmax = max(targinfo.r(:,:),[],2); tmin = min(targinfo.r(:,:),[],2); + wmax = max(abs(tmax-tmin)); + width = max(width,wmax/3); + npxy = chnk.flam.nproxy_square(kerneval,width); + [pr,ptau,pw,pin] = chnk.flam.proxy_square_pts(npxy); + + pxyfun = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ... + ctr,chnkr,kerneval,opdims,pr,ptau,pw,pin); + + optsifmm=[]; optsifmm.Tmax=Inf; + F = ifmm(matfun,targinfo_flam.r,xflam1,200,1e-14,pxyfun,optsifmm); + fints = ifmm_mv(F,dens(:),matfun); + else + wts2 = repmat(wts(:).', opdims(2), 1); + sigma = wts2(:).*dens(:); + fints = kern.fmm(1e-14, chnkr, targinfo.r(:,:), sigma); + end + % delete interactions in flag array (possibly unstable approach) + + + if ~isempty(flag) + for i = 1:nch + densvals = dens(:,:,i); densvals = densvals(:); + dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); + dsdtdt = dsdtdt(:).*w(:); + dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); + densvals = densvals.*(dsdtdt(:)); + srcinfo = []; srcinfo.r = chnkr.r(:,:,i); + srcinfo.n = chnkr.n(:,:,i); + srcinfo.d = chnkr.d(:,:,i); srcinfo.d2 = chnkr.d2(:,:,i); + + delsmooth = find(flag(:,i)); + delsmoothrow = (opdims(1)*(delsmooth(:)-1)).' + (1:opdims(1)).'; + delsmoothrow = delsmoothrow(:); + + targinfo_use = []; + targinfo_use.r = targinfo.r(:,delsmooth); + + if isfield(targinfo, 'd') + targinfo_use.d = targinfo.d(:,delsmooth); + end + + if isfield(targinfo, 'd2') + targinfo_use.d2 = targinfo.d2(:,delsmooth); + end + + if isfield(targinfo, 'n') + targinfo_use.n = targinfo.n(:,delsmooth); + end + + kernmat = kerneval(srcinfo,targinfo_use); + + selfzeroch = selfzero(:, opdims(2)*k*(i-1) + (1:opdims(2)*k)); + [isp,jsp,~] = find(selfzeroch); + linsp = isp + (jsp-1)*length(i(:)); + kernmat(linsp) = 0; + + fints(delsmoothrow) = fints(delsmoothrow) - kernmat*densvals; + end + end +end +% sum(fints) +end diff --git a/chunkie/+chnk/flagself.m b/chunkie/+chnk/flagself.m new file mode 100644 index 0000000..2878f8d --- /dev/null +++ b/chunkie/+chnk/flagself.m @@ -0,0 +1,58 @@ +function flagslf = flagself(srcs, targs, tol) + % identify sources and targets pairs that are within tol (1e-14); + + % todo: make this more robust by searching for all pairs that are close + % not just the closest + if nargin < 3 + tol = 1e-14; + end + + flagslf = []; + + randangle = 2*pi*rand(); + randrot = [cos(randangle), -sin(randangle); ... + sin(randangle), cos(randangle)]; + + srcrot = randrot*srcs(:,:); + targrot = randrot*targs(:,:); + + [srcsortx, jds] = sort(srcrot(1,:)); + [targsortx, ids] = sort(targrot(1,:)); + + binids = cell(length(srcsortx),1); + + idcheck = [1;2]; + for j = 1:length(srcsortx) + while (idcheck(2) < length(targsortx) && ... + targsortx(idcheck(2)) < srcsortx(j)) + idcheck = idcheck+1; + end + [d, kd] = min(abs(targsortx(idcheck) - srcsortx(j))); + idclose = idcheck(kd); + if d < tol + binids{jds(j)} = [binids{jds(j)}, ids(idclose)]; + end + end + + [srcsorty, jds] = sort(srcrot(2,:)); + [targsorty, ids] = sort(targrot(2,:)); + + idcheck = [1;2]; + for j = 1:length(srcsorty) + while (idcheck(2) < length(targsorty) && ... + targsorty(idcheck(2)) < srcsorty(j)) + idcheck = idcheck+1; + end + [d, kd] = min(abs(targsorty(idcheck) - srcsorty(j))); + idclose = idcheck(kd); + if d < tol + if ismember(ids(idclose), binids{jds(j)}) + flagslf = [flagslf, [jds(j);ids(idclose)]]; + end + end + end + if ~isempty(flagslf) + [~, ids] = sort(flagslf(1,:)); + flagslf = flagslf(:,ids); + end +end \ No newline at end of file diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index 46a2ed1..4906252 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -122,7 +122,7 @@ if (dim ~= 2); warning('only dimension two tested'); end if opts_use.forcesmooth - fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ... + fints = chnk.chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ... [],opts_use); return end @@ -137,7 +137,7 @@ if opts_use.forcepquad optsflag = []; optsflag.fac = opts_use.fac; flag = flagnear(chnkr,targinfo.r,optsflag); - fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ... + fints = chnk.chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ... flag,opts_use); fints = fints + chunkerkerneval_ho(chnkr,kern,opdims,dens, ... @@ -157,7 +157,7 @@ nlegnew = max(nlegnew,chnkr.k); [chnkr2,dens2] = upsample(chnkr,nlegnew,dens); -fints = chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targinfo, ... +fints = chnk.chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targinfo, ... flag,opts_use); fints = fints + chunkerkerneval_adap(chnkr,kern,opdims,dens, ... @@ -215,198 +215,6 @@ end -function fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens, ... - targinfo,flag,opts) - -if isa(kern,'kernel') - kerneval = kern.eval; -else - kerneval = kern; -end - -flam = false; -accel = true; -forcefmm = false; - -if nargin < 6 - flag = []; -end -if nargin < 7 - opts = []; -end -if isfield(opts,'flam'); flam = opts.flam; end -if isfield(opts,'accel'); accel = opts.accel; end -if isfield(opts,'forcefmm'); forcefmm = opts.forcefmm; end - -k = chnkr.k; -nch = chnkr.nch; - -assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size') -dens = reshape(dens,opdims(2),k,nch); - -[~,w] = lege.exps(k); -[~,nt] = size(targinfo.r); - -fints = zeros(opdims(1)*nt,1); - -% assume smooth weights are good enough - -% Sequence of checks, first see ifflam is set as it supercedes -% everything, if not flam, then check to see if the fmm -% exists and whether it should be used -% The number of sources set to 200 is currently a hack, -% must be set based on opdims, accuracy, and kernel type -% considerations - -imethod = 'direct'; -if flam - imethod = 'flam'; -elseif isa(kern,'kernel') && ~isempty(kern.fmm) - if forcefmm - imethod = 'fmm'; - elseif accel - if nt > 200 || chnkr.npt > 200 - imethod = 'fmm'; - end - end -end - -if strcmpi(imethod,'fmm') - icheck = exist(['fmm2d.' mexext], 'file'); - if icheck ~=3 - imethod = 'direct'; - fstr = ['CHUNKERKERNEVAL: forcefmm flag used but fmm2d not found\n' ... - 'Using direct computation instead\nMake sure fmm2d.' mexext ... - ' is in MATLAB path']; - warning(sprintf(fstr)) - end -end - -if strcmpi(imethod,'direct') - % do dense version - if isempty(flag) - % nothing to ignore - for i = 1:nch - densvals = dens(:,:,i); densvals = densvals(:); - dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:); - dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); - densvals = densvals.*(dsdtdt(:)); - srcinfo = []; srcinfo.r = chnkr.r(:,:,i); - srcinfo.n = chnkr.n(:,:,i); - srcinfo.d = chnkr.d(:,:,i); srcinfo.d2 = chnkr.d2(:,:,i); - kernmat = kerneval(srcinfo,targinfo); - fints = fints + kernmat*densvals; - end - else - % ignore interactions in flag array - for i = 1:nch - densvals = dens(:,:,i); densvals = densvals(:); - dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:); - dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); - densvals = densvals.*(dsdtdt(:)); - srcinfo = []; srcinfo.r = chnkr.r(:,:,i); - srcinfo.n = chnkr.n(:,:,i); - srcinfo.d = chnkr.d(:,:,i); srcinfo.d2 = chnkr.d2(:,:,i); - kernmat = kerneval(srcinfo,targinfo); - - rowkill = find(flag(:,i)); - rowkill = (opdims(1)*(rowkill(:)-1)).' + (1:opdims(1)).'; - kernmat(rowkill,:) = 0; - - fints = fints + kernmat*densvals; - end - end -else - - wts = chnkr.wts; - wts = wts(:); - - if strcmpi(imethod,'flam') - xflam1 = chnkr.r(:,:); - xflam1 = repmat(xflam1,opdims(2),1); - xflam1 = reshape(xflam1,chnkr.dim,numel(xflam1)/chnkr.dim); - - targinfo_flam = []; - targinfo_flam.r = repelem(targinfo.r(:,:),1,opdims(1)); - if isfield(targinfo, 'd') - targinfo_flam.d = repelem(targinfo.d(:,:),1,opdims(1)); - end - - if isfield(targinfo, 'd2') - targinfo_flam.d2 = repelem(targinfo.d2(:,:),1,opdims(1)); - end - - if isfield(targinfo, 'n') - targinfo_flam.n = repelem(targinfo.n(:,:),1,opdims(1)); - end - -% TODO: Pull through data? - - matfun = @(i,j) chnk.flam.kernbyindexr(i, j, targinfo_flam, ..., - chnkr, kerneval, opdims); - - - width = max(abs(max(chnkr)-min(chnkr)))/3; - tmax = max(targinfo.r(:,:),[],2); tmin = min(targinfo.r(:,:),[],2); - wmax = max(abs(tmax-tmin)); - width = max(width,wmax/3); - npxy = chnk.flam.nproxy_square(kerneval,width); - [pr,ptau,pw,pin] = chnk.flam.proxy_square_pts(npxy); - - pxyfun = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ... - ctr,chnkr,kerneval,opdims,pr,ptau,pw,pin); - - optsifmm=[]; optsifmm.Tmax=Inf; - F = ifmm(matfun,targinfo_flam.r,xflam1,200,1e-14,pxyfun,optsifmm); - fints = ifmm_mv(F,dens(:),matfun); - else - wts2 = repmat(wts(:).', opdims(2), 1); - sigma = wts2(:).*dens(:); - fints = kern.fmm(1e-14, chnkr, targinfo.r(:,:), sigma); - end - % delete interactions in flag array (possibly unstable approach) - - - if ~isempty(flag) - for i = 1:nch - densvals = dens(:,:,i); densvals = densvals(:); - dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:); - dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); - densvals = densvals.*(dsdtdt(:)); - srcinfo = []; srcinfo.r = chnkr.r(:,:,i); - srcinfo.n = chnkr.n(:,:,i); - srcinfo.d = chnkr.d(:,:,i); srcinfo.d2 = chnkr.d2(:,:,i); - - delsmooth = find(flag(:,i)); - delsmoothrow = (opdims(1)*(delsmooth(:)-1)).' + (1:opdims(1)).'; - delsmoothrow = delsmoothrow(:); - - targinfo_use = []; - targinfo_use.r = targinfo.r(:,delsmooth); - - if isfield(targinfo, 'd') - targinfo_use.d = targinfo.d(:,delsmooth); - end - - if isfield(targinfo, 'd2') - targinfo_use.d2 = targinfo.d2(:,delsmooth); - end - - if isfield(targinfo, 'n') - targinfo_use.n = targinfo.n(:,delsmooth); - end - - - kernmat = kerneval(srcinfo,targinfo_use); - fints(delsmoothrow) = fints(delsmoothrow) - kernmat*densvals; - end - end -end - -end function fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ... targinfo,flag,opts) diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 9aa6265..79f2fa5 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -39,6 +39,9 @@ % entries for which a special quadrature is used % (e.g. self and neighbor interactions) and return % in a sparse array. +% opts.corrections = boolean (false), if true, only compute the +% corrections to the smooth quadrature rule and +% return in a sparse array, see opts.nonsmoothonly % opts.l2scale = boolean (false), if true scale rows by % sqrt(whts) and columns by 1/sqrt(whts) % opts.auxquads = struct, struct storing auxilliary nodes @@ -98,9 +101,11 @@ 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) @@ -124,6 +129,7 @@ quad = 'ggq'; nonsmoothonly = false; +corrections = false; l2scale = false; isrcip = true; nsub = 40; @@ -144,6 +150,12 @@ if isfield(opts,'nonsmoothonly') nonsmoothonly = opts.nonsmoothonly; end +if isfield(opts,'corrections') + corrections = opts.corrections; +end +if corrections + nonsmoothonly = true; +end if(isfield(opts,'rcip')) isrcip = opts.rcip; @@ -162,7 +174,6 @@ adaptive_correction = opts.adaptive_correction; end - nchunkers = length(chnkrs); opdims_mat = zeros(2,nchunkers,nchunkers); @@ -199,11 +210,27 @@ irowlocs = zeros(nchunkers+1,1); icollocs = zeros(nchunkers+1,1); +idrowchnk = zeros(2,npttot); +idcolchnk = zeros(2,npttot); + + irowlocs(1) = 1; icollocs(1) = 1; for i=1:nchunkers icollocs(i+1) = icollocs(i) + lchunks(i)*opdims_mat(2,1,i); irowlocs(i+1) = irowlocs(i) + lchunks(i)*opdims_mat(1,i,1); + + % which chunker + idrowchnk(1,irowlocs(i):(irowlocs(i+1)-1)) = i; + % which chunk in the chunker + idrowchnk(2,irowlocs(i):(irowlocs(i+1)-1)) = ceil((1:lchunks(i)*opdims_mat(1,i,1))/ ... + (opdims_mat(1,i,1)*chnkrs(i).k)); + + % which chunker + idcolchnk(1,icollocs(i):(icollocs(i+1)-1)) = i; + % which chunk in the chunker + idcolchnk(2,icollocs(i):(icollocs(i+1)-1)) = ceil((1:lchunks(i)*opdims_mat(2,1,i))/ ... + (opdims_mat(2,1,i)*chnkrs(i).k)); end nrows = irowlocs(end)-1; @@ -470,6 +497,7 @@ chnk.rcip.setup(ngl,ndim,nedge,isstart); optsrcip = opts; optsrcip.nonsmoothonly = false; + optsrcip.corrections = false; R = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ... Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... @@ -504,6 +532,132 @@ end + +if (corrections) + % subtract out the smooth quadrature rule in the points computed using + % a special quadrature rule + sys_loc = sysmat; + iinds = isysmat(idx); + jinds = jsysmat(idx); + icorinds = []; + jcorinds = []; + vcors = []; + + iselfinds = []; + jselfinds = []; + + for ichkr = 1:nchunkers + for i = 1:chnkrs(ichkr).npt + targinfo = []; + targinfo.r = chnkrs(ichkr).r(:,i); + targinfo.d = chnkrs(ichkr).d(:,i); + targinfo.d2 = chnkrs(ichkr).d2(:,i); + targinfo.n = chnkrs(ichkr).n(:,i); + + icortmp = (i-1)*opdims_mat(1,ichkr,1) ... + +(1:opdims_mat(1,ichkr,1)) + irowlocs(ichkr) - 1; + jsrc = []; + for l = 1:opdims_mat(1,ichkr,1) + jsrc = [jsrc;jinds(iinds == icortmp(l))]; + end + + jchnks = unique(idcolchnk(:,jsrc)','rows')'; + jchnkrs = unique(jchnks(1,:)); + if (size(kern) == 1) + srcinfo = []; + srcinfo.r = []; + srcinfo.d = []; + srcinfo.d2 = []; + srcinfo.n = []; + srcw = []; + for jchkr = jchnkrs + jchnklgth = opdims_mat(2,1,jchkr)*chnkrs(jchkr).k; + + jcortmp = jchnklgth*(jchnks(2,jchnks(1,:)==jchkr)-1)... + + ((1:jchnklgth)')+icollocs(jchkr)-1; + jcortmp = jcortmp(:); % global density j indices + + icortmp2 = repmat(icortmp(:),1,length(jcortmp)); + icorinds = [icorinds; icortmp2(:)]; + + jcortmp2 = repmat(jcortmp,1,opdims_mat(1,ichkr,1))'; + jcorinds = [jcorinds; jcortmp2(:)]; + + jloc = chnkrs(jchkr).k*(jchnks(2,jchnks(1,:)==jchkr)-1) + ... + ((1:chnkrs(jchkr).k)'); % local pt j indices + jloc = jloc(:)'; + srcinfo.r = [srcinfo.r, chnkrs(jchkr).r(:,jloc)]; + srcinfo.d = [srcinfo.d, chnkrs(jchkr).d(:,jloc)]; + srcinfo.d2 = [srcinfo.d2,chnkrs(jchkr).d2(:,jloc)]; + srcinfo.n = [srcinfo.n, chnkrs(jchkr).n(:,jloc)]; + + srcwj = repmat(chnkrs(jchkr).wts(jloc), ... + opdims_mat(2,1,jchkr),1); + srcw = [srcw, srcwj(:)']; + end + vcorvals = kern.eval(srcinfo, targinfo).*srcw; + vcors = [vcors; vcorvals(:)]; + else + for jchkr = jchnkrs + jchnklgth = opdims_mat(2,1,jchkr)*chnkrs(jchkr).k; + + jcortmp = jchnklgth*(jchnks(2,jchnks(1,:)==jchkr)-1) + ... + ((1:jchnklgth)')+icollocs(jchkr)-1; + jcortmp = jcortmp(:); % global density j indices + + icortmp2 = repmat(icortmp(:),1,length(jcortmp)); + icorinds = [icorinds; icortmp2(:)]; + + jcortmp2 = repmat(jcortmp,1,opdims_mat(1,ichkr,1))'; + jcorinds = [jcorinds; jcortmp2(:)]; + + jloc = chnkrs(jchkr).k*(jchnks(2,jchnks(1,:)==jchkr)-1) + ... + ((1:chnkrs(jchkr).k)'); % local pt j indices + jloc = jloc(:)'; + + srcinfo = []; + srcinfo.r = chnkrs(jchkr).r(:,jloc); + srcinfo.d = chnkrs(jchkr).d(:,jloc); + srcinfo.d2 = chnkrs(jchkr).d2(:,jloc); + srcinfo.n = chnkrs(jchkr).n(:,jloc); + + srcw = repmat(chnkrs(jchkr).wts(jloc), ... + opdims_mat(2,1,jchkr),1); + srcw = srcw(:)'; + + vcorvals = kern(ichkr,jchkr).eval(srcinfo, targinfo) ... + .*srcw; + vcors = [vcors; vcorvals(:)]; + end + end + end + + % identify self interactions to be set to zero + ids = 1:lchunks(ichkr); + jds = ids; + + opdims = opdims_mat(:,ichkr,ichkr); + islf = repmat((1:opdims(1))',opdims(2),1) + opdims(1)*(ids-1) ... + + irowlocs(ichkr)-1; + tmp = repmat((1:opdims(2)),opdims(1),1); + jslf = tmp(:) + opdims(2)*(jds-1)+ icollocs(ichkr)-1; + + iselfinds = [iselfinds; islf(:)]; + jselfinds = [jselfinds; jslf(:)]; + + end + icorinds = icorinds(:); + jcorinds = jcorinds(:); + vcors = vcors(:); + + vcormat = sparse(icorinds,jcorinds,vcors,nrows,ncols); + + linsp = iselfinds + (jselfinds-1)*nrows; + vcormat(linsp) = 0; % set self interactions to zero to avoid NaNs + + sysmat = sys_loc-vcormat; +end + if (nargout >1) varargout{1} = opts; end diff --git a/chunkie/chunkermatapply.m b/chunkie/chunkermatapply.m new file mode 100644 index 0000000..8b91580 --- /dev/null +++ b/chunkie/chunkermatapply.m @@ -0,0 +1,243 @@ +function u = chunkermatapply(chnkr,kern,dval,dens,cormat,opts) +%CHUNKERMATAPPLY - apply chunkermat system on chunker defined by kern +% +% Syntax: u = chunkermatapply(chnkr,kern,dval,dens,sigma,opts) +% +% Input: +% chnkobj - chunker object describing boundary +% kern - kernel function. By default, this should be a function handle +% accepting input of the form kern(srcinfo,targinfo), where srcinfo +% and targinfo are in the ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.n - unit normals (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% dval - (default 0.0) float or float array. Let A be the matrix +% corresponding to on-curve convolution with the provided kernel. +% If a scalar is provided, the system matrix is +% A + dval*eye(size(A)) +% If a vector is provided, it should be length size(A,1). The +% system matrix is then +% A + diag(dval) +% dens - density on boundary, should have size opdims(2) x k x nch +% where k = chnkr.k, nch = chnkr.nch, where opdims is the +% size of kern for a single src,targ pair +% +% Optional input: +% cormat - sparse matrix of corrections to the smooth quadrature rule +% generated by chunkermat with option corrections = true +% It will be computed if not provided. +% opts - options structure. available options (default settings) +% opts.quad = string ('ggq'), specify quadrature routine to +% use for neighbor and self interactoins. Other +% available options include: +% - 'native' selects standard scaled Gauss-Legendre +% quadrature for smooth integral kernels +% opts.sing = string ('log') +% default type of singularity of kernel in case it +% is not defined by the kernel object. Supported +% types are: +% smooth => smooth kernels +% log => logarithmically singular kernels or +% smooth times log + smooth +% pv => principal value singular kernels + log +% hs => hypersingular kernels + pv +% +% opts.l2scale = boolean (false), if true scale rows by +% sqrt(whts) and columns by 1/sqrt(whts) +% opts.auxquads = struct, struct storing auxilliary nodes +% and weights which might be required for some of +% the quadrature methods like ggq for example. +% There is a different sub structure for each +% quadrature and singularity type which should be named +% as +% +% opts.auxquads. +% +% For example, the structure for logarithmically +% singular kernels integrated using ggq +% quadrature, the relevant struct is +% +% opts.auxquads.ggqlog +% +% The specific precomputed variables and their values +% will depend on the quadrature method used. +% opts.rcip = boolean (true), flag for whether to include rcip +% corrections for near corners if input chnkobj is +% of type chunkergraph +% opts.nsub_or_tol = (40) specify the level of refinements in rcip +% or a tolerance where the number of levels is given by +% ceiling(log_{2}(1/tol^2)); +% opts.adaptive_correction = (false) flag for whether to use +% adaptive quadrature for near touching panels on +% different chunkers +% opts.eps = (1e-14) tolerance for adaptive quadrature +% opts.flam - if = true, use flam utilities. to be replaced by the +% opts.forceflam flag. +% opts.flam supercedes opts.accel, if +% both are true, then flam will be used. (false) +% opts.accel - if = true, use specialized fmm if defined +% for the kernel, if it doesnt exist or if too few +% sources/targets, or if false, +% do direct. (true) +% opts.forcefmm - if = true, use specialized fmm if defined, +% independent of the number of sources/targets. (false) +% +% Outputs: +% u - system matrix applied to sigma +% +% +% +% See also: CHUNKERMAT + + +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 + 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 < 5 + cormat = []; +end +if nargin < 6 + opts = []; +end + + + +% get opts from struct if available + + +eps = 1e-14; +if isfield(opts,'eps') + eps = opts.eps; +end + +% Flag for determining whether input object is a chunkergraph +icgrph = 0; + +if (class(chnkr) == "chunker") + chnkrs = chnkr; +elseif(class(chnkr) == "chunkgraph") + icgrph = 1; + chnkrs = chnkr.echnks; +else + msg = "First input is not a chunker or chunkgraph object"; + error(msg) +end + + +% check for local corrections +if isempty(cormat) + selfopts = opts; + selfopts.corrections = true; + cormat = chunkermat(chnkr,kern,selfopts); +end + +% preproccessing +nchunkers = length(chnkrs); + +opdims_mat = zeros(2,nchunkers,nchunkers); +lchunks = zeros(nchunkers,1); + +%TODO: figure out a way to avoid this nchunkers^2 loop +for i=1:nchunkers + + targinfo = []; + targinfo.r = chnkrs(i).r(:,2); targinfo.d = chnkrs(i).d(:,2); + targinfo.d2 = chnkrs(i).d2(:,2); targinfo.n = chnkrs(i).n(:,2); + lchunks(i) = size(chnkrs(i).r(:,:),2); + + for j=1:nchunkers + + % determine operator dimensions using first two points + + srcinfo = []; + srcinfo.r = chnkrs(j).r(:,1); srcinfo.d = chnkrs(j).d(:,1); + srcinfo.d2 = chnkrs(j).d2(:,1); srcinfo.n = chnkrs(j).n(:,1); + + if (size(kern) == 1) + ftemp = kern.eval(srcinfo,targinfo); + else + ktmp = kern(i,j).eval; + ftemp = ktmp(srcinfo,targinfo); + end + opdims = size(ftemp); + opdims_mat(:,i,j) = opdims; + end +end + +irowlocs = zeros(nchunkers+1,1); +icollocs = zeros(nchunkers+1,1); + +irowlocs(1) = 1; +icollocs(1) = 1; +for i=1:nchunkers + icollocs(i+1) = icollocs(i) + lchunks(i)*opdims_mat(2,1,i); + irowlocs(i+1) = irowlocs(i) + lchunks(i)*opdims_mat(1,i,1); +end + + +% apply local corrections and diagonal scaling +u = dval*dens + cormat*dens; + +% apply smooth quadratures +if size(kern) == 1 + opts_smth = opts; + % opts_smth.flam = true; % todo fix flam on vector valued kernels + + chnkrmerge = merge(chnkrs); + + targinfo = []; + targinfo.r = chnkrmerge.r(:,:); + targinfo.d = chnkrmerge.d(:,:); + targinfo.d2 = chnkrmerge.d2(:,:); + targinfo.n = chnkrmerge.n(:,:); + + u = u + chnk.chunkerkerneval_smooth(chnkrmerge,kern, ... + opdims_mat(:,1,1),dens,targinfo,[],opts_smth); +else + opts_smth = opts; + % opts_smth.flam = true; % todo fix flam on vector valued kernels + + for i = 1:nchunkers + targinfo = []; + targinfo.r = chnkrs(i).r(:,:); + targinfo.d = chnkrs(i).d(:,:); + targinfo.d2 = chnkrs(i).d2(:,:); + targinfo.n = chnkrs(i).n(:,:); + iids = irowlocs(i):(irowlocs(i+1)-1); + for j = 1:nchunkers + jids = icollocs(j):(icollocs(j+1)-1); + u(iids) = u(iids) + chnk.chunkerkerneval_smooth(chnkrs(j), ... + kern(i,j),opdims_mat(:,i,j),dens(jids),targinfo,[], ... + opts_smth); + end + end +end + +end diff --git a/devtools/test/chunkerintegralTest.m b/devtools/test/chunkerintegralTest.m index deeed94..43efdfc 100644 --- a/devtools/test/chunkerintegralTest.m +++ b/devtools/test/chunkerintegralTest.m @@ -18,6 +18,7 @@ amp = 0.5; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams); + % scalar function on boundary fscal = @(xx) cos(xx(1,:)-1.0) + sin(xx(2,:)-0.5); @@ -36,5 +37,4 @@ assert(abs(fscal_int1-fscal_int2)/abs(fscal_int2) < 1e-9); assert(abs(fscal_int3-fscal_int2)/abs(fscal_int2) < 1e-9); -assert(abs(fscal_int4-fscal_int2)/abs(fscal_int2) < 1e-9); - +assert(abs(fscal_int4-fscal_int2)/abs(fscal_int2) < 1e-9); \ No newline at end of file diff --git a/devtools/test/chunkermatapplyTest.m b/devtools/test/chunkermatapplyTest.m new file mode 100644 index 0000000..c488754 --- /dev/null +++ b/devtools/test/chunkermatapplyTest.m @@ -0,0 +1,232 @@ +%CHUNKERMATAPPLYTEST test the routines for a matrix free apply of the system +%matrix +% +% + +clearvars; clear all; +addpaths_loc(); + +seed = 8675309; +rng(seed); + +% chunker geometry parameters and construction + +cparams = []; +cparams.eps = 1.0e-4; +pref = []; +pref.k = 16; +narms = 5; +amp = 0.5; +chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams); + +% scalar chunker test +fkern = kernel('lap','d'); + +start = tic; sysmat = chunkermat(chnkr,fkern); t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = -0.5*eye(chnkr.k*chnkr.nch) + sysmat; + +fkernsrc = kernel('lap','s'); +sources = [1;1]; +strengths = [1]; +% eval u on bdry + +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = chnkr.r(:,:); +targinfo.d = chnkr.d(:,:); +dens = fkernsrc.fmm(1e-12,srcinfo,targinfo,strengths); + +udense = sys*dens; +cormat = chunkermat(chnkr,fkern,struct("corrections",true)); +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) +relerr = norm(udense-u)/norm(udense); +fprintf('relative apply error %5.2e\n',relerr); +assert(relerr < 1e-13) + +sol1 = sys\dens; + +sol2 = gmres(sysapply, dens, [], 1e-14, 100); +relerr = norm(sol1-sol2)/norm(sol1); +fprintf('relative solve error %5.2e\n',relerr); +assert(relerr < 1e-13) + +% vector valued chunker test +nregions = 2; +ks = [1.1;2.1]*30; +coefs = [1.0;1.0]; +ncurve = 1; +cs(1,1:ncurve) = 1; +cs(2,1:ncurve) = 2; +opts = []; +opts.bdry_data_type = 'point sources'; + +sources = cell(1,2); +sources{1} = [0.0;0.1]; +sources{2} = [3.0;-3.2]; +charges{1} = (1.2+1j)*10; +charges{2} = (1+0j)*10; + +opts.sources = cell(1,2); +opts.sources{1} = sources{1}; +opts.sources{2} = sources{2}; +opts.charges = cell(1,2); +opts.charges{1} = charges{1}; +opts.charges{2} = charges{2}; +[kerns,bdry_data] = chnk.helm2d.transmission_helper(chnkr,ks,cs,coefs,opts); + +start = tic; sysmat = chunkermat(chnkr,kerns); +t1 = toc(start); + +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = eye(chnkr.k*chnkr.nch*2) + sysmat; + +udense = sys*bdry_data; +cormat = chunkermat(chnkr,kerns,struct("corrections",true)); +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) + +relerr = norm(udense-u)/norm(udense); +fprintf('relative apply error %5.2e\n',relerr); +assert(relerr < 1e-13) + +sol1 = sys\bdry_data; +sol2 = gmres(sysapply, bdry_data, [], 1e-12, 1000); + +relerr = norm(sol1-sol2)/norm(sol1); +fprintf('relative solve error %5.2e\n',relerr); +assert(relerr < 1e-10) + + +% setup chunkgraph + +nverts = 3; +verts = exp(1i*2*pi*(0:(nverts-1))/nverts); +verts = [real(verts);imag(verts)]; + +iind = 1:nverts; +jind = 1:nverts; + +iind = [iind iind]; +jind = [jind jind + 1]; +jind(jind>nverts) = 1; +svals = [-ones(1,nverts) ones(1,nverts)]; +edge2verts = sparse(iind,jind,svals,nverts,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; +[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams); + +vstruc = procverts(cgrph); +rgns = findregions(cgrph); +cgrph = balance(cgrph); + +% scalar chunkgraph test +fkern = -2*kernel('lap','d'); + +start = tic; sysmat = chunkermat(cgrph,fkern); t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = eye(cgrph.npt) + sysmat; + +fkernsrc = kernel('lap','s'); +sources = [1;1]; +strengths = [1]; +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = cgrph.r(:,:); +targinfo.d = cgrph.d(:,:); +dens = fkernsrc.fmm(1e-12,srcinfo,targinfo,strengths); + +udense = sys*dens; +cormat = chunkermat(cgrph,fkern,struct("corrections",true)); +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) + +relerr = norm(udense-u)/norm(udense); +fprintf('relative apply error %5.2e\n',relerr); +assert(relerr < 1e-13) + +sol1 = sys\dens; + +sol2 = gmres(sysapply, dens, [], 1e-14, 100); +relerr = norm(sol1-sol2)/norm(sol1); +fprintf('relative solve error %5.2e\n',relerr); +assert(relerr < 1e-13) + +% vectorvalued chunkgraph test +nregions = 2; +ks = [1.1;2.1]*30; +coefs = [1.0;1.0]; +ncurve = 3; +cs(1,1:ncurve) = 1; +cs(2,1:ncurve) = 2; +opts = []; +opts.bdry_data_type = 'point sources'; + +sources = cell(1,2); +sources{1} = [0.0;0.1]; +sources{2} = [3.0;-3.2]; +charges{1} = (1.2+1j)*10; +charges{2} = (1+0j)*10; + +opts.sources = cell(1,2); +opts.sources{1} = sources{1}; +opts.sources{2} = sources{2}; +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); + +start = tic; sysmat = chunkermat(cgrph,kerns); t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = eye(size(sysmat,1)) + sysmat; + + +udense = sys*bdry_data; +cormat = chunkermat(cgrph,kerns,struct("corrections",true)); +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) + +relerr = norm(udense-u)/norm(udense); +fprintf('relative apply error %5.2e\n',relerr); +assert(relerr < 1e-13) + +sol1 = sys\bdry_data; +sol2 = gmres(sysapply, bdry_data, [], 1e-12, 200); + +relerr = norm(sol1-sol2)/norm(sol1); +fprintf('relative solve error %5.2e\n',relerr); +assert(relerr < 1e-10) + + + + + +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/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/flagselfTest.m b/devtools/test/flagselfTest.m new file mode 100644 index 0000000..d3981b8 --- /dev/null +++ b/devtools/test/flagselfTest.m @@ -0,0 +1,29 @@ +% FLAGSELFTEST test the routine for flagging when sources and targets +% overlap + +clear all; close all; + +nsrc = 40000; +xtarg = linspace(0,1,floor(sqrt(nsrc))); +ytarg = linspace(0,2,floor(sqrt(nsrc))); +[xxsrc,yysrc] = meshgrid(xtarg,ytarg); +srcs = zeros(2,length(xxsrc(:))); +srcs(1,:) = xxsrc(:); srcs(2,:) = yysrc(:); + +nsrc = size(srcs,2); + + +targs = [srcs, [1;2].*rand(2,100)]; +P = randperm(size(targs,2)); +targs = targs(:,P)+1e-15*(2*rand(size(targs))-1); + +tic +flagslf = chnk.flagself(srcs, targs); +toc +fprintf('number of flagged points: %d\n', size(flagslf,2)) +assert(size(flagslf,2) == nsrc) + + +errs = sum(vecnorm(srcs(:,flagslf(1,:)) - targs(:,flagslf(2,:)))>1e-10); +fprintf('number of incorrectly flagged points: %d\n', errs) +assert(norm(srcs(:,flagslf(1,:)) - srcs(:,P(flagslf(2,:))),1)<1e-6) \ No newline at end of file diff --git a/devtools/test/mixedbcTest.m b/devtools/test/mixedbcTest.m new file mode 100644 index 0000000..89a7b3f --- /dev/null +++ b/devtools/test/mixedbcTest.m @@ -0,0 +1,262 @@ +% MIXEDBCTEST test the code with mixed boundary conditions, first test with +% mixed Dirichlet and Neumann BCs then test with Dirichlet and transmission +% conditions, which has variable opdims + +clearvars; close all; + +nverts = 4; +vertsout = exp(1i*2*pi*(0:(nverts-1))/nverts); +vertsout = [real(vertsout);imag(vertsout)]; + +vertsin = .5*exp(1i*2*pi*(0:(nverts-1))/nverts+1i*pi/4); +vertsin = [real(vertsin);imag(vertsin)]; + +verts = [vertsout, vertsin]; + +iind = 1:nverts; +jind = 1:nverts; + +iind = [iind iind]; +jind = [jind jind + 1]; +jind(jind>nverts) = 1; +svals = [-ones(1,nverts) ones(1,nverts)]; +edge2verts = sparse(iind,jind,svals,nverts,nverts); +edge2verts = [edge2verts, 0*edge2verts; 0*edge2verts, edge2verts]; + +edir = 1:4; % indices of edges with Dirichlet conditions +eneu = 5:8; % indices of edges with Neumann conditions + +fchnks = cell(1,size(edge2verts,1)); + +cparams = []; +cparams.nover = 2; +[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams); + +vstruc = procverts(cgrph); +rgns = findregions(cgrph); +cgrph = balance(cgrph); + +% dirichlet and neumann test +zk = 30; + +kerns(length([edir,eneu]),length([edir,eneu])) = kernel(); +kerns(edir,edir) = -2*kernel('helm', 'd', zk); % Dirichlet conditions +kerns(eneu,edir) = -2*kernel('helm', 'dp', zk); % Neumann conditions + +kerns(edir,eneu) = -2*kernel('helm', 's', zk); +kerns(eneu,eneu) = -2*kernel('helm', 'sp', zk); + + +start = tic; sysmat = chunkermat(cgrph,kerns); t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = eye(size(sysmat,1)) + sysmat; + +fkernsrc = kernel('helm','s',zk); +sources = [1;1]; +charges = [1]; +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = merge(cgrph.echnks(edir)).r(:,:); +targinfo.d = merge(cgrph.echnks(edir)).d(:,:); +bdrydatad = fkernsrc.fmm(1e-12,srcinfo,targinfo,charges); + +fkernsrc = kernel('helm','sp',zk); +targinfo = []; targinfo.r = merge(cgrph.echnks(eneu)).r(:,:); +targinfo.d = merge(cgrph.echnks(eneu)).d(:,:); +targinfo.n = merge(cgrph.echnks(eneu)).n(:,:); +bdrydatan = fkernsrc.fmm(1e-12,srcinfo,targinfo,charges); +bdrydata = [bdrydatad;bdrydatan]; + +sol1 = sys\bdrydata; + +cormat = chunkermat(cgrph,kerns,struct("corrections",true)); +sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat); + +xapply1 = sys*bdrydata; +xapply2 = sysapply(bdrydata); +relerr = norm(xapply1-xapply2)/norm(xapply1); +fprintf('relative matrix free apply error %5.2e\n',relerr); +assert(relerr < 1e-10) + + +rmin = min(cgrph.r(:,:)'); rmax = max(cgrph.r(:,:)'); +xl = rmax(1)-rmin(1); +yl = rmax(2)-rmin(2); +nplot = 100; +xtarg = linspace(rmin(1)-0.1*xl,rmax(1)+0.1*xl,nplot); +ytarg = linspace(rmin(2)-0.1*yl,rmax(2)+0.1*yl,nplot); +[xxtarg,yytarg] = meshgrid(xtarg,ytarg); +targets = zeros(2,length(xxtarg(:))); +targets(1,:) = xxtarg(:); targets(2,:) = yytarg(:); + +start = tic; in1 = chunkerinterior(merge(cgrph.echnks(edir)),{xtarg,ytarg}); +in2 = chunkerinterior(merge(cgrph.echnks(eneu)),{xtarg,ytarg}); +t1 = toc(start); +in = in1 & ~in2; +out = ~in; + +fprintf('%5.2e s : time to find points in domain\n',t1) + +% compute layer potential based on oversample boundary + +start = tic; +fkernd = 2*kernel('helm', 'd', zk); +iddir = 1:merge(cgrph.echnks(edir)).npt; +uscat = chunkerkerneval(merge(cgrph.echnks(edir)),fkernd,sol1(iddir), ... + targets(:,in)); + +fkern = 2*kernel('helm', 's', zk); +idneu = (1:merge(cgrph.echnks(eneu)).npt) + merge(cgrph.echnks(edir)).npt; +uscat = uscat + chunkerkerneval(merge(cgrph.echnks(eneu)),fkern, ... + sol1(idneu),targets(:,in)); +t1 = toc(start); + +fprintf('%5.2e s : time for kernel eval (for plotting)\n',t1) +fkernsrc = kernel('helm','s',zk); +targinfo = []; targinfo.r = targets(:,in); +uin = fkernsrc.fmm(1e-12,srcinfo,targinfo,charges); +utot = uscat(:)+uin(:); + +figure(1) +t = tiledlayout(1,2,'TileSpacing','compact'); + +nexttile +zztarg = nan(size(xxtarg)); +zztarg(in) = uscat; +h=pcolor(xxtarg,yytarg,imag(zztarg)); +set(h,'EdgeColor','none') +hold on +plot(cgrph,'k','LineWidth',2) +axis equal tight +set(gca, "box","off","Xtick",[],"Ytick",[]); +title('$u$','Interpreter','latex','FontSize',12) + +nexttile +zztarg = nan(size(xxtarg)); +zztarg(in) = utot; +h=pcolor(xxtarg,yytarg,log10(abs((zztarg)))); +set(h,'EdgeColor','none') +hold on +plot(cgrph,'k','LineWidth',2) +axis equal tight +set(gca, "box","off","Xtick",[],"Ytick",[]); +colorbar() + +title('$\log10($error$)$','Interpreter','latex','FontSize',12) +relerr = max(abs(utot)); +fprintf('relative field error %5.2e\n',relerr); +assert(relerr < 1e-4) + +% dirichlet and transmission test +nregions = 2; +ks = [1.1;2.1]*30; + + +kerns(length([edir,eneu]),length([edir,eneu])) = kernel(); +kerns(edir,edir) = -2*kernel('helm', 'd', ks(1)); + +kerntmp = @(s,t) -2*[chnk.helm2d.kern(ks(1),s,t,'d'); ... + -chnk.helm2d.kern(ks(1),s,t,'dprime')]; +kerns(eneu,edir) = (kerntmp); + +kerntmp = @(s,t) -chnk.helm2d.kern(ks(1),s,t,'trans_rep',[-1,-1]); +kerns(edir,eneu) = (kerntmp); + +cc1 = [-1,-1;1,1]; +cc2 = cc1; +kerntmp = @(s,t) -(chnk.helm2d.kern(ks(1),s,t,'all',cc1)- ... + chnk.helm2d.kern(ks(2),s,t,'all',cc2)); +kerns(eneu,eneu) = (kerntmp); + + +start = tic; sysmat = chunkermat(cgrph,kerns); t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sys = eye(size(sysmat,1)) + sysmat; + +fkernsrc = kernel('helm','s',ks(1)); +sources = [.1;.1]; +charges = [1]; +srcinfo = []; srcinfo.r = sources; +targinfo = []; targinfo.r = merge(cgrph.echnks(edir)).r(:,:); +targinfo.d = merge(cgrph.echnks(edir)).d(:,:); +bdrydatad = fkernsrc.fmm(1e-12,srcinfo,targinfo,charges); + +srcinfo.sources = sources; +srcinfo.charges = charges; +targinfo = []; +targinfo.sources = merge(cgrph.echnks(eneu)).r(:,:); +targinfo.n = merge(cgrph.echnks(eneu)).n(:,:); +U = hfmm2d(1e-12,ks(1),srcinfo,0,targinfo.sources,2); +bdrydatat = [U.pottarg; -sum(targinfo.n.*U.gradtarg,1)]; +bdrydata = [bdrydatad;bdrydatat(:)]; + +sol1 = sys\bdrydata; + +cormat = chunkermat(cgrph,kerns,struct("corrections",true)); +sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat); + +xapply1 = sys*bdrydata; +xapply2 = sysapply(bdrydata); +relerr = norm(xapply1-xapply2)/norm(xapply1); +fprintf('relative matrix free apply error %5.2e\n',relerr); +assert(relerr < 1e-10) + +start = tic; +fkernd = 2*kernel('helm', 'd', ks(1)); +iddir = 1:merge(cgrph.echnks(edir)).npt; +uscat1 = chunkerkerneval(merge(cgrph.echnks(edir)),fkernd,sol1(iddir), ... + targets(:,in)); + +kerntmp = @(s,t) -chnk.helm2d.kern(ks(1),s,t,'trans_rep',[1,1]); +fkern = kernel(kerntmp); +idneu = (1:2*merge(cgrph.echnks(eneu)).npt) + merge(cgrph.echnks(edir)).npt; +uscat1 = uscat1 - chunkerkerneval(merge(cgrph.echnks(eneu)),fkern, ... + sol1(idneu),targets(:,in)); + +kerntmp = @(s,t) -chnk.helm2d.kern(ks(2),s,t,'trans_rep',[1,1]); +fkern = kernel(kerntmp); +uscat2 = chunkerkerneval(merge(cgrph.echnks(eneu)),fkern,sol1(idneu), ... + targets(:,in2)); +t1 = toc(start); +fprintf('%5.2e s : time for kernel eval (for plotting)\n',t1) + +fkernsrc = kernel('helm','s',ks(1)); +targinfo = []; targinfo.r = targets(:,in); +uin1 = fkernsrc.fmm(1e-12,srcinfo,targinfo,charges); +utot1 = uscat1(:)-uin1(:); +utot2 = uscat2(:); + +figure(3) +t = tiledlayout(1,2,'TileSpacing','compact'); + +nexttile +zztarg = nan(size(xxtarg)); +zztarg(in) = uscat1; +zztarg(in2) = uscat2; +h=pcolor(xxtarg,yytarg,imag(zztarg)); +set(h,'EdgeColor','none') +hold on +plot(cgrph,'k','LineWidth',2) +axis equal tight +set(gca, "box","off","Xtick",[],"Ytick",[]); +title('$u$','Interpreter','latex','FontSize',12) + +nexttile +zztarg = nan(size(xxtarg)); +zztarg(in) = utot1; +zztarg(in2) = uscat2; +h=pcolor(xxtarg,yytarg,log10(abs((zztarg)))); +set(h,'EdgeColor','none') +hold on +plot(cgrph,'k','LineWidth',2) +axis equal tight +set(gca, "box","off","Xtick",[],"Ytick",[]); +colorbar() + +title('$\log10($error$)$','Interpreter','latex','FontSize',12) +relerr1 = max(abs(utot1)); +relerr2 = max(abs(uscat2)); +relerr = max(relerr1, relerr2); +fprintf('relative field error %5.2e\n',relerr); +assert(relerr < 1e-4)