function pe=PenaltyCreatePVector(pengrid,grids,materials,inputs,analysis_type,kmax);
%PENALTYCREATEPVECTOR - creates a penalty vector for a particular grid
%
%   This routine applies for thermal, melting, and Stokes solution.
%
%   Usage:
%      pe=PenaltyCreatePVector(pengrid,grids,materials,inputs,analysis_type,kmax)
% 
%   See also CREATEPVECTOR, PENALTYCREATEKMATRIX, PENALTYCONTRAIN

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

		pe=PenaltyCreatePVectorThermal(pengrid,grids,materials,inputs,kmax);

	elseif strcmpi(analysis_type,'diagnostic_stokes')

		pe=PenaltyCreatePVectorDiagnosticStokes(pengrid,grids,materials,inputs,kmax);

	elseif strcmpi(analysis_type,'melting')

		pe=PenaltyCreatePVectorMelting(pengrid,grids,materials,inputs,kmax);

	end
end %end function

function pe=PenaltyCreatePVectorDiagnosticStokes(pengrid,grids,materials,inputs,kmax);
	
	pe=elemvector(0);

end %end function

function pe=PenaltyCreatePVectorThermal(pengrid,grids,materials,inputs,kmax);
	
	if ~pengrid.active,
		pe=elemvector(0);
	else
		
		%initialize load vector
		pe=elemvector(1);

		%First recover pressure
		[pressure_param pressure_is_present]=recover_input(inputs,'pressure');

		if (~pressure_is_present),
			error('CreatePVectorThermal error message: missing pressure as input');
		end

		%Recover pressure  at the grid and initialize elementary vector dof.
		doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
		dof=doflist(1);
		pe.row_indices=dof;
		pressure=pressure_param(dof);

		%Compute pressure melting point
		matpar=materials(end).constants;
		meltingpoint=matpar.meltingpoint;
		beta=matpar.beta;
		t_pmp=meltingpoint-beta*pressure;

		%Add penalty load
		pe.terms(1)=kmax*10^(pengrid.penalty_offset)*t_pmp;
	end
end %end function


function pe=PenaltyCreatePVectorMelting(pengrid,grids,materials,inputs,kmax);
			
	%initialize load vector
	pe=elemvector(1);

	%First recover beta, pressure and temperature vectors;
	[pressure_param pressure_is_present]=recover_input(inputs,'pressure');
	[temperature_param temperature_is_present]=recover_input(inputs,'temperature');
	[melting_offset melting_offset_is_present]=recover_input(inputs,'melting_offset');
	if ~pengrid.thermal_steadystate,
		[dt dt_is_present]=recover_input(inputs,'dt');
		if ~dt_is_present
			error('CreatePVectorMelting error message: missing dt');
		end
	end

	if (~temperature_is_present | ~pressure_is_present | ~melting_offset_is_present),
		error('CreatePVectorMelting error message: missing input parameters');
	end

	%Recover pressure  at the grid and initialize elementary vector dof.
	doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
	dof=doflist(1);
	pe.row_indices=dof;
	pressure=pressure_param(dof);
	temperature=temperature_param(dof);

	%Recover parameters
	matpar=materials(end).constants;
	meltingpoint=matpar.meltingpoint;
	beta=matpar.beta;
	heatcapacity=matpar.heatcapacity;
	latentheat=matpar.latentheat;

	%Compute pressure melting point
	t_pmp=meltingpoint-beta*pressure;

	%Add penalty load
	%This time, the penalty must have the same value as the one used for the thermal computation
	%so that the corresponding melting can be computed correctly
	%In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
	if temperature<t_pmp %no melting
		pe.terms(1)=0;
	else
		if pengrid.thermal_steadystate,
			pe.terms(1)=melting_offset*10^(pengrid.penalty_offset)*(temperature-t_pmp);
		else
			pe.terms(1)=melting_offset*10^(pengrid.penalty_offset)*(temperature-t_pmp)/dt;
		end
	end
end %end function
