function pe=CreatePVector(element,grids,materials,inputs,analysis_type);
%CREATEPVECTOR - create the load vector for an acceleratedtriaelem
%
%   global finite element load vector of the whole model
%   works only for MaAyeal's model
%
%   Usage:
%      Pe=CreatePVector(element,grids,materials,inputs,analysis_type);
% 
%   See also CREATEKMATRIX

global aire alpha beta

%some variables
numgrids=3;
NDOF2=2;
numdof=numgrids*NDOF2;

%recover extra inputs from users, at current convergence iteration.
[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
[surface_param surface_is_present]=recover_input(inputs,'surface');


if thickness_is_present,
	thickness=thickness_param(1:6:6*element.nods);
else
	thickness=element.thickness;
end

if surface_is_present,
	surface=surface_param(1:6:6*element.nods);
else
	surface=element.surface;
end


%retrieve rho_ice,g and thickness
matpar=materials(end).constants;
rho_ice=matpar.rho_ice;
rho_water=matpar.rho_water;
g=matpar.g;

%Build length_icefront and normal_icefront:
%[length_icefront,normal_icefront]=buildicefrontnormal(element.x,element.y,element.index_icefront);

Rhs_x=zeros(element.nel*27,1);
Rhs_y=zeros(element.nel*27,1);
Rhs_y=ones(element.nel*27,1);
row_rhs=zeros(element.nel*27,1);

count=-element.nel+1;
for i=1:3
	for n=1:3
		for m=1:3
			count=count+element.nel;

			Rhs_x(count:count+element.nel-1)= -rho_ice * g  * ...
			 thickness(element.index(:,n)).*surface(element.index(:,m)) ...
			.* aire(:) .* alpha(:,m) * ( (n==i)/6 + (n~=i)/12 );

			Rhs_y(count:count+element.nel-1)= -rho_ice * g  * ...
			 thickness(element.index(:,n)).*surface(element.index(:,m)) ...
			.* aire(:) .* beta(:,m) .* ( (n==i)/6 + (n~=i)/12 );

			row_rhs(count:count+element.nel-1)=element.index(:,i);

		end
	end
end

Rhs=full([sparse2(row_rhs,ones(element.nel*27,1),Rhs_x,element.nods,1)
   sparse2(row_rhs,ones(element.nel*27,1),Rhs_y,element.nods,1)]);

%Index icefront is skipped: 
%for k=1:2
%	for l=1:2
%		for j=1:2
%				Rhs(element.index_icefront(:,k))=Rhs(element.index_icefront(:,k)) + ...
%(rho_ice*g/2*thickness(element.index_icefront(:,l)).*thickness(element.index_icefront(:,j)) ...
%+ rho_water*g/2*(element.gridoniceshelf(element.index_icefront(:,k))) ...
%.*(-min(0,element.bed(element.index_icefront(:,l))).*min(0,element.bed(element.index_icefront(:,j))) ...
%+min(0,thickness(element.index_icefront(:,l))+element.bed(element.index_icefront(:,l)))...
%.*min(0,element.thickness(element.index_icefront(:,j))+element.bed(element.index_icefront(:,j)))))...
%.*normal_icefront(:,1).*length_icefront(:) ...
%/(4*(l==k & j==k) + 12*(l~=k | j~=k) );
%
%				Rhs(element.index_icefront(:,k)+element.nods)=Rhs(element.index_icefront(:,k)+element.nods) + ...
%(rho_ice*g/2*thickness(element.index_icefront(:,l)).*thickness(element.index_icefront(:,j)) ...
%+ rho_water*g/2*(element.gridoniceshelf(element.index_icefront(:,k))) ...
%.*(-min(0,element.bed(element.index_icefront(:,l))).*min(0,element.bed(element.index_icefront(:,j))) ...
%+min(0,thickness(element.index_icefront(:,l))+element.bed(element.index_icefront(:,l)))...
%.*min(0,thickness(element.index_icefront(:,j))+element.bed(element.index_icefront(:,j)))))...
%.*normal_icefront(:,2).*length_icefront(:) ...
%/(4*(l==k & j==k) + 12*(l~=k | j~=k) );
%
%		end
%	end
%end
%
clear Rhs_x Rhs_y row_rhs


%Create global vector.
pe=globalvector;
pe.nrows=numdof*element.nods;
pe.vector=spalloc(6*element.nods,1,2*element.nods);
pe.vector(1:6:6*element.nods)=Rhs(1:element.nods);
pe.vector(2:6:6*element.nods)=Rhs(element.nods+1:2*element.nods);
