md=triangle(model(),'../Exp/Square.exp',150000.);
md=setmask(md,'all','');
md.materials=matdamageice();
md=parameterize(md,'../Par/SquareShelf.par');
md.damage.isdamage=1;
md.damage.D=zeros(md.mesh.numberofvertices,1);
md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
md.damage.law=1;

%boundary conditions for damage, to be put in SquareShelf.par
boundary=zeros(md.mesh.numberofvertices,1);
boundary(md.mesh.segments(:,1:2))=1;
pos=find(~md.mask.ice_levelset==0 & boundary==1);
md.damage.spcdamage(pos)=0.5;
md.damage.stabilization=1;

md.damage.requested_outputs={'default','DamageF'};

md=setflowequation(md,'SSA','all');
md.cluster=generic('name',oshostname(),'np',3);
md=solve(md,DamageEvolutionSolutionEnum());

%Fields and tolerances to track changes
field_names={'D','F'};
field_tolerances={1e-13,1e-13};
field_values={...
		(md.results.DamageEvolutionSolution.DamageDbar),...
		(md.results.DamageEvolutionSolution.DamageF),...
	};
