function viscousheating=ViscousHeating(triaelem,grids,materials,inputs)
%VISCOUSHEATING - compute the viscous heating on each triaelem
%
%   Usage:
%      viscousheating=ViscousHeating(triaelem,grids,materials,inputs)
%
%   See also GETSTRAINRATE, GETVISCOSITY3D


	%initialize
	viscousheating=0;
	volume=0;

	%some variables
	NDOF1=1;
	numgrids=3;

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

	%Get all element grid data:
	xyz_list=getgriddata(triaelem,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('ViscousHeating error message: input velocity not present!');
	end

	%initialize vxvyvz_list
	vxvy_list=zeros(numgrids,2);

	%Build row indices for elementary vector.
	for i=1:numgrids,
		doflist=grids(triaelem.g(i)).grid.doflist;
		for j=1:2,
			dof=doflist(j);
			vxvy_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(triaelem,B_list,gauss_coord);
		matice.B=B_param; clear B_param.
	end

	% Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
	%get tria gaussian points as well as segment gaussian points. For tria gaussian 
	%points, order of integration is 2, because we need to integrate the product tB*D*B' 
	%which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
	%points, same deal, which yields 3 gaussian points.
	[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)];
	    
		%Build Deformational heating
		epsilon=GetStrainRate(triaelem,vxvy_list,xyz_list,gauss_coord);
		viscosity=GetViscosity2d(matice,epsilon); 

		%Build the strain rate matrix
		epsilon_matrix=[epsilon(1) epsilon(3)
				epsilon(3) epsilon(2)];

		epsilon_e=effective_value(epsilon_matrix);

		phi=2*viscosity*epsilon_e^2;

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

		viscousheating=viscousheating+phi*Jdet*gauss_weight;

		volume=volume+Jdet*gauss_weight;

	end %for ig=1:num_area_gauss,

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


end %end function
