%INITIALIZATION class definition
%
%   Usage:
%      initialization=initialization();

classdef initialization
	properties (SetAccess=public) 
		vx            = NaN;
		vy            = NaN;
		vz            = NaN;
		vel           = NaN;
		pressure      = NaN;
		temperature   = NaN;
		waterfraction = NaN;
		sediment_head = NaN;
		epl_head      = NaN;
		epl_thickness = NaN;
		watercolumn   = NaN;
	end
	methods
		function obj = initialization(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{
			if ismember(StressbalanceAnalysisEnum(),analyses)
				if ~(isnan(md.initialization.vx) | isnan(md.initialization.vy)),
					md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
					md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
			end
			if ismember(MasstransportAnalysisEnum(),analyses),
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
			end
			if ismember(BalancethicknessAnalysisEnum(),analyses) & solution==BalancethicknessSolutionEnum,
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
				%Triangle with zero velocity
				if any(sum(abs(md.initialization.vx(md.mesh.elements)),2)==0 & sum(abs(md.initialization.vy(md.mesh.elements)),2)==0)
					md = checkmessage(md,'at least one triangle has all its vertices with a zero velocity');
				end
			end
			if ismember(ThermalAnalysisEnum(),analyses),
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
				if strcmp(meshtype(md.mesh),'3D')
					md = checkfield(md,'fieldname','initialization.vz','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
				md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]);
			end
			if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy)
				md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices 1]);
			end
			if ismember(HydrologyShreveAnalysisEnum(),analyses),
				if isa(md.hydrology,'hydrologyshreve'),
					md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
			end
			if ismember(HydrologyDCInefficientAnalysisEnum(),analyses),
				if isa(md.hydrology,'hydrologydc'),
					md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
			end
			if ismember(HydrologyDCEfficientAnalysisEnum(),analyses),
				if isa(md.hydrology,'hydrologydc'),
					if md.hydrology.isefficientlayer==1,
						md = checkfield(md,'fieldname', ...
														'initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
						md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'size',[md.mesh.numberofvertices 1]);
					end
				end
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   initial field values:'));

			fielddisplay(obj,'vx','x component of velocity [m/yr]');
			fielddisplay(obj,'vy','y component of velocity [m/yr]');
			fielddisplay(obj,'vz','z component of velocity [m/yr]');
			fielddisplay(obj,'vel','velocity norm [m/yr]');
			fielddisplay(obj,'pressure','pressure field [Pa]');
			fielddisplay(obj,'temperature','temperature [K]');
			fielddisplay(obj,'waterfraction','fraction of water in the ice');
			fielddisplay(obj,'sediment_head','sediment water head of subglacial system [m]');
			fielddisplay(obj,'epl_head','epl water head of subglacial system [m]');
			fielddisplay(obj,'epl_thickness','epl layer thickness [m]');
			fielddisplay(obj,'watercolumn','thickness of subglacial water [m]');

		end % }}}
		function marshall(obj,md,fid) % {{{

			yts=365.0*24.0*3600.0;

			WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum,'scale',1./yts);
			WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum,'scale',1./yts);
			WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum,'scale',1./yts);
			WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum);
			WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum);
			WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum);
			WriteData(fid,'data',obj.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum);
			WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum);
			WriteData(fid,'data',obj.epl_thickness,'format','DoubleMat','mattype',1,'enum',HydrologydcEplThicknessEnum);
			WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);

			if md.thermal.isenthalpy,
				tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
				pos  = find(md.initialization.temperature>tpmp);
				enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
				enthalpy(pos) = md.materials.heatcapacity*tpmp(pos) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction(pos);
				WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum());
			end
		end % }}}
	end
end
