1% Example 10: Stokes flow in a multiply-connected pipe-like domain, using
2% the combined-layer Stokes representation u = (D - S)[sigma] with the
3% nullspace correction W = normonesmat.
5mip load --install magland/magland/chunkie;
7tic;
9% outer boundary: peanut shape
10modes = [2.5,0,0,1];
11ctr = [0;0];
12chnkrouter = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
14% inner boundaries: small circles, reversed so orientations are consistent
15chnkrcirc = chunkerfunc(@(t) chnk.curves.bymode(t,0.25,[0;0]));
16chnkrcirc = reverse(chnkrcirc);
17centers = [ [-2:2, -2:2]; [(0.7 + 0.2*(-1).^(-2:2)) , ...
18 (-0.7 + 0.2*(-1).^(-2:2))]];
19centers = centers + 0.1*randn(size(centers));
21% assemble outer + shifted circles into one chunker
22chnkrlist = [chnkrouter];
23for j = 1:size(centers,2)
24 chnkr1 = chnkrcirc;
25 chnkr1 = chnkr1.move([0;0],centers(:,j));
26 chnkrlist = [chnkrlist chnkr1];
27end
28chnkr = merge(chnkrlist);
30% boundary velocity (Gaussian profile on the outer boundary, zero on holes)
31wid = 0.3;
32f = @(r) [exp(-r(2,:).^2/(2*wid^2)); zeros(size(r(2,:)))];
33rhsout = f(chnkrouter.r(:,:)); rhsout = rhsout(:);
34rhs = [rhsout; zeros(2*10*chnkrcirc.npt,1)];
36% combined-layer Stokes kernel D - S
37c = -1;
38mu = 1;
39kerncvel = kernel('stok','dvel',mu) + c*kernel('stok','svel',mu);
41% matrix discretization of the BIE
42cmat = chunkermat(chnkr,kerncvel);
44% identity term and nullspace correction
45W = normonesmat(chnkr);
46sysmat = cmat - 0.5*eye(2*chnkr.npt) + W;
48% solve
49sigma = gmres(sysmat,rhs,[],1e-10,100);
51% grid for plotting (interior)
52x1 = linspace(-3.75,3.75,200);
53y1 = linspace(-2,2,100);
54[xx,yy] = meshgrid(x1,y1);
55targs = [xx(:).'; yy(:).'];
56in = chunkerinterior(chnkr,{x1,y1});
57uu = nan([2,size(xx)]);
58pres = nan(size(xx));
60% velocity from layer potential
61uu(:,in) = reshape(chunkerkerneval(chnkr,kerncvel,sigma,targs(:,in)),2,nnz(in));
63% pressure from the corresponding pressure kernel
64kerncpres = kernel('stok','dpres',mu) + c*kernel('stok','spres',mu);
65opts = []; opts.eps = 1e-3;
66pres(in) = chunkerkerneval(chnkr,kerncpres,sigma,targs(:,in),opts);
68toc;
70figure(1); clf
71h = pcolor(xx,yy,pres); set(h,'EdgeColor','none'); colorbar
72colormap("parula");
73hold on
74plot(chnkr,'k','LineWidth',2); axis equal
76u = reshape(uu(1,:,:),size(xx)); v = reshape(uu(2,:,:),size(xx));
77startt = linspace(pi-pi/12,pi+pi/12,20);
78startr = 0.99*chnk.curves.bymode(startt,modes,[0;0]);
79startx = startr(1,:); starty = startr(2,:);
81sl = streamline(xx,yy,u,v,startx,starty);
82set(sl,'LineWidth',2,'Color','w')