%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      = {};
	end
	methods
        function createxml(obj,fid) % {{{
            fprintf(fid, '\n\n');
            fprintf(fid, '%s\n', '<!-- StressBalance solution parameters: -->');
            
            % Convergence criteria            
            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Convergence criteria">','<section name="stressbalance" />'); 
			 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="restol" type="',      class(obj.restol),'" default="',        convert2str(obj.restol),'">',	'     <section name="stressbalance" />','     <help> mechanical equilibrium residual convergence criterion </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="reltol" type="',          class(obj.reltol),'" default="',            convert2str(obj.reltol),'">',     '     <section name="stressbalance" />','     <help> velocity relative convergence criterion, NaN: not applied </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',	'<parameter key ="abstol" type="',        class(obj.abstol),'" default="',          convert2str(obj.abstol),'">',	'     <section name="stressbalance" />','     <help> velocity absolute convergence criterion, NaN: not applied </help>','</parameter>');
             
             %is newton drop-down (0,1,or 2)
             fprintf(fid,'%s\n%s\n%s\n%s\n',    '<parameter key ="isnewton" type="alternative" optional="false">',	'     <section name="stressbalance" />','     <help> 0: Picards fixed point, 1: Newtons method, 2: hybrid </help>');
             fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
             fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
             fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
    
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="maxiter" type="',  	class(obj.maxiter),'" default="',       convert2str(obj.maxiter),'">',     '     <section name="stressbalance" />','     <help> maximum number of nonlinear iterations </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="viscosity_overshoot" type="',         class(obj.viscosity_overshoot),'" default="',           convert2str(obj.viscosity_overshoot),'">',	'     <section name="stressbalance" />','     <help> over-shooting constant new=new+C*(new-old) </help>','</parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');
            
            % boundary conditions    
            fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Boundary conditions">','<section name="stressbalance" />'); 			 
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcvx" type="',     	class(obj.spcvx),'" default="',    	convert2str(obj.spcvx),'">',     '     <section name="stressbalance" />','     <help> x-axis velocity constraint (NaN means no constraint) [m/yr] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcvy" type="',     	class(obj.spcvy),'" default="',    	convert2str(obj.spcvy),'">',	'     <section name="stressbalance" />','     <help> y-axis velocity constraint (NaN means no constraint) [m/yr] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcvz" type="',    	class(obj.spcvz),'" default="',    	convert2str(obj.spcvz),'">',     '     <section name="stressbalance" />','     <help> z-axis velocity constraint (NaN means no constraint) [m/yr] </help>','</parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');
            
            % Rift options   
            fprintf(fid,'%s\n%s\n%s\n','<frame key="3" label="Rift options">','<section name="stressbalance" />'); 	
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="rift_penalty_threshold" type="',    	class(obj.rift_penalty_threshold),'" default="',       convert2str(obj.rift_penalty_threshold),'">',     '     <section name="stressbalance" />','     <help> threshold for instability of mechanical constraints </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    	'<parameter key ="rift_penalty_lock" type="',         class(obj.rift_penalty_lock),'" default="',           convert2str(obj.rift_penalty_lock),'">',	'     <section name="stressbalance" />','     <help> number of iterations before rift penalties are locked </help>','</parameter>');
            fprintf(fid,'%s\n%s\n','</frame>'); 
            
            %others
            fprintf(fid,'%s\n%s\n%s\n','<frame key="4" label="Other">','<section name="stressbalance" />'); 	
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="shelf_dampening" type="',  	class(obj.shelf_dampening),'" default="',       convert2str(obj.shelf_dampening),'">',     '     <section name="stressbalance" />','     <help> use dampening for floating ice ? Only for FS model </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',  	'<parameter key ="FSreconditioning" type="', 	class(obj.FSreconditioning),'" default="',           convert2str(obj.FSreconditioning),'">',	'     <section name="stressbalance" />','     <help> multiplier for incompressibility equation. Only for FS model </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="referential" type="',  	class(obj.referential),'" default="',       convert2str(obj.referential),'">',     '     <section name="stressbalance" />','     <help> local referential </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',  	'<parameter key ="loadingforce" type="', 	class(obj.loadingforce),'" default="',           convert2str(obj.loadingforce),'">',	'     <section name="stressbalance" />','     <help> loading force applied on each point [N/m^3] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="requested_outputs" type="',  	class(obj.requested_outputs),'" default="',       convert2str(obj.requested_outputs),'">',     '     <section name="stressbalance" />','     <help> additional outputs requested </help>','</parameter>');
             
            
            fprintf(fid,'%s\n%s\n','</frame>'); 
        end % }}}
		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;

			 %output default:
			 obj.requested_outputs={'default'};

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

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

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

			%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 strcmp(domaintype(md.mesh),'3D') & md.flowequation.isFS,
				pos=find(md.mask.groundedice_levelset>0. & md.mesh.vertexonbase);
				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,'fieldname','stressbalance.FSreconditioning','>',0);
			end
		end % }}}
		function list=defaultoutputs(self,md) % {{{

			if dimension(md.mesh)==3,
				list = {'Vx','Vy','Vz','Vel','Pressure'};
			elseif dimension(md.mesh)==2,
				list = {'Vx','Vy','Vel','Pressure'};
			else
				error('mesh type not supported yet');
			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,'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);

			%process requested outputs
			outputs = obj.requested_outputs;
			pos  = find(ismember(outputs,'default'));
			if ~isempty(pos),
				outputs(pos) = [];                         %remove 'default' from outputs
				outputs      = [outputs defaultoutputs(obj,md)]; %add defaults
			end
			WriteData(fid,'data',outputs,'enum',StressbalanceRequestedOutputsEnum,'format','StringArray');
		end % }}}
	end
end
