function stress_bed=StressBedCompute(m,inputs,type);
%STRESSBEDCOMPUTE - compute the stress on the basal boundary condition
%
%   return a vector of size (numberofelements,1), holding the stress for 
%   every element on bed
%
%   Usage:
%      stress_bed=StressBedCompute(m,inputs,type)

%global variables
global cluster gridset

%recover fem model fields
elements=m.elements;
grids=m.grids;
materials=m.materials;
loads=m.loads;
gridset=m.gridset;

%figure out active elements that will take part in the stiffness and load generation
[n1,n2]=GetNumberOfActiveElements(elements);

%initialization
stress_bed=zeros((n2-n1)+1,1);

if strcmpi(type,'2d')
	disp('stress_bed not computed for 2d meshes')
	return
end

%initialize vectors
stress_bed=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'stress_n','stress_nn','normal_x',[],'normal_y',[],'normal_z',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[],'effectivevalue',[]);
stress_bed1=zeros((n2-n1)+1,6);
normal1=zeros((n2-n1)+1,3);
stress_n1=zeros((n2-n1)+1,3);
stress_nn1=zeros((n2-n1)+1,1);
A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);

%Go through all elements and call the stress_bed routine, then compute eigen values and vector
for n=n1:n2,
	if ~isempty(elements(n).element),
		[stress_vector,normal]=StressBed(elements(n).element,grids,materials,inputs);

		if ~isnan(stress_vector)
			stress_matrix=[stress_vector(1) stress_vector(4) stress_vector(5)
				      stress_vector(4)  stress_vector(2)  stress_vector(6)
				      stress_vector(5)  stress_vector(6)  stress_vector(3)];
			stress_n=stress_matrix*normal;
			stress_nn=normal'*stress_matrix*normal;

			%eigen values and vectors
			[directions,value]=eig(stress_matrix);

			%Plug into global vectors
			stress_bed1(n,:)=stress_vector;
			stress_n1(n,:)=stress_n;
			stress_nn1(n)=stress_nn;
			normal1(n,:)=normal';
			A1(n,1)=value(1,1); A2(n,1)=value(2,2); A3(n,1)=value(3,3);
			Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2); Vx3(n,1)=directions(1,3);
			Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2); Vy3(n,1)=directions(2,3);
			Vz1(n,1)=directions(3,1); Vz2(n,1)=directions(3,2); Vz3(n,1)=directions(3,3);
		else
			stress_bed1(n,:)=NaN*ones(1,6);
			stress_n1(n,:)=NaN*ones(1,3);
			stress_nn1(n)=NaN;
			normal1(n,:)=NaN*ones(1,3);
			A1(n,1)=NaN; A2(n,1)=NaN; A3(n,1)=NaN;
			Vx1(n,1)=NaN; Vx2(n,1)=NaN; Vx3(n,1)=NaN;
			Vy1(n,1)=NaN; Vy2(n,1)=NaN; Vy3(n,1)=NaN;
			Vz1(n,1)=NaN; Vz2(n,1)=NaN; Vz3(n,1)=NaN;
		end
	end
end

%plug results into outputs
%NB: Matlab sorts the eigen value in increasing order, we want the reverse
%components
stress_bed.xx=stress_bed1(:,1);
stress_bed.yy=stress_bed1(:,2);
stress_bed.zz=stress_bed1(:,3);
stress_bed.xy=stress_bed1(:,4);
stress_bed.xz=stress_bed1(:,5);
stress_bed.yz=stress_bed1(:,6);
%principal axis
stress_bed.principalvalue3=A1;
stress_bed.principalaxis3=[Vx1 Vy1 Vz1];
stress_bed.principalvalue2=A2;
stress_bed.principalaxis2=[Vx2 Vy2 Vz2];
stress_bed.principalvalue1=A3;
stress_bed.principalaxis1=[Vx3 Vy3 Vz3];
%norm or effective value
stress_bed.effectivevalue=1/sqrt(2)*sqrt(stress_bed.xx.^2+stress_bed.yy.^2+stress_bed.zz.^2+2*stress_bed.xy.^2+2*stress_bed.xz.^2+2*stress_bed.yz.^2);
%force
stress_bed.normal_x=normal1(:,1);
stress_bed.normal_y=normal1(:,2);
stress_bed.normal_z=normal1(:,3);
stress_bed.stress_n=stress_n1;
stress_bed.stress_nn=stress_nn1;
