# Marker-and-Cell (MAC) Finite Volume Solver for the steady incompressible
# Stokes equations in 2 dimensions. Domain boundaries (inflow,outflow,walls) are
# defined using a mask vector.
#
# Author: Gustavo Buscaglia
# Last Modification: 16/9/2016

#Problem definition
global N1 N2
Lx=15.0;
Ly=10.0;
N1=40; #number of cells in the domain along x
N2=40;  #number of cells in the domain along y
#Material properties
mu = 1/10;
#Default mesh
aux=Lx/N1;
Nx=N1+2; #total cells along x
X=[-aux, 0:aux:Lx, Lx+aux]; 
aux=Ly/N2;
Ny=N2+2; #total cells along x
Y=[-aux, 0:aux:Ly, Ly+aux];
#Cell sizes
dx=X(2:Nx+1)-X(1:Nx);
dy=Y(2:Ny+1)-Y(1:Ny);
#Cell center coordinates
Xh=0.5*(X(2:Nx+1)+X(1:Nx));
Yh=0.5*(Y(2:Ny+1)+Y(1:Ny));
dxh=Xh(2:Nx)-Xh(1:Nx-1);
dyh=Yh(2:Ny)-Yh(1:Ny-1);
#Postprocessing
NPRIN = 10;
NSAVE = 20;
#Arrays
dimU = (Nx+1)*Ny;
dimV = Nx*(Ny+1);
dimP = Nx*Ny;
#Global matrix and rhs
nunk = dimU+dimV+dimP;
Ag = spalloc(nunk,nunk,27*nunk);
rhs= zeros(nunk,1);
#Boundary conditions and Mask definition
ibcx=zeros(10);vbcx=ibcx;
ibcy=zeros(10);vbcy=ibcy;
FLUID=0;
OPEN1=1;
ibcx(OPEN1)=2; #fx imposed
vbcx(OPEN1)=10.0; #imposed value
ibcy(OPEN1)=1; #vy imposed
vbcy(OPEN1)=0.0; #imposed value
OPEN2=2;
ibcx(OPEN2)=1; #vx imposed
vbcx(OPEN2)=0.0; #imposed value
ibcy(OPEN2)=2; #fy imposed
vbcy(OPEN2)=-10.0; #imposed value
OPEN3=3;
ibcx(OPEN3)=2; #fx imposed
vbcx(OPEN3)=0.0; #imposed value
ibcy(OPEN3)=1; #vy imposed
vbcy(OPEN3)=0.0; #imposed value
OPEN4=4;
ibcx(OPEN4)=1; #vx imposed
vbcx(OPEN4)=0.0; #imposed value
ibcy(OPEN4)=2; #fy imposed
vbcy(OPEN4)=-5.0; #imposed value
WALL =5;
ibcx(WALL)=1; #vx imposed
vbcx(WALL)=0.0; #imposed value
ibcy(WALL)=1; #vy imposed
vbcy(WALL)=0.0; #imposed value
mask=FLUID*ones(Nx,Ny);
disp("creating mask ..."); fflush(stdout);
for i=1:Nx
   for j=1:Ny
     ###DOMAIN BOUNDARY
     if (i == 1 )
       mask(i,j) = WALL;
       if ((Yh(j)>0.4*Ly)&&(Yh(j)<0.6*Ly))
         mask(i,j)=OPEN1;
       endif
     endif
     if (i == Nx)
       mask(i,j) = WALL;
       if ((Yh(j)>0.2*Ly)&&(Yh(j)<0.6*Ly))
         mask(i,j)=OPEN3;
       endif
     endif
     if (j == 1)
       mask(i,j) = WALL;
       if ((Xh(i)>0.4*Lx)&&(Xh(i)<0.6*Lx))
         mask(i,j)=OPEN4;
       endif
     endif
     if (j == Ny)
      mask(i,j) = WALL;
      if ((Xh(i)>0.2*Lx)&&(Xh(i)<0.4*Lx))
       mask(i,j) = OPEN2;
      endif
     endif
   endfor
endfor
##########################################################
###X-momentum equation ((Nx+1)*Ny control volumes)
##The leftmost and rightmost unknowns are always dummy
for j=1:Ny
   nguP=ij2ng(Nx+1,1,j);
   Ag(nguP,nguP)=1; rhs(nguP)=0;
   nguP=ij2ng(Nx+1,Nx+1,j);
   Ag(nguP,nguP)=1; rhs(nguP)=0;
endfor
##The upper and lower rows are always dummy
for i=1:Nx+1
   nguP=ij2ng(Nx+1,i,1);
   Ag(nguP,nguP)=1; rhs(nguP)=0;
   nguP=ij2ng(Nx+1,i,Ny);
   Ag(nguP,nguP)=1; rhs(nguP)=0;
endfor
##Now the inner rows and columns
for i=2:Nx
 for j=2:Ny-1
  nguP=ij2ng(Nx+1,i,j);
  mskl=mask(i-1,j);mskr=mask(i,j); ## masks of left and right half cells
  if ((mskl*mskr)!=0) ## cell inside boundary, dummy
   Ag(nguP,nguP)=1; rhs(nguP)=0;
  elseif ((mskl!=0) && (ibcx(mskl)==1)) ## velocity imposed from the left
   Ag(nguP,nguP)=1; rhs(nguP)=vbcx(mskl);
  elseif ((mskr!=0) && (ibcx(mskr)==1)) ## velocity imposed from the right
   Ag(nguP,nguP)=1; rhs(nguP)=vbcx(mskr);
  else                               ## unknown is alive and not imposed
   if ((mskl==0)&&(mskr==0))                      ## fully fluid cell
    ### W boundary (integral of -p+2*muW*du/dx)
    muW=mu;
    nguWW=ij2ng(Nx+1,i-1,j); ngpW=ij2ng(Nx,i-1,j)+dimU+dimV;
    Ag(nguP,nguP)=Ag(nguP,nguP)+2*muW*dy(j)/dx(i-1);
    Ag(nguP,nguWW)=Ag(nguP,nguWW)-2*muW*dy(j)/dx(i-1);
    Ag(nguP,ngpW)=-dy(j);
    ### E boundary (integral of p-2*muE*du/dx)
    muE=mu;
    nguEE=ij2ng(Nx+1,i+1,j); ngpE=ij2ng(Nx,i,j)+dimU+dimV;
    Ag(nguP,nguP)=Ag(nguP,nguP)+2*muE*dy(j)/dx(i);
    Ag(nguP,nguEE)=Ag(nguP,nguEE)-2*muE*dy(j)/dx(i);
    Ag(nguP,ngpE)=dy(j);
    ### N boundary (integral of -mu(dv/dx))
    muN=mu;
    nguNN=ij2ng(Nx+1,i,j+1);ngvNW=ij2ng(Nx,i-1,j+1)+dimU;ngvNE=ij2ng(Nx,i,j+1)+dimU;
    Ag(nguP,ngvNE)=Ag(nguP,ngvNE)-muN*dxh(i-1)/dxh(i-1);
    Ag(nguP,ngvNW)=Ag(nguP,ngvNW)+muN*dxh(i-1)/dxh(i-1);
    ### N boundary (integral of -mu(du/dy))
    if ((mask(i-1,j+1)==0)&&(mask(i,j+1)==0)) ## north cell boundary fully fluid
     Ag(nguP,nguP) =Ag(nguP,nguP) +muN*dxh(i-1)/dyh(j);
     Ag(nguP,nguNN)=Ag(nguP,nguNN)-muN*dxh(i-1)/dyh(j);
    elseif ((mask(i-1,j+1)!=0)&&(mask(i,j+1))!=0) ##N cells both boundary
     Ag(nguP,nguP) =Ag(nguP,nguP) +muN*dxh(i-1)/(dy(j)/2);
     rhs(nguP)     =rhs(nguP)     +...
      muN*dxh(i-1)/(dy(j)/2)*(vbcx(mask(i-1,j+1))+vbcx(mask(i,j+1)))/2;
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
    ### S boundary (integral of mu(dv/dx))
    muS=mu;
    nguSS=ij2ng(Nx+1,i,j-1);ngvSW=ij2ng(Nx,i-1,j)+dimU;ngvSE=ij2ng(Nx,i,j)+dimU;
    Ag(nguP,ngvSE)=Ag(nguP,ngvSE)+muS*dxh(i-1)/dxh(i-1);
    Ag(nguP,ngvSW)=Ag(nguP,ngvSW)-muS*dxh(i-1)/dxh(i-1);
    ### S boundary (integral of mu(du/dy))
    if ((mask(i-1,j-1)==0)&&(mask(i,j-1)==0)) ## south cell boundary fully fluid
     Ag(nguP,nguP) =Ag(nguP,nguP) +muS*dxh(i-1)/dyh(j-1);
     Ag(nguP,nguSS)=Ag(nguP,nguSS)-muS*dxh(i-1)/dyh(j-1);
    elseif ((mask(i-1,j-1)!=0)&&(mask(i,j-1))!=0) ##S cells both boundary
                             ## assumed that both have tangential velocity imposed
     Ag(nguP,nguP) =Ag(nguP,nguP) +muS*dxh(i-1)/(dy(j)/2);
     rhs(nguP)     =rhs(nguP)     +...
      muS*dxh(i-1)/(dy(j)/2)*(vbcx(mask(i-1,j-1))+vbcx(mask(i,j-1)))/2;
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
   elseif (mskl!=0) ## left half of cell is boundary (right half is fluid)
                     ## boundary imposes x force and y velocity
    ### W boundary (rhs+=integral of force)
    rhs(nguP)=rhs(nguP)+dy(j)*vbcx(mskl);
    ### E boundary (integral of p-2*muE*du/dx)
    muE=mu;
    nguEE=ij2ng(Nx+1,i+1,j); ngpE=ij2ng(Nx,i,j)+dimU+dimV;
    Ag(nguP,nguP)=Ag(nguP,nguP)+2*muE*dy(j)/dx(i);
    Ag(nguP,nguEE)=Ag(nguP,nguEE)-2*muE*dy(j)/dx(i);
    Ag(nguP,ngpE)=dy(j);
    ### N boundary (integral of -mu(dv/dx), E half of upper cell bound.)
    muN=mu;
    if (mask(i,j+1)!=0)
     vN=vbcy(mskl);
    else
     vN=0.5*(vbcy(mskl)+vbcy(mask(i-1,j+1)));
    endif
    nguNN=ij2ng(Nx+1,i,j+1);ngvNE=ij2ng(Nx,i,j+1)+dimU;
    Ag(nguP,ngvNE)=Ag(nguP,ngvNE)-muN*(dx(i)/2)/(dx(i)/2);
    rhs(nguP)=rhs(nguP)-muN*(dx(i)/2)/(dx(i)/2)*vN;
    ### N boundary (integral of -mu(du/dy), E half)
    if (mask(i,j+1)==0) ## NE cell fluid, uNN is an unknown
     Ag(nguP,nguP) =Ag(nguP,nguP) +muN*0.5*dx(i)/dyh(j);
     Ag(nguP,nguNN)=Ag(nguP,nguNN)-muN*0.5*dx(i)/dyh(j);
    elseif ((mask(i-1,j+1)!=0)&&(mask(i,j+1))!=0) ##N cells both boundary
               ## in this case we take the value of u from cell (i,j+1)
     Ag(nguP,nguP) =Ag(nguP,nguP) +muN*(dx(i)/2)/(dy(j)/2);
     rhs(nguP)     =rhs(nguP)     +...
      muN*(dx(i)/2)/(dy(j)/2)*vbcx(mask(i,j+1));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
    ### S boundary (integral of mu(dv/dx), E half)
    muS=mu;
    if (mask(i,j-1)!=0) 
     vS=vbcy(mskl);
    else
     vS=0.5*(vbcy(mskl)+vbcy(mask(i-1,j-1)));
    endif
    nguSS=ij2ng(Nx+1,i,j-1);ngvSE=ij2ng(Nx,i,j)+dimU;
    Ag(nguP,ngvSE)=Ag(nguP,ngvSE)+muS*(dx(i)/2)/(dx(i)/2);
    rhs(nguP)=rhs(nguP)+muS*(dx(i)/2)/(dx(i)/2)*vS;
    ### S boundary (integral of mu(du/dy), E half)
    if (mask(i,j-1)==0) ## SE cell fluid, uSS is an unknown
     Ag(nguP,nguP) =Ag(nguP,nguP) +muS*0.5*dx(i)/dyh(j-1);
     Ag(nguP,nguSS)=Ag(nguP,nguSS)-muS*0.5*dx(i)/dyh(j-1);
    elseif ((mask(i-1,j-1)!=0)&&(mask(i,j-1))!=0) ##S cells both boundary
               ## in this case we take the value of u from cell (i,j-1)
     Ag(nguP,nguP) =Ag(nguP,nguP) +muS*(dx(i)/2)/(dy(j)/2);
     rhs(nguP)     =rhs(nguP)     +...
      muS*(dx(i)/2)/(dy(j)/2)*vbcx(mask(i,j-1));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
   elseif (mskr!=0) ## right half of cell is boundary (left half is fluid)
                     ## boundary imposes x force and y velocity
    ### E boundary (rhs+=integral of force)
    rhs(nguP)=rhs(nguP)+dy(j)*vbcx(mskr);
    ### W boundary (integral of -p+2*muW*du/dx)
    muW=mu;
    nguWW=ij2ng(Nx+1,i-1,j); ngpW=ij2ng(Nx,i-1,j)+dimU+dimV;
    Ag(nguP,nguP)=Ag(nguP,nguP)+2*muW*dy(j)/dx(i-1);
    Ag(nguP,nguWW)=Ag(nguP,nguWW)-2*muW*dy(j)/dx(i-1);
    Ag(nguP,ngpE)=-dy(j);
    ### N boundary (integral of -mu(dv/dx), W half)
    muN=mu;
    if (mask(i-1,j+1)!=0)
     vN=vbcy(mskr);
    else
     vN=0.5*(vbcy(mskr)+vbcy(mask(i,j+1)));
    endif
    nguNN=ij2ng(Nx+1,i,j+1);ngvNW=ij2ng(Nx,i-1,j+1)+dimU;
    Ag(nguP,ngvNW)=Ag(nguP,ngvNW)+muN*(dx(i-1)/2)/(dx(i-1)/2);
    rhs(nguP)=rhs(nguP)+muN*(dx(i-1)/2)/(dx(i-1)/2)*vN;
    ### N boundary (integral of -mu(du/dy), W half)
    if (mask(i-1,j+1)==0) ## NW cell fluid, uNN is an unknown
     Ag(nguP,nguP) =Ag(nguP,nguP) +muN*0.5*dx(i-1)/dyh(j);
     Ag(nguP,nguNN)=Ag(nguP,nguNN)-muN*0.5*dx(i-1)/dyh(j);
    elseif ((mask(i-1,j+1)!=0)&&(mask(i,j+1))!=0) ##N cells both boundary
               ## in this case we take the value of u from cell (i-1,j+1)
     Ag(nguP,nguP) =Ag(nguP,nguP) +muN*(dx(i-1)/2)/(dy(j)/2);
     rhs(nguP)     =rhs(nguP)     +...
      muN*(dx(i-1)/2)/(dy(j)/2)*vbcx(mask(i,j+1));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
    ### S boundary (integral of mu(dv/dx), W half)
    muS=mu;
    if (mask(i-1,j-1)!=0) 
     vS=vbcy(mskr);
    else
     vS=0.5*(vbcy(mskr)+vbcy(mask(i,j-1)));
    endif
    nguSS=ij2ng(Nx+1,i,j-1);ngvSW=ij2ng(Nx,i-1,j)+dimU;
    Ag(nguP,ngvSW)=Ag(nguP,ngvSW)-muS*(dx(i-1)/2)/(dx(i-1)/2);
    rhs(nguP)=rhs(nguP)-muS*(dx(i-1)/2)/(dx(i-1)/2)*vS;
    ### S boundary (integral of mu(du/dy), E half)
    if (mask(i-1,j-1)==0) ## SE cell fluid, uSS is an unknown
     Ag(nguP,nguP) =Ag(nguP,nguP) +muS*0.5*dx(i-1)/dyh(j-1);
     Ag(nguP,nguSS)=Ag(nguP,nguSS)-muS*0.5*dx(i-1)/dyh(j-1);
    elseif ((mask(i-1,j-1)!=0)&&(mask(i,j-1))!=0) ##S cells both boundary
               ## in this case we take the value of u from cell (i-1,j-1)
     Ag(nguP,nguP) =Ag(nguP,nguP) +muS*(dx(i-1)/2)/(dy(j)/2);
     rhs(nguP)     =rhs(nguP)     +...
      muS*(dx(i-1)/2)/(dy(j)/2)*vbcx(mask(i-1,j-1));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
   else
    disp("BC crash..."); fflush(stdout);
   endif
  endif
 endfor
endfor
###############End of X-Momentum equations
###############Now Y-Momentum
###Y-momentum equation (Nx*(Ny+1) control volumes)
##The northmost and southmost unknowns are always dummy
for i=1:Nx
   ngvP=ij2ng(Nx,i,1)+dimU;
   Ag(ngvP,ngvP)=1; rhs(ngvP)=0;
   ngvP=ij2ng(Nx,i,Ny+1)+dimU;
   Ag(ngvP,ngvP)=1; rhs(ngvP)=0;
endfor
##The east and west columns are always dummy
for j=1:Ny+1
   ngvP=ij2ng(Nx,1,j)+dimU;
   Ag(ngvP,ngvP)=1; rhs(ngvP)=0;
   ngvP=ij2ng(Nx,Nx,j)+dimU;
   Ag(ngvP,ngvP)=1; rhs(ngvP)=0;
endfor
##Now the inner rows and columns
for i=2:Nx-1
 for j=2:Ny
  ngvP=ij2ng(Nx,i,j)+dimU;
  msks=mask(i,j-1);mskn=mask(i,j); ## masks of south and north half cells
  if (msks*mskn!=0) ## cell inside boundary, dummy
   Ag(ngvP,ngvP)=1; rhs(ngvP)=0;
  elseif (msks!=0 && ibcy(msks)==1) ## velocity imposed from the south
   Ag(ngvP,ngvP)=1; rhs(ngvP)=vbcy(msks);
  elseif (mskn!=0 && ibcy(mskn)==1) ## velocity imposed from the north
   Ag(ngvP,ngvP)=1; rhs(ngvP)=vbcy(mskn);
  else                               ## unknown is alive and not imposed
   if ((msks==0)&&(mskn==0))                      ## fully fluid cell
    ### S boundary (integral of -p+2*muS*dv/dy)
    muS=mu;
    ngvSS=ij2ng(Nx,i,j-1)+dimU; ngpS=ij2ng(Nx,i,j-1)+dimU+dimV;
    Ag(ngvP,ngvP)=Ag(ngvP,ngvP)+2*muS*dx(i)/dy(j-1);
    Ag(ngvP,ngvSS)=Ag(ngvP,ngvSS)-2*muS*dx(i)/dy(j-1);
    Ag(ngvP,ngpS)=-dx(i);
    ### N boundary (integral of p-2*muN*dv/dy)
    muN=mu;
    ngvNN=ij2ng(Nx,i,j+1)+dimU; ngpN=ij2ng(Nx,i,j)+dimU+dimV;
    Ag(ngvP,ngvP)=Ag(ngvP,ngvP)+2*muN*dx(i)/dy(j);
    Ag(ngvP,ngvNN)=Ag(ngvP,ngvNN)-2*muN*dx(i)/dy(j);
    Ag(ngvP,ngpN)=dx(i);
    ### E boundary (integral of -mu(du/dy))
    muE=mu;
    ngvEE=ij2ng(Nx,i+1,j)+dimU;nguNE=ij2ng(Nx+1,i+1,j);nguSE=ij2ng(Nx+1,i+1,j-1);
    Ag(ngvP,nguNE)=Ag(ngvP,nguNE)-muE*dyh(j-1)/dyh(j-1);
    Ag(ngvP,nguSE)=Ag(ngvP,ngvSE)+muE*dyh(j-1)/dyh(j-1);
    ### E boundary (integral of -mu(dv/dx))
    if ((mask(i+1,j)==0)&&(mask(i+1,j-1)==0)) ## east cell boundary fully fluid
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muE*dyh(j-1)/dxh(i);
     Ag(ngvP,ngvEE)=Ag(ngvP,ngvEE)-muE*dyh(j-1)/dxh(i);
    elseif ((mask(i+1,j)!=0)&&(mask(i+1,j-1))!=0) ## E cells both boundary
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muE*dyh(j-1)/(dx(i)/2);
     rhs(ngvP)     =rhs(ngvP)     +...
      muE*dyh(j-1)/(dx(i)/2)*(vbcy(mask(i+1,j))+vbcy(mask(i+1,j-1)))/2;
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
    ### W boundary (integral of mu(du/dy))
    muW=mu;
    ngvWW=ij2ng(Nx,i-1,j)+dimU;nguNW=ij2ng(Nx+1,i,j);nguSW=ij2ng(Nx+1,i,j-1);
    Ag(ngvP,nguNW)=Ag(ngvP,nguNW)+muW*dyh(j-1)/dyh(j-1);
    Ag(ngvP,nguSW)=Ag(ngvP,nguSW)-muW*dyh(j-1)/dyh(j-1);
    ### W boundary (integral of mu(dv/dx))
    if ((mask(i-1,j-1)==0)&&(mask(i-1,j)==0)) ## W cell boundary fully fluid
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muW*dyh(j-1)/dxh(i-1);
     Ag(ngvP,ngvWW)=Ag(ngvP,ngvWW)-muW*dyh(j-1)/dxh(i-1);
    elseif ((mask(i-1,j-1)!=0)&&(mask(i-1,j))!=0) ## W cells both boundary
                             ## assumed that both have tangential velocity imposed
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muW*dyh(j-1)/(dx(i)/2);
     rhs(ngvP)     =rhs(ngvP)     +...
      muW*dyh(j-1)/(dx(i)/2)*(vbcy(mask(i-1,j-1))+vbcy(mask(i-1,j)))/2;
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
   elseif (msks!=0) ## south half of cell is boundary (north half is fluid)
                     ## boundary imposes y force and x velocity
    ### S boundary (rhs+=integral of force)
    rhs(ngvP)=rhs(ngvP)+dx(i)*vbcy(msks);
    ### N boundary (integral of p-2*muN*dv/dy)
    muN=mu;
    ngvNN=ij2ng(Nx,i,j+1)+dimU; ngpN=ij2ng(Nx,i,j)+dimU+dimV;
    Ag(ngvP,ngvP)=Ag(ngvP,ngvP)+2*muN*dx(i)/dy(j);
    Ag(ngvP,ngvNN)=Ag(ngvP,ngvNN)-2*muN*dx(i)/dy(j);
    Ag(ngvP,ngpN)=dx(i);
    ### E boundary (integral of -mu(du/dy), N half of E cell bound.)
    muE=mu;
    if (mask(i+1,j)!=0)
     uE=vbcx(msks);
    else
     uE=0.5*(vbcx(msks)+vbcx(mask(i+1,j-1)));
    endif
    ngvEE=ij2ng(Nx,i+1,j)+dimU;nguNE=ij2ng(Nx+1,i+1,j);
    Ag(ngvP,nguNE)=Ag(ngvP,nguNE)-muE*(dy(j)/2)/(dy(j)/2);
    rhs(ngvP)=rhs(ngvP)-muE*(dy(j)/2)/(dy(j)/2)*uE;
    ### E boundary (integral of -mu(dv/dx), N half)
    if (mask(i+1,j)==0) ## NE cell fluid, vEE is an unknown
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muE*0.5*dy(j)/dxh(i);
     Ag(ngvP,ngvEE)=Ag(ngvP,ngvEE)-muE*0.5*dy(j)/dxh(i);
    elseif ((mask(i+1,j)!=0)&&(mask(i+1,j-1))!=0) ##E cells both boundary
               ## in this case we take the value of v from cell (i+1,j)
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muE*(dy(j)/2)/(dx(i)/2);
     rhs(ngvP)     =rhs(ngvP)     +...
      muE*(dy(j)/2)/(dx(i)/2)*vbcy(mask(i+1,j));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
    ### W boundary (integral of mu(du/dy), N half)
    muW=mu;
    if (mask(i-1,j)!=0) 
     uW=vbcx(msks);
    else
     uW=0.5*(vbcx(msks)+vbcx(mask(i-1,j-1)));
    endif
    ngvWW=ij2ng(Nx,i-1,j)+dimU;nguNW=ij2ng(Nx+1,i,j);
    Ag(ngvP,nguNW)=Ag(ngvP,nguNW)+muW*(dy(j)/2)/(dy(j)/2);
    rhs(ngvP)=rhs(ngvP)+muW*(dy(j)/2)/(dy(j)/2)*uW;
    ### W boundary (integral of mu(dv/dx), N half)
    if (mask(i-1,j)==0) ## NW cell fluid, vWW is an unknown
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muW*0.5*dy(j)/dxh(i-1);
     Ag(ngvP,ngvWW)=Ag(ngvP,ngvWW)-muW*0.5*dy(j)/dxh(i-1);
    elseif ((mask(i-1,j-1)!=0)&&(mask(i-1,j))!=0) ## W cells both boundary
               ## in this case we take the value of v from cell (i-1,j)
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muW*(dy(j)/2)/(dx(i)/2);
     rhs(ngvP)     =rhs(ngvP)     +...
      muW*(dy(j)/2)/(dx(i)/2)*vbcy(mask(i-1,j));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
   elseif (mskn!=0) ## north half of cell is boundary (south half is fluid)
                     ## boundary imposes y force and x velocity
    ### N boundary (rhs+=integral of force)
    rhs(ngvP)=rhs(ngvP)+dx(i)*vbcy(mskn);
    ### S boundary (integral of -p+2*muS*dv/dy)
    muS=mu;
    ngvSS=ij2ng(Nx,i,j-1)+dimU; ngpS=ij2ng(Nx,i,j-1)+dimU+dimV;
    Ag(ngvP,ngvP)=Ag(ngvP,ngvP)+2*muN*dx(i)/dy(j-1);
    Ag(ngvP,ngvSS)=Ag(ngvP,ngvSS)-2*muN*dx(i)/dy(j-1);
    Ag(ngvP,ngpS)=-dx(i);
    ### E boundary (integral of -mu(du/dy), S half of E cell bound.)
    muE=mu;
    if (mask(i+1,j-1)!=0)
     uE=vbcx(mskn);
    else
     uE=0.5*(vbcx(mskn)+vbcx(mask(i+1,j)));
    endif
    ngvEE=ij2ng(Nx,i+1,j)+dimU;nguSE=ij2ng(Nx+1,i+1,j-1);
    Ag(ngvP,nguSE)=Ag(ngvP,nguSE)+muE*(dy(j-1)/2)/(dy(j-1)/2);
    rhs(ngvP)=rhs(ngvP)+muE*(dy(j-1)/2)/(dy(j-1)/2)*uE;
    ### E boundary (integral of -mu(dv/dx), S half)
    if (mask(i+1,j-1)==0) ## SE cell fluid, vEE is an unknown
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muE*0.5*dy(j-1)/dxh(i);
     Ag(ngvP,ngvEE)=Ag(ngvP,ngvEE)-muE*0.5*dy(j-1)/dxh(i);
    elseif ((mask(i+1,j)!=0)&&(mask(i+1,j-1))!=0) ##E cells both boundary
               ## in this case we take the value of v from cell (i+1,j-1)
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muE*(dy(j-1)/2)/(dx(i)/2);
     rhs(ngvP)     =rhs(ngvP)     +...
      muE*(dy(j-1)/2)/(dx(i)/2)*vbcy(mask(i+1,j-1));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
    ### W boundary (integral of mu(du/dy), S half)
    muW=mu;
    if (mask(i-1,j-1)!=0) 
     uW=vbcx(mskn);
    else
     uW=0.5*(vbcx(mskn)+vbcx(mask(i-1,j)));
    endif
    ngvWW=ij2ng(Nx,i-1,j)+dimU;nguSW=ij2ng(Nx+1,i,j-1);
    Ag(ngvP,nguSW)=Ag(ngvP,nguSW)-muW*(dy(j-1)/2)/(dy(j-1)/2);
    rhs(ngvP)=rhs(ngvP)-muW*(dy(j-1)/2)/(dy(j-1)/2)*uW;
    ### W boundary (integral of mu(dv/dx), S half)
    if (mask(i-1,j-1)==0) ## SW cell fluid, vWW is an unknown
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muW*0.5*dy(j-1)/dxh(i-1);
     Ag(ngvP,ngvWW)=Ag(ngvP,ngvWW)-muW*0.5*dy(j-1)/dxh(i-1);
    elseif ((mask(i-1,j-1)!=0)&&(mask(i-1,j))!=0) ## W cells both boundary
               ## in this case we take the value of v from cell (i-1,j-1)
     Ag(ngvP,ngvP) =Ag(ngvP,ngvP) +muW*(dy(j-1)/2)/(dx(i)/2);
     rhs(ngvP)     =rhs(ngvP)     +...
      muW*(dy(j-1)/2)/(dx(i)/2)*vbcy(mask(i-1,j-1));
    else
     disp("BC not implemented! ..."); fflush(stdout);
    endif 
   else
    disp("BC crash..."); fflush(stdout);
   endif
  endif
 endfor
endfor
###############End of Y-Momentum equations
###############Beginning of Mass equations
###Mass equation (Nx*Ny control volumes)
##The northmost (j=Ny) and southmost (j=1) lines are always dummy
for i=1:Nx
   ngpP=ij2ng(Nx,i,1)+dimU+dimV;
   Ag(ngpP,ngpP)=1; rhs(ngpP)=0;
   ngpP=ij2ng(Nx,i,Ny)+dimU+dimV;
   Ag(ngpP,ngpP)=1; rhs(ngpP)=0;
endfor
##The east (i=1) and west (i=Nx) columns are always dummy
for j=1:Ny
   ngpP=ij2ng(Nx,1,j)+dimU+dimV;
   Ag(ngpP,ngpP)=1; rhs(ngpP)=0;
   ngpP=ij2ng(Nx,Nx,j)+dimU+dimV;
   Ag(ngpP,ngpP)=1; rhs(ngpP)=0;
endfor
##Now the inner rows and columns
for i=2:Nx-1
 for j=2:Ny-1
  ngpP=ij2ng(Nx,i,j)+dimU+dimV;
  if (mask(i,j)!=0) ## cell in boundary, dummy
   Ag(ngpP,ngpP)=1; rhs(ngpP)=0;
  else                               ## pressure unknown is alive (fluid cell)
   ### S boundary (integral of -v)
   ngvS=ij2ng(Nx,i,j)+dimU;
   Ag(ngpP,ngvS)=Ag(ngpP,ngvS)-dx(i);
   ### N boundary (integral of v)
   ngvN=ij2ng(Nx,i,j+1)+dimU;
   Ag(ngpP,ngvN)=Ag(ngpP,ngvN)+dx(i);
   ### E boundary (integral of u)
   nguE=ij2ng(Nx+1,i+1,j);
   Ag(ngpP,nguE)=Ag(ngpP,nguE)+dy(j);
   ### W boundary (integral of -u)
   nguW=ij2ng(Nx+1,i,j);
   Ag(ngpP,nguW)=Ag(ngpP,nguW)-dy(j);
  endif
 endfor
endfor
###############End of Mass equations
xxx = Ag \ rhs;
##########Cell-average velocities
Uav=zeros(Nx,Ny);
Vav=zeros(Nx,Ny);
P=zeros(Nx,Ny);
U=zeros(Nx+1,Ny);
V=zeros(Nx,Ny+1);
for i=1:Nx+1
 for j=1:Ny
  U(i,j)=xxx(ij2ng(Nx+1,i,j));
 endfor
endfor
for i=1:Nx
 for j=1:Ny+1
  V(i,j)=xxx(ij2ng(Nx,i,j)+dimU);
 endfor
endfor
for i=1:Nx
 for j=1:Ny
  if (mask(i,j)==0)
   Uav(i,j)=0.5*(xxx(ij2ng(Nx+1,i,j))+xxx(ij2ng(Nx+1,i+1,j)));
   Vav(i,j)=0.5*(xxx(ij2ng(Nx,i,j)+dimU)+xxx(ij2ng(Nx,i,j+1)+dimU));
   P(i,j)=xxx(ij2ng(Nx,i,j)+dimU+dimV);
  endif
 endfor
endfor
###############
figure(1);
h=quiver(Xh,Yh,Uav',Vav',3);
set(h,"linewidth",2);
figure(2);
##contourf(X,Yh,U');
##contourf(Xh,Y,V');
##contourf(Xh(2:Nx-1),Y(2:Ny),V(2:Nx-1,2:Ny)')   ##only active ones
contourf(Xh(2:Nx-1),Yh(2:Ny-1),P(2:Nx-1,2:Ny-1)');