function pe=CreatePVector(singelem,grids,materials,inputs,analysis_type);
%CREATEPVECTOR - create the load vector for a singelem
%
%   works only for Hutter's model
%
%   Usage:
%      Pe=CreatePVector(singelem,grids,materials,inputs,analysis_type);
% 
%   See also CREATEKMATRIX


	if strcmpi(analysis_type,'diagnostic_hutter'),

		pe=CreatePVectorHutter(singelem,grids,materials,inputs);

	else
		error('singelem/CreatePVector error message: analysis type not implemented yet')
	end
end %end function

function pe=CreatePVectorHutter(singelem,grids,materials,inputs);

	%some variables
	numgrids=1;
	NDOF2=2;

	%Create elementary vector.
	pe=elemvector(numgrids*NDOF2);

	%recover material
	matice=materials(singelem.matid).material;
	matpar=materials(end).constants;

	%recover material parameters
	gravity=matpar.g;
	rho_ice=matpar.rho_ice;
	n=matice.n;

	%recover extra inputs
	[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
	[B_param B_is_present]=recover_input(inputs,'B');
	[temperature_param temperature_is_present]=recover_input(inputs,'temperature');
	[slopesurface_param slopesurface_is_present]=recover_input(inputs,'slopesurface');

	if ~slopesurface_is_present,
		error('CreatePVector/singelem error message: slope_surface is missing in inputs');
	end

	%Initialize inputs
	B_list=zeros(numgrids,1);
	temperature_list=zeros(numgrids,1);
	slopesurface_list=zeros(numgrids,2);

	%Build row indices for elementary vector.
	for i=1:numgrids,
		doflist=grids(singelem.g(i)).grid.doflist;
		for j=1:NDOF2,
			pe.row_indices((i-1)*NDOF2+j)=doflist(j);
			slopesurface_list(i,j)=slopesurface_param(doflist(j));
		end
		dof=doflist(1);
		if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
		if(B_is_present) B_list(i)=B_param(dof);end;
		if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
	end

	%compute slope2 slopex and slopy
	slope=zeros(2,1);
	slope(1)=slopesurface_list(1,1); %2 and 1 have the same slopes (mesh has been verticaly extruded)
	slope(2)=slopesurface_list(1,2);
	slope2=slope(1)^2+slope(2)^2;

	%compute thickness
	if thickness_is_present
		thickness=thickness_list(1);
	else
		thickness=singelem.h(1);
	end

	%compute ub
	ub=-1.58*10^-10*rho_ice*gravity*thickness*slope(1);
	vb=-1.58*10^-10*rho_ice*gravity*thickness*slope(2);

	%compute constant_part
	constant_part=-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2);

	%Update material if temperature is provided.
	if B_is_present,
		B=GetParameterValue(singelem,B_list,gauss_coord);
	elseif temperature_is_present,
		temperature=GetParameterValue(singelem,temperature_list,gauss_coord);
		B=paterson(temperature);
	else
		B=matice.B;
	end

	pe.terms(1)=ub-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2)*thickness^(n)/(B^n*(n+1))*slope(1);
	pe.terms(2)=vb-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2)*thickness^(n)/(B^n*(n+1))*slope(2);

end %end function
