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

%TestD 
L=40000; %in m
nx=45; %numberof nodes in x direction
ny=45;
md=model;
md=squaremesh(md,L,L,nx,ny);
md.z=zeros(md.numberofgrids,1);
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'); 

%No dirichlet in this model, just periodic boundary conditions
md.gridondirichlet_diag=zeros(md.numberofgrids,1);
pos=find((md.x==0 |md.x==80000 ) & (md.y==0 | md.y==80000));
md.gridondirichlet_diag(pos)=1;

%We need one grd on dirichlet (in the middle)
%md.gridondirichlet_diag(33)=1;
%md.dirichletvalues_diag(1,:)=[18/md.yts -3/md.yts];

%%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];

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

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