#Test Name: SquareSheetConstrainedGia2d
from socket import gethostname

import numpy as np

from giaivins import giaivins
from model import *
from parameterize import *
from setflowequation import *
from setmask import *
from solve import *
from triangle import *


md = triangle(model(), '../Exp/Square.exp', 100000.)
md = setmask(md, '', '')
md = parameterize(md, '../Par/SquareSheetConstrained.py')

#GIA
md.gia = giaivins()
md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km
md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s
md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa
md.materials.lithosphere_density = 3.32 # in g/cm^3
md.materials.mantle_shear_modulus = 1.45e11 # in Pa
md.materials.mantle_density = 3.34 # in g/cm^3

#Indicate what you want to compute
md.gia.cross_section_shape = 1 # for square-edged x-section

#Define loading history (see test2001.m for the description)
md.timestepping.start_time = 2400000  # 2, 400 kyr
md.timestepping.final_time = 2500000  # 2, 500 kyr
md.geometry.thickness = np.array([
    np.append(md.geometry.thickness * 0.0, 0.0),
    np.append(md.geometry.thickness / 2.0, 0.1),
    np.append(md.geometry.thickness, 0.2),
    np.append(md.geometry.thickness, 1.0),
    np.append(md.geometry.thickness, md.timestepping.start_time)
    ]).T

#Solve for GIA deflection
md.cluster = generic('name', gethostname(), 'np', 3)
md.verbose = verbose('1111111')
md = solve(md, 'Gia')

#Fields and tolerances to track changes
field_names = ['UGia', 'UGiaRate']
field_tolerances = [1e-13, 1e-13]
field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]
