%STRESSBALANC Eclass definition
%
%   Usage:
%      stressbalance=stressbalance();

classdef stressbalance
	properties (SetAccess=public) 
		spcvx                  = NaN;
		spcvy                  = NaN;
		spcvz                  = NaN;
		restol                 = 0;
		reltol                 = 0;
		abstol                 = 0;
		isnewton               = 0;
		FSreconditioning       = 0;
		viscosity_overshoot    = 0;
		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 = stressbalance(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					inputstruct=varargin{1};
					list1 = properties('stressbalance');
					list2 = fieldnames(inputstruct);
					for i=1:length(list1)
						fieldname = list1{i};
						if ismember(fieldname,list2),
							obj.(fieldname) = inputstruct.(fieldname);
						end
					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.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 stressbalance 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 md = checkconsistency(obj,md,solution,analyses) % {{{

			%Early return
			if ~ismember(StressbalanceAnalysisEnum(),analyses), return; end

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

			%singular solution
			if ~(any(~isnan(md.stressbalance.spcvx)) & any(~isnan(md.stressbalance.spcvy))),
				md = checkmessage(md,['model 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.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6),
				md = checkmessage(md,['Each line of stressbalance.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.stressbalance.referential),2)==0),
				pos=find(sum(isnan(md.stressbalance.referential),2)==0);
				if any(abs(dot(md.stressbalance.referential(pos,1:3),md.stressbalance.referential(pos,4:6),2))>eps),
					md = checkmessage(md,['Vectors in stressbalance.referential (columns 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.isFS,
				pos=find(md.mask.groundedice_levelset>0. & md.mesh.vertexonbed);
				if any(~isnan(md.stressbalance.referential(pos,:))),
					md = checkmessage(md,['no referential should be specified for basal vertices of grounded ice']);
				end
				md = checkfield(md,'stressbalance.FSreconditioning','>',0);
			end
		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]');

			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) % {{{

			yts=365.0*24.0*3600.0;

			WriteData(fid,'object',obj,'class','stressbalance','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','restol','format','Double');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','reltol','format','Double');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','abstol','format','Double');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','isnewton','format','Integer');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','FSreconditioning','format','Double');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','viscosity_overshoot','format','Double');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','maxiter','format','Integer');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','shelf_dampening','format','Integer');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','penalty_factor','format','Double');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer');
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'class','stressbalance','fieldname','requested_outputs','format','DoubleMat','mattype',3);
			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);
		end % }}}
	end
end
