%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;
		requested_outputs = {};
	end
	methods
		function obj = thermal(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function list = defaultoutputs(self,md) % {{{

			if self.isenthalpy,
				list = {'Enthalpy','Temperature','Waterfraction','Watercolumn','BasalforcingsMeltingRate'};
			else
				list = {'Temperature','BasalforcingsMeltingRate'};
			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;

			%default output
			obj.requested_outputs={'default'};
		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,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2]);
			md = checkfield(md,'fieldname','thermal.spctemperature','forcing',1);
			if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy & strcmp(meshtype(md.mesh),'3D')),
				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,'fieldname','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,'fieldname','thermal.isenthalpy','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[0 1]);
				if(md.thermal.isenthalpy)
					md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[1], 'message',['for enthalpy run thermal.isdynamicbasalspc should be 1']);
				end
	    end

		 md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1);
    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,'penalty_factor','scaling exponent (default is 3)');
			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)']);
			fielddisplay(obj,'requested_outputs','additional outputs requested');

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

			%process requested outputs
			outputs = obj.requested_outputs;
			pos  = find(ismember(outputs,'default'));
			if ~isempty(pos),
				outputs(pos) = [];                         %remove 'default' from outputs
				outputs      = [outputs defaultoutputs(obj,md)]; %add defaults
			end
			WriteData(fid,'data',outputs,'enum',ThermalRequestedOutputsEnum(),'format','StringArray');
        	end % }}}
	end
end
