"""
== == == == == == == == == == == == == == == == == == ==
Auto generated python script for ISSM:   test218.m
Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
== == == == == == == == == == == == == == == == == == ==

Matlab script conversion into python
translateToPy.py Author: Michael Pellegrin
translateToPy.py Date: 09/24/12
== == == == == == == == == == == == == == == == == == ==
"""

from MatlabFuncs import *
from model import *
from EnumDefinitions import *
from numpy import *
from setmask import *
from parameterize import *
from setflowequation import *
from solve import *

md=squaremesh(model(),1000000,1000000,5,5)
md=setmask(md,'all','')
md=parameterize(md,'../Par/SquareShelf.py')
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)


# redo the parameter file for this special shelf. 

# constant thickness, constrained (vy=0) flow into an icefront, 

# from 0 m/yr at the grounding line.



# needed later

ymin=min(md.mesh.y)
ymax=max(md.mesh.y)
xmin=min(md.mesh.x)
xmax=max(md.mesh.x)


di=md.materials.rho_ice/md.materials.rho_water


h=1000
md.geometry.thickness=h*ones(md.mesh.numberofvertices,1)
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

md.initialization.vx=zeros(md.mesh.numberofvertices,1)
md.initialization.vy=zeros(md.mesh.numberofvertices,1)
md.initialization.vz=zeros(md.mesh.numberofvertices,1)
md.initialization.pressure=zeros(md.mesh.numberofvertices,1)


# 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)


# Boundary conditions:

md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1)
md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1)
md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1)


# constrain flanks to 0 normal velocity

pos=numpy.nonzero(md.mesh.x==xmin | md.mesh.x==xmax)
md.diagnostic.spcvx(pos)=0
md.diagnostic.spcvz(pos)=NaN


# constrain grounding line to 0 velocity

pos=numpy.nonzero(md.mesh.y==ymin)
md.diagnostic.spcvx(pos)=0
md.diagnostic.spcvy(pos)=0


# icefront

nodeonicefront=zeros(md.mesh.numberofvertices,1)
pos=numpy.nonzero(md.mesh.y==ymax)
nodeonicefront(pos)=1
pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)))
diagnostic.icefront=md.mesh.segments(pos,:)
diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]
md.diagnostic.icefront=diagnostic.icefront


# partitioning

md.qmu.numberofpartitions=md.mesh.numberofvertices
md=partitioner(md,'package','linear','npart',md.qmu.numberofpartitions)
md.qmu.pytition=md.qmu.pytition-1


# Dakota options

# variables

md.qmu.variables.rheology_B=normal_uncertain('scaled_MaterialsRheologyB',1,.05)


# responses

md.qmu.responses.MaxVel=response_function('MaxVel',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999])


# method

md.qmu.method     =dakota_method('nond_l')


# parameters

md.qmu.pyams.direct=true
md.qmu.pyams.analysis_driver='diagnostic'
md.qmu.pyams.evaluation_concurrency=1
md.qmu.pyams.interval_type='forward'


# imperative! 

md.diagnostic.reltol=10^-10#tighten for qmu analysese
md.qmu.isdakota=1


# solve

md=solve(md,DiagnosticSolutionEnum(),'overwrite','y')


# Fields and tolerances to track changes

md.qmu.results=md.results.dakota
md.results.dakota.importancefactors=importancefactors(md,'scaled_MaterialsRheologyB','MaxVel')
field_names     =['importancefactors']
field_tolerances=[1e-10]
field_values=[\
md.results.dakota.importancefactors,\
]
