function data = initialize()

%% THIS DATA MAY BE SET BY THE USER*********************************

%% Constant for splitting to processors
% Number of MATLAB labs available to the program
P = 4;


%% Data input for initial conditions

% 2^N is the number of grid lines in each direction in the original
% square domain
N = 5;

% Coefficient of rotation: analogous to a coriolis force
b = 1.5;

% Maximum time before computation stops
Tmax = 10;

% Bottom friction coefficient: a=0 corresponds to no friction
a = 0;

% Initial x-coordinate of the fluid center of mass
xc = 1/sqrt(2);

% Initial y-coordinate of the fluid center of mass
yc = 0;

% Eddy viscosity
mu= .05;

% Ratio of dt/dx
a1 = .01;
 
%% END USER INPUT SECTION****************************************
       
%% Setting up domain and initial values of X, Y, u, and v

% P = 2^(p-1) therefore p = lg(P) + 1
p = log2(P) + 1;

m=2^N;

x=-.5:1/m:.5;
y=ones(size(1:m+1));

X=(y'*x)';
Y=(x'*y)';

dx=1/m ;
H=1;

dt=a1*dx;

d=0;
s=0;
w=0; 
T=0;
u=((d+w)*X+(w+T)*Y)/2;
v=((T-w)*X+(d-s)*Y)/2;
X=X+xc;
Y=Y+yc;

%% Calculating area of the fluid elements

DxX1(1:m,1:2:m-1)=X(2:m+1,1:2:m-1)-X(1:m,1:2:m-1);
DxY1(1:m,1:2:m-1)=Y(2:m+1,1:2:m-1)-Y(1:m,1:2:m-1);
DxX1(1:m,2:2:m)=X(2:m+1,3:2:m+1)-X(1:m,3:2:m+1);
DxY1(1:m,2:2:m)=Y(2:m+1,3:2:m+1)-Y(1:m,3:2:m+1);

DxX2(1:m,1:2:m-1)=X(2:m+1,2:2:m)-X(1:m,2:2:m);
DxY2(1:m,1:2:m-1)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m);
DxX2(1:m,2:2:m)=X(2:m+1,2:2:m)-X(1:m,2:2:m);
DxY2(1:m,2:2:m)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m);

DyX1(2:2:m,1:m)=X(2:2:m,2:m+1)-X(2:2:m,1:m); 
DyY1(2:2:m,1:m)=Y(2:2:m,2:m+1)-Y(2:2:m,1:m); 
DyX1(1:2:m-1,1:m)=X(2:2:m,2:m+1)-X(2:2:m,1:m);
DyY1(1:2:m-1,1:m)=Y(2:2:m,2:m+1)-Y(2:2:m,1:m);

DyX2(1:2:m-1,1:m)=X(1:2:m-1,2:m+1)-X(1:2:m-1,1:m);
DyY2(1:2:m-1,1:m)=Y(1:2:m-1,2:m+1)-Y(1:2:m-1,1:m);
DyX2(2:2:m,1:m)=X(3:2:m+1,2:m+1)-X(3:2:m+1,1:m);
DyY2(2:2:m,1:m)=Y(3:2:m+1,2:m+1)-Y(3:2:m+1,1:m);


% area of type 1 and type 2 triangles assigned to each node
A1=(DxX1.*DyY1-DyX1.*DxY1)/2;
A2=(DxX2.*DyY2-DyX2.*DxY2)/2;


% mass in each fluid element
m1=H*A1;
m2=H*A2;


%%  Mass Matrix Computation

mass=sum(sum(sum(A1))+sum(sum(A2)))/m/m;
M=zeros(m+1);
M(1:m+1,1:m+1)=mass;

% Each processor calculates its own mass matrix
M(1,2:m)=mass/2;
M(m+1,2:m)=mass/2;
M(2:m,1)=mass/2;
M(2:m,m+1)=mass/2;
M(1,1)=mass/4;
M(m+1,m+1)=mass/4;
M(m+1,1)=mass/4;
M(1,m+1)=mass/4;

data.MM=sum(sum(M));

n = 2^(N-p+1);


data.m = m;
data.n = n;
data.b = b;
data.a = a;
data.mu = mu;
data.dt = dt;
data.dx = dx;
data.Tmax = Tmax;
data.P = P;
data.xc = xc;
data.yc = yc;

data.X = X;
data.Y = Y;
data.u = u;
data.v = v;
data.M = M;
data.m1 = m1;
data.m2 = m2;


    
end