#Test Name: Esa2Dsurface
#Northern hemisphere example for north-south, east-west components of horiz motion 
#Same as test2111.m except that AIS is assumed to have located in Northern Hemisphere 

import numpy as np
from model import *
from socket import gethostname
from solve import *
from roundmesh import *
from love_numbers import *
from paterson import *

#mesh ais: {{{
md = model()
md = triangle(md,'../Exp/Ais.exp',200000); # max element size
# }}}
#define load: {{{
md.esa.deltathickness = np.zeros((md.mesh.numberofelements,))
disc_radius = 500 # km
index = md.mesh.elements
x_element = np.mean(md.mesh.x[index - 1], 1) - 1.0e6
y_element = np.mean(md.mesh.y[index - 1], 1) - 1.0e6
rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1   # 1 m water withdrawl
# }}}
#love numbers: {{{
nlov = 10000	# horizontal displacements do not work for low degree truncation, e.g., 101
md.esa.love_h = np.array(love_numbers('h','CF'))
md.esa.love_h = np.resize(md.esa.love_h, nlov + 1)
md.esa.love_l = np.array(love_numbers('l','CF'))
md.esa.love_l = np.resize(md.esa.love_l, nlov + 1)
# }}}
#mask:  {{{
#make sure wherever there is an ice load, that the mask is set to ice: 
md.mask.ice_levelset = np.ones((md.mesh.numberofvertices,))
pos = np.where(md.esa.deltathickness)
md.mask.ice_levelset[md.mesh.elements[pos,:]] = -1

#is ice grounded? 
md.mask.groundedice_levelset = -np.ones((md.mesh.numberofvertices,))
pos = np.where(md.mask.ice_levelset <= 0)
md.mask.groundedice_levelset[pos] = 1
# }}}
#geometry:  {{{
di = md.materials.rho_ice / md.materials.rho_water
md.geometry.thickness = np.ones((md.mesh.numberofvertices,))
md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices,))
md.geometry.base = md.geometry.surface - md.geometry.thickness
md.geometry.bed = md.geometry.base
# }}}
#materials:  {{{
md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices,))
md.materials.rheology_B = paterson(md.initialization.temperature)
md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
# }}}
#additional parameters, miscellaneous: {{{
md.miscellaneous.name='test2113';
md.esa.degacc=0.01;
md.esa.hemisphere = 1; # AIS is placed in Northern Hemisphere 
# }}}

#solve esa: {{{
md.esa.requested_outputs = ['EsaUmotion','EsaNmotion','EsaEmotion']
md.cluster = generic('name',gethostname(),'np',3)
md.verbose = verbose('111111111')
md = solve(md,'Esa')
# }}}
#Fields and tolerances to track changes: {{{
field_names     = ['EsaUmotion','EsaNmotion','EsaEmotion']
field_tolerances = [1e-13,1e-13,1e-13]
field_values = [
	md.results.EsaSolution.EsaUmotion,
	md.results.EsaSolution.EsaNmotion,
	md.results.EsaSolution.EsaEmotion,
	]
# }}}