%This test is a test from the ISMP-HOM Intercomparison project
%Pattyn and Payne 2006

%TestE 
Lx=10; %in m
Ly=5000; %in m
nx=3; %number of nodes in x direction
ny=51;
md=model;
md=squaremesh(md,Lx,Ly,nx,ny);
md=geography(md,'',''); %ice sheet test
md=parameterize(md,'Square.par');
md=extrude(md,10,1);

md.gridondirichlet_diag(find(md.gridonbed))=1; %case 1: no slip
%md.gridondirichlet_diag(find(md.gridonbed & (md.y<2200 | md.y> 2500)))=1; %case 2: slip between 2200 and 2500m

%md=setelementstype(md,'pattyn','all');
md=setelementstype(md,'pattyn','all','stokes','all'); 

md.gridonstokes(find(md.x==0 | md.x==Lx))=0;
%Create MPCs to have periodic boundary conditions
posx=find(md.x==0);
posx2=find(md.x==max(md.x));

md.penalties=[posx,posx2];

%Compute the diagnostic
md=solve(md,'diagnostic','ice');

plotmodel(md,'data','vx','data','vy','data','vz','layer#all',md.numlayers)
