source: issm/trunk/test/Par/ISMIPE.py@ 21341

import numpy
from arch import *
from SetIceSheetBC import SetIceSheetBC

#Ok, start defining model parameters here

print "      creating thickness"
data = numpy.array(archread('../Data/ISMIPE.arch','data'));
md.geometry.surface=numpy.zeros((md.mesh.numberofvertices,1))
md.geometry.base=numpy.zeros((md.mesh.numberofvertices,1))
for i in xrange(0,md.mesh.numberofvertices):
        y=md.mesh.y[i]
        point1=numpy.floor(y/100.)
        point2=numpy.minimum(point1+1,50)
        coeff=(y-(point1-1.)*100.)/100.
        md.geometry.base[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1]
        md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2]
md.geometry.thickness=md.geometry.surface-md.geometry.base
md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01
md.geometry.base=md.geometry.surface-md.geometry.thickness

print "      creating drag"
md.friction.coefficient=numpy.zeros((md.mesh.numberofvertices,1))
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 stressbalance model:"
#Create node on boundary first (because we can not use mesh)
md=SetIceSheetBC(md)