%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                          = '';
		augmented_lagrangian_r         = 1.;
		augmented_lagrangian_rhop      = 1.;
		augmented_lagrangian_rlambda   = 1.;
		augmented_lagrangian_rholambda = 1.;
		XTH_theta                      = 0.;
		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 createxml(obj,fid) % {{{
			fprintf(fid, '\n\n');
			fprintf(fid, '%s\n', '<!-- flowequation -->');
			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Flow equation parameters">','<section name="flowequation" />');                    

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="isSIA" type="',              class(obj.isSIA),'" default="',                  convert2str(obj.isSIA),'">',              '     <section name="flowequation" />','     <help> is the Shallow Ice Approximation (SIA) used ? </help>','</parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="isSSA" type="',            class(obj.isSSA),'" default="',                convert2str(obj.isSSA),'">',   '     <section name="flowequation" />','     <help> is the Shelfy-Stream Approximation (SSA) used ? </help>','</parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="isL1L2" type="',      class(obj.isL1L2),'" default="',          convert2str(obj.isL1L2),'">',            '     <section name="flowequation" />','     <help> is the L1L2 approximation used ? </help>','</parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="isHO" type="',     class(obj.isHO),'" default="',         convert2str(obj.isHO),'">',            '     <section name="flowequation" />','     <help> is the Higher-Order (HO) approximation used ? </help>','</parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="isFS" type="',       class(obj.isFS),'" default="',              convert2str(obj.isFS),'">',              '     <section name="flowequation" />','     <help> are the Full-FS (FS) equations used ? </help>','</parameter>');

			% fe_SSA drop-down (P1, P1bubble, P1bubblecondensed, P2)
			fprintf(fid,'%s\n%s\n%s\n%s\n',    '<parameter key ="fe_SSA" type="alternative" optional="false">','     <section name="flowequation" />','     <help> Finite Element for SSA  "P1", "P1bubble" "P1bubblecondensed" "P2" </help>');
			fprintf(fid,'%s\n','       <option value="P1" type="string" default="true"> </option>');
			fprintf(fid,'%s\n','       <option value="P1bubble" type="string" default="false"> </option>');
			fprintf(fid,'%s\n','       <option value="P1bubblecondensed" type="string" default="false"> </option>');
			fprintf(fid,'%s\n%s\n','       <option value="P2" type="string" default="false"> </option>','</parameter>');

			%fe_HO drop-down (P1, P1bubble, P1bubblecondensed, P1xP2, P2xP1, P2)
			fprintf(fid,'%s\n%s\n%s\n%s\n',        '<parameter key ="fe_HO" type="alternative" optional="false">','     <section name="flowequation" />','     <help> Finite Element for HO   "P1" "P1bubble" "P1bubblecondensed" "P1xP2" "P2xP1" "P2" </help>');
			fprintf(fid,'%s\n','       <option value="P1" type="string" default="true"> </option>');
			fprintf(fid,'%s\n','       <option value="P1bubble" type="string" default="false"> </option>');
			fprintf(fid,'%s\n','       <option value="P1bubblecondensed" type="string" default="false"> </option>');
			fprintf(fid,'%s\n','       <option value="P1xP2" type="string" default="false"> </option>');
			fprintf(fid,'%s\n','       <option value="P2xP1" type="string" default="false"> </option>');
			fprintf(fid,'%s\n%s\n','       <option value="P2" type="string" default="false"> </option>','</parameter>');

			% fe_FS drop-down (P1P1, P1P1GLS, MINIcondensed, MINI, TaylowHood)
			fprintf(fid,'%s\n%s\n%s\n%s\n',        '<parameter key ="fe_FS" type="alternative" optional="false">','     <section name="flowequation" />','     <help> Finite Element for FS   "P1P1" (debugging only) "P1P1GLS" "MINIcondensed" "MINI" "TaylorHood" </help>');
			fprintf(fid,'%s\n','       <option value="P1P1" type="string" default="true"> </option>');
			fprintf(fid,'%s\n','       <option value="P1P1GLS" type="string" default="false"> </option>');
			fprintf(fid,'%s\n','       <option value="MINIcondensed" type="string" default="false"> </option>');
			fprintf(fid,'%s\n','       <option value="MINI" type="string" default="false"> </option>');
			fprintf(fid,'%s\n%s\n','       <option value="TaylorHood" type="string" default="false"> </option>','</parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="vertex_equation" type="',      class(obj.vertex_equation),'" default="',            convert2str(obj.vertex_equation),'">',            '     <section name="flowequation" />','     <help> flow equation for each vertex </help>','</parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="element_equation" type="',            class(obj.element_equation),'" default="',              convert2str(obj.element_equation),'">',              '     <section name="flowequation" />','     <help> flow equation for each element </help>','</parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="borderSSA" type="', class(obj.borderSSA),'" default="',   convert2str(obj.borderSSA),'">',   '     <section name="flowequation" />','     <help> vertices on SSAs border (for tiling) </help>','</parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="borderHO" type="',          class(obj.borderHO),'" default="',            convert2str(obj.borderHO),'">',            '     <section name="flowequation" />','     <help> vertices on HOs border (for tiling) </help>','</parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="borderFS" type="',          class(obj.borderFS),'" default="',            convert2str(obj.borderFS),'">',            '     <section name="flowequation" />','     <help> vertices on FS border (for tiling) </help>','</parameter>');

			fprintf(fid,'%s\n%s\n','</frame>');
		end % }}}
		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','P2bubble'});
				md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'});
				md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'});
				md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>',0.);
				md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rlambda','numel',[1],'>',0.);
				md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rhop','numel',[1],'>',0.);
				md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rholambda','numel',[1],'>',0.);
				md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5);
				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(domaintype(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(domaintype(md.mesh),'2Dvertical')
					md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[2,4,5]);
					md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[2,4,5]);
				elseif strcmp(domaintype(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,'enum',AugmentedLagrangianREnum(),'data',obj.augmented_lagrangian_r ,'format','Double');
			WriteData(fid,'enum',AugmentedLagrangianRhopEnum(),'data',obj.augmented_lagrangian_rhop ,'format','Double');
			WriteData(fid,'enum',AugmentedLagrangianRlambdaEnum(),'data',obj.augmented_lagrangian_rlambda ,'format','Double');
			WriteData(fid,'enum',AugmentedLagrangianRholambdaEnum(),'data',obj.augmented_lagrangian_rholambda ,'format','Double');
			WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',obj.XTH_theta ,'format','Double');
			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)=L1L2ApproximationEnum();
			pos=find(data==4); data(pos,end)=HOApproximationEnum();
			pos=find(data==5); data(pos,end)=FSApproximationEnum();

			pos=find(data==6); data(pos,end)=SSAHOApproximationEnum();
			pos=find(data==7); data(pos,end)=HOFSApproximationEnum();
			pos=find(data==8); data(pos,end)=SSAFSApproximationEnum();
			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)=L1L2ApproximationEnum();
			pos=find(data==4); data(pos,end)=HOApproximationEnum();
			pos=find(data==5); data(pos,end)=FSApproximationEnum();

			pos=find(data==6); data(pos,end)=SSAHOApproximationEnum();
			pos=find(data==7); data(pos,end)=SSAFSApproximationEnum();
			pos=find(data==8); data(pos,end)=HOFSApproximationEnum();
			WriteData(fid,'data',data,'enum',FlowequationElementEquationEnum(),'format','DoubleMat','mattype',2);
		end % }}}
	end
end
