function [pengrid unstable]=PenaltyConstrain(pengrid,grids,materials,inputs)
%PENALTYCONSTRAIN - indicate if a penalty on a grid is stable
%
%   The penalty is stable if it doesn't change during to successive iterations.   
%
%   Usage:
%      [pengrid unstable]=PenaltyConstrain(pengrid,grids,materials,inputs)
% 
%   See also PENALTYCREATEKMATRIX, PENALTYCREATEKMATRIX

%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');

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

%Recover pressure and temperature at the grid
doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
dof=doflist(1);

pressure=pressure_param(dof);
temperature=temperature_param(dof);

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

t_pmp=meltingpoint-beta*pressure;

%Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
if temperature>t_pmp,
	active=1;
else
	active=0;
end

%Figure out stability of this penalty
if pengrid.active==active,
	unstable=0;
else
	unstable=1;
end
	
%Set penalty flag
pengrid.active=active;
