%Test Name: SquareSheetConstrainedGia2d
%GIA test, inspired on test101. Running default GIA Ivins class.
md=triangle(model(),'../Exp/Square.exp',100000.);
md=setmask(md,'','');
md=parameterize(md,'../Par/SquareSheetConstrained.par');

%GIA: 
md.gia=giaivins();
md.gia.lithosphere_thickness=100.*ones(md.mesh.numberofvertices,1); % in km
md.gia.mantle_viscosity=1.0*10^21*ones(md.mesh.numberofvertices,1); % in Pa.s
md.materials.lithosphere_shear_modulus=6.7*10^10;                   % in Pa
md.materials.lithosphere_density=3.32;                              % in g/cm^-3
md.materials.mantle_shear_modulus=1.45*10^11;                       % 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 
md.timestepping.start_time=2400000; %2,400 kyr :: EVALUATION TIME
% to get rid of default final_time: make sure final_time>start_time
md.timestepping.final_time=2500000; %2,500 kyr
md.geometry.thickness=[...
	[md.geometry.thickness; 0.0],...
	[md.geometry.thickness; md.timestepping.start_time],...
	];

%% solve for GIA deflection 
md.cluster=generic('name',oshostname(),'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),...
	};
