%DAMAGEICE class definition
%
%   Usage:
%      damage=damage();

classdef damage
	properties (SetAccess=public)  
		%damage 
		D                   = NaN;
		law                 = '';
		spcdamage           = NaN; 
		max_damage          = NaN;
	
		%numerical
		stabilization       = NaN;
		penalty_threshold   = NaN;
		maxiter             = NaN;
		penalty_lock        = NaN;
		penalty_factor      = NaN;
		
		%general parameters for evolution law: 
		stress_threshold    = NaN;
		c1                  = NaN;
		c2                  = NaN;
		c3                  = NaN;
		c4                  = NaN;
		healing             = NaN;
		equiv_stress		  = NaN;
		requested_outputs   = {};
	end
	methods
		function obj = damage(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					inputstruct=varargin{1};
					list1 = properties('damage');
					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) % {{{

			%damage parameters: 
			obj.D=0;
			obj.law='undamaged';
			
			obj.max_damage=1-1e-5; %if damage reaches 1, solve becomes singular, as viscosity becomes nil
		
			%Type of stabilization used
			obj.stabilization=2;
			
			%Maximum number of iterations
			obj.maxiter=100;

			%factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
			obj.penalty_factor=3;
			
			%stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)
			obj.penalty_lock=0;
			
			%threshold to declare convergence of thermal solution (default is 0)
			obj.penalty_threshold=0;
		
			%damage evolution parameters 
			obj.stress_threshold=0;
			obj.healing=0;
			obj.c1=0;
			obj.c2=0;
			obj.c3=0;
			obj.c4=0;
			obj.equiv_stress=0;

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

		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{
			
			md = checkfield(md,'fieldname','damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]);
			md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0);
			md = checkfield(md,'fieldname','damage.law','values',{'undamaged','pralong'});
			md = checkfield(md,'fieldname','damage.spcdamage','forcing',1);
			
			md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]);
			md = checkfield(md,'fieldname','damage.maxiter','>=0',0);
			md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0);
			md = checkfield(md,'fieldname','damage.penalty_lock','>=0',0);
			md = checkfield(md,'fieldname','damage.penalty_threshold','>=0',0);

			if strcmpi(obj.law,'pralong'),
				md = checkfield(md,'fieldname','damage.healing','>=',0);
				md = checkfield(md,'fieldname','damage.c1','>=',0);
				md = checkfield(md,'fieldname','damage.c2','>=',0);
				md = checkfield(md,'fieldname','damage.c3','>=',0);
				md = checkfield(md,'fieldname','damage.c4','>=',0);
				md = checkfield(md,'fieldname','damage.stress_threshold','>=',0);
				md = checkfield(md,'fieldname','damage.equiv_stress','numel',[1],'values',[0 1]);
				md = checkfield(md,'fieldname','damage.requested_outputs','stringrow',1);
			elseif strcmpi(obj.law,'undamaged'),
				if (solution==DamageEvolutionSolutionEnum),
					error('Invalid evolution law (md.damage.law) for a damage solution');
				end
			else 
				error('invalid damage evolution law');
			end

		end % }}}
		function list=defaultoutputs(self,md) % {{{

			if strcmp(meshtype(md.mesh),'3D'),
				list = {'DamageD',};
			elseif strcmp(meshtype(md.mesh),'2Dhorizontal'),
				list = {'DamageD'};
			elseif strcmp(meshtype(md.mesh),'2Dvertical'),
				list = {'DamageD'};
			else
				error('mesh type not supported yet');
			end

		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   Damage:\n'));

			fielddisplay(obj,'D','damage tensor (scalar)');
			fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}');
			fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint)');
			fielddisplay(obj,'max_damage','maximum possible damage (0<=max_damage<1)');
			
			fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
			fielddisplay(obj,'maxiter','maximum number of non linear iterations');
			fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
			fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)');
			fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');

			if strcmpi(obj.law,'pralong'),
				fielddisplay(obj,'c1','damage parameter 1');
				fielddisplay(obj,'c2','damage parameter 2');
				fielddisplay(obj,'c3','damage parameter 3');
				fielddisplay(obj,'c4','damage parameter 4');
				fielddisplay(obj,'healing','damage healing parameter');
				fielddisplay(obj,'stress_threshold','damage stress threshold [Pa]');
				fielddisplay(obj,'equiv_stress','0: von Mises, 1: max principal');
				fielddisplay(obj,'requested_outputs','additional outputs requested');
			end

		end % }}}
		function marshall(obj,md,fid) % {{{
		
			WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','law','format','String');
			WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'fieldname','max_damage','format','Double');

			WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
			WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
			WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
			WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
			WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
	
			if strcmpi(obj.law,'pralong'),
				WriteData(fid,'object',obj,'fieldname','c1','format','Double');
				WriteData(fid,'object',obj,'fieldname','c2','format','Double');
				WriteData(fid,'object',obj,'fieldname','c3','format','Double');
				WriteData(fid,'object',obj,'fieldname','c4','format','Double');
				WriteData(fid,'object',obj,'fieldname','stress_threshold','format','Double');
				WriteData(fid,'object',obj,'fieldname','healing','format','Double');
				WriteData(fid,'object',obj,'fieldname','equiv_stress','format','Integer');
					
				%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',DamageEvolutionRequestedOutputsEnum,'format','StringArray');
			end

		end % }}}
	end
end
