import os.path
from arch import *
import numpy
import inspect
from verbose import verbose
from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
from paterson import paterson
from SetIceShelfBC import SetIceShelfBC

#Start defining model parameters here
#Geometry
hmin=300.
hmax=1000.
ymin=numpy.min(md.mesh.y)
ymax=numpy.max(md.mesh.y)
xmin=min(md.mesh.x)
xmax=max(md.mesh.x)
md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x.reshape(-1,1)-xmin)/(xmax-xmin)
md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
md.geometry.surface=md.geometry.base+md.geometry.thickness
md.geometry.bed=md.geometry.base-10;

#Initial velocity 
#x         = numpy.reshape(numpy.array(archread('../Data/SquareShelfConstrained.arch','x')),(-1))
#y         = numpy.reshape(numpy.array(archread('../Data/SquareShelfConstrained.arch','y')),(-1))
x         = numpy.array(archread('../Data/SquareShelfConstrained.arch','x'))
y         = numpy.array(archread('../Data/SquareShelfConstrained.arch','y'))
vx        = numpy.array(archread('../Data/SquareShelfConstrained.arch','vx'))
vy        = numpy.array(archread('../Data/SquareShelfConstrained.arch','vy'))
index     = numpy.array(archread('../Data/SquareShelfConstrained.arch','index').astype(int))

[md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
[md.initialization.vy]=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)
md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))

#Materials
md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
md.materials.rheology_B=paterson(md.initialization.temperature)
md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))

#Surface mass balance and basal melting
md.smb.mass_balance=10.*numpy.ones((md.mesh.numberofvertices,1))
md.basalforcings.groundedice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))
md.basalforcings.floatingice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))

#Friction
md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
md.friction.p=numpy.ones((md.mesh.numberofelements,1))
md.friction.q=numpy.ones((md.mesh.numberofelements,1))

#Numerical parameters
md.stressbalance.viscosity_overshoot=0.0
md.masstransport.stabilization=1.
md.thermal.stabilization=1.
md.verbose = verbose(0)
md.settings.waitonlock=30
md.stressbalance.restol=0.05
md.stressbalance.reltol=0.05
md.steadystate.reltol=0.05
md.stressbalance.abstol=float('nan')
md.timestepping.time_step=1.
md.timestepping.final_time=3.

#Deal with boundary conditions:
md = SetIceShelfBC(md)

#Change name so that no tests have the same name
if len(inspect.stack()) > 2:
	md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
