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

	if strcmpi(analysis_type,'diagnostic_hutter'),

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

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

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

	%some variables
	numgrids=2;
	NDOF2=2;

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

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

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

	%recover extra inputs
	[surface_param surface_is_present]=recover_input(inputs,'surface');
	[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/beamelem error message: slope_surface is missing in inputs');
	end

	%Get all element grid data:
	z_list=getgriddata(beamelem,grids);
	
	%Initialize inputs
	surface_list=zeros(numgrids,1);
	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(beamelem.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(surface_is_present) surface_list(i)=surface_param(dof);end;
		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 surface
	if surface_is_present
		surface=surface_list(1);
	else
		surface=beamelem.s(1);
	end

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

	% Get gaussian points and weights.
	num_gauss=3;
	[gauss_coords,gauss_weights]=GaussSegment(num_gauss);

	%Start  looping on the number of gaussian points:
	for ig=1:num_gauss,
		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=gauss_coords(ig);
	
		%Update material if temperature is provided.
		if B_is_present,
			B=GetParameterValue(beamelem,B_list,gauss_coord);
		elseif temperature_is_present,
			temperature=GetParameterValue(beamelem,temperature_list,gauss_coord);
			B=paterson(temperature);
		else
			B=matice.B;
		end

		%compute constant_part
		z_g=GetParameterValue(beamelem,z_list,gauss_coord);

		%Get Jacobian determinant:
		Jdet=GetJacobianDeterminant(beamelem,z_list,gauss_coord);
		%disp(sprintf('Jacobian determinant: %f\n', Jdet));

		%Get nodal functions
		l1l2=GetNodalFunctions(beamelem,gauss_coord);

		%Build pe_gaussian vector:
		pe_gaussian=zeros(numgrids*NDOF2,1);

		for j=1:NDOF2,
			pe_gaussian(NDOF2+j)=constant_part*((surface-z_g)/B)^n*slope(j)*Jdet*gauss_weight;
		end

		%Add pe_gaussian vector to pe:
		for i=1:pe.nrows,
			pe.terms(i)=pe.terms(i)+pe_gaussian(i);
		end
	end %for ig=1:num_area_gauss,

	%Deal with lower surface
	if (beamelem.onbed==1),

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

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

		%Add pe_gg_drag_gaussian to pe
		pe.terms(1)=pe.terms(1)+ub;
		pe.terms(2)=pe.terms(2)+vb;
	end
end %end function
