md=triangle(model(),'../Exp/Square.exp',40000.);
md=setmask(md,'all','');
md=parameterize(md,'../Par/SquareShelfConstrained.par');
md=setflowequation(md,'SSA','all');
md.cluster=generic('name',oshostname(),'np',3);
md.verbose=verbose('mprocessor',true,'module',true,'solution',true,'solver',true,'convergence',true,'control',true,'qmu',true,'autodiff',true);
md.stressbalance.requested_outputs={'StressTensor'};
md=solve(md,StressbalanceSolutionEnum());

	
index=md.mesh.elements;
nel=length(index);
x=md.mesh.x;
y=md.mesh.y;
z=md.mesh.z;
vel=md.results.StressbalanceSolution.Vel;

patch=zeros(nel,12);
for i=1:nel,
	patch(i,1)=x(index(i,1));
	patch(i,2)=y(index(i,1));
	patch(i,3)=z(index(i,1));
	patch(i,4)=x(index(i,2));
	patch(i,5)=y(index(i,2));
	patch(i,6)=z(index(i,2));
	patch(i,7)=x(index(i,3));
	patch(i,8)=y(index(i,3));
	patch(i,9)=z(index(i,3));
	patch(i,10)=vel(index(i,1));
	patch(i,11)=vel(index(i,2));
	patch(i,12)=vel(index(i,3));
end
