diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m index 0b535c4..401e0e4 100644 --- a/chunkie/+chnk/+helm2d/kern.m +++ b/chunkie/+chnk/+helm2d/kern.m @@ -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] @@ -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(:,:); @@ -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 diff --git a/chunkie/+chnk/+helm2dquas/green.m b/chunkie/+chnk/+helm2dquas/green.m index 781b5ff..2aaca02 100644 --- a/chunkie/+chnk/+helm2dquas/green.m +++ b/chunkie/+chnk/+helm2dquas/green.m @@ -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: @@ -13,7 +13,6 @@ [~,nsrc] = size(src); [~,ntarg] = size(targ); - xs = repmat(src(1,:),ntarg,1); ys = repmat(src(2,:),ntarg,1); @@ -40,8 +39,7 @@ npt = size(r,1); -ythresh = d/2; -% ythresh = d/10; +ythresh = 2*d/2; iclose = abs(ry) < ythresh; ifar = ~iclose; @@ -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)); @@ -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.']); @@ -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); @@ -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; @@ -202,7 +193,5 @@ if nargout>2 hess = reshape(quasi_phase.*hess,nkappa*ntarg,nsrc,3); end -% end - end diff --git a/chunkie/+chnk/+helm2dquas/kern.m b/chunkie/+chnk/+helm2dquas/kern.m index ecb4568..bb5f387 100644 --- a/chunkie/+chnk/+helm2dquas/kern.m +++ b/chunkie/+chnk/+helm2dquas/kern.m @@ -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] @@ -80,6 +80,7 @@ src = srcinfo.r(:,:); targ = targinfo.r(:,:); + kappa = quas_param.kappa; d = quas_param.d; l = quas_param.l; @@ -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); @@ -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); @@ -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(:,:); @@ -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(:,:); @@ -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 @@ -269,26 +267,27 @@ % 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 @@ -296,7 +295,7 @@ % 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(:,:); @@ -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 ... @@ -324,7 +325,7 @@ % 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(:,:); @@ -332,18 +333,19 @@ % 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 diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m index 3e3aa6c..0a644cb 100644 --- a/chunkie/@chunker/chunker.m +++ b/chunkie/@chunker/chunker.m @@ -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 @@ -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) @@ -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) diff --git a/chunkie/@chunker/tochunkgraph.m b/chunkie/@chunker/tochunkgraph.m new file mode 100644 index 0000000..ea49e94 --- /dev/null +++ b/chunkie/@chunker/tochunkgraph.m @@ -0,0 +1,50 @@ +function cgrph = tochunkgraph(chnkr) +% TOCHUNKGRAPH converts a chunker into a chunkgraph with one edge per +% connected component + +chnkr = sort(chnkr); +[~,~,info] = sortinfo(chnkr); + +ncomp = info.ncomp; +nchs = info.nchs; + +istart = 1; +nverts = 0; +verts = []; +edge2verts = []; +cparams = cell(1,ncomp); +fchnks = cell(1,ncomp); + +for i = 1:ncomp + nch = nchs(i); + vs = chunkends(chnkr,[istart,istart+nch-1]); + + ifclosed = info.ifclosed(i); + if ifclosed + vs = vs(:,1); + verts = [verts, vs]; + edge2verts = [edge2verts, [nverts+1;nverts+1]]; + nverts = nverts + 1; + else + vs = vs(:, [1,4]); + verts = [verts, vs]; + edge2verts = [edge2verts, [nverts+1;nverts+2]]; + nverts = nverts + 2; + end + + rchnkr = []; + rchnkr.r = chnkr.r(:,:,istart:(istart+nch-1)); + rchnkr.d = chnkr.d(:,:,istart:(istart+nch-1)); + rchnkr.d2 = chnkr.d2(:,:,istart:(istart+nch-1)); + fchnks{i} = chunkerpoints(rchnkr,struct('ifclosed',ifclosed)); + cparams{i}.ifclosed = ifclosed; + istart = istart + nch; +end + +pref = []; + +pref.dim = chnkr.dim; +pref.k = chnkr.k; +cgrph = chunkgraph(verts,edge2verts,fchnks,cparams,pref); + +end \ No newline at end of file diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m index fe5809b..d719bd1 100644 --- a/chunkie/@chunkgraph/chunkgraph.m +++ b/chunkie/@chunkgraph/chunkgraph.m @@ -67,7 +67,6 @@ % using bounding rectangles % flag = flagnear_rectangle_grid(obj,x,y) - find points defined via % a meshgrid of points that are close to the chunkgraph -% obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc % rmin = min(obj) - minimum of coordinate values % rmax = max(obj) - maximum of coordinate values % onesmat = onesmat(obj) - matrix that corresponds to integration of a @@ -94,14 +93,14 @@ % edgesendverts - (2 x nedges) array specifying starting and ending % vertices of an edge. % Optional input: -% fchnks - (nedges x 1) cell array of function handles, specifying a -% smooth curve to connect the given vertices. fchnk{j} should be a -% function in the format expected by CHUNKERFUNC but with the default -% that fchnk{j} is a function from [0,1] to the curve. If the specified -% curve does not actually connect the given vertices, the curve will be -% translated, rotated, and scaled to connect them. If the specified -% curve should be a loop, it is only translated to start at the correct -% vertex. +% fchnks - (1 x nedges) cell array of function handles or chunkers, +% specifying a smooth curve to connect the given vertices. If fchnk{j} +% is a function handle, it should be a function in the format expected +% by CHUNKERFUNC but with the default that fchnk{j} is a function from +% [0,1] to the curve. If the specified curve does not actually connect +% the given vertices, the curve will be translated, rotated, and scaled +% to connect them. If the specified curve should be a loop, it is only +% translated to start at the correct vertex. % cparams - struct or (nedges x 1) cell array of structs specifying curve % parameters in the format expected by CHUNKERFUNC. % pref - struct specifying CHUNKER preferences. @@ -179,7 +178,7 @@ fchnks = []; end - if isa(fchnks,"function_handle") + if isa(fchnks,"function_handle") || isa(fchnks,"chunker") fchnks0 = fchnks; fchnks = cell(nedge,1); for j = 1:nedge @@ -242,7 +241,7 @@ chnkr = sort(chnkr); %chnkr.vert = [v1,v2]; echnks(i) = chnkr; - elseif (~isempty(fchnks{i}) && isa(fchnks{i},'function_handle')) + else if iscell(cparams) cploc = cparams{i}; else @@ -250,46 +249,71 @@ end % set cploc.ifclosed in a way that makes sense cploc.ifclosed = false; + i1 = edgesendverts(1,i); i2 = edgesendverts(2,i); - if isnan(i1) || isnan(i2) cploc.ifclosed = true; end - % chunkgraph edges need at least 4 chunks - if isfield(cploc,'nchmin') - cploc.nchmin = max(4,cploc.nchmin); - else - cploc.nchmin = 4; - end - - ta = 0; tb = 1; - if isfield(cploc,'ta') - ta = cploc.ta; - else - cploc.ta = ta; - end - if isfield(cploc,'tb') - tb = cploc.tb; - else - cploc.tb = tb; - end - vs =fchnks{i}([ta,tb]); - if isfield(cploc,'maxchunklen') - if ~isnan(i1) && ~isnan(i2) - if i1 ~= i2 - vfin0 = verts(:,i1); - vfin1 = verts(:,i2); - scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro'); - cploc.maxchunklen = cploc.maxchunklen/scale; + if isa(fchnks{i},'function_handle') + + + % chunkgraph edges need at least 4 chunks + if isfield(cploc,'nchmin') + cploc.nchmin = max(4,cploc.nchmin); + else + cploc.nchmin = 4; + end + + ta = 0; tb = 1; + if isfield(cploc,'ta') + ta = cploc.ta; + else + cploc.ta = ta; + end + if isfield(cploc,'tb') + tb = cploc.tb; + else + cploc.tb = tb; + end + + vs =fchnks{i}([ta,tb]); + if isfield(cploc,'maxchunklen') + if ~isnan(i1) && ~isnan(i2) + if i1 ~= i2 + vfin0 = verts(:,i1); + vfin1 = verts(:,i2); + scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro'); + cploc.maxchunklen = cploc.maxchunklen/scale; + end end end + + chnkr = chunkerfunc(fchnks{i}, cploc, pref); + + elseif isa(fchnks{i},'chunker') + chnkr = fchnks{i}; + [~,~,info] = sortinfo(chnkr); + if info.ncomp > 1 + msg = "CHUNKGRAPH:CONSTRUCTOR: invalid edge descriptor." + ... + " fchnks{" + num2str(i) + "} must be a " + ... + "single connected component."; + warning(msg); + end + + vs = chunkends(chnkr,[1,chnkr.nch]); + vs = vs(:, [1,4]); + + else + msg = "CHUNKGRAPH:CONSTRUCTOR: invalid edge descriptor." + ... + " fchnks{" + num2str(i) + "} must be empty, a" + ... + " function handle, or a chunker."; + warning(msg); end - - chnkr = chunkerfunc(fchnks{i}, cploc, pref); chnkr = sort(chnkr); + if ~isnan(i1) && ~isnan(i2) if i1 ~= i2 vfin0 = verts(:,i1); @@ -317,7 +341,6 @@ chnkr = move(chnkr,zeros(size(vs(:,1))),vfin-vs(:,1),0,1); end end - echnks(i) = chnkr; end end diff --git a/chunkie/@kernel/helm2dquas.m b/chunkie/@kernel/helm2dquas.m index d71db42..a469177 100644 --- a/chunkie/@kernel/helm2dquas.m +++ b/chunkie/@kernel/helm2dquas.m @@ -60,7 +60,7 @@ l=2; N = 40; a = 15; M = 1e4; if nargin == 6 - if isefield(quad_opts,'l') + if isfield(quad_opts,'l') l = quad_opts.l; end if isfield(quad_opts,'N') @@ -83,11 +83,6 @@ quas_param.l = l; quas_param.sn = sn; -% obj.params.kappa = kappa; -% obj.params.d = d; -% obj.params.l = l; -% obj.params.Sn = Sn; - obj.params.quas_param = quas_param; switch lower(type) @@ -160,7 +155,6 @@ error('Unknown Helmholtz kernel type ''%s''.', type); end - end function f = helm2dquas_s_split(zk,s,t,quas_param) diff --git a/devtools/test/quasiperiodicTest.m b/devtools/test/quasiperiodicTest.m index e30fedd..30302aa 100644 --- a/devtools/test/quasiperiodicTest.m +++ b/devtools/test/quasiperiodicTest.m @@ -68,10 +68,10 @@ function quasiperiodicTest0() % test all all_assemb = zeros(size(allval)); -all_assemb(1:2:end, 1:2:end) = coefa(1,1)*dval; -all_assemb(1:2:end, 2:2:end) = coefa(1,2)*sval; -all_assemb(2:2:end, 1:2:end) = coefa(2,1)*dpval; -all_assemb(2:2:end, 2:2:end) = coefa(2,2)*spval; +all_assemb(1:2, 1:2:end) = coefa(1,1)*dval; +all_assemb(1:2, 2:2:end) = coefa(1,2)*sval; +all_assemb(3:end, 1:2:end) = coefa(2,1)*dpval; +all_assemb(3:end, 2:2:end) = coefa(2,2)*spval; assert(norm(all_assemb - allval) < 1e-12) % test transmission representiatoon @@ -84,8 +84,8 @@ function quasiperiodicTest0() % test c2trans c2trval_assemb = zeros(size(c2trval)); -c2trval_assemb(1:2:end, :) = coefs(1)*dval + coefs(2)*sval; -c2trval_assemb(2:2:end, :) = coefs(1)*dpval + coefs(2)*spval; +c2trval_assemb(1:2, :) = coefs(1)*dval + coefs(2)*sval; +c2trval_assemb(3:end, :) = coefs(1)*dpval + coefs(2)*spval; assert(norm(c2trval_assemb-c2trval) < 1e-12) end @@ -97,7 +97,7 @@ function quasiperiodicTest1() % problem parameters d= 8; zk = 1; -kappa = .05+0.1i; +kappa = .05+.1i; % setup geometry nch = 2^3; diff --git a/devtools/test/tochunkgraphTest.m b/devtools/test/tochunkgraphTest.m new file mode 100644 index 0000000..5cb31b7 --- /dev/null +++ b/devtools/test/tochunkgraphTest.m @@ -0,0 +1,62 @@ +tochunkgraphTest0(); + +function tochunkgraphTest0() +%TOCHUNKGRAPHTEST tests the routine for converting chunkers to chunkgraphs + +seed = 8675309; +rng(seed); + +% geometry parameters and construction + + +cparams = []; +cparams.eps = 1.0e-9; +pref = []; +pref.k = 16; +narms = 0; +amp = 0.0; +% make a circle +chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); + +% make an open arc +cparams.ifclosed = 0; +chnkr2 = chunkerfunc(@(t) cos_func(t,pi,1),cparams,pref); +chnkr2 = move(chnkr2,[0;3]); + +% merge arcs +chnkrs(3) = chunker(); +chnkrs(1) = chnkr; +rfac = 1.1; +chnkrs(2) = move(chnkr, [0;0], [3;0], 0, rfac); +chnkrs(3) = chnkr2; +chnkrtotal = merge(chnkrs); + +% form chunkgraph +cgrph = tochunkgraph(chnkrtotal); + +assert(size(cgrph.verts,2) == 4) +assert(length(cgrph.echnks) == 3) +assert(cgrph.npt == chnkrtotal.npt) +assert(norm(cgrph.echnks(1).r(:) - chnkr2.r(:))<1e-14) + +verts = [[0;0],[1;1],[3;0]]; edge2verts = [[1;2],[3;3]]; +cgrph = chunkgraph(verts,edge2verts,{chnkr2,chnkr}); +assert(isequal(cgrph.verts,verts)) +assert(isequal(cgrph.edgesendverts,edge2verts)) + +% verify chunker shifted and scaled correctly +vend = cgrph.echnks(1).r(:,1); +assert(norm(vend-verts(:,1)) < 1e-2); + +vend = cgrph.echnks(1).r(:,end); +assert(norm(vend-verts(:,2)) < 1e-2); + +end + +function [r,d,d2] = cos_func(t,d,A) +% parameterization of sinusoidal boundary with period d and amplitude A +omega = 2*pi/d; +r = [t(:), A*cos(omega*t(:))].'; +d = [ones(length(t),1), -omega*A*sin(omega*t(:))].'; +d2 = [zeros(length(t),1), -omega^2*A*cos(omega*t(:))].'; +end \ No newline at end of file