function interstitial_flow_discrete

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%              Moffitt PhD Program in Cancer Biology               %%
%%            lecture in Integrated Mathematical Oncology I         %%
%% This is a companion code for a lecture in Tumor Microenvironment %%
%% K.A. Rejniak, PhD, Moffitt Cancer Center & Research Institute,   %%
%% October 30, 2018                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% cell boundary coordinates
[cell,hb]=DefineCellBoundaryPoints; % cell coordinates & point separation
CellNum=size(cell,1);               % number of cells

% computational domain
xeps=15;                            % marigin from cell boundaries
xmin=min(min(cell(:,1)))-2*xeps; xmax=max(max(cell(:,1)))+xeps;  % x-axis boundaries
ymin=min(min(cell(:,2)))-xeps; ymax=max(max(cell(:,2)))+xeps;    % y-axis boundaries
x0=(xmax+xmin)/2; y0=(ymax+ymin)/2;                              % center of the domain

% time step
dt = 0.05;  

% fluid properties
mu = 0.04;   % fluid viscosity
hg = 10;     % fluid grid width

% fluid grid
Ngx=floor((xmax-xmin)/hg);
Ngy=floor((ymax-ymin)/hg);
for ii=1:Ngx+1
  for jj=1:Ngy+1
    xyg(ii,jj,1)=xmin+(ii-1)*hg;
    xyg(ii,jj,2)=ymin+(jj-1)*hg;
  end
end

% fluid velocities
uxg=zeros(Ngx+1,Ngy+1);  %% fluid velocity x-coordinate
uyg=zeros(Ngx+1,Ngy+1);  %% fluid velocity y-coordinate
    
% influx bundary velocities
Nz=floor((ymax-ymin)/hb);
xz(1:Nz,1)=xmin*ones(Nz,1);     % left edge coordinates
xz(1:Nz,2)=[ymin+hb:hb:ymax];
uz=zeros(Nz,2);
uz(1:Nz,1)=60*ones(Nz,1);       % parallel fluid flow values
uz1=uz;

% zero bundary velocities (top & bottom edge)
Nx=1+floor((xmax-xmin)/hb);
xx=ones(2*Nx,2);
xx(1:Nx,1)=[xmin:hb:xmax];  % bottom edge coordinates
xx(1:Nx,2)=ymin;
xx(Nx+1:2*Nx,1)=[xmin:hb:xmax];
xx(Nx+1:2*Nx,2)=ymax;       % top edge coordinates
ux=zeros(2*Nx,2);           % zero velocities on both edges 

% zero cell boundary velocities
Npts=CellNum;
xb(1:Npts,1:2)=cell(1:Npts,1:2); % cell boundaries coordinates
ub=zeros(Npts,2);                % zero velocity on cell boundaries


% draw a figure
h_fig=figure('position',[300,300,750,500]);
 set(h_fig,'PaperPositionMode','auto')  
 set(h_fig,'InvertHardcopy','off')
 plot(cell(:,1),cell(:,2),'r.','MarkerSize',10)   % cells
 axis([xmin,xmax,ymin,ymax])
 axis equal
 hold on
 plot(xx(:,1),xx(:,2),'c.')   % boundaries
 plot(xz(:,1),xz(:,2),'m.')
 pause(0.1)
 
 
% solve forces ff and fb from velocities: uf and ub
  xtot=zeros(Npts+Nz+2*Nx,2); % create one long vector of cell and edge pts
  xtot(1:Npts,1:2)=xb(1:Npts,1:2);
  xtot(Npts+1:Npts+Nz,1:2)=xz(1:Nz,1:2);
  xtot(Npts+Nz+1:end,1:2)=xx(1:2*Nx,1:2);
  
  utot=zeros(Npts+Nz+2*Nx,2); % create one long vector of cell and edge velo
  utot(1:Npts,1:2)=ub(1:Npts,1:2);
  utot(Npts+1:Npts+Nz,1:2)=uz(1:Nz,1:2);
  utot(Npts+Nz+1:end,1:2)=ux(1:2*Nx,1:2);

  disp(['total number of points: ',num2str(size(xtot,1))]);
  disp('  ');
    
  % solve regularized Stokelets equations for forces from known velocities  
  disp('calculate forces for a given velocity: '); 
  ftot=ForceFromVelo(xtot,Npts+Nz+2*Nx,xtot,utot,Npts+Nz+2*Nx,mu,hb);
  
  % split forces into separate vectors
  fb(1:Npts,1:2)=ftot(1:Npts,1:2); 
  fz(1:Nz  ,1:2)=ftot(Npts+1:Npts+Nz,1:2); 
  fx(1:2*Nx,1:2)=ftot(Npts+Nz+1:end,1:2); 
    
  % draw all forces
  quiver(xb(:,1),xb(:,2),fb(:,1),fb(:,2),1,'c')
  quiver(xz(1:end,1),xz(1:end,2),fz(1:end,1),fz(1:end,2),1.5,'c')
  quiver(xx(1:end,1),xx(1:end,2),fx(1:end,1),fx(1:end,2),5,'c')
  axis([xmin,xmax,ymin,ymax])
  pause(0.1)
  

  % solve velocities on the grid from the known forces
  disp('calculate velocities for given forces: '); 
  [uxg,uyg]=VeloFromForce2D(xyg,Ngx+1,Ngy+1,xtot,ftot,Npts+Nz+2*Nx,mu,hb);
    
  % solve velocities on the cell boundaries  
  ub=VeloFromForce(xb,Npts,xtot,ftot,Npts+Nz+2*Nx,mu,hb);

  % solve velocities on the domain edge boundaries
  ux=VeloFromForce(xx,2*Nx,xtot,ftot,Npts+Nz+2*Nx,mu,hb);
  
  % solve velocities on the influx domain 
  uz=VeloFromForce(xz,Nz,xtot,ftot,Npts+Nz+2*Nx,mu,hb);



  % draw the fluid velocity field
  quiver(xyg(:,:,1),xyg(:,:,2),uxg(:,:),uyg(:,:),1.5,'b') 
 
  % draw the velocities on domain influx  
  quiver(xz(:,1),xz(:,2),uz(:,1),uz(:,2),0.25,'m')
  
  axis([xmin,xmax,ymin,ymax])
  pause(1) 
  
  
  
  disp('draw drug particles moving with the fluid flow: ')
  Ndrug=0; Nmolecule=size(ymin:10:ymax,2);
  for kk=1:175  % move the drug particles
    if (kk<30)
      drug(Ndrug+1:Ndrug+Nmolecule,1)=60*ones(Nmolecule,1); % x coordinates
      drug(Ndrug+1:Ndrug+Nmolecule,2)=ymin:10:ymax;         % y coordinates
      Ndrug=Ndrug+Nmolecule;
    end
    
    % calculate advection
    udrug=VeloFromForce(drug,Ndrug,xtot,ftot,Npts+Nz+2*Nx,mu,hb);
    drug=drug+udrug*dt;   % calculate new drug particle positions

    clf  
    axis([xmin,xmax,ymin,ymax])
    axis equal
    axis([xmin,xmax,ymin,ymax])
    hold on
    plot(drug(:,1),drug(:,2),'k.','Markersize',10)   
    quiver(xyg(:,:,1),xyg(:,:,2),uxg(:,:),uyg(:,:),1.5,'b') % fluid flow
    plot(cell(:,1),cell(:,2),'r.','MarkerSize',10)   % cells
    plot(xx(:,1),xx(:,2),'c.')   % boundaries
    plot(xz(:,1),xz(:,2),'m.')
        
    axis([xmin,xmax,ymin,ymax])
    axis equal
    axis([xmin,xmax,ymin,ymax])
    hold on
    pause(0.1)
  end

end    % end program


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function ForceFromVelo                                              %
% calculates forces at points yy of size Ny that will generate given  %
% velocities ux at points xx; calculates the regularized Stokeslets   %
% transition matrix MM, uses Matlab routines for the LU decomposition % 
% and the gmres method.                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fy=ForceFromVelo(yy,Ny,xx,ux,Nx,mu,eps)
  fy=zeros(Ny,2);  fyy=zeros(2*Ny,1);
  MM=zeros(2*Ny,2*Ny);
  uy=zeros(2*Ny,1); uy(1:2:2*Ny,1)=ux(1:Nx,1); uy(2:2:2*Ny,1)=ux(1:Nx,2);
  
  W=1/(4.0*pi*mu); 
  ep2=eps*eps;
  for ii=1:Ny  
    for jj=1:Nx
      dx=yy(ii,1)-xx(jj,1); dy=yy(ii,2)-xx(jj,2);
      r2=dx*dx+dy*dy;
      Hep=0.5*log(r2+ep2)-(ep2/(r2+ep2)); Jep=1/(r2+ep2);
      AA=W*[Jep*dx*dx-Hep,Jep*dx*dy;Jep*dx*dy,Jep*dy*dy-Hep];
      MM(2*ii-1:2*ii,2*jj-1:2*jj)=AA;
    end
  end

  [L1,U1] = lu(MM);
  fyy=gmres(MM,uy,5,1e-5,10,L1,U1);  
  fy(1:end,1)=fyy(1:2:end,1); fy(1:end,2)=fyy(2:2:end,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function vx=VeloFromForce                                            %
% calculates velocities at points xx of size Nx generated by the given %
% forces fy at points yy; calculates the regularized Stokeslets tran-  %
% sition matrix MM                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vx=VeloFromForce(xx,Nx,yy,fy,Ny,mu,eps)
  vx=zeros(Nx,2);

  W=1/(4.0*pi*mu);
  ep2=eps*eps;
  for ii=1:Nx
    for jj=1:Ny
      dx=xx(ii,1)-yy(jj,1); dy=xx(ii,2)-yy(jj,2);
      r2=dx*dx+dy*dy;
      Hep=0.5*log(r2+ep2)-(ep2/(r2+ep2)); Jep=1/(r2+ep2);

      vx(ii,1)=vx(ii,1)-W*Hep*fy(jj,1)+W*dx*Jep*(fy(jj,1)*dx+fy(jj,2)*dy);
      vx(ii,2)=vx(ii,2)-W*Hep*fy(jj,2)+W*dy*Jep*(fy(jj,1)*dx+fy(jj,2)*dy);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function vx=VeloFromForce2D                                          %
% calculates velocities at points xx of sizes Nx1*Nx2 generated by the %
% given forces fy at points yy; calculates the regularized Stokeslets  %
% transition matrix MM                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [vx,vy]=VeloFromForce2D(xx,Nx1,Nx2,yy,fy,Ny,mu,eps)
  vx=zeros(Nx1,Nx2); vy=zeros(Nx1,Nx2);
  
  W=1/(4.0*pi*mu);
  ep2=eps*eps;
  for jj=1:Ny
    for ii=1:Nx1
      for kk=1:Nx2  
        dx=xx(ii,kk,1)-yy(jj,1); dy=xx(ii,kk,2)-yy(jj,2);
        r2=dx*dx+dy*dy;
        Hep=0.5*log(r2+ep2)-(ep2/(r2+ep2)); Jep=1/(r2+ep2);

        vx(ii,kk)=vx(ii,kk)-W*Hep*fy(jj,1)+W*dx*Jep*(fy(jj,1)*dx+fy(jj,2)*dy);
        vy(ii,kk)=vy(ii,kk)-W*Hep*fy(jj,2)+W*dy*Jep*(fy(jj,1)*dx+fy(jj,2)*dy);
      end
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function DefineCellBoundaryPoints                           %
% defines the boundary points of all cells                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pts,hb]=DefineCellBoundaryPoints
  cells=1.0e+03*[...
    0.1037    0.3544; 0.3166    0.2106; 0.1756    0.1083; ...
    0.1147   -0.0659; 0.3083   -0.2096; 0.6429   -0.3783; ...
    0.5682   -0.1930; 0.8558   -0.3147; 1.0272   -0.3894; ...
    1.1433   -0.1654; 0.9276   -0.0686; 0.6124   -0.0023; ...
    0.3608    0.1249; 0.4991    0.2936; 0.6235    0.4291; ... 
    0.8917    0.3931; 1.1157    0.2715; 1.0382   -0.0133; ...
    0.9304    0.1664; 0.7092    0.1111; 0.2344   -0.3501];

  Ncell=size(cells,1);
  rad=25; hb=2*pi*rad/25; thn=hb/rad; th=0:thn:2*pi-thn; 
  Npts=0;
  for ii=1:Ncell 
    pts(Npts+1:Npts+length(th),1:2)=[cells(ii,1)+rad*cos(th);...
       cells(ii,2)+rad*sin(th)]'; 
    Npts=Npts+length(th); 
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%