% This file can be run to compare an icesheet in hydrostatic equilibrium with an iceshelf.
% This test deals with a  mesh includind a iceshelf/icesheet with no ice front. The geometry 
% is square. Just run this file in Matlab, with a properly setup Ice code. 
% The result should be exactly the same since the ice sheet is in hydrostatic equilibrium
% so the friction at the base should be zero.

md=model;
md=mesh(md,'DomainOutline.exp',50000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'Macayeal','all');

%Compute solution for a 2d model
md=solve(md,'analysis_type','diagnostic');
vel_shelf=md.results.diagnostic.vel;

%Compute solution for a 3d model
md=geography(md,'','');
md=solve(md,'analysis_type','diagnostic');
vel_sheet=md.results.diagnostic.vel;

%Plot of the velocity from the ice sheet and he ice shelf model
plotmodel(md,'data',vel_shelf,'title','Velocity of ice shelf','data',vel_sheet,'title','Velocity of ice sheet','data',vel_sheet-vel_shelf,'title','Absolute difference between ice sheet and ice shelf velocities','data',abs((vel_sheet-vel_shelf)./vel_shelf)*100,'title','Relative difference between ice sheet and ice shelf velocities [in %]')