%FLOWEQUATION class definition
%
%   Usage:
%      flowequation=flowequation();

classdef flowequation
	properties (SetAccess=public) 
		isSIA            = 0;
		isSSA            = 0;
		isL1L2           = 0;
		isHO             = 0;
		isFS             = 0;
		fe_SSA           = '';
		fe_HO            = '';
		fe_FS            = '';
		vertex_equation  = NaN;
		element_equation = NaN;
		borderSSA        = NaN;
		borderHO         = NaN;
		borderFS         = NaN;
	end
	methods (Static)
		function obj = loadobj(obj) % {{{
			% This function is directly called by matlab when a model object is
			% loaded. If the input is a struct it is an old version of this class and
			% old fields must be recovered (make sure they are in the deprecated
			% model properties)

			if verLessThan('matlab','7.9'),
				disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
				disp('         if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');

				% This is a Matlab bug: all the fields of md have their default value
				% Example of error message:
				% Warning: Error loading an object of class 'model':
				% Undefined function or method 'exist' for input arguments of type 'cell'
				%
				% This has been fixed in MATLAB 7.9 (R2009b) and later versions
			end

			if isstruct(obj)
				disp('Recovering flowequation from older version');
				objstruct = obj;
				obj = structtoobj(flowequation(),objstruct);

				%2013 July 23rd
				if isfield(objstruct,'ishutter'),      obj.isSIA     = objstruct.ishutter;       end; 
				if isfield(objstruct,'ismacayeal'),    obj.isSSA     = objstruct.ismacayeal;     end; 
				if isfield(objstruct,'ispattyn'),      obj.isHO      = objstruct.ispattyn;       end; 
				if isfield(objstruct,'isstokes'),      obj.isFS      = objstruct.isstokes;       end; 
				if isfield(objstruct,'bordermacayeal'),obj.borderSSA = objstruct.bordermacayeal; end; 
				if isfield(objstruct,'borderpattyn'),  obj.borderHO  = objstruct.borderpattyn;   end; 
				if isfield(objstruct,'borderstokes'),  obj.borderFS  = objstruct.borderstokes;   end; 
			end

		end% }}}
	end
	methods
		function obj = flowequation(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			%P1 for SSA
			obj.fe_SSA= 'P1';

			%P1 for HO
			obj.fe_HO= 'P1';

			%MINI condensed element for FS by default
			obj.fe_FS = 'MINIcondensed';
		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{

			if ismember(StressbalanceAnalysisEnum(),analyses),

				md = checkfield(md,'fieldname','flowequation.isSIA','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.isSSA','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.isL1L2','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2'});
				md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4'});
				md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'});
				md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices 1],'values',[0 1]);
				md = checkfield(md,'fieldname','flowequation.borderFS','size',[md.mesh.numberofvertices 1],'values',[0 1]);
				if strcmp(meshtype(md.mesh),'2Dhorizontal')
					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[1:2]);
					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[1:2]);
				elseif strcmp(meshtype(md.mesh),'2Dvertical')
					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[2:4]);
					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[2:4]);
				elseif strcmp(meshtype(md.mesh),'3D'),
					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[0:8]);
					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[0:8]);
				else
					error('Case not supported yet');
				end
				if ~(obj.isSIA || obj.isSSA || obj.isL1L2 || obj.isHO || obj.isFS),
					md = checkmessage(md,['no element types set for this model']);
				end
			end
			if ismember(StressbalanceSIAAnalysisEnum(),analyses),
				if any(obj.element_equation==1),
					if(obj.vertex_equation & md.mask.groundedice_levelset<0.),
						disp(sprintf('\n !!! Warning: SIA''s model is not consistent on ice shelves !!!\n'));
					end
				end
			end

		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   flow equation parameters:'));

			fielddisplay(obj,'isSIA','is the Shallow Ice Approximation (SIA) used ?');
			fielddisplay(obj,'isSSA','is the Shelfy-Stream Approximation (SSA) used ?');
			fielddisplay(obj,'isL1L2','is the L1L2 approximation used ?');
			fielddisplay(obj,'isHO','is the Higher-Order (HO) approximation used ?');
			fielddisplay(obj,'isFS','are the Full-FS (FS) equations used ?');
			fielddisplay(obj,'fe_SSA','Finite Element for SSA  ''P1'', ''P1bubble'' ''P1bubblecondensed'' ''P2''');
			fielddisplay(obj,'fe_HO', 'Finite Element for HO   ''P1'' ''P1bubble'' ''P1bubblecondensed'' ''P1xP2'' ''P2xP1'' ''P2''');
			fielddisplay(obj,'fe_FS', 'Finite Element for FS   ''P1P1'' (debugging only) ''P1P1GLS'' ''MINIcondensed'' ''MINI'' ''TaylorHood'' ''XTaylorHood''');
			fielddisplay(obj,'vertex_equation','flow equation for each vertex');
			fielddisplay(obj,'element_equation','flow equation for each element');
			fielddisplay(obj,'borderSSA','vertices on SSA''s border (for tiling)');
			fielddisplay(obj,'borderHO','vertices on HO''s border (for tiling)');
			fielddisplay(obj,'borderFS','vertices on FS'' border (for tiling)');

		end % }}}
		function marshall(obj,md,fid) % {{{
			WriteData(fid,'object',obj,'fieldname','isSIA','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','isSSA','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','isL1L2','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','isHO','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','isFS','format','Boolean');
			WriteData(fid,'enum',FlowequationFeSSAEnum(),'data',StringToEnum(obj.fe_SSA),'format','Integer');
			WriteData(fid,'enum',FlowequationFeHOEnum() ,'data',StringToEnum(obj.fe_HO) ,'format','Integer');
			WriteData(fid,'enum',FlowequationFeFSEnum() ,'data',StringToEnum(obj.fe_FS) ,'format','Integer');
			WriteData(fid,'object',obj,'fieldname','borderSSA','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','borderHO','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','borderFS','format','DoubleMat','mattype',1);
			%convert approximations to enums
			data=obj.vertex_equation;
			pos=find(data==0); data(pos,end)=NoneApproximationEnum();
			pos=find(data==1); data(pos,end)=SIAApproximationEnum();
			pos=find(data==2); data(pos,end)=SSAApproximationEnum();
			pos=find(data==3); data(pos,end)=HOApproximationEnum();
			pos=find(data==4); data(pos,end)=FSApproximationEnum();
			pos=find(data==5); data(pos,end)=SSAHOApproximationEnum();
			pos=find(data==6); data(pos,end)=SSAFSApproximationEnum();
			pos=find(data==7); data(pos,end)=HOFSApproximationEnum();
			pos=find(data==8); data(pos,end)=L1L2ApproximationEnum();
			WriteData(fid,'data',data,'enum',FlowequationVertexEquationEnum(),'format','DoubleMat','mattype',1);
			data=obj.element_equation;
			pos=find(data==0); data(pos,end)=NoneApproximationEnum();
			pos=find(data==1); data(pos,end)=SIAApproximationEnum();
			pos=find(data==2); data(pos,end)=SSAApproximationEnum();
			pos=find(data==3); data(pos,end)=HOApproximationEnum();
			pos=find(data==4); data(pos,end)=FSApproximationEnum();
			pos=find(data==5); data(pos,end)=SSAHOApproximationEnum();
			pos=find(data==6); data(pos,end)=SSAFSApproximationEnum();
			pos=find(data==7); data(pos,end)=HOFSApproximationEnum();
			pos=find(data==8); data(pos,end)=L1L2ApproximationEnum();
			WriteData(fid,'data',data,'enum',FlowequationElementEquationEnum(),'format','DoubleMat','mattype',2);
		end % }}}
	end
end
