Docs

chunkie/+chnk/+helm2d/transmission_helper.m

@@ -84,7 +84,9 @@
     ncurve = length(chnkrs);
     sources = cell(1,ncurve);
     charges = cell(1,ncurve);
-    exposed_curves = ones(1,ncurve);
+
+    exposed_curves = (cs(1,:)==1) -(cs(2,:)==1);
+
     for i=1:ncurve
         sources{i} = zeros(2,1);
         charges{i} = 0 + 1j*0;

chunkie/@chunker/merge.m

@@ -28,7 +28,7 @@
 w = chnkrs(1).wstor;
 chnkrout = chunker(pref,t,w);
 
-for i = 1:length(chnkrs)
+for i = 1:numel(chnkrs)
   chnkrtemp = chnkrs(i);
   assert(chnkrtemp.dim == chnkrout.dim && chnkrtemp.k == chnkrout.k,...
       'chunkers to merge must be in same dimension and same order');

chunkie/@chunkgraph/chunkgraph.m

@@ -58,6 +58,11 @@
 %   obj = refine(obj,varargin) - refine the curve
 %   wts = weights(obj) - scaled integration weights on curve
 %   rn = normals(obj) - recompute normal vectors
+%   flag = flagnear(obj,pts) - find points close to the chunkgraph
+%   flag = flagnear_rectangle(obj,pts) - find points close to chunkgraph
+%           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
@@ -70,9 +75,8 @@
 %   kappa = signed_curvature(obj) - get signed curvature along curve
 %   obj = makedatarows(obj,nrows) - add nrows rows to the data storage.
 %   rflag = datares(obj,opts) - check if data in data rows is resolved
-%
-%   To add:
-%     flagnear
+%   edge_reg = find_regions_of_edges(obj) - determine regions on either 
+%                side of each edge
 %
 % Syntax:
 %
@@ -418,6 +422,7 @@
         obj = makedatarows(obj, nrows)
         scatter(obj, varargin)
         rres = datares(obj, opts)
+        edge_regs = find_edge_regions(obj)
 
 
     end

chunkie/@chunkgraph/find_edge_regions.m

@@ -0,0 +1,28 @@
+function edge_regs = find_edge_regions(obj)
+%FIND_EDGE_REGIONS find regions on either side of the edges
+% in a chunkgraph.
+%
+% Syntax: edge_regs = find_edge_regions(obj)
+%
+% Input:
+%   obj - chunkgraph object describing curve
+%   
+% Output:
+%   edge_regs - (2,nedge) array where nedges
+%     is the number of edges in the chunkgraph. 
+%     edge_reg(1,i) is index of the region in which the 
+%     normal to edge i is pointing and edge_reg(2,i)
+%     is the index corresponding to the negative of
+%     the normal.
+%
+% author: Tristan Goodwill
+
+
+    edge_regs = zeros(2,size(obj.edgesendverts,2));
+    nreg = length(obj.regions);
+    for i = 1:nreg
+        id_edge = obj.regions{i}{1};
+        edge_regs(1,id_edge(id_edge>0)) = i;
+        edge_regs(2,-id_edge(id_edge<0)) = i;
+    end
+end

chunkie/@kernel/elast2d.m

@@ -30,6 +30,8 @@
 %
 % See also CHNK.ELAST2D.KERN
 
+% author : Dan Fortunato, Travis Askham, Manas Rachh
+
 if ( nargin < 1 )
     error('Missing elasticity kernel type.');
 end

chunkie/@kernel/helm2d.m

@@ -20,6 +20,8 @@
 %   COEFS(2)*KERNEL.HELM2D('sp', ZK).
 % See also CHNK.HELM2D.KERN.
 
+% author: Dan Fortunato
+
 if ( nargin < 1 )
     error('Missing Helmholtz kernel type.');
 end

chunkie/@kernel/kernel.m

@@ -48,6 +48,15 @@
 %        single source and target, respectively. For scalar kernels,
 %        opdims = [1 1].
 %
+%      - K.sing: a string/character array specifying the singularity  
+%        strength of the kernel for sources and targets which are on
+%        the same curve. Currently recognized singularities are considered
+%        as part of a hierarchy
+%        - 'smooth' a smooth integral kernel
+%        - 'log' sum of above type kernel and phi(s)log(s-t)
+%        - 'pv' sum of above type kernels and phi(s)/(s-t)
+%        - 'hs' sum of above type kernels and phi(s)/(s-t)^2
+%
 %      - K.fmm: A function handle which calls the FMM for the corresponding
 %        kernel. K.fmm(eps, s, t, sigma) evaluates the kernel with density
 %        sigma from sources s to targets t with accuracy eps.

chunkie/@kernel/stok2d.m

@@ -55,6 +55,8 @@
 
 % See also CHNK.STOK2D.KERN.
 
+% author: Dan Fortunato  
+
 if ( nargin < 1 )
     error('Missing Stokes kernel type.');
 end

chunkie/demo/demo_Neumann_combined.m

@@ -18,7 +18,7 @@
 rado = 5;
 for i = 1:narms
     verts(:,2*i-1) = radi*[real(rots(i));imag(rots(i))];
-    verts(:,2*i  ) = rado*[real(rots(i));imag(rots(i))];
+    verts(:,2*i) = rado*[real(rots(i));imag(rots(i))];
 end
 nverts = size(verts,2);
 edge2verts = [1:nverts;circshift(1:nverts,-1)];
@@ -28,27 +28,18 @@
 frq = 5;
 
 % define curve for each edge
-fchnks = {};
+fchnks = cell(2*narms,1);
 for i = 1:narms
     % odd edges are straight
     fchnks{2*i-1} = [];
     % even edges are curved
-    fchnks{2*i  } = @(t) sinearc(t,amp,frq);
+    fchnks{2*i} = @(t) sinearc(t,amp,frq);
 end
 
 cparams = [];
 cparams.maxchunklen = min(4.0/zk,.5);
 [cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams);
 
-
-% plot
-figure(1);clf
-plot(cgrph)
-% hold on
-% quiver(cgrph)
-% hold off
-axis equal
-
 %% Define system
 
 % define kernels
@@ -74,6 +65,8 @@
    c3*Sik  Z        Z        ;
    c3*Sikp Z        Z        ];
 K = kernel(K);
+
+
 Keval = c1*kernel([Sk 1i*alpha*Dk Z]);
 
 npts = cgrph.npt;
@@ -132,6 +125,8 @@
 plot(cgrph,'k')
 axis equal
 title('$u^{\textrm{tot}}$','Interpreter','latex','FontSize',12)
+% END DEMO NEUMANN COMBINED
+saveas(figure(2),"demo_neumann_combined_plot.png")
 
 
 function [r,d,d2] = sinearc(t,amp,frq)

chunkie/demo/demo_concentric_transmission.m

@@ -21,27 +21,31 @@
 % geometry will be a series of loops through subsets of vertices
 id_vert_out = [5,1,6,2,7,3,8,4];
 % call helper function
-e2v_out = build_loop_edges(id_vert_out);
+e2v_out = [id_vert_out;circshift(id_vert_out,1)];
 
-id_vert_in = [5:8];
-e2v_in = build_loop_edges(id_vert_in);
+id_vert_in = 5:8;
+e2v_in = [id_vert_in;circshift(id_vert_in,1)];
 
 id_vert_in2 = [5,9,7,10];
-e2v_in2 =  build_loop_edges(id_vert_in2);
+e2v_in2 =  [id_vert_in2;circshift(id_vert_in2,1)];
 
 id_vert_in3 = [5,11,7,12];
-e2v_in3 =  build_loop_edges(id_vert_in3);
+e2v_in3 =  [id_vert_in3;circshift(id_vert_in3,1)];
 
 % combine lists of edges
 edge_2_verts = [e2v_out,e2v_in,e2v_in2,e2v_in3];
 
 cparams = [];
 cparams.maxchunklen = 4.0/max(zks);
 cgrph = chunkgraph(verts,edge_2_verts,[],cparams);
-figure(1);clf
-% plot(cgrph)
-% quiver(cgrph)
+figure(1);
+clf
 plot_regions(cgrph)
+xlim([-2,2])
+axis equal 
+% END DOMAIN
+
+saveas(figure(1),"demo_concentric_domain.png")
 
 %% Build system
 nreg = length(cgrph.regions);
@@ -52,18 +56,11 @@
 % jump condition coefficients
 coefs = ones(1,nreg);
 
-% identify region on each side of each edge
-edge_regs = zeros(2,size(edge_2_verts,2));
-for i = 1:nreg
-    id_edge = cgrph.regions{i}{1};
-    edge_regs(1,id_edge(id_edge>0)) = i;
-    edge_regs(2,-id_edge(id_edge<0)) = i;
-end
 
+edge_regs = find_edge_regions(cgrph);
 % set up parameters for planewave data
 opts = [];
 opts.bdry_data_type = 'pw';
-opts.exposed_curves = (edge_regs(1,:)==1) -(edge_regs(2,:)==1);
 
 % build system and get boundary data using the transmission helper
 [kerns, bdry_data, kerns_eval] = chnk.helm2d.transmission_helper(cgrph, ...
@@ -104,18 +101,15 @@
 utot = uin + uscat;
 
 umax = max(abs(utot(:))); 
-figure(2);clf
+figure(2); clf
 h = pcolor(xx,yy,imag(utot)); set(h,'EdgeColor','none'); colorbar
 colormap(redblue); clim([-umax,umax]);
 hold on 
 plot(cgrph,'k')
 axis equal
 title('$u^{\textrm{tot}}$','Interpreter','latex','FontSize',12)
 
+% END TRANSMISSION PROBLEM
+saveas(figure(2),"demo_concentric_fields.png")
 
 
-
-
-function edge_2_verts = build_loop_edges(id_verts)
-edge_2_verts = [id_verts;circshift(id_verts,1)];
-end

chunkie/demo/demo_many_scatterers.m

@@ -5,7 +5,6 @@
 % Demonstrate transmission problem with interleaved kernel
 % Demonstrate derived quantity
 
-
 % planewave direction
 phi = 0;
 % ambient wavenumber
@@ -64,30 +63,25 @@
 figure(3); clf;
 plot(chnkr)
 axis equal
-% hold on
-% quiver(chnkr)
-% hold off
 
 %% Make system
 % define system kernel
 
-
 coefs = ones(2,2,2);
-% we multiply normal derivative by negative one
-coefs(2,:,:) = -1;
-fkern = kernel('helmdiff','all',zks,coefs);
+coefs(:,2,:) = -1;
+fkern = kernel('helmdiff', 'all', zks, coefs);
 
 % define eval kernels
 fkern_eval(2,1) = kernel();
 for i=1:2
     Dk = kernel('helm', 'd', zks(i));
     Sk = kernel('helm', 's', zks(i));
-    fkern_eval(i) = kernel([Dk, Sk]);
+    fkern_eval(i) = kernel([Dk, (-1)*Sk]);
 end
 
 % define boundary data
 rhs_val = -planewave(kvec,chnkr.r(:,:));
-rhs_grad = -1i*sum(kvec.*chnkr.n(:,:),1).*rhs_val;
+rhs_grad = 1i*sum(kvec.*chnkr.n(:,:),1).*rhs_val;
 
 % interleave data
 rhs = zeros(2*chnkr.npt,1);

chunkie/demo/demo_mixed_bc.m

@@ -0,0 +1,77 @@
+%DEMO MIXED BOUNDARY VALUE PROBLEM
+%
+% This code illustrates the solution of the Laplace equation
+% with mixed boundary conditions. The domain is a star-shaped
+% domain where the radius has 5 Fourier modes, and the Dirichlet
+% part of the boundary corresponds to y>0, and the Neumann
+% part of the boundary corresponds to y<0. The boundary
+% data itself is given by a Fourier series in x
+% BEGIN DEMO MIXED BC
+%% Define geometry
+narms = 5;
+amp = 0.3;
+
+verts = [1+amp, -1+amp; 0, 0];
+fchnks = @(t) starfish(t, narms, amp);
+cparams = cell(2,1);
+cparams{1}.ta = 0;
+cparams{1}.tb = pi;
+cparams{1}.maxchunklen = 0.1;
+cparams{2}.ta = -pi;
+cparams{2}.tb = 0;
+cparams{2}.maxchunklen = 0.1;
+edgesendverts = [1 2; 2 1];
+
+cgrph = chunkgraph(verts, edgesendverts, fchnks, cparams);
+
+%% Setup kernels
+S = 2*kernel('lap', 's');
+D = (-2)*kernel('lap', 'd');
+Sp = 2*kernel('lap', 'sprime');
+Dp = (-2)*kernel('lap', 'dprime');
+
+K(2,2) = kernel();
+K(1,1) = D;
+K(1,2) = S;
+K(2,1) = Dp;
+K(2,2) = Sp;
+
+Keval(1,2) = kernel();
+Keval(1,1) = D;
+Keval(1,2) = S;
+
+%% Build matrix
+A = chunkermat(cgrph, K);
+A = A + eye(size(A,1));
+
+%% Build test boundary data;
+rhs = zeros(cgrph.npt,1);   
+npt1 = cgrph.echnks(1).npt;
+
+rng(1);
+rhs(1:npt1) = real(exp(1j*40*cgrph.echnks(1).r(1,:))*exp(1j*2*pi*rand)).';
+rhs(npt1+1:end) = real(exp(1j*43*cgrph.echnks(2).r(1,:))*exp(1j*2*pi*rand)).';
+
+sig = A \ rhs;
+
+%% Postprocess
+% define targets
+x1 = linspace(-1-amp,1+amp,300);
+[xx,yy] = meshgrid(x1,x1);
+targs = [xx(:).'; yy(:).'];
+
+% identify points in computational domain
+in = chunkerinterior(cgrph,{x1,x1});
+
+uu = chunkerkerneval(cgrph, Keval, sig, targs);
+u = nan(size(xx(:)));
+u(in) = uu(in);
+
+figure(1)
+clf
+plot(cgrph, 'LineWidth', 2);
+hold on;
+pcolor(xx, yy, reshape(u, size(xx))); shading interp
+axis equal tight
+% END DEMO MIXED BC
+saveas(figure(1),"mixed_bc_plot.png")