function deviatoricstress=DeviatoricStressCompute(m,inputs,type);
%DEVIATORICSTRESSCOMPUTE - compute the deviatoric stress components
%
%   Usage:
%      deviatoricstress=DeviatoricStressCompute(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);

if strcmpi(type,'2d')
	%initialize vectors
	deviatoricstress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[]);
	deviatoricstress1=zeros((n2-n1)+1,3);
	A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); 
	A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1);

	%Go through all elements and call the deviatoricstress routine, then compute eigen values and vector
	for n=n1:n2,
		if ~isempty(elements(n).element),
			deviatoricstressvector=DeviatoricStress(elements(n).element,grids,materials,inputs)';

			deviatoricstressmatrix=[deviatoricstressvector(1) deviatoricstressvector(3)
					      deviatoricstressvector(3)  deviatoricstressvector(2)]; 

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

                        %Plug into global vectors
			deviatoricstress1(n,:)=deviatoricstressvector;
                        A1(n,1)=value(1,1); A2(n,1)=value(2,2);
                        Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2);
                        Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2);
		end
	end

	%plug results into outputs
	%NB: Matlab sorts the eigen value in increasing order, we want the reverse
	deviatoricstress.xx=deviatoricstress1(:,1);
	deviatoricstress.yy=deviatoricstress1(:,2);
	deviatoricstress.xy=deviatoricstress1(:,3);
	deviatoricstress.principalvalue2=A1;
	deviatoricstress.principalaxis2=[Vx1 Vy1];
	deviatoricstress.principalvalue1=A2;
	deviatoricstress.principalaxis1=[Vx2 Vy2];
	%norm or effective value
	deviatoricstress.effectivevalue=1/sqrt(2)*sqrt(deviatoricstress.xx.^2+deviatoricstress.yy.^2+2*deviatoricstress.xy.^2);
else
	%initialize vectors
	deviatoricstress=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[]);
	deviatoricstress1=zeros((n2-n1)+1,6);
	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 deviatoricstress routine, then compute eigen values and vector
	for n=n1:n2,
		if ~isempty(elements(n).element),
			deviatoricstressvector=DeviatoricStress(elements(n).element,grids,materials,inputs)';

			deviatoricstressmatrix=[deviatoricstressvector(1) deviatoricstressvector(4) deviatoricstressvector(5)
				      deviatoricstressvector(4)  deviatoricstressvector(2)  deviatoricstressvector(6)
				      deviatoricstressvector(5)  deviatoricstressvector(6)  deviatoricstressvector(3)];

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

                        %Plug into global vectors
			deviatoricstress1(n,:)=deviatoricstressvector;
                        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);
		end
	end

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