% with mesh refinement it is rquired to manually change the nodal
% location at which source is applied
% M=10x10,111; M=20X20,421;M=40x40, 1641; M=80x80,6481;

clear all;clc;

global X1 X2 elcon;

%Create geometry file
el_X = 10;
el_Y = 10;
left_X = 0;
right_X = 1;
bottom_Y = 0;
top_Y = 1;
%filename='mesh_10x10.geom';
GenRectMesh('filename',el_X,el_Y,left_X,right_X,bottom_Y,top_Y);
%*********** Geometry ***********************************
%Read geometry from file.
fid = fopen('filename');
var1 = fscanf(fid,'%s',1);
var2 = fscanf(fid,'%s',1);
if( var1 == 'X1')
	size_X1X2 = fscanf(fid,'%i',[1 2]);
	X1X2 = (fscanf(fid,'%f',[size_X1X2(2) size_X1X2(1)]))';
	X1 = X1X2(:,1);
	X2 = X1X2(:,2);
else
fprinf('Error:The domain variable, X1, has not been identified.')
end
var = fscanf(fid,'%s',1);
if( var == 'elcon')
	size_elcon = fscanf(fid,'%i',[1 2]);
	elcon = (fscanf(fid,'%i',[size_elcon(2) size_elcon(1)]))';
else
fprinf('Error:The element connectivity, elcon, has not been identified.')
end
while(~feof(fid))
var = fscanf(fid,'%s',1);
if(strcmp(var,'XL'))
size_bc = fscanf(fid,'%i',1);
XL = (fscanf(fid,'%i',size_bc))';
end
if(strcmp(var,'YT'))
size_bc = fscanf(fid,'%i',1);
YT = (fscanf(fid,'%i',size_bc))';
end
if(strcmp(var,'YB'))
size_bc = fscanf(fid,'%i',1);
YB = (fscanf(fid,'%i',size_bc))';
end
if(strcmp(var,'XR'))
size_bc = fscanf(fid,'%i',1);
XR = (fscanf(fid,'%i',size_bc))';
end
end
fclose(fid);

%*********** Boundary nodes ***********************
SRC = [111];

%*********** Physical Parameters ***********************
Pa = 0.01;
%inlet = 1;
walls = 0;
u = 1;
v = 1;
src = 1;

%*********** Initial Conditions ***********************
[Xnodes Ynodes] = size(X1);
q = ones(Xnodes.*Ynodes,1)*0;

%** Enforce Boundary Conditions on Initial Condition for Set 0 

for ii=1:1:length(XL)
q(XL(ii)) = +walls;
end
for ii=1:1:length(YT)
q(YT(ii)) = +walls;
end
for ii=1:1:length(YB)
q(YB(ii)) = +walls;
end
for ii=1:1:length(XR)
q(XR(ii)) = +walls;
end
for ii=1:1:length(SRC)
q(SRC(ii)) = +src;
end

%********************* Equations for Set 0*****************
res_0 = +asres2D(u,[],[],0,'B201L',q)...
+asres2D(v,[],[],0,'B202L',q)...
+asres2D(Pa,[],[],-1,'B211L',q)...
+asres2D(Pa,[],[],-1,'B212L',q)...
+asres2D(Pa,[],[],-1,'B221L',q)...
+asres2D(Pa,[],[],-1,'B222L',q);

%Apply Dirichlet boundary condition to boundary XL.
for ii=1:1:length(XL)
res_0(XL(ii)) = 0;
end

%Apply Dirichlet boundary condition to boundary YB.
for ii=1:1:length(YB)
res_0(YB(ii)) = 0;
end

%Apply Dirichlet boundary condition to boundary SRC.
for ii=1:1:length(SRC)
res_0(SRC(ii)) = 0;
end

%***************** Jacobians for Set 0**********
%Equation q by q
jac_0 = +asjac2D(u,[],[],0,'B201L',[])...
+asjac2D(v,[],[],0,'B202L',[])...
+asjac2D(Pa,[],[],-1,'B211L',[])...
+asjac2D(Pa,[],[],-1,'B212L',[])...
+asjac2D(Pa,[],[],-1,'B221L',[])...
+asjac2D(Pa,[],[],-1,'B222L',[]);

%Apply Dirichlet boundary condition to boundary XL.
for ii=1:1:length(XL)
jac_0(XL(ii),:) = 0;
jac_0(:,XL(ii)) = 0;
jac_0(XL(ii),XL(ii)) = 1;
end

%Apply Dirichlet boundary condition to boundary YB.
for ii=1:1:length(YB)
jac_0(YB(ii),:) = 0;
jac_0(:,YB(ii)) = 0;
jac_0(YB(ii),YB(ii)) = 1;
end

%Apply Dirichlet boundary condition to boundary SRC.
for ii=1:1:length(SRC)
jac_0(SRC(ii),:) = 0;
jac_0(:,SRC(ii)) = 0;
jac_0(SRC(ii),SRC(ii)) = 1;
end

%***************** Solver ********************************************
dq = jac_0\(-res_0);
q = q + dq
format long
enorm=q'*jac_0*q
%
%Post Processing
eltall = 10
elwide = 10;
X_new = zeros(eltall+1,elwide+1);
Y_new = zeros(eltall+1,elwide+1);
Q_new = zeros(eltall+1,elwide+1);
X_new(:) = X1;
Y_new(:) = X2;
Q_new(:) = q;
figure(1)
surf(X_new,Y_new,Q_new)
xlabel('x', 'fontsize',12), ylabel('y','fontsize',12), zlabel('q','fontsize',12)
title('GWS, Pe=10, Mesh=', 'Fontsize',14)
colorbar('vert')
