%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;
        	isdynamicbasalspc = 0;
	end
	methods
		function obj = thermal(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				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;
            
            		%will basal boundary conditions be set dynamically
            		obj.isdynamicbasalspc=0;
		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{

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

			md = checkfield(md,'thermal.stabilization','numel',[1],'values',[0 1 2]);
			md = checkfield(md,'thermal.spctemperature','forcing',1);
			if (ismember(EnthalpyAnalysisEnum(),analyses) & (md.thermal.isenthalpy | solution==EnthalpySolutionEnum()) & md.mesh.dimension==3),
				pos=find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN);
				replicate=repmat(md.geometry.surface-md.mesh.z,1,size(md.thermal.spctemperature,2));
				md = checkfield(md,'thermal.spctemperature(find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN))','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate(pos),'message','spctemperature should be below the adjusted melting point');
				md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0 1]);
                		md = checkfield(md,'thermal.isdynamicbasalspc','numel',[1],'values',[0 1]);
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   Thermal solution parameters:'));

			fielddisplay(obj,'spctemperature','temperature constraints (NaN means no constraint) [K]');
			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)');
            		fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
            
		end % }}}
		function marshall(obj,md,fid) % {{{
			WriteData(fid,'object',obj,'fieldname','spctemperature','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+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');
            		WriteData(fid,'object',obj,'fieldname','isdynamicbasalspc','format','Boolean');
        	end % }}}
	end
end
