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

%TestB 
L=80000; %in m
nx=20; %numberof 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);

%Create dirichlet on the bed on top of the boundary
pos=find(md.gridonbed);
md.gridondirichlet_diag(pos)=1;

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

%Create dirichlet on the bed only
md.gridondirichlet_diag=zeros(md.numberofgrids,1);
pos=find(md.gridonbed);
md.gridondirichlet_diag(pos)=1;

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

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

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)
