Added demos

chunkie/demo/demo_concentric_transmission.m

@@ -0,0 +1,121 @@
+%DEMO_CONCENTRIC_TRANSMISSION
+%
+% demonstrate Helmholtz transmission helper
+% demonstrate chunkgraph regions
+%
+% code takes about 2 minutes
+
+
+% wavenumbers
+zks = 2.^(1+(1:5));
+
+%% Build geometry
+% define verts
+verts_out = [[1;1],[1;-1],[-1;-1],[-1;1]];
+verts_in = [[0;1],[1;0],[0;-1],[-1;0]];
+verts_in2 = [[1/2;0],[-1/2;0]];
+verts_in3 = verts_in2/2;
+
+verts = [verts_out, verts_in, verts_in2, verts_in3];
+
+% 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);
+
+id_vert_in = [5:8];
+e2v_in = build_loop_edges(id_vert_in);
+
+id_vert_in2 = [5,9,7,10];
+e2v_in2 =  build_loop_edges(id_vert_in2);
+
+id_vert_in3 = [5,11,7,12];
+e2v_in3 =  build_loop_edges(id_vert_in3);
+
+% 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)
+plot_regions(cgrph)
+
+%% Build system
+nreg = length(cgrph.regions);
+
+% assign region wavenumbers
+ks = [zks(1),repmat(zks(2),1,4),repmat(zks(3),1,2),repmat(zks(4),1,2),zks(5)];
+
+% 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
+
+% 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, ...
+                                   ks, edge_regs, coefs, opts);
+% build system matrix
+tic;
+[sysmat] = chunkermat(cgrph, kerns);
+sysmat = sysmat + eye(size(sysmat,2));
+tbuild = toc
+
+% solve
+tic; 
+dens = sysmat\bdry_data;
+tsolve = toc
+
+%% Compute field
+
+L = 1.5*max(vecnorm(cgrph.r(:,:)));
+x1 = linspace(-L,L,300);
+[xx,yy] = meshgrid(x1,x1);
+targs = [xx(:).'; yy(:).'];
+ntargs = size(targs,2);
+
+% evaluate potentials and return region labels
+tic;
+[uscat, targdomain] = chunkerkerneval(cgrph, kerns_eval, dens, targs);
+uscat = reshape(uscat,size(xx));
+tplot = toc
+
+% identify points in computational domain
+in = chunkerinterior(cgrph,{x1,x1});
+out = ~in;
+
+% get incoming solution
+uin = zeros(size(xx));
+uin(targdomain==1) = planewave([zks(1);0],targs(:,targdomain==1));
+
+utot = uin + uscat;
+
+umax = max(abs(utot(:))); 
+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)
+
+
+
+
+
+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

@@ -0,0 +1,182 @@
+%DEMO_MANY_SCATTERERS
+%
+% Solve transmission scattering problems with many inclusions
+%
+% Demonstrate transmission problem with interleaved kernel
+% Demonstrate derived quantity
+
+
+% planewave direction
+phi = 0;
+% ambient wavenumber
+zk1 = 20;
+kvec = zk1*[cos(phi);sin(phi)];
+% inclusion wavenumber
+zk2 = 1.5*zk1;
+zks = [zk1, zk2];
+
+%% Make geometry
+
+% geometry preferences
+cparams = [];
+cparams.maxchunklength = min(4.0/max(zks),0.125);
+
+pref = []; 
+pref.k = 16;
+narms =5;
+amp = 0.25;
+scale = .6;
+rad = scale*(amp+1);
+ntry = 1000;
+
+% make interior boundaries with random locations
+chnkr = [];
+L = 3;
+theta = 2*pi*rand();
+ctrs = L*rand()*[cos(theta);sin(theta)];
+n_pts = [];
+nscat = 3;
+for i  = 1:nscat
+    % each boundary is rotated starfish with a random number of arms
+    phi = 2*pi*rand();
+    narms = randi([3,6]);
+    chnkr_i = chunkerfunc(@(t) starfish(t,narms,amp,ctrs(:,i),phi,scale), ...
+        cparams,pref); 
+    % track number of points
+    n_pts(i) = chnkr_i.npt;
+    chnkr = [chnkr,chnkr_i];
+
+    % try to find location for next chunker
+    for j = 1:ntry
+        theta = 2*pi*rand();
+        tmp = L*rand()*[cos(theta);sin(theta)];
+        rmin = min(vecnorm(tmp - ctrs));
+        if (rmin > rad * 3); break; end
+    end
+    if j == ntry; error('Could not place next boundary'); end
+    ctrs = [ctrs,tmp];
+end
+ctrs = ctrs(:,1:end-1);
+chnkr = merge(chnkr);
+npt_int = chnkr.npt;
+
+fprintf('Geometry generated\n')
+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);
+
+% 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]);
+end
+
+% define boundary data
+rhs_val = -planewave(kvec,chnkr.r(:,:));
+rhs_grad = -1i*sum(kvec.*chnkr.n(:,:),1).*rhs_val;
+
+% interleave data
+rhs = zeros(2*chnkr.npt,1);
+rhs(1:2:end) = rhs_val;
+rhs(2:2:end) = rhs_grad;
+
+tic;
+% build fast direct solver
+F = chunkerflam(chnkr,fkern,1.0);
+t_build_solver = toc
+tic;
+% solve
+sol = rskelf_sv(F,rhs);
+tsolve = toc
+
+
+%% Compute example field
+
+L = 1.1*max(vecnorm(chnkr.r(:,:)));
+x1 = linspace(-L,L,500);
+[xx,yy] = meshgrid(x1,x1);
+targs = [xx(:).'; yy(:).'];
+ntargs = size(targs,2);
+
+% identify points in computational domain
+in = chunkerinterior(chnkr,{x1,x1});
+out = ~in;
+
+% get incoming solution
+uin = zeros(size(xx));
+uin(out) = planewave(kvec(:),targs(:,out));
+
+% get solution
+tic;
+uscat = zeros(size(xx));
+uscat(out) = chunkerkerneval(chnkr,fkern_eval(1),sol,targs(:,out));
+uscat(in) = chunkerkerneval(chnkr,fkern_eval(2),sol,targs(:,in));
+tplot = toc
+
+utot = uin + uscat;
+%% make plots
+umax = max(abs(utot(:))); 
+figure(2);clf
+h = pcolor(xx,yy,imag(utot)); set(h,'EdgeColor','none'); colorbar
+colormap(redblue); clim([-umax,umax]);
+hold on 
+plot(chnkr,'k')
+axis equal
+title('$u^{\textrm{tot}}$','Interpreter','latex','FontSize',12)
+
+%% Compute monostatic radar cross section
+
+nphi = 500;
+phis = 2*pi*(1:nphi)/nphi;
+
+tic;
+cross = zeros(nphi,1);
+for i = 1:nphi
+    phi = phis(i);
+    kvec = zk1*[cos(phi);sin(phi)];
+
+    % define boundary data
+    rhs_val = -planewave(kvec,chnkr.r(:,:));
+    rhs_grad = -1i*sum(kvec.*chnkr.n(:,:),1).*rhs_val;
+
+    % interleave data
+    rhs = zeros(2*chnkr.npt,1);
+    rhs(1:2:end) = rhs_val;
+    rhs(2:2:end) = rhs_grad;
+
+    % solve
+    sol = rskelf_sv(F,rhs);
+
+    % extract individual densities
+    ddens = sol(1:2:end);
+    sdens = sol(2:2:end);
+
+    % compute intermediate quantities
+    d_dot_n = sum(kvec/zk1.*chnkr.n(:,:),1).';
+    d_fac = exp(-1i*pi/2)*zk1;
+
+    %compute cross section as integral against density times farfield
+    %signature of kernels
+    cross(i) = chnkr.wts(:).' * (sdens + ddens.*d_fac.*d_dot_n);
+    cross(i) = cross(i)/4i * sqrt(zk1) * exp(-1i*pi/4) * sqrt(2/pi);
+end
+tcross = toc
+
+figure(3);clf
+plot(phis,abs(cross),'-','Linewidth',2)
+xlabel('$\phi$','Interpreter','latex')
+ylabel('abs of Monostatic cross section','Interpreter','latex')

chunkie/demo/demo_mult_connect.m

@@ -140,10 +140,10 @@
 % identify points in computational domain
 in = chunkerinterior(chnkr,{x1,x1});
 
-% evaluate electric field
+% evaluate electric field, the negative gradient of the potential
 uu = nan(size(xx));
 uugrad = nan(2,size(xx(:),1));
-uugrad(:,in(:)) = reshape(chunkerkerneval(chnkr,ckern_grad,sol,targs(:,in(:))),2,[]);
+uugrad(:,in(:)) = -reshape(chunkerkerneval(chnkr,ckern_grad,sol,targs(:,in(:))),2,[]);
 tstream = toc(tstart)
 
 % plot fieldlines