from numpy import *
from verbose import *
import scipy.io as matio
import InterpFromMeshToMesh2d as im
from   paterson import  *
from SetIceShelfBC import *
import inspect

#Start defining model parameters here
#Geometry
hmin = 300.
hmax = 1000.
ymin = min(md.mesh.y)
ymax = max(md.mesh.y)

md.geometry.thickness = hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)
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 
mat=matio.loadmat('../Data/SquareShelfConstrained.data')
#Reshape as Rank-1 arrays
x=reshape(mat['x'],(-1))
y=reshape(mat['y'],(-1))
vx=mat['vx']
vy=mat['vy']
#deal with 'F' oriented matlab matrices!
index=mat['index'].astype(float)
index=reshape(index.T,(len(index),3),order='F')

md.initialization.vx = im.InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
md.initialization.vy = im.InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
md.initialization.vz = zeros(md.mesh.numberofvertices)
md.initialization.pressure = zeros(md.mesh.numberofvertices)
#Materials
md.initialization.temperature = (273.-20)*ones(md.mesh.numberofvertices)
md.materials.rheology_B = paterson(md.initialization.temperature)
md.materials.rheology_n = 3.*ones(md.mesh.numberofelements)
#Surface mass balance and basal melting
md.surfaceforcings.mass_balance = 10.*ones(md.mesh.numberofvertices)
md.basalforcings.melting_rate = 5.*ones(md.mesh.numberofvertices)
#Friction
pos = nonzero(md.mask.elementonfloatingice)
md.friction.coefficient = 20.*ones(md.mesh.numberofvertices)
md.friction.coefficient[md.mesh.elements[pos,:].astype(int)-1] =0.
md.friction.p = ones(md.mesh.numberofelements)
md.friction.q = ones(md.mesh.numberofelements)
#Numerical parameters
md.diagnostic.viscosity_overshoot = 0.0
md.prognostic.stabilization = 1.
md.thermal.stabilization = 1.
md.verbose = verbose()
md.settings.waitonlock = 30.
md.diagnostic.restol = 0.05
md.diagnostic.reltol = 0.05
md.steadystate.reltol = 0.05
md.diagnostic.abstol = 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 = inspect.stack()[2][1].split('.')[0]

