% This file can be run to compare MacAyeal and Pattyn's model ie a 2d and a 3d model.
% This test deals with a  mesh includind an iceshelf ans an icesheet. The geometry 
% is square. Just run this file in Matlab, with a properly setup Ice code. 

%Compute solution for a Pattyn model
md=model;
md=mesh(md,'DomainOutline.exp',50000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,5,3);
md=extrude(md,'Pattyn','all');

md=solve(md,'diagnostic','ice');
vel_pattyn=md.vel;
vx_pattyn=md.vx;
vy_pattyn=md.vy;
vz_pattyn=md.vz;
pressure_pattyn=md.pressure;

mdpattyn=md;
save modelpattyn mdpattyn;

%Compute solution for a Stokes model
md=model;
md=mesh(md,'DomainOutline.exp',50000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,5,3);
md=extrude(md,'Pattyn','all','Stokes','all');

md=solve(md,'diagnostic','ice');
vel_stokes=md.vel;
vx_stokes=md.vx;
vy_stokes=md.vy;
vz_stokes=md.vz;
pressure_stokes=md.pressure;

mdstokes=md;
save modelstokes mdstokes;

error('debug')



%Calculate the average velocity on each grid
vel_3d=zeros(md.numberofgrids2d,1);
grid_vel=0;

for i=1:md.numberofgrids2d
	for j=1:(md.numlayers-1)
		grid_vel=grid_vel+1/(2*(md.numlayers-1))*(md.vel(i+j*md.numberofgrids2d,1)+md.vel(i+(j-1)*md.numberofgrids2d,1));
	end
	vel_3d(i,1)=grid_vel;
	grid_vel=0;
end

vel_diff=(vel_2d-vel_3d)./vel_2d;
vel_diff(find(vel_2d==vel_3d))=0;

%Plot of the velocity from the stokes and pattyn model 
figure(1)
subplot(2,2,1)
p1=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_pattyn,'FaceColor','interp','EdgeColor','none');
title('MacAyeal model [m/yr]','FontSize',14,'FontWeight','bold')
colorbar; 
   
subplot(2,2,2)
p2=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_stokes,'FaceColor','interp','EdgeColor','none');
title('Pattyn model [m/yr]','FontSize',14,'FontWeight','bold')
colorbar; 

subplot(2,2,3)
p3=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
(vel_pattyn-vel_stokes),'FaceColor','interp','EdgeColor','none');
title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
colorbar;

subplot(2,2,4)
p3=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
(vel_pattyn-vel_stokes)*100,'FaceColor','interp','EdgeColor','none');
title('Absolute misfit [m/yr]','FontSize',14,'FontWeight','bold')
colorbar;

