Skip to content

Update shape of vector valued quasiperiodic kernels #143

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Update shape of vector valued quasiperiodic kernels #143

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a5dd01f
update shape of vector valued quasiperiodic kernels
tristangdwl Jun 30, 2025
458dd2e
update quasi green ls
tristangdwl Jun 30, 2025
db745cf
add tochunkgraph
tristangdwl Jul 11, 2025
8ee0e79
repair documentation
tristangdwl Jul 11, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions chunkie/+chnk/+helm2d/kern.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
% [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
% type == 'c2trans' returns the combined field, and the
% normal derivative of the combined field
% [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
% [coef(1,1)*D + coef(1,2)*S; coef(2,1)*D' + coef(2,2)*S']
% type == 'trans_rep' returns the potential corresponding
% to the transmission representation
% [coef(1)*D coef(2)*S]
Expand Down Expand Up @@ -250,6 +250,7 @@
if strcmpi(type,'c2trans')
coef = ones(2,1);
if(nargin == 5); coef = varargin{1}; end
if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
targnorm = targinfo.n(:,:);
srcnorm = srcinfo.n(:,:);

Expand All @@ -272,8 +273,8 @@

submat = zeros(2*nt,ns);

submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
submat(1:2:2*nt,:) = coef(1,1)*submatd + coef(1,2)*submats;
submat(2:2:2*nt,:) = coef(2,1)*submatdp + coef(2,2)*submatsp;

end

Expand Down
21 changes: 5 additions & 16 deletions chunkie/+chnk/+helm2dquas/green.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l)
%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic helmholtz green's function
%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic Helmholtz Green's function
% for the given sources and targets
%
% Input:
Expand All @@ -13,7 +13,6 @@
[~,nsrc] = size(src);
[~,ntarg] = size(targ);


xs = repmat(src(1,:),ntarg,1);
ys = repmat(src(2,:),ntarg,1);

Expand All @@ -40,8 +39,7 @@

npt = size(r,1);

ythresh = d/2;
% ythresh = d/10;
ythresh = 2*d/2;
iclose = abs(ry) < ythresh;
ifar = ~iclose;

Expand Down Expand Up @@ -76,23 +74,19 @@
% beta = sqrt((xi_m.^2-zk^2));
beta = sqrt(1i*(xi_m-zk)).*sqrt(-1i*(xi_m+zk));

% fhat = exp(-beta.*sqrt(ryfar.^2))./beta.*exp(1i*xi_m.*rxfar)/2;
fhat = exp(-beta.*sqrt(ryfar.^2) + 1i*xi_m.*rxfar)./(2*beta);
val(:,ifar,:) = sum(fhat,3)/(d);
if nargout > 1
grad(:,ifar,1) = sum(1i*xi_m.*fhat,3)/d;
grad(:,ifar,2) = sum(-beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
% grad(:,ifar,:) =[gx,gy];
end

if nargout >2
hess(:,ifar,1) = sum(-xi_m.^2.*fhat,3)/d;
hess(:,ifar,2) = sum(-1i*xi_m.*beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
hess(:,ifar,3) = sum((beta.*(sqrt(ryfar.^2)./ryfar)).^2.*fhat,3)/d;
% hess(:,ifar,:) = [hxx,hxy,hyy];
end


end

alpha = (exp(1i*kappa(:)*d));
Expand All @@ -101,7 +95,8 @@
grad_near = zeros(1,1,2);
hess_near = zeros(1,1,3);
if ~isempty(rxclose)
for i = -l:l
ls = -l:l;
for i = ls
rxi = rxclose - i*d;
if nargout>2
[vali,gradi,hessi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
Expand All @@ -123,8 +118,7 @@
val_near = val_near + vali.*alpha.^i;
end
end

% sn = sn.';

N = size(sn,2)-1;
ns = (0:N);
ns_use = (0:N+2);
Expand All @@ -149,9 +143,6 @@
val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
val(:,iclose) = val_near+val_far;




if nargout >1
DJs = cat(3,-Js(:,:,2),.5*(Js(:,:,1:end-3)-Js(:,:,3:end-1)))*zk;
ss = (eipn-1./eipn)/2i;
Expand Down Expand Up @@ -202,7 +193,5 @@
if nargout>2
hess = reshape(quasi_phase.*hess,nkappa*ntarg,nsrc,3);
end
% end

end

90 changes: 46 additions & 44 deletions chunkie/+chnk/+helm2dquas/kern.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@
% [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
% type == 'c2trans' returns the combined field, and the
% normal derivative of the combined field
% [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
% [coef(1,1)*D + coef(1,2)*S; coef(2,1)*D' + coef(2,2)*S']
% type == 'trans_rep' returns the potential corresponding
% to the transmission representation
% [coef(1)*D coef(2)*S]
Expand Down Expand Up @@ -80,6 +80,7 @@
src = srcinfo.r(:,:);
targ = targinfo.r(:,:);


kappa = quas_param.kappa;
d = quas_param.d;
l = quas_param.l;
Expand Down Expand Up @@ -111,13 +112,6 @@
submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
end

% x derivative of single layer
if strcmpi(type,'sx')
[~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);

submat = (grad(:,:,1));
end

% single later
if strcmpi(type,'s')
submat = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
Expand Down Expand Up @@ -145,21 +139,21 @@
% gradient of single layer
if strcmpi(type,'sgrad')
[~,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
submat = reshape(permute(grad,[3,1,2]),[],ns);
submat = reshape(permute(grad,[1,3,2]),[],ns);
end

% gradient of double layer
if strcmpi(type,'dgrad')
[~,~,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
submat = -(hess(:,:,1:2).*reshape(srcinfo.n(1,:),1,[],1)+hess(:,:,2:3).*reshape(srcinfo.n(2,:),1,[],1));
submat = reshape(permute(submat,[3,1,2]),[],ns);
submat = reshape(permute(submat,[1,3,2]),[],ns);
end

% Combined field
if strcmpi(type,'c')
srcnorm = srcinfo.n(:,:);
coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end
if(nargin >= 6); coef = varargin{1}; end
[submats,grad] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);
nx = repmat(srcnorm(1,:),nkappa*nt,1);
ny = repmat(srcnorm(2,:),nkappa*nt,1);
Expand All @@ -170,7 +164,7 @@
% normal derivative of combined field
if strcmpi(type,'cprime')
coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end
if(nargin >= 6); coef = varargin{1}; end
targnorm = targinfo.n(:,:);
srcnorm = srcinfo.n(:,:);

Expand All @@ -196,8 +190,9 @@

% Dirichlet and neumann data corresponding to combined field
if strcmpi(type,'c2trans')
coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end
coef = ones(1,2);
if(nargin >= 6); coef = varargin{1}; end
if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
targnorm = targinfo.n(:,:);
srcnorm = srcinfo.n(:,:);

Expand All @@ -220,23 +215,26 @@
% D
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);

submat = zeros(nkappa*2*nt,ns);

submat(1:2:end,:) = coef(1)*submatd + coef(2)*reshape(submats,[],ns);
submat(2:2:end,:) = coef(1)*submatdp + coef(2)*submatsp;
% submat = zeros(nkappa*2*nt,ns);
% submat(1:2:end,:) = coef(1)*submatd + coef(2)*reshape(submats,[],ns);
% submat(2:2:end,:) = coef(1)*submatdp + coef(2)*submatsp;

submat = zeros(nkappa,2,nt,ns);
submat(:,1,:,:) = reshape(coef(1,1)*submatd + coef(1,2)*reshape(submats,[],ns),nkappa,1,nt,[]);
submat(:,2,:,:) = reshape(coef(2,1)*submatdp + coef(2,2)*submatsp,nkappa,1,nt,[]);
submat = reshape(submat, [],ns);
end

% gradient of combined field
if strcmpi(type,'cgrad')
coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end
if(nargin >= 6); coef = varargin{1}; end
[~,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);

submats = reshape(permute(grad,[3,1,2]),[],ns);
submats = reshape(permute(grad,[1,3,2]),[],ns);

submatd = -(hess(:,:,1:2).*srcinfo.n(1,:)+hess(:,:,2:3).*srcinfo.n(2,:));
submatd = reshape(permute(submatd,[3,1,2]),[],ns);
submatd = reshape(permute(submatd,[1,3,2]),[],ns);

submat = coef(1)*submatd + coef(2)*submats;
end
Expand Down Expand Up @@ -269,34 +267,35 @@
% D
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);

submat = zeros(nkappa*2*nt,2*ns);
submat(1:2:end,1:2:2*ns) = submatd*cc(1,1);
submat(1:2:end,2:2:2*ns) = submats*cc(1,2);
submat(2:2:end,1:2:2*ns) = submatdp*cc(2,1);
submat(2:2:end,2:2:2*ns) = submatsp*cc(2,2);
submat = zeros(nkappa,2,nt,2*ns);
submat(:,1,:,1:2:2*ns) = reshape(submatd,nkappa,1,nt,[])*cc(1,1);
submat(:,1,:,2:2:2*ns) = reshape(submats,nkappa,1,nt,[])*cc(1,2);
submat(:,2,:,1:2:2*ns) = reshape(submatdp,nkappa,1,nt,[])*cc(2,1);
submat(:,2,:,2:2:2*ns) = reshape(submatsp,nkappa,1,nt,[])*cc(2,2);
submat = reshape(submat, [],2*ns);
end
% Dirichlet data/potential correpsonding to transmission rep
if strcmpi(type,'trans_rep')

coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end;
if(nargin >= 6); coef = varargin{1}; end;
srcnorm = srcinfo.n(:,:);
[submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);

submats = reshape(submats,[],ns);
nx = repmat(srcnorm(1,:),nt,1);
ny = repmat(srcnorm(2,:),nt,1);
submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);

submat = zeros(nkappa*nt,2*ns);
submat = zeros(nkappa*nt,2*ns,'like',1i);
submat(1:1:end,1:2:2*ns) = coef(1)*submatd;
submat(1:1:end,2:2:2*ns) = coef(2)*submats;
end

% Neumann data corresponding to transmission rep
if strcmpi(type,'trans_rep_prime')
coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end;
if(nargin >= 6); coef = varargin{1}; end;
targnorm = targinfo.n(:,:);
srcnorm = srcinfo.n(:,:);

Expand All @@ -305,10 +304,12 @@
% Get gradient and hessian info
[submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);

nxsrc = repmat(srcnorm(1,:),nt,1);
nysrc = repmat(srcnorm(2,:),nt,1);
nxtarg = repmat((targnorm(1,:)).',1,ns);
nytarg = repmat((targnorm(2,:)).',1,ns);
nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
nxtarg = repmat(reshape(targnorm(1,:),1,nt,1),nkappa,1,ns);
nxtarg = reshape(nxtarg,[],ns);
nytarg = repmat(reshape(targnorm(2,:),1,nt,1),nkappa,1,ns);
nytarg = reshape(nytarg,[],ns);

% D'
submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
Expand All @@ -324,26 +325,27 @@
% Gradient correpsonding to transmission rep
if strcmpi(type,'trans_rep_grad')
coef = ones(2,1);
if(nargin == 6); coef = varargin{1}; end;
if(nargin >= 6); coef = varargin{1}; end;

srcnorm = srcinfo.n(:,:);

% submat = zeros(nt,ns,6);
% S
[submats,grad,hess] = chnk.helm2dquas.green(src,targ,zk,kappa,d,sn,l);

nxsrc = repmat(srcnorm(1,:),nt,1);
nysrc = repmat(srcnorm(2,:),nt,1);
nxsrc = repmat(srcnorm(1,:),nkappa*nt,1);
nysrc = repmat(srcnorm(2,:),nkappa*nt,1);
% D
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);

submat = zeros(nkappa*2*nt,2*ns);
submat = zeros(nkappa,2*nt,2*ns);

submat(1:2:end,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc);
submat(1:2:end,2:2:2*ns) = coef(2)*grad(:,:,1);
submat(:,1,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc);
submat(:,1,2:2:2*ns) = coef(2)*grad(:,:,1);

submat(2:2:end,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
submat(2:2:end,2:2:2*ns) = coef(2)*grad(:,:,2);
submat(:,2,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
submat(:,2,2:2:2*ns) = coef(2)*grad(:,:,2);
submat = reshape(submat,nkappa*2*nt,2*ns);
end


Expand Down
3 changes: 2 additions & 1 deletion chunkie/@chunker/chunker.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@
% quiver(obj,varargin) - quiver plot of the chnkr points and normals
% scatter(obj,varargin) - scatter plot of the chnkr nodes
% tau = taus(obj) - unit tangents to curve
% cgrph = tochunkgraph(obj) - convert to chunkgraph
% obj = refine(obj,varargin) - refine the curve
% a = area(obj) - for a closed curve, area inside
% s = arclength(obj) - get values of arclength along curve
Expand Down Expand Up @@ -363,6 +364,7 @@
quiver(obj,varargin)
scatter(obj,varargin)
tau = taus(obj)
cgrph = tochunkgraph(obj)
obj = refine(obj,varargin)
a = area(obj)
s = arclength(obj)
Expand All @@ -377,7 +379,6 @@
obj = mtimes(A,obj)
obj = rotate(obj,theta,r0,r1)
obj = reflect(obj,theta,r0,r1)
du = arclengthder(obj,u)
end
methods(Static)
obj = chunkerfunc(fcurve,varargin)
Expand Down
Loading
Loading