% Stokes Flow forcing function test case

%% Forcing Function
app.nu = 1;
Fu = @(time,p) ... 
              pi.*cos(time+1).*sin(2.*pi.*p(:,2)).*sin(pi.*p(:,1)).^2-...
        2.*pi.^3.*sin(time+1).*sin(2.*pi.*p(:,2)).*cos(pi.*p(:,1)).^2+...
        6.*pi.^3.*sin(time+1).*sin(2.*pi.*p(:,2)).*sin(pi.*p(:,1)).^2-...
                  sin(time+1).*sin(pi.*p(:,1)).*pi.*sin(pi.*p(:,2));
Fv = @(time,p) ...
             -pi.*cos(time+1).*sin(2.*pi.*p(:,1)).*sin(pi.*p(:,2)).^2-...
        6.*pi.^3.*sin(time+1).*sin(2.*pi.*p(:,1)).*sin(pi.*p(:,2)).^2+...
        2.*pi.^3.*sin(time+1).*sin(2.*pi.*p(:,1)).*cos(pi.*p(:,2)).^2+...
                  sin(time+1).*cos(pi.*p(:,1)).*cos(pi.*p(:,2)).*pi;

%% Exact Solution (used in plotting function)
uexact = @(time,x)  pi*sin(time+1).*(sin(2*pi*x(:,2)).*(sin(pi*x(:,1))).^2);
vexact = @(time,x) -pi*sin(time+1).*(sin(2*pi*x(:,1)).*(sin(pi*x(:,2))).^2);
pexact = @(time,x)     sin(time+1).* cos(  pi*x(:,1)).* sin(pi*x(:,2)) ;

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

%% Set boundary conditions
bcsDP= [];         bcsNP= [0 0 0 0];        
bcsDu= [0 0 0 0];  bcsNu= [];        
bcsDv= [0 0 0 0];  bcsNv= [];        

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

%% create geometry mask
Nbcs = length(bcsDP) + length(bcsNP);
Mask = zeros(app.Ny, app.Nx);
%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,[1 1 1 1],1);


%% Initialize fields
%create vectors of coordinates for p, u, and v control volumes
xp = linspace(dx/2/a-0.5, 0.5-dx/2/a, app.Nx  )*a;
yp = fliplr(linspace(dy/2/b-0.5, 0.5-dy/2/b, app.Ny  ))*b;
xu = linspace( dx/a-0.5   ,0.5-dx/a , app.Nx-1)*a;
yu = fliplr(linspace(dy/2/b-0.5, 0.5-dy/2/b, app.Ny  ))*b;
xv = linspace(dx/2/a-0.5, 0.5-dx/2/a, app.Nx  )*a;
yv = fliplr(linspace(dy/b-0.5    , 0.5-dy/b, app.Ny-1))*b;

%create matrix of x and y coordinates
[XP YP] = meshgrid(xp,yp);
[XU YU] = meshgrid(xu,yu);
[XV YV] = meshgrid(xv,yv);

%Initialize vectors
P = [bcsDP';bcsNP';pexact(0,[XP(idsP) YP(idsP)])];
u = [bcsDu';bcsNu';uexact(0,[XU(idsu) YU(idsu)])];
v = [bcsDv';bcsNv';vexact(0,[XV(idsv) YV(idsv)])];
Pex = [bcsDP';bcsNP';pexact(0,[XP(idsP) YP(idsP)])];
uex = [bcsDu';bcsNu';uexact(0,[XU(idsu) YU(idsu)])];
vex = [bcsDv';bcsNv';vexact(0,[XV(idsv) YV(idsv)])];
