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

#Start defining model parameters here
#Geometry
hmin=300.
hmax=1000.
ymin=min(md.mesh.y)
ymax=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.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
md.geometry.surface=md.geometry.bed+md.geometry.thickness

#Initial velocity and pressure
iVelF = netCDF4.Dataset('../Data/SquareShelf.nc','r')

x=reshape(iVelF.variables['x'][:],(-1))
y=reshape(iVelF.variables['y'][:],(-1))
vx=iVelF.variables['vx'][:]
vy=iVelF.variables['vy'][:]
index=iVelF.variables['index'][:].astype(int)

#dbg - begin
# #print 'vars in SquareShelf.nc:'
# #for v in iVelF.variables:
# #	print v
#dbg - end 

iVelF.close()

[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=zeros((md.mesh.numberofvertices,1))
md.initialization.pressure=zeros((md.mesh.numberofvertices,1))

#dbg - begin
#print '...vx:'
#print md.initialization.vx
#print '...vy:'
#print md.initialization.vy
##print '...vz:'
##print md.initialization.vz
##print '...pressure:'
##print md.initialization.pressure
#dbg - end 


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

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

#Numerical parameters
md.stressbalance.viscosity_overshoot = 0.3
md.masstransport.stabilization = 1.
md.thermal.stabilization = 1.
md.settings.waitonlock = 30
md.verbose=verbose()
md.stressbalance.restol = 0.10
md.steadystate.reltol = 0.02
md.stressbalance.reltol = 0.02
md.stressbalance.abstol = nan
md.timestepping.time_step = 1.
md.timestepping.final_time = 3.

#Boundary conditions:
# #md=SetIceShelfBC(md)
md=SetIceShelfBC(md,'../Exp/SquareFront.exp')

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