function deviatoricstress=DeviatoricStress(pentaelem,grids,materials,inputs)
%DEVIATORICSTRESS - computes the deviatoric stress of a pentaelem
% 
%   Usage:
%      deviatoricstress=DeviatoricStress(pentaelem,grids,materials,inputs)

	%initialize
	deviatoricstress=zeros(6,1);
	volume=0;

	%some variables
	NDOF1=1;
	numgrids=6;

	%recover material parameters
	matice=materials(pentaelem.matid).material;
	B=matice.B;

	%Get all element grid data:
	xyz_list=getgriddata(pentaelem,grids);

	%recover extra inputs
	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
	[flow_law_param flow_law_is_present]=recover_input(inputs,'B');

	%we need velocities to compute thermal profiles (even if it is a zero 
	%vector). 
	if ~velocity_is_present,
		error('DeviatoricStress error message: input velocity not present!');
	end

	%initialize vxvyvz_list
	vxvyvz_list=zeros(numgrids,3);

	%Build row indices for elementary vector.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist;
		for j=1:3,
			dof=doflist(j);
			vxvyvz_list(i,j)=velocity_param(dof);
		end
		dof=doflist(1);
		if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
	end
	
	%Update material parameter that deals with ice rigidity: 
	if flow_law_is_present,
		B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
		matice.B=B_param; clear B_param.
	end
	
	% Get gaussian points and weights 
	area_order=2;
	num_vert_gauss=2;
	[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

	%Start  looping on the number of gaussian points:
	for igarea=1:num_area_gauss,
		for igvert=1:num_vert_gauss,

			%Pick up the gaussian point and its weight:
			gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
			gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
		    
			%Build Stress
			strainrate_g=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
			viscosity=GetViscosity3d(matice,strainrate_g);
			deviatoricstress_g=viscosity*strainrate_g;

			%Get Jacobian determinant:
			Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);

			deviatoricstress=deviatoricstress+deviatoricstress_g*Jdet*gauss_weight;
			volume=volume+Jdet*gauss_weight;

		end %for ig=1:vert_gauss,
	end %for ig=1:num_area_gauss,

%Divide deviatoricstress, integrated over volume, by volume.
deviatoricstress=deviatoricstress/volume;
