% Cylinder Cross Flow setup file

%Define Forcing functions
Fu = @(time,p) 0;
Fv = @(time,p) 0;

%Set size of domain
a    = 60;          b    = 3;
dx   = a./(app.Nx);    dy   = b./(app.Ny); 

%Set boundary conditions
bcsDrho = [0 0 0 0];bcsNrho = [0];        
bcsDP= [0];         bcsNP = [0 0 0 0];        
bcsOP= [];
bcsDu= [0 1];  bcsNu= [0 0 0];        
bcsOu= [];
bcsDv= [0 0 0 0];  bcsNv=  [0];        

%Set number of timesteps
Nt = ceil(app.T/app.dt);

%% create geometry mask
Nbcs = length(bcsDP) + length(bcsNP);
bnd = [1 1 1 1];
Mask = zeros(app.Ny, app.Nx);

% specify center of the disk
xcdisk = 4 + 0.5;
ycdisk = b/2;
rdisk = 1/2;

xp = linspace(dx/2, a-dx/2, app.Nx  );yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xu = linspace( dx   ,a-dx , app.Nx-1);yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xv = linspace(dx/2, a-dx/2, app.Nx  );yv = fliplr(linspace(dy    , b-dy, app.Ny-1));
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);

if ~isfield(app,'sq')
    %Cylinder
    Mask((XP - xcdisk).^2 + (YP - ycdisk).^2 < rdisk.^2) = 1 ;
else
    % Do a square instead
    Mask(abs(XP - xcdisk) <= (rdisk) ...
        & abs(YP - ycdisk) <= (rdisk)) = 1 ;
end
%Make sure Mask is symmetric
Mask=(Mask+flipud(Mask));
Mask=Mask.*(max(1,Mask)).^-1;
xcdisk=app.Nx/5*dx;
ycdisk=b/2;
rdisk=1/2;

Mask=(Mask+flipud(Mask));
Mask=Mask.*(max(1,Mask)).^-1;
clear xtemp ytemp 
pcolor(Mask)

%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);
NodeRho=NodeP;
NodeP(NodeP==1)=-1;NodeP(NodeP==3)=1;NodeP(NodeP==-1)=3;
Nodev(Nodev==3)=-1;Nodev(Nodev==5)=3;Nodev(Nodev==-1)=5;
NodeRho(NodeRho==3)=-1;NodeRho(NodeRho==5)=3;NodeRho(NodeRho==-1)=5;

%% Initialize fields
if isfield(app,'S') || isfield(app,'nopert')
    rhoinit=@(X,Y) 0*(X>5).*(X<6).*(Y>2);
else
    rhoinit=@(X,Y) (X>5).*(X<6).*(Y>2);
end


pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) zeros(size(X));
vinit = @(X,Y) zeros(size(X));

base = size(Nodeu,1)-2;
UScale = base./repmat(base-sum((Nodeu==1)),base,1);
UScale(:,1)=[];
rho=[bcsDrho' ;bcsNrho' ;rhoinit(XP(idsP), YP(idsP))];
P = [bcsDP';bcsNP';  pinit(XP(idsP), YP(idsP))];
u = [bcsDu';bcsNu';  UScale(idsu)];
v = [bcsDv';bcsNv';  vinit(XV(idsv), YV(idsv))];
%Correct u-velocity to be divergence-free
% u(Nodeu(2:end-1,round(app.Nx/5)+2:end-1))=1/3;
