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

%TestC 
L=80000; %in m
nx=20; %number of nodes in x direction
ny=20;
md=model;
md=squaremesh(md,L,L,nx,ny);
md=geography(md,'',''); %ice sheet test
md=parameterize(md,'Square.par');
md=extrude(md,10,1);

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

%Add dirichlet values on the corners (we need dirichlet somewhere
md.gridondirichlet_diag=zeros(md.numberofgrids,1);
pos=find((md.x==0 |md.x==L) & (md.y==0 | md.y==L));
%pos=find((md.x==0 |md.x==L | md.y==0 | md.y==L)); %dirichlet on all the boundaries
md.gridondirichlet_diag(pos)=1;

%Create MPCs to have periodic boundary conditions
posx=find(md.x==0);
posx2=find(md.x==L);

posy=find(md.y==0 & md.x~=0 & md.x~=L); %Don't take the same grids two times
posy2=find(md.y==L & md.x~=0 & md.x~=L);

md.penalties=[posx,posx2;posy,posy2];

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

plotmodel(md,'data','vx','data','vy','data','vz','layer#all',md.numlayers,'contourlevels',{10.878,16.317,21.756,43.51,59.829},'contourlevels',{-2.421,-1.614,-0.807,3.228},'contourlevels',{-2.48,-1.86,0,0.62,1.24,1.86,2.48})
