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

radius=1.e6
shelfextent=2.e5

md=roundmesh(model(),radius,50000.)
#fix center node to 0,0
rad=numpy.sqrt(md.mesh.x**2+md.mesh.y**2)
pos=numpy.argmin(rad)
md.mesh.x[pos]=0.
md.mesh.y[pos]=0.    #the closest node to the center is changed to be exactly at the center
xelem=numpy.mean(md.mesh.x[md.mesh.elements.astype(int)-1],axis=1)
yelem=numpy.mean(md.mesh.y[md.mesh.elements.astype(int)-1],axis=1)
rad=numpy.sqrt(xelem**2+yelem**2)
flags=numpy.zeros(md.mesh.numberofelements)
pos=numpy.nonzero(rad>=(radius-shelfextent))
flags[pos]=1
md=setmask(md,flags,'') 
md=parameterize(md,'../Par/RoundSheetShelf.py')
md=setflowequation(md,'SSA','all')
md.cluster=generic('name',oshostname(),'np',3)

md.transient.isthermal=False
md.transient.ismasstransport=False
md.transient.isstressbalance=True
md.transient.isgroundingline=True

#test different grounding line dynamics.
md.groundingline.migration='AggressiveMigration'
md=solve(md,TransientSolutionEnum())
element_on_iceshelf_agressive=md.results.TransientSolution[0].MaskGroundediceLevelset
vel_agressive=md.results.TransientSolution[0].Vel

md.groundingline.migration='SoftMigration'
md=solve(md,TransientSolutionEnum())
element_on_iceshelf_soft=md.results.TransientSolution[0].MaskGroundediceLevelset
vel_soft=md.results.TransientSolution[0].Vel

md.mask.groundedice_levelset=md.geometry.thickness + md.materials.rho_water/md.materials.rho_ice*md.geometry.bed
md.groundingline.migration='SubelementMigration'
md=solve(md,TransientSolutionEnum())
element_on_iceshelf_subelement=md.results.TransientSolution[0].MaskGroundediceLevelset
vel_subelement=md.results.TransientSolution[0].Vel

md.groundingline.migration='SubelementMigration2'
md=solve(md,TransientSolutionEnum())
element_on_iceshelf_subelement2=md.results.TransientSolution[0].MaskGroundediceLevelset
vel_subelement2=md.results.TransientSolution[0].Vel

#Fields and tolerances to track changes
field_names     =['ElementOnIceShelfAggressive','VelAggressive','ElementOnIceShelfSoft','VelSoft','ElementOnIceShelfSubelement','VelSubelement','ElementOnIceShelfSubelement2','VelSubelement2']
field_tolerances=[1e-13,2e-13,1e-13,2e-13,1e-13,2e-13,1e-13,2e-13]
field_values=[element_on_iceshelf_agressive,vel_agressive,element_on_iceshelf_soft,vel_soft,element_on_iceshelf_subelement,vel_subelement,element_on_iceshelf_subelement2,vel_subelement2]
