% Example 10: Stokes flow in a multiply-connected pipe-like domain, using
% the combined-layer Stokes representation u = (D - S)[sigma] with the
% nullspace correction W = normonesmat.

mip load --install magland/magland/chunkie;

tic;

% outer boundary: peanut shape
modes = [2.5,0,0,1];
ctr = [0;0];
chnkrouter = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));

% inner boundaries: small circles, reversed so orientations are consistent
chnkrcirc = chunkerfunc(@(t) chnk.curves.bymode(t,0.25,[0;0]));
chnkrcirc = reverse(chnkrcirc);
centers = [ [-2:2, -2:2]; [(0.7 + 0.2*(-1).^(-2:2)) , ...
    (-0.7 + 0.2*(-1).^(-2:2))]];
centers = centers + 0.1*randn(size(centers));

% assemble outer + shifted circles into one chunker
chnkrlist = [chnkrouter];
for j = 1:size(centers,2)
    chnkr1 = chnkrcirc;
    chnkr1 = chnkr1.move([0;0],centers(:,j));
    chnkrlist = [chnkrlist chnkr1];
end
chnkr = merge(chnkrlist);

% boundary velocity (Gaussian profile on the outer boundary, zero on holes)
wid = 0.3;
f = @(r) [exp(-r(2,:).^2/(2*wid^2)); zeros(size(r(2,:)))];
rhsout = f(chnkrouter.r(:,:)); rhsout = rhsout(:);
rhs = [rhsout; zeros(2*10*chnkrcirc.npt,1)];

% combined-layer Stokes kernel D - S
c = -1;
mu = 1;
kerncvel = kernel('stok','dvel',mu) + c*kernel('stok','svel',mu);

% matrix discretization of the BIE
cmat = chunkermat(chnkr,kerncvel);

% identity term and nullspace correction
W = normonesmat(chnkr);
sysmat = cmat - 0.5*eye(2*chnkr.npt) + W;

% solve
sigma = gmres(sysmat,rhs,[],1e-10,100);

% grid for plotting (interior)
x1 = linspace(-3.75,3.75,200);
y1 = linspace(-2,2,100);
[xx,yy] = meshgrid(x1,y1);
targs = [xx(:).'; yy(:).'];
in = chunkerinterior(chnkr,{x1,y1});
uu = nan([2,size(xx)]);
pres = nan(size(xx));

% velocity from layer potential
uu(:,in) = reshape(chunkerkerneval(chnkr,kerncvel,sigma,targs(:,in)),2,nnz(in));

% pressure from the corresponding pressure kernel
kerncpres = kernel('stok','dpres',mu) + c*kernel('stok','spres',mu);
opts = []; opts.eps = 1e-3;
pres(in) = chunkerkerneval(chnkr,kerncpres,sigma,targs(:,in),opts);

toc;

figure(1); clf
h = pcolor(xx,yy,pres); set(h,'EdgeColor','none'); colorbar
colormap("parula");
hold on
plot(chnkr,'k','LineWidth',2); axis equal

u = reshape(uu(1,:,:),size(xx)); v = reshape(uu(2,:,:),size(xx));
startt = linspace(pi-pi/12,pi+pi/12,20);
startr = 0.99*chnk.curves.bymode(startt,modes,[0;0]);
startx = startr(1,:); starty = startr(2,:);

sl = streamline(xx,yy,u,v,startx,starty);
set(sl,'LineWidth',2,'Color','w')