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

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

		Ke=CreateKMatrixThermal(pengrid,grids,materials,inputs,kmax);

	elseif strcmpi(analysis_type,'diagnostic_stokes') 

		Ke=CreateKMatrixStokes(pengrid,grids,materials,inputs,kmax);

	elseif strcmpi(analysis_type,'melting')

		Ke=CreateKMatrixMelting(pengrid,grids,materials,inputs,kmax);

	end
end %end function

function Ke=CreateKMatrixThermal(pengrid,grids,materials,inputs,kmax)

	if ~pengrid.active,
		Ke=elemmatrix(0);
		return;
	else
		
		%initialize stiffness matrix
		Ke=elemmatrix(1);

		%Build linear indices for elementary stiffness matrix.
		doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
		dof=doflist(1);
		Ke.row_indices=dof;
			
		Ke.terms(1,1)=kmax*10^(pengrid.penalty_offset);
	end

end %function

function Ke=CreateKMatrixMelting(pengrid,grids,materials,inputs,kmax)
	
	%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 (~temperature_is_present | ~pressure_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);
	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;

	%Add penalty load
	if temperature>t_pmp %melting
		Ke=elemmatrix(0);

	else                 %no melting
		%If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp
		%instead of using spcs, one uses penalties
		Ke=elemmatrix(1);
		Ke.row_indices=dof;
		Ke.terms(1,1)=kmax*10^pengrid.penalty_offset;
	end	

end %function

function Ke=CreateKMatrixStokes(pengrid,grids,materials,inputs,kmax)

	%First recover beta, pressure and temperature vectors;
	[slope_param slope_is_present]=recover_input(inputs,'slopebed');
	
	if ~slope_is_present,
		error('CreatePVectorStokes error message: missing input parameters');
	end
	
	%Recover slope at the grid and initialize elementary vector dof.
	slope=zeros(1,2);
	doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
	slope(1,1)=slope_param(doflist(1));
	slope(1,2)=slope_param(doflist(2));

	%Create elementary matrix
	Ke=elemmatrix(3);
	for i=1:3,
		Ke.row_indices(i)=doflist(i);
	end

	%Add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
	Ke.terms(3,1)=-slope(1,1)*kmax*10^pengrid.penalty_offset;
	Ke.terms(3,2)=-slope(1,2)*kmax*10^pengrid.penalty_offset;
	Ke.terms(3,3)=1*kmax*10^pengrid.penalty_offset;

end %function
