import numpy
from SetIceSheetBC import *

#Ok, start defining model parameters here

print "      creating thickness"
md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(0.5*numpy.pi/180.)
md.geometry.bed=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))
md.geometry.thickness=md.geometry.surface-md.geometry.bed

print "      creating drag"
md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))    #q=1.
#Take care of iceshelves: no basal drag
pos=numpy.nonzero(md.mask.elementonfloatingice)[0]
md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
md.friction.p=numpy.ones((md.mesh.numberofelements,1))
md.friction.q=numpy.ones((md.mesh.numberofelements,1))

print "      creating flow law parameter"
md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))

print "      boundary conditions for diagnostic model"
#Create node on boundary first (because we cannot use mesh)
md=SetIceSheetBC(md)
