%DIAGNOSTIC class definition
%
%   Usage:
%      diagnostic=diagnostic();

classdef diagnostic
	properties (SetAccess=public) 
		spcvx                    = NaN;
		spcvy                    = NaN;
		spcvz                    = NaN;
		restol                   = 0;
		reltol                   = 0;
		abstol                   = 0;
		isnewton                 = 0;
		stokesreconditioning     = 0;
		viscosity_overshoot      = 0;
		icefront                 = NaN;
		maxiter                  = 0;
		shelf_dampening          = 0;
		vertex_pairing           = NaN;
		penalty_factor           = NaN;
		rift_penalty_lock        = NaN;
		rift_penalty_threshold   = 0;
		referential              = NaN;
		requested_outputs        = NaN;
	end
	methods
		function obj = diagnostic(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					in=varargin{1};
					if (isa(in,'numeric') & in==0),
						% requesting templates do nothing
					else
						error('constructor not supported');
					end 
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			 %maximum of non-linear iterations.
			 obj.maxiter=100;
			 
			 %Convergence criterion: absolute, relative and residual
			 obj.restol=10^-4; 
			 obj.reltol=0.01;
			 obj.abstol=10;
			 
			 obj.stokesreconditioning=10^13;
			 obj.shelf_dampening=0;

			 %Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
			 obj.penalty_factor=3;

			 %coefficient to update the viscosity between each iteration of
			 %a diagnostic according to the following formula
			 %viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
			 obj.viscosity_overshoot=0;

			 %Stop the iterations of rift if below a threshold
			 obj.rift_penalty_threshold=0;

			 %in some solutions, it might be needed to stop a run when only
			 %a few constraints remain unstable. For thermal computation, this
			 %parameter is often used.
			 obj.rift_penalty_lock=10;

		end % }}}
		function checkconsistency(obj,md,solution,analyses) % {{{

			%Early return
			if ~ismember(DiagnosticHorizAnalysisEnum,analyses), return; end
			%if ~ismember(DiagnosticHorizAnalysisEnum,analyses) |  (solution==TransientSolutionEnum & md.transient.isdiagnostic==0), return; end

			checkfield(md,'diagnostic.spcvx','forcing',1);
			checkfield(md,'diagnostic.spcvy','forcing',1);
			if md.mesh.dimension==3, checkfield(md,'diagnostic.spcvz','forcing',1); end
			checkfield(md,'diagnostic.restol','size',[1 1],'>',0);
			checkfield(md,'diagnostic.reltol','size',[1 1]);
			checkfield(md,'diagnostic.abstol','size',[1 1]);
			checkfield(md,'diagnostic.isnewton','numel',1,'values',[0 1]);
			checkfield(md,'diagnostic.stokesreconditioning','size',[1 1],'NaN',1);
			checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1);
			if md.mesh.dimension==2,
				checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1);
			else
				checkfield(md,'diagnostic.icefront','size',[NaN 6],'NaN',1);
			end
			checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]);
			checkfield(md,'diagnostic.maxiter','size',[1 1],'>=',1);
			checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices 6]);
			if ~isempty(md.diagnostic.requested_outputs),
				checkfield(md,'diagnostic.requested_outputs','size',[NaN 1]);
			end

			%singular solution
			if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2),
				checkmessage(['model ' md.miscellaneous.name ' is not well posed (singular). You need at least one node with fixed velocity!'])
			end
			%CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
			if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6),
				checkmessage(['model ' md.miscellaneous.name ' has problem with rotated spc. Each line of diagnostic.referential should contain either only NaN values or no NaN values']);
			end
			%CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
			if any(sum(isnan(md.diagnostic.referential),2)==0),
				pos=find(sum(isnan(md.diagnostic.referential),2)==0);
				if any(abs(dot(md.diagnostic.referential(pos,1:3)',md.diagnostic.referential(pos,4:6)'))>eps),
					checkmessage(['model ' md.miscellaneous.name ' has problem with referential. Vectors in diagnostic.referential (colums 1 to 3 and 4 to 6) must be orthogonal']);
				end
			end
			%CHECK THAT NO rotation specified for FS Grounded ice at base
			if md.mesh.dimension==3 & md.flowequation.isstokes,
				pos=find(md.mask.vertexongroundedice & md.mesh.vertexonbed);
				if any(~isnan(md.diagnostic.referential(pos,:))),
					checkmessage(['model ' md.miscellaneous.name ' not consistent. No referential should be specified for basal vertices of grounded ice']);
				end
			end
		end % }}}
		function disp(obj) % {{{

			disp(sprintf('   Diagnostic solution parameters:'));

			disp(sprintf('\n      %s','Convergence criteria:'));
			fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion');
			fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN -> not applied');
			fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN -> not applied');
			fielddisplay(obj,'isnewton','Apply Newton''s method instead of a Picard fixed point method');
			fielddisplay(obj,'maxiter','maximum number of nonlinear iterations');
			fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');

			disp(sprintf('\n      %s','boundary conditions:'));
			fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint)');
			fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint)');
			fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint)');
			fielddisplay(obj,'icefront','segments on ice front list (last column 0-> Air, 1-> Water, 2->Ice');

			disp(sprintf('\n      %s','Rift options:'));
			fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints');
			fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked');

			disp(sprintf('\n      %s','Penalty options:'));
			fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
			fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized');

			disp(sprintf('\n      %s','Other:'));
			fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for Stokes model');
			fielddisplay(obj,'stokesreconditioning','multiplier for incompressibility equation. Only for Stokes model');
			fielddisplay(obj,'referential','local referential');
			fielddisplay(obj,'requested_outputs','additional outputs requested');

		end % }}}
		function marshall(obj,fid) % {{{
			WriteData(fid,'object',obj,'fieldname','spcvx','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','spcvy','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','spcvz','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','restol','format','Double');
			WriteData(fid,'object',obj,'fieldname','reltol','format','Double');
			WriteData(fid,'object',obj,'fieldname','abstol','format','Double');
			WriteData(fid,'object',obj,'fieldname','isnewton','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','stokesreconditioning','format','Double');
			WriteData(fid,'object',obj,'fieldname','viscosity_overshoot','format','Double');
			WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
			WriteData(fid,'object',obj,'fieldname','shelf_dampening','format','Integer');
			WriteData(fid,'object',obj,'fieldname','vertex_pairing','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
			WriteData(fid,'object',obj,'fieldname','rift_penalty_lock','format','Integer');
			WriteData(fid,'object',obj,'fieldname','rift_penalty_threshold','format','Integer');
			WriteData(fid,'object',obj,'fieldname','referential','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','requested_outputs','format','DoubleMat','mattype',3);

			%marshall ice front
			data=obj.icefront;
			pos=find(data(:,end)==0); data(pos,end)=AirEnum;
			pos=find(data(:,end)==1); data(pos,end)=WaterEnum;
			pos=find(data(:,end)==2); data(pos,end)=IceEnum;
			WriteData(fid,'data',data,'enum',DiagnosticIcefrontEnum,'format','DoubleMat','mattype',3);
		end % }}}
	end
end
