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 *

"""
EISMINT benchmark experiment A
"""

numlayers=8
resolution=50000.

#To begin with the numerical model
md=triangle(model(),'../Exp/SquareEISMINT750000.exp',resolution)
md=setmask(md,'','')
md=parameterize(md,'../Par/RoundSheetEISMINT.py')

#We extrude the model to have a 3d model
md.extrude(numlayers,1.)
md=setflowequation(md,'SIA','all')

#Spc the nodes on the bed
pos=numpy.nonzero(md.mesh.vertexonbed)
md.stressbalance.spcvx[pos]=0.
md.stressbalance.spcvy[pos]=0.
md.stressbalance.spcvz[pos]=0.

#Adapt the time steps to the resolution
md.timestepping.time_step=15.
md.settings.output_frequency=500
md.timestepping.final_time=30000.
md.masstransport.stabilization=1
md.thermal.stabilization=1

#Now we can solve the problem 
md.cluster=generic('name',oshostname(),'np',8)
md=solve(md,TransientSolutionEnum())

#Fields and tolerances to track changes
field_names     =['Vx','Vy','Vz','Vel','Pressure','Thickness','Base','Surface','Temperature','BasalforcingsMeltingRate']
field_tolerances=[1e-08,1e-08,1e-07,1e-08,1e-08,1e-08,1e-08,1e-08,1e-07,1e-07]
field_values=[\
	md.results.TransientSolution[-1].Vx,\
	md.results.TransientSolution[-1].Vy,\
	md.results.TransientSolution[-1].Vz,\
	md.results.TransientSolution[-1].Vel,\
	md.results.TransientSolution[-1].Pressure,\
	md.results.TransientSolution[-1].Thickness,\
	md.results.TransientSolution[-1].Base,\
	md.results.TransientSolution[-1].Surface,\
	md.results.TransientSolution[-1].Temperature,\
	md.results.TransientSolution[-1].BasalforcingsMeltingRate,\
	]
