function pe=CreatePVector(pentaelem,grids,materials,inputs,analysis_type);
%CREATEPVECTOR - create the load vector for pentaelem
%
%   this load vector works for MacAyeal, Pattyn's  and Stokes' model ,
%   thermal model (transient ans steady), prognostic, melting and 
%   for the computation of the bed slope and surface slope
%
%   Usage:
%      pe=CreatePVector(pentaelem,grids,materials,inputs,analysis_type)
%
%   See also CREATEKMATRIX


if strcmpi(analysis_type,'diagnostic_horiz'),

	pe=CreatePVectorHoriz(pentaelem,grids,materials,inputs,analysis_type);

elseif strcmpi(analysis_type,'diagnostic_vert'),

	pe=CreatePVectorVert(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'diagnostic_basevert'),

	pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'prognostic'),

	pe=CreatePVectorPrognostic(pentaelem,grids,materials,inputs);

elseif (strcmpi(analysis_type,'thermalsteady') | strcmpi(analysis_type,'thermaltransient')),

	pe=CreatePVectorThermal(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'melting')

	pe=elemvector(0);

elseif strcmpi(analysis_type,'diagnostic_stokes'),

	pe=CreatePVectorStokes(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') | strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),

	pe=CreatePVectorSlopeCompute(pentaelem,grids,materials,inputs,analysis_type);


end
end %end function

function pe=CreatePVectorStokes(pentaelem,grids,materials,inputs,analysis_type);

%some variables
numgrids=6;
numgrids2d=3;
NDOF4=4;

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

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

%recover material parameters
gravity=matpar.g;
viscosity_overshoot=matpar.viscosity_overshoot;
rho_ice=matpar.rho_ice;
rho_water=matpar.rho_water;

%recover conditioning number for the matrix
conditioning=pentaelem.reconditioning_number;

%recover extra inputs
[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
[temperature_param temperature_is_present]=recover_input(inputs,'temperature');
[surface_param surface_is_present]=recover_input(inputs,'surface');
[bed_param bed_is_present]=recover_input(inputs,'bed');
[melting_param melting_is_present]=recover_input(inputs,'melting');
[basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
[flow_law_param flow_law_is_present]=recover_input(inputs,'B');
[velocity_param velocity_is_present]=recover_input(inputs,'velocity');

%Get all element grid data:
xyz_list=getgriddata(pentaelem,grids);

%Initialize inputs
if velocity_is_present,
	vxvyvz_list=zeros(numgrids,3);
end
B_list=zeros(numgrids,1);
basal_drag_list=zeros(numgrids,1);
melting_list=zeros(numgrids,1);
thickness_list=zeros(numgrids,1);
surface_list=zeros(numgrids,1);
bed_list=zeros(numgrids,1);
temperature_list=zeros(numgrids,1);

%Build row indices for elementary vector.
for i=1:numgrids,
	
	doflist=grids(pentaelem.g(i)).grid.doflist;
	for j=1:NDOF4,
		dof=doflist(j);
		pe.row_indices((i-1)*NDOF4+j)=dof;
		if j<4, %the fourth dof corresponds to pressure and not velocity
			if velocity_is_present, vxvyvz_list(i,j)=velocity_param(dof); end;
		end
	end
	dof=doflist(1);
	if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
	if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
	if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
	if(surface_is_present) surface_list(i)=surface_param(dof);end;
	if(bed_is_present) bed_list(i)=bed_param(dof);end;
	if(melting_is_present) melting_list(i)=melting_param(dof);end;
	if(basal_drag_is_present) basal_drag_list(i)=basal_drag_param(dof);end;

end

% Get gaussian points and weights 
area_order=5;
num_vert_gauss=5;
[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

%for i=1:num_gauss,
%	disp(sprintf('Gauss coord %i: %f %f %f Weight: %f\n',i,first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
%end
%
pe_temp=zeros(27,1);
Kebubble_temp=zeros(27,3);

%Start  looping on the number of gaussian points:
for igarea=1:num_area_gauss,
	for igvert=1:num_vert_gauss,
		%Pick up the gaussian point and its weight:
		gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
		gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];

		%In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 
		%element itself: 
		if (~pentaelem.shelf) & (pentaelem.friction_type==1),
			if basal_drag_is_present,
				plastic_stress=GetParameterValueStokes(pentaelem,basal_drag_list,gauss_coord);
			else
				plastic_stress=GetParameterValueStokes(pentaelem,pentaelem.k,gauss_coord);
			end
		end

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

		%Get nodal functions
		l1l7=GetNodalFunctionsStokes(pentaelem,gauss_coord);

		%Build pe_gaussian vector:
		pe_gaussian=zeros(numgrids*NDOF4+3,1);

		if(~pentaelem.shelf) & (pentaelem.friction_type==1),
			for i=1:numgrids+1,
				for j=1:2,
					pe_gaussian((i-1)*NDOF4+j)=-plastic_stress*Jdet*gauss_weight*l1l7(i);
				end
				pe_gaussian((i-1)*NDOF4+3)=-rho_ice*gravity*Jdet*gauss_weight*l1l7(i);
			end
		else
			for i=1:numgrids+1,
				pe_gaussian((i-1)*NDOF4+3)=-rho_ice*gravity*Jdet*gauss_weight*l1l7(i);
			end
		end

		%Add pe_gaussian vector to pe:
		pe_temp=pe_temp+pe_gaussian;

		%Compute the bubble part of the elementary stiffness matrix
		
		%Get strain rate, if velocity has been supplied: 
		if velocity_is_present,
			epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
		else
			epsilon=zeros(6,1);
		end
		%disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3),epsilon(4),epsilon(5),epsilon(6)));
		
		%Update material parameter that deals with ice rigidity. The update depends either on the 
		%temperature (if we are running thermal, or transient), or on B (if we are running inverse 
		%control methods on B). The latter prevails. 
		
		%Update material if temperature is provided.
		if temperature_is_present,
			temperature=GetParameterValueStokes(pentaelem,temperature_list,gauss_coord);
			matice.B=paterson(temperature);
		end

		%Update material if flow law is specified. This will erase the previous change 
		%on B when temperature is provided. 
		if flow_law_is_present,
			B_param=GetParameterValueStokes(pentaelem,B_list,gauss_coord);
			matice.B=B_param; clear B_param.
		end

		%Get viscosity at last two iterations: 
		viscosity=GetViscosity3d(matice,epsilon);
		%disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));

		%Get B and Bprime matrices:
		B=GetBStokes(pentaelem,xyz_list,gauss_coord);
		Bprime=GetBprimeStokes(pentaelem,xyz_list,gauss_coord);
		Bprime_bubble=Bprime(:,25:27);

		% Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 
		% onto this scalar matrix, so that we win some computational time: */
		D=gauss_weight*Jdet*diag([viscosity*ones(6,1);-conditioning*ones(2,1)]);

		% Do the triple product tB*D*Bprime: 
		Ke_gg_gaussian=B'*D*Bprime_bubble;

		%Add Ke_gg_gaussian to Ke_temp: 
		Kebubble_temp=Kebubble_temp+Ke_gg_gaussian;

	end %for igarea=1:num_area_gauss,
end %for igvert=1:num_vert_gauss,

%Deal with water pressure integration: 

if (pentaelem.onbed==1 & pentaelem.shelf==1),

	%Plug the grids on surface into the triaelem
	if thickness_is_present,
		thickness=thickness_list(1:3);
	else			
		thickness=pentaelem.h(1:3);
	end
	if bed_is_present,
		bed=bed_list(1:3);
	else			
		bed=pentaelem.b(1:3);
	end
	if basal_drag_is_present,
		basal_drag=basal_drag_list(1:3);
	else			
		basal_drag=pentaelem.k(1:3);
	end

	%Get the coordinates of the three grids
	xyz_list_tria=xyz_list(1:3,:);

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	%Build normal vector to the surface pointing outward the glacier system
	n=zeros(3,1);
	normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(3,:)-xyz_list_tria(2,:)));
	n(1)=1/norm(normal)*normal(1);
	n(2)=1/norm(normal)*normal(2);
	n(3)=1/norm(normal)*normal(3);

	pe_gg_water=zeros(27,1);

	for ig=1:num_gauss2D,

		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
	
		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);

		%Get bed at gaussian point
		bed_gg=GetParameterValue(triaelem,bed,gauss_coord);
		
		%Get L matrix 
		L=GetL(triaelem,gauss_coord,1);
   
		%Compute pressure
		pressure=gravity*rho_water*bed_gg;

		%Compute elementary load vector
		for i=1:numgrids2d,
			for j=1:3,
				pe_gg_water((i-1)*NDOF4+j)=pe_gg_water((i-1)*NDOF4+j)+pressure*gauss_weight*Jdettria*L(i)*n(j);
			end
		end
	end

	%Add pe_gaussian vector to pe:
	pe_temp=pe_temp+pe_gg_water;

end%ends of water pressure

%Create reduced vector 
pe_reduced=ReduceVectorStokes(pentaelem,Kebubble_temp, pe_temp);

pe.terms=pe.terms+pe_reduced;

end %end function


function pe=CreatePVectorSlopeCompute(pentaelem,grids,materials,inputs,analysis_type);

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.
	if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
		if pentaelem.onbed~=1,
			pe=elemvector(0);
			return;
		end
	elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
		if pentaelem.onsurface~=1,
			pe=elemvector(0);
			return;
		end
	end
	
	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%Create elementary stiffness matrix 
	pe=elemvector(numdof);

	%Get all element grid data:
	xyz_list=getgriddata(pentaelem,grids); 

	%Just keep the first 3 grids and recover extra inputs
	if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
		xyz_list=xyz_list(1:3,:);
		[param param_is_present]=recover_input(inputs,'bed');
	elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
		xyz_list=xyz_list(4:6,:);
		[param param_is_present]=recover_input(inputs,'surface');
	end

	param_list=zeros(numgrids,1);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
			doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
			doflist=grids(pentaelem.g(i+3)).grid.doflist; %list of dofs in the g-set: last three grids
		end	
		if(param_is_present) param_list(i)=param(doflist(1));end;
		dof=doflist(1);
		pe.row_indices(i)=dof;
	end

	if param_is_present, 
		param=param_list(1:3); 
	else 
		if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
			param=pentaelem.b(1:3); 
		elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
			param=pentaelem.s(4:6); 
		end
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	%for icesheets
	for ig=1:num_gauss2D,
	
		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		
		%Get bed slope
		dparam=GetParameterDerivativeValue(triaelem,param,xyz_list,gauss_coord);
		dparamdx=dparam(1);
		dparamdy=dparam(2);

		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L: 
		L=GetL(triaelem,gauss_coord,DOFPERGRID);

		%Build gaussian vector 
		if strcmp(analysis_type,'bed_slope_compute_x') | strcmp(analysis_type,'surface_slope_compute_x'), 
			pe_g_gaussian=Jdettria*gauss_weight*dparamdx*L';
		elseif strcmp(analysis_type,'bed_slope_compute_y') | strcmp(analysis_type,'surface_slope_compute_y'),
			pe_g_gaussian=Jdettria*gauss_weight*dparamdy*L';
		end

		%Add Ke_gg_drag_gaussian to Ke
		pe.terms=pe.terms+pe_g_gaussian;
	end

end %end function



function pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.
	if pentaelem.onbed~=1,
		pe=elemvector(0);
		return;
	end

	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%recover material parameters
	matice=materials(pentaelem.matid).material;
	matpar=materials(end).constants;
	rho_ice=matpar.rho_ice;
	rho_water=matpar.rho_water;
	di=rho_ice/rho_water;

	%Create elementary stiffness matrix 
	pe=elemvector(numdof);

	%Get all element grid data:
	xyz_list=getgriddata(pentaelem,grids); 

	%Just keep the first 3 grids
	xyz_list=xyz_list(1:3,:);

	%recover extra inputs
	[bed_param bed_is_present]=recover_input(inputs,'bed');
	[melting_param melting_is_present]=recover_input(inputs,'melting');
	[accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
	[velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
	[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
	if (~velocity_is_present | ~velocity_average_is_present),
		error('CreatePVector error message: no input velocity!');
	end

	%initialize velocity array
	vxvy_list=zeros(numgrids,2);
	vxvy_average_list=zeros(numgrids,2);
	bed_list=zeros(numgrids,1);
	melting_list=zeros(numgrids,1);
	accumulation_list=zeros(numgrids,1);
	thickness_list=zeros(numgrids,1);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		for j=1:2,
			dof=doflist(j);
			vxvy_list(i,j)=velocity_param(dof);
			vxvy_average_list(i,j)=velocity_average_param(dof);
		end
		if(bed_is_present) bed_list(i)=bed_param(doflist(1));end;
		if(melting_is_present) melting_list(i)=melting_param(doflist(1));end;
		if(accumulation_is_present) accumulation_list(i)=accumulation_param(doflist(1));end;
		if(thickness_is_present) thickness_list(i)=thickness_param(doflist(1));end;
		dof=doflist(3); %z degree of freedom
		pe.row_indices(i)=dof;
	end

	if bed_is_present, bed=bed_list(1:3); else bed=pentaelem.b(1:3); end
	if thickness_is_present, thickness=thickness_list(1:3); else thickness=pentaelem.h(1:3); end
	if melting_is_present, melting=melting_list(1:3); else melting=pentaelem.melting(1:3); end
	if accumulation_is_present, accumulation=accumulation_list(1:3); else accumulation=pentaelem.accumulation(1:3); end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	%for icesheets
	for ig=1:num_gauss2D,
	
		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		
		%Get melting at gaussian point.
		melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
		
		%Get velocity at gaussian point.
		vx=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
		vy=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);

		%Get bed slope
		db=GetParameterDerivativeValue(triaelem,bed,xyz_list,gauss_coord);
		dbdx=db(1);
		dbdy=db(2);

		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L: 
		L=GetL(triaelem,gauss_coord,DOFPERGRID);

		%Build gaussian vector 
		pe_g_basevert_gaussian=Jdettria*gauss_weight*(vx*dbdx+vy*dbdy-melting_gauss)*L';

		%Add Ke_gg_drag_gaussian to Ke
		pe.terms=pe.terms+pe_g_basevert_gaussian;
	end

	%for iceshelves, vertical velocity at the base is determined using the hydrostatic equilibrium -> we add a corrective term 
	%to the icesheet term.
	if pentaelem.shelf,

		%build the vector thickness*velocity_average
		HU=thickness.*vxvy_average_list(:,1);
		HV=thickness.*vxvy_average_list(:,2);

		for ig=1:num_gauss2D,
	
			%Pick up the gaussian point and its weight:
			gauss_weight=gauss_weights(ig);
			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
			
			%get div(H*Vel)
			dHU=GetParameterDerivativeValue(triaelem,HU,xyz_list,gauss_coord);
			dHV=GetParameterDerivativeValue(triaelem,HV,xyz_list,gauss_coord);
			divHVel=dHU(1)+dHV(2);

			%Get melting and accumulation at gaussian point.
			melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
			accumulation_gauss=GetParameterValue(triaelem,accumulation,gauss_coord);

			%Get L: 
			L=GetL(triaelem,gauss_coord,DOFPERGRID);

			%Get the Jacobian determinant
			Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

			%Build gaussian vector 
			pe_g_basevert_iceshelf_gaussian=Jdettria*gauss_weight*(-di*(-divHVel+accumulation_gauss-melting_gauss))*L';
			
			%Add Ke_gg_drag_gaussian to Ke
			pe.terms=pe.terms+pe_g_basevert_iceshelf_gaussian;
		end
	end
end %end function

function pe=CreatePVectorVert(pentaelem,grids,materials,inputs);

	%some variables
	numgrids=6;
	NDOF1=1;
	NDOF2=2;

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

	%recover extra inputs
	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
	if ~velocity_is_present,
		error('CreatePVector error message: no input velocity!');
	end

	%initialize velocity array
	vxvy_list=zeros(numgrids,2);

	%Get all element grid data:
	xyz_list=getgriddata(pentaelem,grids);

	%Build row indices for elementary vector.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist;
		for j=1:NDOF2,
			dof=doflist(j);
			vxvy_list(i,j)=velocity_param(dof);
		end
		dof=doflist(3);
		pe.row_indices(i)=dof;
	end

	% Get gaussian points and weights 
	area_order=2;
	num_vert_gauss=2;
	[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

	%Start  looping on the number of gaussian points:
	for igarea=1:num_area_gauss,
	    for igvert=1:num_vert_gauss,
		%Pick up the gaussian point and its weight:
		gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
		gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];

		%Get velocity derivative, with respect to x and y: 
		dudx=GetParameterDerivativeValue(pentaelem,vxvy_list(:,1),xyz_list,gauss_coord); dudx=dudx(1);
		dvdy=GetParameterDerivativeValue(pentaelem,vxvy_list(:,2),xyz_list,gauss_coord); dvdy=dvdy(2);

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

		%Get nodal functions
		l1l6=GetNodalFunctions(pentaelem,gauss_coord);

	    pe.terms=pe.terms+l1l6'*(dudx+dvdy)*Jdet*gauss_weight;

	    end %for ig=1:num_area_gauss,
	end %for ig=1:num_vert_gauss,

end %end function

function pe=CreatePVectorThermal(pentaelem,grids,materials,inputs)

	%some variables
	NDOF1=1;
	numgrids=6;

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

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

	rho_ice=matpar.rho_ice;
	rho_water=matpar.rho_water;
	gravity=matpar.g;
	heatcapacity=matpar.heatcapacity;
	beta=matpar.beta;
	meltingpoint=matpar.meltingpoint;

	%Get all element grid data:
	xyz_list=getgriddata(pentaelem,grids);

	%recover extra inputs
	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
	[flow_law_param flow_law_is_present]=recover_input(inputs,'B');
	[basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
	[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
	[surface_param surface_is_present]=recover_input(inputs,'surface');
	[bed_param bed_is_present]=recover_input(inputs,'bed');
	[temperature_param temperature_is_present]=recover_input(inputs,'temperature');
	[dt dt_is_present]=recover_input(inputs,'dt');
	[pressure_param pressure_is_present]=recover_input(inputs,'pressure');
	
	%we need velocities to compute thermal profiles (even if it is a zero 
	%vector). 
	if ~velocity_is_present,
		error('CreatePVectorThermal error message: input velocity not present!');
	end
	if ~dt_is_present & ~pentaelem.thermal_steadystate,
		error('CreatePVectorThermal error message: transient thermal solution requires present of time step as parameter!');
	end
	if ~temperature_is_present & ~pentaelem.thermal_steadystate,
		error('CreatePVectorThermal error message: transient thermal solution requires temperature at previous iteration!');
	end
	if ~pressure_is_present,
		error('CreatePVectorThermal error message: thermal solution requires pressure');
	end



	%initialize vxvyvz_list
	vxvyvz_list=zeros(numgrids,3); 
	B_list=zeros(numgrids,1);
	K_list=zeros(numgrids,1);
	thickness_list=zeros(numgrids,1);
	surface_list=zeros(numgrids,1);
	bed_list=zeros(numgrids,1);
	temperature_list=zeros(numgrids,1);
	pressure_list=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	
		dof=doflist(1);
		pe.row_indices(i)=dof;
		if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
		if(basal_drag_is_present), K_list(i)=basal_drag_param(dof);end;
		if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
		if(surface_is_present) surface_list(i)=surface_param(dof);end;
		if(bed_is_present) bed_list(i)=bed_param(dof);end;
		if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
		if (pressure_is_present) pressure_list(i)=pressure_param(dof); end;
	end
	vxvy_list= vxvyvz_list(:,1:2);

	% 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.
	area_order=2;
	num_vert_gauss=3;
	[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

	%Start  looping on the number of gaussian points:
	for igarea=1:num_area_gauss,
		for igvert=1:num_vert_gauss,
			%Pick up the gaussian point and its weight:
			gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
			gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
		   
			%Update material if temperature is provided
			if temperature_is_present,
				temperature=GetParameterValue(pentaelem,temperature_list,gauss_coord);
				matice.B=paterson(temperature);
			end
 
			%Build Deformational heating
			epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
			viscosity=GetViscosity3d(matice,epsilon); 

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

			epsilon_e=effective_value(epsilon_matrix);

			phi=2*viscosity*epsilon_e^2;

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

			%Get nodal functions
			l1l6=GetNodalFunctions(pentaelem,gauss_coord);

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

			if pentaelem.thermal_steadystate,
				pe_gaussian_deformational=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight*l1l6';
			else
				pe_gaussian_deformational=dt*phi/(rho_ice*heatcapacity)*Jdet*gauss_weight*l1l6';
			end

			%transient
			if pentaelem.thermal_steadystate,
				pe_gaussian_transient=0;
			else
				pe_gaussian_transient=temperature*Jdet*gauss_weight*l1l6';
			end

			%*Add pe_gaussian vector to pe:
			pe.terms=pe.terms+pe_gaussian_deformational+pe_gaussian_transient;

		end %for ig=1:vert_gauss,
	end %for ig=1:num_area_gauss,

	%Ice/ocean heat exchange flux on ice shelf base 
	if (pentaelem.onbed & pentaelem.shelf),

		%recover material parameters
		mixed_layer_capacity=matpar.mixed_layer_capacity;
		thermal_exchange_velocity=matpar.thermal_exchange_velocity;

		%pressure of the 3 active grids
		pressure_tria=pressure_list(1:3);

		% Get gaussian points and weights.
		[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

		for ig=1:num_gauss2D,

			%Pick up the gaussian point and its weight:
			gauss_weight=gauss_weights(ig);
			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
			
			%Get the Jacobian determinant
			xyz_list_tria=xyz_list(1:3,:);
			Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
			
			%Get nodal functions value
			l1l2l3=GetNodalFunctions(triaelem,gauss_coord);

			%Calculate pressure melting point at gaussian point
			pressure_g=GetParameterValue(triaelem,pressure_tria,gauss_coord);
			t_pmp=meltingpoint-beta*pressure_g;

			%Compute the elementary load vector
			if pentaelem.thermal_steadystate,
				pe_ocean_gaussian=gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity*t_pmp/(heatcapacity*rho_ice)*l1l2l3';
			else
				pe_ocean_gaussian=dt*gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity*t_pmp/(heatcapacity*rho_ice)*l1l2l3';
			end

			%Add Ke_gg_drag_gaussian to Ke
			pe.terms(1:3)=pe.terms(1:3)+pe_ocean_gaussian;
		end
	end%ends ice/ocean exchange flux 

	%Geothermal flux on ice sheet base and basal friction G+tau_b.ub
	if (pentaelem.onbed & ~pentaelem.shelf),

		%Plug the grids on bed into the triaelem
		geothermalflux=pentaelem.geothermalflux(1:3,1);

		%Get the basal friction tau_b.u_b=alpha2*ub^2
		basal_friction=zeros(3,1);
		if (pentaelem.friction_type==2),  %friction_type=2 implies viscous, friction_type=1 implies plastic

			%Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
			frictionparameters=struct();
			frictionparameters.element_type='3d';
			frictionparameters.rho_ice=rho_ice;
			frictionparameters.rho_water=rho_water;
			frictionparameters.g=gravity;	
			frictionparameters.p=pentaelem.p;
			frictionparameters.q=pentaelem.q;
			frictionparameters.velocities=vxvy_list(1:3,:);

			if thickness_is_present,
				frictionparameters.h=thickness_list(1:3,:);
			else
				frictionparameters.h=pentaelem.h(1:3,:);
			end
			if(bed_is_present),
				frictionparameters.b=bed_list(1:3,:);
			else
				frictionparameters.b=pentaelem.b(1:3,:);
			end
			if basal_drag_is_present,
				frictionparameters.k=K_list(1:3,:);
			else
				frictionparameters.k=pentaelem.k(1:3,:);
			end

			alpha2=Getalpha2(frictionparameters);
	
			%We need the velocity magnitude to evaluate the basal stress:
			u_b2=vxvy_list(1:3,1).^2+vxvy_list(1:3,2).^2;

			%Compute basal friction
			basal_friction=alpha2.*u_b2;
		end

		% Get gaussian points and weights.
		[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

		for ig=1:num_gauss2D,

			%Pick up the gaussian point and its weight:
			gauss_weight=gauss_weights(ig);
			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
			
			%Get the Jacobian determinant
			xyz_list_tria=xyz_list(1:3,:);
			Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
			
			%Get geothermal flux and basal friction
			G_g=GetParameterValue(triaelem,geothermalflux,gauss_coord);
			basal_friction_g=GetParameterValue(triaelem,basal_friction,gauss_coord);
			
			%Get nodal functions value
			l1l2l3=GetNodalFunctions(triaelem,gauss_coord);

			%Compute the elementary load vector
			if pentaelem.thermal_steadystate,
				pe_geo_gaussian=gauss_weight*Jdettria*(basal_friction_g+G_g)/(heatcapacity*rho_ice)*l1l2l3';
			else
				pe_geo_gaussian=dt*gauss_weight*Jdettria*(basal_friction_g+G_g)/(heatcapacity*rho_ice)*l1l2l3';
			end

			%Add Ke_gg_drag_gaussian to Ke
			pe.terms(1:3)=pe.terms(1:3)+pe_geo_gaussian;
		end
	end%ends geothermal 


end %end function


function pe=CreatePVectorHoriz(pentaelem,grids,materials,inputs,analysis_type);

	%Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 
	%bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
	%the load vector

	if pentaelem.collapse & ~pentaelem.onbed

		pe=elemvector(0);
		return;

	elseif pentaelem.collapse & pentaelem.onbed

		%we collapse the penta into its base tria, and use the tria to build the stiffness matrix.
		element=triaelem;
		element.type='triaelem';
		element.id=NaN; %not needed here, we are going to destroy this triaelem very soon.
		element.matid=pentaelem.matid; %same materials, not dependent on element.
		element.g=pentaelem.g(1:3); %first three grids correspond to the base of the penta element.
		element.h=pentaelem.h(1:3);
		element.s=pentaelem.s(1:3);
		element.b=pentaelem.b(1:3);
		element.friction_type=pentaelem.friction_type;
		element.k=pentaelem.k(1:3);
		element.p=pentaelem.p;
		element.q=pentaelem.q;
		element.shelf=pentaelem.shelf;
		element.meanvel=pentaelem.meanvel;
		element.epsvel=pentaelem.epsvel;
		element.acceleration=pentaelem.acceleration;

		%Call CreateKMatrix for this new element.
		pe=CreatePVector(element,grids,materials,inputs,analysis_type);

		%That's it, element will get destroyed upon return, and we have the correct 
		%collapsed stiffness matrix. 
		return;

	else 

		%implement standard penta element.

			
		%some variables
		numgrids=6;
		numgrids2d=3;
		NDOF2=2;

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

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

		%recover material parameters
		gravity=matpar.g;
		rho_ice=matpar.rho_ice;
		rho_water=matpar.rho_water;

		%recover extra inputs
		[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
		[surface_param surface_is_present]=recover_input(inputs,'surface');
		[bed_param bed_is_present]=recover_input(inputs,'bed');
		[melting_param melting_is_present]=recover_input(inputs,'melting');
		[basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
	
		%Get all element grid data:
		xyz_list=getgriddata(pentaelem,grids);
		
		%Initialize inputs
		basal_drag_list=zeros(numgrids,1);
		melting_list=zeros(numgrids,1);
		thickness_list=zeros(numgrids,1);
		surface_list=zeros(numgrids,1);
		bed_list=zeros(numgrids,1);
		temperature_list=zeros(numgrids,1);

		%Build row indices for elementary vector.
		for i=1:numgrids,
			
			doflist=grids(pentaelem.g(i)).grid.doflist;
			for j=1:NDOF2,
				pe.row_indices((i-1)*NDOF2+j)=doflist(j);
			end
			dof=doflist(1);
			if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
			if(surface_is_present) surface_list(i)=surface_param(dof);end;
			if(bed_is_present) bed_list(i)=bed_param(dof);end;
			if(melting_is_present) melting_list(i)=melting_param(dof);end;
			if(basal_drag_is_present) basal_drag_list(i)=basal_drag_param(dof);end;

		end

%We consider that we have hydrostatic equilibrium for Pattyn's model, so no reason to have a water pressure ... to decide !!!
%		%Deal with water pressure integration: 
%
%		if (pentaelem.onbed==1 & pentaelem.shelf==1),
%			
%			%Create the element at the upper base of pentaelem, which is a triangle
%			triabed=triaelem;
%
%			%Plug the grids on surface into the triaelem
%			triabed.g=pentaelem.g(1:3);
%			if thickness_is_present,
%				triabed.h=thickness_list(1:3);
%			else			
%				triabed.h=pentaelem.h(1:3);
%			end
%			if bed_is_present,
%				triabed.b=bed_list(1:3);
%			else			
%				triabed.b=pentaelem.b(1:3);
%			end
%			if basal_drag_is_present,
%				triabed.k=basal_drag_list(1:3);
%			else			
%				triabe.k=pentaelem.k(1:3);
%			end
%
%
%			%Get the coordinates of the three grids
%			xyz_list_tria=xyz_list(1:3,:);
%
%			% Get gaussian points and weights.
%			[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
%
%			%Build normal vector to the surface
%			n=zeros(2,1);
%			normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(2,:)-xyz_list_tria(3,:)));
%			n(1)=1/norm(normal)*normal(1);
%			n(2)=1/norm(normal)*normal(2);
%
%			for ig=1:num_gauss2D,
%
%				%Pick up the gaussian point and its weight:
%				gauss_weight=gauss_weights(ig);
%				gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
%				
%				%Get the Jacobian determinant
%				Jdettria=GetJacobianDeterminant(triabed,xyz_list_tria,gauss_coord);
%
%				%Get thickness and bed at gaussian point
%				thickness_gg=GetParameterValue(triabed,triabed.h,gauss_coord);
%				bed_gg=GetParameterValue(triabed,triabed.b,gauss_coord);
%				
%				%Get L matrix 
%				L=GetL(triabed,gauss_coord,1);
%		   
%				%Calculate pressure
%				pressure=gravity*(rho_ice*thickness_gg+rho_water*bed_gg)*gauss_weight*Jdettria;
%
%				%Calculate elementary load vector
%				for i=1:numgrids2d,
%				    for j=1:NDOF2,
%					pe_gg_water((i-1)*NDOF2+j)=pressure*L(i)*n(j);
%				    end
%				end
%
%				%Add pe_gaussian vector to pe:
%				for i=1:numgrids2d*NDOF2,
%					pe.terms(i)=pe.terms(i)+pe_gg_water(i);
%				end
%			end
%		end%ends of water pressure

		% Get gaussian points and weights 
		area_order=2;
		num_vert_gauss=3;
		[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

		%for i=1:num_gauss,
		%	disp(sprintf('Gauss coord %i: %f %f %f Weight: %f\n',i,first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
		%end

		%Start  looping on the number of gaussian points:
		for igarea=1:num_area_gauss,
			for igvert=1:num_vert_gauss,
				%Pick up the gaussian point and its weight:
				gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
				gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
			    
				%Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
				if thickness_is_present,
					thickness=GetParameterValue(pentaelem,thickness_list,gauss_coord);
				else
					thickness=GetParameterValue(pentaelem,pentaelem.h,gauss_coord);
				end
				%disp(sprintf('Thickness: %f\n', thickness));

				if surface_is_present,
					slope=GetParameterDerivativeValue(pentaelem,surface_list,xyz_list,gauss_coord);
				else
					slope=GetParameterDerivativeValue(pentaelem,pentaelem.s,xyz_list,gauss_coord);
				end
				%disp(sprintf('Slope: sx %f sy %f\n',slope(1),slope(2)));

				%In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 
				%element itself: 
				if (~pentaelem.shelf) & (pentaelem.friction_type==1),
					if basal_drag_is_present,
						plastic_stress=GetParameterValue(pentaelem,basal_drag_list,gauss_coord);
					else			
						plastic_stress=GetParameterValue(pentaelem,pentaelem.k,gauss_coord);
					end
				end

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

				%Get nodal functions
				l1l6=GetNodalFunctions(pentaelem,gauss_coord);

				%Compute driving stress: 
				driving_stress_baseline=rho_ice*gravity;

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

				if(~pentaelem.shelf) & (pentaelem.friction_type==1),
					for i=1:numgrids,
						for j=1:NDOF2,
							pe_gaussian((i-1)*NDOF2+j)=(-driving_stress_baseline*slope(j)-plastic_stress)*Jdet*gauss_weight*l1l6(i);
						end
					end
				else
					for i=1:numgrids,
						for j=1:NDOF2,
							pe_gaussian((i-1)*NDOF2+j)=-driving_stress_baseline*slope(j)*Jdet*gauss_weight*l1l6(i);
						end
					end
				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,
		end %for ig=1:num_vert_gauss,
	end % elseif pentaelem.collapse & pentaelem.onbed
end %end function




function pe=CreatePVectorPrognostic(pentaelem,grids,materials,inputs);

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.
	if pentaelem.onbed~=1,
		pe=elemvector(0);
		return;
	end

	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%Create elementary stiffness matrix 
	pe=elemvector(numdof);

	%Get all element grid data:
	xyz_list=getgriddata(pentaelem,grids); 

	%Just keep the first 3 grids
	xyz_list=xyz_list(1:3,:);

	%recover extra inputs from users, at current convergence iteration.
	[accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
	[melting_param melting_is_present]=recover_input(inputs,'melting');
	[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
	[dt dt_is_present]=recover_input(inputs,'dt');

	%check on all parameters
	if ( ~accumulation_is_present | ~melting_is_present | ~thickness_is_present | ~dt_is_present),
		error('CreateKMatrixPrognostic error message: missing input parameters!');
	end
		
	%initialize parameter lists (for our 3 grids)
	thickness_list=zeros(numgrids,1);
	accumulation_list=zeros(numgrids,1);
	melting_list=zeros(numgrids,1);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		dof=doflist(1); %first degree of freedom
		
		pe.row_indices(i)=dof;
	
		%recover thickness_list on first dof of every grid.
		thickness_list(i)=thickness_param(dof);
		melting_list(i)=melting_param(dof);
		accumulation_list(i)=accumulation_param(dof);
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	for ig=1:num_gauss2D,
	
		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		
		%Get melting at gaussian point
		melting=GetParameterValue(triaelem,melting_list,gauss_coord);
	
		%Get accumulation at gaussian point
		accumulation=GetParameterValue(triaelem,accumulation_list,gauss_coord);

		%Get thickness at gaussian point
		thickness=GetParameterValue(triaelem,thickness_list,gauss_coord);

		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L: 
		L=GetL(triaelem,gauss_coord,DOFPERGRID);

		%Build gaussian vector 
		pe_g_gaussian=Jdettria*gauss_weight*(thickness+dt*(accumulation-melting))*L';

		%Add Ke_gg_drag_gaussian to Ke
		pe.terms=pe.terms+pe_g_gaussian;
	end

end %end function
