#Test Name: ISMIPCFS
import numpy
from model import *
from triangle 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
"""

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

for L in L_list:
	md=triangle(model(),"../Exp/Square_%d.exp" % L,L/10.)    #size 3*L 
	md=setmask(md,'','')    #ice sheet test
	md=parameterize(md,'../Par/ISMIPC.py')
	md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/L)*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/L)))
	md.extrude(10,1.)

	#Add spc on the borders
	pos=numpy.nonzero(logical_or_n(md.mesh.x==0.,md.mesh.x==numpy.max(md.mesh.x),md.mesh.y==0.,md.mesh.y==numpy.max(md.mesh.y)))
	md.stressbalance.spcvx[pos]=0.
	md.stressbalance.spcvy[pos]=0.
	if   (L==5000.):
		md.stressbalance.spcvx[pos]=15.66
		md.stressbalance.spcvy[pos]=-0.1967
	elif (L==10000.):
		md.stressbalance.spcvx[pos]=16.04
		md.stressbalance.spcvy[pos]=-0.1977
	elif (L==20000.):
		md.stressbalance.spcvx[pos]=16.53
		md.stressbalance.spcvy[pos]=-1.27
	elif (L==40000.):
		md.stressbalance.spcvx[pos]=17.23
		md.stressbalance.spcvy[pos]=-3.17
	elif (L==80000.):
		md.stressbalance.spcvx[pos]=16.68
		md.stressbalance.spcvy[pos]=-2.69
	elif (L==160000.):
		md.stressbalance.spcvx[pos]=16.03
		md.stressbalance.spcvy[pos]=-1.27

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

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

	#Plot the results and save them
	vx=md.results.StressbalanceSolution.Vx
	vy=md.results.StressbalanceSolution.Vy
	vz=md.results.StressbalanceSolution.Vz
	results.append(md.results.StressbalanceSolution)

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

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