%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.groundedice_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.damage.D=zeros(md.mesh.numberofvertices,1);
md.damage.spcdamage=NaN(md.mesh.numberofvertices,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 = 'flowline';

%Go solve
md.cluster=generic('np',2);
md=solve(md,StressbalanceSolutionEnum());

%Fields and tolerances to track changes
field_names     ={'Vx','Vy','Vel','Pressure'};
field_tolerances={2e-08,2e-08,2e-08,2e-08};
field_values={...
	(md.results.StressbalanceSolution.Vx),...
	(md.results.StressbalanceSolution.Vy),...
	(md.results.StressbalanceSolution.Vel),...
	(md.results.StressbalanceSolution.Pressure)};
