import numpy
import shutil
from model import *
from squaremesh import *
from setmask import *
from parameterize import *
from setflowequation import *
from EnumDefinitions import *
from solve import *
from MatlabFuncs import *
from PythonFuncs import *

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

printingflag=False

L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
results=[]
minvx=[]
maxvx=[]

for L in L_list:    #in m (3 times the desired length for BC problems)  
	nx=30    #number of nodes in x direction
	ny=30
	md=model()
	md=squaremesh(md,L,L,nx,ny)
	md=setmask(md,'','')    #ice sheet test
	md=parameterize(md,'../Par/ISMIPC.py')
	md.extrude(10,1.)

	md=setflowequation(md,'pattyn','all')

	#Create MPCs to have periodic boundary conditions
	md.diagnostic.spcvx=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
	md.diagnostic.spcvy=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
	md.diagnostic.spcvz=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))

	posx=numpy.nonzero(logical_and_n(md.mesh.x==0.,md.mesh.y!=0.,md.mesh.y!=L))[0]
	posx2=numpy.nonzero(logical_and_n(md.mesh.x==L,md.mesh.y!=0.,md.mesh.y!=L))[0]

	posy=numpy.nonzero(logical_and_n(md.mesh.y==0.,md.mesh.x!=0.,md.mesh.x!=L))[0]    #Don't take the same nodes two times
	posy2=numpy.nonzero(logical_and_n(md.mesh.y==L,md.mesh.x!=0.,md.mesh.x!=L))[0]

	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))

	#Add spc on the corners
	pos=numpy.nonzero(logical_and_n(numpy.logical_or(md.mesh.x==0.,md.mesh.x==L),numpy.logical_or(md.mesh.y==0.,md.mesh.y==L),md.mesh.vertexonbed))
	md.diagnostic.spcvx[pos]=0.
	md.diagnostic.spcvy[pos]=0.
	if   (L==5000.):
		md.diagnostic.spcvx[pos]=15.66
		md.diagnostic.spcvy[pos]=-0.1967
	elif (L==10000.):
		md.diagnostic.spcvx[pos]=16.04
		md.diagnostic.spcvy[pos]=-0.1977
	elif (L==20000.):
		md.diagnostic.spcvx[pos]=16.53
		md.diagnostic.spcvy[pos]=-1.27
	elif (L==40000.):
		md.diagnostic.spcvx[pos]=17.23
		md.diagnostic.spcvy[pos]=-3.17
	elif (L==80000.):
		md.diagnostic.spcvx[pos]=16.68
		md.diagnostic.spcvy[pos]=-2.69
	elif (L==160000.):
		md.diagnostic.spcvx[pos]=16.03
		md.diagnostic.spcvy[pos]=-1.27
	
	#Spc the bed at zero for vz
	pos=numpy.nonzero(md.mesh.vertexonbed)
	md.diagnostic.spcvz[pos]=0.

	#Compute the diagnostic
	md.cluster=generic('name',oshostname(),'np',8)
	md=solve(md,DiagnosticSolutionEnum())

	#Plot the results and save them
	vx=md.results.DiagnosticSolution.Vx
	vy=md.results.DiagnosticSolution.Vy
	vz=md.results.DiagnosticSolution.Vz
	results.append(md.results.DiagnosticSolution)
	minvx.append(numpy.min(vx[-md.mesh.numberofvertices2d:]))
	maxvx.append(numpy.max(vx[-md.mesh.numberofvertices2d:]))

	#Now plot vx, vy, vz and vx on a cross section
#	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
	if printingflag:
		pass
#		set(gcf,'Color','w')
#		printmodel(['ismipcpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
#		shutil.move("ismipcpattynvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
#	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
	if printingflag:
		pass
#		set(gcf,'Color','w')
#		printmodel(['ismipcpattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
#		shutil.move("ismipcpattynvy%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
#	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
	if printingflag:
		pass
#		set(gcf,'Color','w')
#		printmodel(['ismipcpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
#		shutil.move("ismipcpattynvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')

	if   (L==5000.):
		pass
#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
#			'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',5)
	elif (L==10000.):
		pass
#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
#			'resolution',[10 10],'ylim',[13 18],'xlim',[0 10000],'title','','xlabel','')
	elif (L==20000.):
		pass
#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
#			'resolution',[10 10],'ylim',[14 22],'xlim',[0 20000],'title','','xlabel','')
	elif (L==40000.):
		pass
#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
#			'resolution',[10 10],'ylim',[10 40],'xlim',[0 40000],'title','','xlabel','')
	elif (L==80000.):
		pass
#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
#			'resolution',[10 10],'ylim',[0 80],'xlim',[0 80000],'title','','xlabel','')
	elif (L==160000.):
		pass
#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
#			'resolution',[10 10],'ylim',[0 200],'xlim',[0 160000],'title','','xlabel','')
	if printingflag:
		pass
#		set(gcf,'Color','w')
#		printmodel(['ismipcpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
#		shutil.move("ismipcpattynvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')

#Now plot the min and max values of vx for each size of the square
#plot([5 10 20 40 80 160],minvx);ylim([4 18]);xlim([0 160])
if printingflag:
	pass
#	set(gcf,'Color','w')
#	printmodel('ismipcpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
#	shutil.move('ismipcpattynminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
#plot([5 10 20 40 80 160],maxvx);ylim([0 200]); xlim([0 160])
if printingflag:
	pass
#	set(gcf,'Color','w')
#	printmodel('ismipcpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
#	shutil.move('ismipcpattynmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')

#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-08,1e-07,1e-07,\
	1e-09,1e-07,1e-07,\
	1e-09,1e-09,1e-07,\
	1e-09,1e-09,1e-08,\
	1e-09,1e-08,1e-08,\
	1e-09,1e-08,1e-08,\
]
field_values=[]
for result in results:
	field_values=field_values+[\
		result.Vx,\
		result.Vy,\
		result.Vz,\
		]
