source: issm/trunk/test/Validation/HydrostaticIceSheet/NoFront/runme.m@ 4868

% 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 of ice shelf
md=solve(md,'analysis_type',DiagnosticSolutionEnum);
vel_shelf=zeros(md.numberofgrids,1);
vel_shelf(md.results.DiagnosticSolution.Vel.index)=md.results.DiagnosticSolution.Vel.value;

%Compute solution for a 2d model of ice sheet
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=solve(md,'analysis_type',DiagnosticSolutionEnum);
vel_sheet=zeros(md.numberofgrids,1);
vel_sheet(md.results.DiagnosticSolution.Vel.index)=md.results.DiagnosticSolution.Vel.value;

%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 %]')