function [stress_bed,normal]=StressBed(pentaelem,grids,materials,inputs)
%STRESSBED - compute the stress on the bed
%
%   This routine compute the stress_bed of each element, return NaN if not on bed
%
%   Usage:
%      [stress_bed,normal]=StressBed(pentaelem,grids,materials,inputs)
%
%   See also STRESS, STRESSSURFACE, GETSTRAINRATE, GETVISCOSITY3D

if ~pentaelem.onbed
	stress_bed=NaN;
	normal=NaN;
	return
end

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

%some variables
NDOF1=1;
numgrids=6;

%recover material parameters
matice=materials(pentaelem.matid).material;
matpar=materials(end).constants;
gravity=matpar.g;
rho_ice=matpar.rho_ice;
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');
[surface_param surface_is_present]=recover_input(inputs,'surface');
[flow_law_param flow_law_is_present]=recover_input(inputs,'B');
[elementonstokes_param elementonstokes_is_present]=recover_input(inputs,'elementonstokes');

%we need velocities to compute thermal profiles (even if it is a zero 
%vector). 
if ~velocity_is_present &~elementonstokes_is_present
	error('Stress error message: missing inputs!');
end

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

%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
	pressure(i)=velocity_param(doflist(4));	
	dof=doflist(1);
	if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
	if(surface_is_present) surface_list(i)=surface_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 
[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights]=GaussTria(2);

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

	%Pick up the gaussian point and its weight:
	gauss_weight=area_gauss_weights(igarea);
	gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea), -1];
    
	%Build deviatoric stress_bed 
	strainrate_g=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
	viscosity=GetViscosity3d(matice,strainrate_g);
	deviatoricstress_bed=viscosity*strainrate_g;

	%Build Pressure
	if elementonstokes_param(pentaelem.id)==stokesenum(),
		pressure_g=GetParameterValue(pentaelem,pressure,gauss_coord);
	else %use pattyn's and MacAyeal's assumption: P=sigma_zz'+rho_ice*g*(s-z)
		z_av=mean(xyz_list(:,3));
		if surface_is_present,
			s_av=mean(surface_list);
		else
			s_av=mean(pentaelem.s);
		end
		pressure_g=deviatoricstress_bed(3)+rho_ice*gravity*(s_av-z_av);
	end

	%Build stress_bed
	stress_bed_g=deviatoricstress_bed;
	stress_bed_g([1;2;3],1)=stress_bed_g(1:3,1)-pressure_g;         %stress_bed=deviatoricstress_bed-pressure

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

	stress_bed=stress_bed+stress_bed_g*Jdet*gauss_weight;
	volume=volume+Jdet*gauss_weight;

end %for ig=1:num_area_gauss,

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

%Build normal
normal=zeros(3,1);
n=cross((xyz_list(1,:)-xyz_list(3,:)),(xyz_list(3,:)-xyz_list(2,:)));
normal(1)=1/norm(n)*n(1);
normal(2)=1/norm(n)*n(2);
normal(3)=1/norm(n)*n(3);
