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

if strcmpi(analysis_type,'diagnostic_horiz'),

	pe=CreatePVectorHoriz(triaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'diagnostic_vert'),

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

elseif strcmpi(analysis_type,'temperature'),

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

elseif strcmpi(analysis_type,'prognostic'),

	pe=CreatePVectorPrognostic(triaelem,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(triaelem,grids,materials,inputs,analysis_type);

end
end %end function

function pe=CreatePVectorVert(triaelem,grids,materials,inputs)
error('CreatePVector error message: Vertical analysis not implemented yet');
end %end function

function pe=CreatePVectorThermal(triaelem,grids,materials,inputs)
error('CreatePVector error message: Thermal analysis not implemented yet');
end %end function

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

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

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

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

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

	%recover extra inputs
	if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') ,
		[param param_is_present]=recover_input(inputs,'bed');
	elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
		[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,
		doflist=grids(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
		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=triaelem.b(1:3); 
		elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
			param=triaelem.s(1:3); 
		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=CreatePVectorHoriz(triaelem,grids,materials,inputs);

	global element_debug  element_debugid

	%Check if MacAyeal's acceleration is not on? 
	if triaelem.acceleration==1,
		pe={};
	end

	%some variables
	numgrids=3;
	NDOF2=2;

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

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

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

	%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');
	[basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
		
	%Get all element grid data:
	xyz_list=getgriddata(triaelem,grids);

	%initialize extra inputs
	K_list=zeros(numgrids,1);
	thickness_list=zeros(numgrids,1);
	surface_list=zeros(numgrids,1);
	bed_list=zeros(numgrids,1);

	%Build row indices for elementary vector.
	for i=1:numgrids,
		doflist=grids(triaelem.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);
		else
			thickness_list(i)=triaelem.h(i);
		end
		if(surface_is_present),
			surface_list(i)=surface_param(dof);
		else
			surface_list(i)=triaelem.s(i);
		end;
		if(bed_is_present),
			bed_list(i)=bed_param(dof);
		else
			bed_list(i)=triaelem.b(i);
		end
		if(basal_drag_is_present),
			K_list(i)=basal_drag_param(dof);
		else
			K_list(i)=triaelem.k(i);
		end
	end

	if (element_debug & triaelem.id==element_debugid),
		disp(sprintf('gravity %g\n',gravity));
		disp(sprintf('rho_ice %g\n',rho_ice));
		disp(sprintf('thickness_list (%g,%g,%g)\n',thickness_list(1),thickness_list(2),thickness_list(3)));
		disp(sprintf('surface_list (%g,%g,%g)\n',surface_list(1),surface_list(2),surface_list(3)));
		disp(sprintf('bed_list (%g,%g,%g)\n',bed_list(1),bed_list(2),bed_list(3)));
		disp(sprintf('K_list (%g,%g,%g)\n',K_list(1),K_list(2),K_list(3)));
	end

	% Get gaussian points and weights 
	[num_gauss,first_gauss_area_coord,second_gauss_area_coord,third_gauss_area_coord,gauss_weights]=GaussTria(2);

	if (element_debug & triaelem.id==element_debugid),
		disp(sprintf('   gaussian points: '));
		for i=1:num_gauss,
			disp(sprintf('    %g %g %g : %g',first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
		end
	end
 
	%Start  looping on the number of gaussian points:
	for ig=1:num_gauss,
		%Pick up the gaussian point:
		gauss_weight=gauss_weights(ig);
		gauss_l1l2l3=[first_gauss_area_coord(ig) second_gauss_area_coord(ig) third_gauss_area_coord(ig)];

		%Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
		thickness=GetParameterValue(triaelem,thickness_list,gauss_l1l2l3);

		slopevector=GetParameterDerivativeValue(triaelem,surface_list,xyz_list,gauss_l1l2l3);

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

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

		%Get nodal functions
		l1l2l3=GetNodalFunctions(triaelem,gauss_l1l2l3);

		%Compute driving stress: 
		driving_stress_baseline=rho_ice*gravity*thickness;
	
		if (element_debug & triaelem.id==element_debugid),
			disp(sprintf('      gaussian %i\n',ig));
			disp(sprintf('      thickness %g\n',thickness));
			disp(sprintf('      slope(%g,%g)\n',slopevector(1),slopevector(2)));
			disp(sprintf('      Jdet %g\n',Jdet));
			disp(sprintf('      gaussweigth %g\n',gauss_weight));
			disp(sprintf('      l1l2l3 (%g,%g,%g)\n',l1l2l3(1),l1l2l3(2),l1l2l3(3)));
			if(triaelem.friction_type==1), disp(sprintf('      plastic_stress(%g)\n',plastic_stress));end;
		end

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

		if(~triaelem.shelf) & (triaelem.friction_type==1),
			for i=1:numgrids,
				for j=1:NDOF2,
					pe_gaussian((i-1)*NDOF2+j)=(-driving_stress_baseline*slopevector(j)-plastic_stress)*Jdet*gauss_weight*l1l2l3(i); 
				end
			end
		else
			for i=1:numgrids,
				for j=1:NDOF2,
					pe_gaussian((i-1)*NDOF2+j)=-driving_stress_baseline*slopevector(j)*Jdet*gauss_weight*l1l2l3(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_gauss,
		
	if (element_debug & triaelem.id==element_debugid),
		disp(sprintf('      pe_g->terms\n',ig));
		pe.terms
		disp(sprintf('      pe_g->row_indices\n',ig));
		pe.row_indices
	end

end %end function

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

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

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

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

	%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(triaelem.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
