From 549e436b6472a09e3bfcbbe441efcaad00fa378b Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Wed, 28 Feb 2024 22:20:13 -0500 Subject: [PATCH] more stuff in fixing chunker mat, adding non-legendre points option to proxyfun, and derivative of c in kernel in helmholtz, fixing small typo of abs vs sqrt(z2) --- chunkie/+chnk/+flam/proxy_square_pts.m | 25 +++++++++-- chunkie/+chnk/+helm2d/fmm.m | 6 ++- chunkie/+chnk/+helm2d/kern.m | 8 ++-- chunkie/+chnk/+quadnative/buildmat.m | 7 +-- chunkie/@kernel/helm2d.m | 15 +++++++ chunkie/@kernel/kernel.m | 1 + chunkie/chunkerkernevalmat.m | 62 ++++++++++++++++++++++---- chunkie/chunkermat.m | 2 +- 8 files changed, 100 insertions(+), 26 deletions(-) diff --git a/chunkie/+chnk/+flam/proxy_square_pts.m b/chunkie/+chnk/+flam/proxy_square_pts.m index 1d7bb7f..3f1d362 100644 --- a/chunkie/+chnk/+flam/proxy_square_pts.m +++ b/chunkie/+chnk/+flam/proxy_square_pts.m @@ -1,4 +1,4 @@ -function [pr,ptau,pw,pin] = proxy_square_pts(porder) +function [pr,ptau,pw,pin] = proxy_square_pts(porder, opts) %PROXY_SQUARE_PTS return the proxy points, unit tangents, weights, and % function handle for determining if points are within proxy surface. % This function is for the proxy surface around a unit box centered at @@ -16,12 +16,31 @@ porder = 64; end +if nargin < 2 + opts = []; +end + +iflege = 0; +if isfield(opts, 'iflege') + iflege = opts.iflege; +end + + assert(mod(porder,4) == 0, ... 'number of proxy points on square should be multiple of 4') po4 = porder/4; -pts = -1.5+3*(0:(po4-1))*1.0/po4; +if ~iflege + pts = -1.5+3*(0:(po4-1))*1.0/po4; + pw = ones(porder,1)*(4.0*3.0/porder); +else + [xleg, wleg] = lege.exps(po4); + pts = xleg.'*1.5; + wleg = wleg*1.5; + pw = [wleg; wleg; wleg; wleg]; +end + one = ones(size(pts)); pr = [pts, one*1.5, -pts, -1.5*one; @@ -30,8 +49,6 @@ ptau = [one, one*0, -one, one*0; one*0, one, one*0, -one]; -pw = ones(porder,1)*(4.0*3.0/porder); - pin = @(x) max(abs(x),[],1) < 1.5; end \ No newline at end of file diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m index 4d50789..d7eceb0 100644 --- a/chunkie/+chnk/+helm2d/fmm.m +++ b/chunkie/+chnk/+helm2d/fmm.m @@ -38,6 +38,8 @@ % note that targinfo must contain targ.n % type = 'dprime' normal derivative of double layer, % note that targinfo must contain targ.n +% type = 'cprime' normal derivative of combined layer, +% note that targinfo must contain targ.n % sigma - density % varargin{1} - coef in the combined layer formula, otherwise % does nothing @@ -61,7 +63,7 @@ case {'d', 'dprime'} srcuse.dipstr = sigma(:).'; srcuse.dipvec = srcinfo.n(1:2,:); - case 'c' + case {'c', 'cprime'} coefs = varargin{1}; srcuse.charges = coefs(2)*sigma(:).'; srcuse.dipstr = coefs(1)*sigma(:).'; @@ -83,7 +85,7 @@ switch lower(type) case {'s', 'd', 'c'} varargout{1} = U.pottarg.'; - case {'sprime', 'dprime'} + case {'sprime', 'dprime', 'cprime'} if ( ~isfield(targinfo, 'n') ) error('CHUNKIE:helm2d:fmm:normals', ... 'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type); diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m index a90e636..d697efd 100644 --- a/chunkie/+chnk/+helm2d/kern.m +++ b/chunkie/+chnk/+helm2d/kern.m @@ -133,14 +133,14 @@ % normal derivative of combined field if strcmpi(type,'cprime') coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; + if(nargin == 5); coef = varargin{1}; end targnorm = targinfo.n; srcnorm = srcinfo.n; - submat = zeros(nt,ns); + % Get gradient and hessian info - [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + [~,grad,hess] = chnk.helm2d.green(zk,src,targ); nxsrc = repmat(srcnorm(1,:),nt,1); nysrc = repmat(srcnorm(2,:),nt,1); @@ -161,7 +161,7 @@ % Dirichlet and neumann data corresponding to combined field if strcmpi(type,'c2trans') coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; + if(nargin == 5); coef = varargin{1}; end targnorm = targinfo.n; srcnorm = srcinfo.n; diff --git a/chunkie/+chnk/+quadnative/buildmat.m b/chunkie/+chnk/+quadnative/buildmat.m index 262180b..7dbbd01 100644 --- a/chunkie/+chnk/+quadnative/buildmat.m +++ b/chunkie/+chnk/+quadnative/buildmat.m @@ -36,12 +36,7 @@ targinfo = []; targinfo.r = rt; targinfo.d = dt; targinfo.d2 = d2t; targinfo.n = nt; hs = h(j); ht = h(i); -dsnrms = sqrt(sum(abs(ds).^2,1)); -%taus = bsxfun(@rdivide,ds,dsnrms); - -%dtnrms = sqrt(sum(abs(dt).^2,1)); -%taut = bsxfun(@rdivide,dt,dtnrms); - +dsnrms = sqrt(sum(ds.^2,1)); ws = kron(hs(:),wts(:)); dsdt = dsnrms(:).*ws; diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m index 64bed00..73cc8ce 100644 --- a/chunkie/@kernel/helm2d.m +++ b/chunkie/@kernel/helm2d.m @@ -14,6 +14,10 @@ % parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('d', ZK) + % COEFS(2)*KERNEL.HELM2D('s', ZK). % +% KERNEL.HELM2D('cp', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS) +% constructs the derivative of the combined-layer Helmholtz kernel with +% wavenumber ZK and parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + +% COEFS(2)*KERNEL.HELM2D('sp', ZK). % See also CHNK.HELM2D.KERN. if ( nargin < 1 ) @@ -66,6 +70,17 @@ obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs); obj.sing = 'log'; + case {'cp', 'cprime'} + if ( nargin < 3 ) + warning('Missing combined layer coefficients. Defaulting to [1,1i].'); + coefs = [1,1i]; + end + obj.type = 'cprime'; + obj.params.coefs = coefs; + obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'cprime', coefs); + obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs); + obj.sing = 'hs'; + otherwise error('Unknown Helmholtz kernel type ''%s''.', type); diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index f824811..30803a6 100644 --- a/chunkie/@kernel/kernel.m +++ b/chunkie/@kernel/kernel.m @@ -9,6 +9,7 @@ % ---- ---- % 'laplace' ('lap', 'l') 's', 'd', 'sp', 'c' % 'helmholtz' ('helm', 'h') 's', 'd', 'sp', 'dp', 'c' +% 'cp' % 'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp' % 'elasticity' ('elast', 'e') 's', 'strac', 'd', 'dalt' % 'stokes' ('stok', 's') 'svel', 'spres', 'strac', diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m index 31b9526..f3efb12 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/chunkie/chunkerkernevalmat.m @@ -6,8 +6,17 @@ % Syntax: mat = chunkerkernevalmat(chnkr,kern,targs,opts) % % Input: -% chnkr - chunker object description of curve -% kern - integral kernel taking inputs kern(srcinfo,targinfo) +% chnkr - chunker object describing boundary, currently +% only supports chunkers, and not chunkgraphs +% 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,:) % targobj - object describing the target points, can be specified as % * array of points % * chunker object @@ -41,6 +50,37 @@ % author: Travis Askham (askhamwhat@gmail.com) +% convert kernel to kernel object, put in singularity info +% 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 + end + kern = kern2; + +elseif ~isa(kern,'kernel') + msg = "Second input is not a kernel object, function handle, " ... + + "or cell array"; + error(msg); +end + % determine operator dimensions using first two points @@ -51,7 +91,9 @@ targinfo.r = chnkr.r(:,2); targinfo.d = chnkr.d(:,2); targinfo.d2 = chnkr.d2(:,2); targinfo.n = chnkr.n(:,2); -ftemp = kern(srcinfo,targinfo); +ftmp = kern.eval; + +ftemp = ftmp(srcinfo,targinfo); opdims = size(ftemp); if nargin < 4 @@ -96,14 +138,16 @@ optsadap = []; optsadap.eps = eps; + + if forcesmooth - mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targinfo, ... + mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... [],optssmooth); return end if forceadap - mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ... + mat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ... targinfo,[],optsadap); return end @@ -111,9 +155,9 @@ if forcepqud optsflag = []; optsflag.fac = fac; flag = flagnear(chnkr,targinfo.r,optsflag); - spmat = chunkerkernevalmat_ho(chnkr,kern,opdims, ... + spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ... targinfo,flag,optsadap); - mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targinfo, ... + mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... flag,opts); mat = mat + spmat; return @@ -123,7 +167,7 @@ optsflag = []; optsflag.fac = fac; flag = flagnear(chnkr,targinfo.r,optsflag); -spmat = chunkerkernevalmat_adap(chnkr,kern,opdims, ... +spmat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ... targinfo,flag,optsadap); if nonsmoothonly @@ -131,7 +175,7 @@ return; end -mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targinfo, ... +mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... flag,opts); mat = mat + spmat; diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 5e1897f..f374acd 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -7,7 +7,7 @@ % Syntax: sysmat = chunkermat(chnkr,kern,opts) % % Input: -% chnkr - chunker object describing boundary +% 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.