%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 createxml(obj,fid) % {{{
            fprintf(fid, '<!-- thermal -->\n');            
                    
            % thermal solution parameters
            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Thermal solution parameters">','<section name="thermal" />');                    
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="spctemperature" type="',class(obj.spctemperature),'" default="',convert2str(obj.spctemperature),'">','     <section name="thermal" />','     <help> temperature constraints (NaN means no constraint) [K] </help>','  </parameter>');
                
            % penalty_threshold drop-down (0, 1, or 2)
            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="penalty_threshold" type="alternative"  optional="false">','     <section name="thermal" />','     <help> 0: no, 1: artificial_diffusivity, 2: SUPG </help>');
            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
            fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
            fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');

            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="',class(obj.stabilization),'" default="',convert2str(obj.stabilization),'">','     <section name="thermal" />','     <help> maximum number of non linear iterations </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter" type="',class(obj.maxiter),'" default="',convert2str(obj.maxiter),'">','     <section name="thermal" />','     <help> stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_lock" type="',class(obj.penalty_lock),'" default="',convert2str(obj.penalty_lock),'">','     <section name="thermal" />','     <help> threshold to declare convergence of thermal solution (default is 0)  </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_factor" type="',class(obj.penalty_factor),'" default="',convert2str(obj.penalty_factor),'">','     <section name="thermal" />','     <help> scaling exponent (default is 3) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isenthalpy" type="',class(obj.isenthalpy),'" default="',convert2str(obj.isenthalpy),'">','     <section name="thermal" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdynamicbasalspc" type="',class(obj.isdynamicbasalspc),'" default="',convert2str(obj.isdynamicbasalspc),'">','     <section name="thermal" />','     <help> enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)  </help>','  </parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(obj.requested_outputs),'" default="',convert2str(obj.requested_outputs),'">','     <section name="thermal" />','     <help> additional outputs requested </help>','  </parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');    
        
        end % }}}
		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','BasalforcingsGroundediceMeltingRate'};
			else
				list = {'Temperature','BasalforcingsGroundediceMeltingRate'};
			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 & dimension(md.mesh)==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,'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
