%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;
		FSreconditioning     = 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;
		loadingforce             = NaN;
		requested_outputs        = NaN;
	end
	methods
		function obj = diagnostic(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				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.FSreconditioning=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 disp(obj) % {{{

			disp(sprintf('   Stressbalance 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','0: Picard''s fixed point, 1: Newton''s method, 2: hybrid');
			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) [m/yr]');
			fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]');
			fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]');
			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 FS model');
			fielddisplay(obj,'FSreconditioning','multiplier for incompressibility equation. Only for FS model');
			fielddisplay(obj,'referential','local referential');
			fielddisplay(obj,'loadingforce','loading force applied on each point [N/m^3]');
			fielddisplay(obj,'requested_outputs','additional outputs requested');

		end % }}}

		function marshall(obj,md,fid) % {{{
			WriteData(fid,'object',obj,'fieldname','spcvx','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'fieldname','spcvy','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'fieldname','spcvz','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+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','Integer');
			WriteData(fid,'object',obj,'fieldname','FSreconditioning','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,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum);
			WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum);
			WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum);
			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',StressbalanceIcefrontEnum(),'format','DoubleMat','mattype',3);
		end % }}}
	end
end
