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

L_list={5000,10000,20000,40000,80000,160000};
results={};

for i=1:length(L_list),
	L=3*L_list{i};
	nx=60; %numberof nodes in x direction
	ny=60;
	md=model;
	md=squaremesh(md,L,L,nx,ny);
	md=geography(md,'',''); %ice sheet test
	md=parameterize(md,'../Par/ISMIPD.par');
	md=extrude(md,6,1);

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

	%We need one grd on dirichlet: the 4 corners are set to zero
	%md.spcvelocity=zeros(md.numberofgrids,6);
	%pos=find((md.x==0 | md.x==max(md.x)) & (md.y==0 | md.y==max(md.y)));
	%md.spcvelocity(pos,1:3)=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,'analysis_type',DiagnosticSolutionEnum);

	%Plot the results and save them
	vx=PatchToVec(md.results.DiagnosticSolution.Vx);
	vy=PatchToVec(md.results.DiagnosticSolution.Vy);
	vz=PatchToVec(md.results.DiagnosticSolution.Vz);
	results{i}=md.results.DiagnosticSolution;

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

%Fields and tolerances to track changes
field_names     ={...
	'Vx5km','Vy5km','Vz5km',...
	'Vx10km','Vy10km','Vz10km',...
	'Vx20km','Vy20km','Vz20km',...
	'Vx40km','Vy40km','Vz40km',...
	'Vx80km','Vy80km','Vz80km',...
	'Vx160km','Vy160km','Vz160km'
}
field_tolerances={...
	1e-13,1e-13,1e-13,...
	1e-13,1e-13,1e-13,...
	1e-13,1e-13,1e-13,...
	1e-13,1e-13,1e-13,...
	1e-13,1e-13,1e-13,...
	1e-13,1e-13,1e-13,...
};
field_values={};
for i=1:6,
	result=results{i};
	field_values={field_values{:},...
		PatchToVec(result.Vx),...
		PatchToVec(result.Vy),...
		PatchToVec(result.Vz),...
		};
end
