%Test Name: FlowbandFSsheetshelfTrans
%mesh parameters
x =[-5:.5:5]';
[b h sea]=NowickiProfile(x);
x = x*10^3;
h = h*10^3;
b = (b-sea)*10^3;

%mesh domain
md=bamgflowband(model(),x,b+h,b,'hmax',150);

%parameterize
md.geometry.surface   = interp1(x,b+h,md.mesh.x);
md.geometry.base       = interp1(x,b,md.mesh.x);
md.geometry.thickness = md.geometry.surface-md.geometry.base;
md.mask.ice_levelset  = - ones(md.mesh.numberofvertices,1);
md.mask.ice_levelset(find(vertexflags(md.mesh,2))) = 0;
md.mask.ocean_levelset = double(md.mesh.x<0)-.5;

md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1);
md.materials.rheology_B=paterson(md.initialization.temperature);
md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
md.friction.coefficient=zeros(md.mesh.numberofvertices,1);
md.friction.coefficient(find(vertexflags(md.mesh,1)))=20;
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=ones(md.mesh.numberofelements,1);
md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
md.stressbalance.spcvx(find(vertexflags(md.mesh,4)))=800;
md.stressbalance.spcvy(find(vertexflags(md.mesh,4)))=0;

%Misc
md=setflowequation(md,'FS','all');
md.flowequation.fe_FS='TaylorHood';
md.stressbalance.abstol=NaN;
md.miscellaneous.name = 'test703';

%Transient settings
md.timestepping.time_step=0.000001;
md.timestepping.final_time=0.000005;
md.stressbalance.shelf_dampening=1;
md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
posb=find((md.mesh.x>0) & md.mesh.vertexonbase);
md.basalforcings.groundedice_melting_rate(posb)=18;
md.basalforcings.floatingice_melting_rate(posb)=18;
md.initialization.vx=zeros(md.mesh.numberofvertices,1);
md.initialization.vy=zeros(md.mesh.numberofvertices,1);
md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
md.transient.isthermal=0;
md.masstransport.isfreesurface=1;
md.groundingline.migration='None';

%Go solve
md.cluster=generic('np',3);
md.stressbalance.shelf_dampening=1;
md1=solve(md,'Transient');

md.stressbalance.shelf_dampening=0;
md=solve(md,'Transient');

%Fields and tolerances to track changes
field_names     ={...
	'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',...
	'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',...
	'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3',...
	'Vx1_damp','Vy1_damp','Vel1_damp','Pressure1_damp','Bed1_damp','Surface1_damp','Thickness1_damp',...
	'Vx2_damp','Vy2_damp','Vel2_damp','Pressure2_damp','Bed2_damp','Surface2_damp','Thickness2_damp',...
	'Vx3_damp','Vy3_damp','Vel3_damp','Pressure3_damp','Bed3_damp','Surface3_damp','Thickness3_damp'};
field_tolerances={...
	2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,...
	2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,...
	2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,...
	5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,...
	5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,...
	5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10};
field_values={...
	(md.results.TransientSolution(1).Vx),...
	(md.results.TransientSolution(1).Vy),...
	(md.results.TransientSolution(1).Vel),...
	(md.results.TransientSolution(1).Pressure),...
	(md.results.TransientSolution(1).Base),...
	(md.results.TransientSolution(1).Surface),...
	(md.results.TransientSolution(1).Thickness),...
	(md.results.TransientSolution(2).Vx),...
	(md.results.TransientSolution(2).Vy),...
	(md.results.TransientSolution(2).Vel),...
	(md.results.TransientSolution(2).Pressure),...
	(md.results.TransientSolution(2).Base),...
	(md.results.TransientSolution(2).Surface),...
	(md.results.TransientSolution(2).Thickness),...
	(md.results.TransientSolution(3).Vx),...
	(md.results.TransientSolution(3).Vy),...
	(md.results.TransientSolution(3).Vel),...
	(md.results.TransientSolution(3).Pressure),...
	(md.results.TransientSolution(3).Base),...
	(md.results.TransientSolution(3).Surface),...
	(md.results.TransientSolution(3).Thickness),...
	(md1.results.TransientSolution(1).Vx),...
	(md1.results.TransientSolution(1).Vy),...
	(md1.results.TransientSolution(1).Vel),...
	(md1.results.TransientSolution(1).Pressure),...
	(md1.results.TransientSolution(1).Base),...
	(md1.results.TransientSolution(1).Surface),...
	(md1.results.TransientSolution(1).Thickness),...
	(md1.results.TransientSolution(2).Vx),...
	(md1.results.TransientSolution(2).Vy),...
	(md1.results.TransientSolution(2).Vel),...
	(md1.results.TransientSolution(2).Pressure),...
	(md1.results.TransientSolution(2).Base),...
	(md1.results.TransientSolution(2).Surface),...
	(md1.results.TransientSolution(2).Thickness),...
	(md1.results.TransientSolution(3).Vx),...
	(md1.results.TransientSolution(3).Vy),...
	(md1.results.TransientSolution(3).Vel),...
	(md1.results.TransientSolution(3).Pressure),...
	(md1.results.TransientSolution(3).Base),...
	(md1.results.TransientSolution(3).Surface),...
	(md1.results.TransientSolution(3).Thickness),...
	};
