import numpy
from model import *
from triangle import *
from setmask import *
from parameterize import *
from matdamageice import *
from paterson import *
from setflowequation import *
from EnumDefinitions import *
from solve import *
from MatlabFuncs import *

md=triangle(model(),'../Exp/Square.exp',150000.)
md=setmask(md,'all','')
md=parameterize(md,'../Par/SquareShelf.py')
md.materials=matdamageice()
md.materials.rheology_B=paterson(md.initialization.temperature)
md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
md.materials.rheology_Z=0.5*numpy.ones((md.mesh.numberofvertices,1))
md=setflowequation(md,'SSA','all')

#control parameters
md.inversion.iscontrol=1
md.inversion.control_parameters=['MaterialsRheologyZbar']
md.inversion.min_parameters=10**-13*numpy.ones((md.mesh.numberofvertices,1))
md.inversion.max_parameters=numpy.ones((md.mesh.numberofvertices,1))
md.inversion.nsteps=2
md.inversion.cost_functions=101.*numpy.ones((md.inversion.nsteps,1))
md.inversion.cost_functions_coefficients=numpy.ones((md.mesh.numberofvertices,1))
md.inversion.gradient_scaling=0.9*numpy.ones((md.inversion.nsteps,1))
md.inversion.maxiter_per_step=2.*numpy.ones((md.inversion.nsteps,1))
md.inversion.step_threshold=0.99*numpy.ones((md.inversion.nsteps,1))
md.inversion.vx_obs=md.initialization.vx 
md.inversion.vy_obs=md.initialization.vy

md.cluster=generic('name',oshostname(),'np',3)
md=solve(md,StressbalanceSolutionEnum())

#Fields and tolerances to track changes
field_names     =['Gradient','Misfits','MaterialsRheologyZbar','Pressure','Vel','Vx','Vy']
field_tolerances=[1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12]
field_values=[\
   md.results.StressbalanceSolution.Gradient1,\
   md.results.StressbalanceSolution.J,\
   md.results.StressbalanceSolution.MaterialsRheologyZbar,\
   md.results.StressbalanceSolution.Pressure,\
   md.results.StressbalanceSolution.Vel,\
   md.results.StressbalanceSolution.Vx,\
   md.results.StressbalanceSolution.Vy,\
]
