%THERMAL class definition
%
%   Usage:
%      thermal=thermal();

classdef thermal
	properties (SetAccess=public) 
		spctemperature    = NaN;
		penalty_threshold = 0;
		stabilization     = 0;
		maxiter           = 0;
		penalty_lock      = 0;
		penalty_factor    = 0;
		isenthalpy        = 0;
	end
	methods
		function obj = thermal(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					in=varargin{1};
					if (isa(in,'numeric') & in==0),
						% requesting templates do nothing
					else
						error('constructor not supported');
					end 
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			%Number of unstable constraints acceptable
			obj.penalty_threshold=0;

			%Type of stabilization used
			obj.stabilization=1;

			%Maximum number of iterations
			obj.maxiter=100;

			%factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
			obj.penalty_factor=3;

			%Should we use cold ice (default) or enthalpy formulation
			obj.isenthalpy=0;
		end % }}}
		function checkconsistency(obj,md,solution,analyses) % {{{

			%Early return
			if (~ismember(ThermalAnalysisEnum,analyses) & ~ismember(EnthalpyAnalysisEnum,analyses)) | (solution==TransientSolutionEnum & md.transient.isthermal==0), return; end

			checkfield(md,'thermal.stabilization','numel',1,'values',[0 1 2]);
			checkfield(md,'thermal.spctemperature','forcing',1);
			if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy),
			checkfield(md,'thermal.spctemperature','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*md.geometry.thickness,'message','spctemperature should be below the adjusted melting point');
			checkfield(md,'thermal.isenthalpy','numel',1,'values',[0 1]);
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   Thermal solution parameters:'));

			fielddisplay(obj,'spctemperature','temperature constraints (NaN means no constraint)');
			fielddisplay(obj,'stabilization','0->no, 1->artificial_diffusivity, 2->SUPG');
			fielddisplay(obj,'maxiter','maximum number of non linear iterations');
			fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
			fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)');
			fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');

		end % }}}
		function marshall(obj,fid) % {{{
			WriteData(fid,'object',obj,'fieldname','spctemperature','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
			WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
			WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
			WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
			WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
			WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean');
		end % }}}
	end
end
