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

classdef damage
	properties (SetAccess=public)  
		%damage 
		isdamage            = NaN;
		D                   = NaN;
		law                 = '';
		spcdamage           = NaN; 
		max_damage          = NaN;
	
		%numerical
		stabilization       = NaN;
		maxiter             = NaN;
		elementinterp       = '';
		penalty_threshold   = 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 createxml(obj,fid) % {{{
            fprintf(fid, '\n\n');
            fprintf(fid, '%s\n', '<!-- damage -->');
            fprintf(fid, '%s\n', '<!-- Note: this class depends on different input of law -->');
            
            %fprintf(fid,'%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="law" type="logical"',           '" default="',                num2str(obj.law),'">',   '     <section name="damage" />','     <help> damage law (string) from {"undamaged","pralong"} </help>','</parameter>');
            % drop-down 
            fprintf(fid,'%s%s%s%s%s\n\t%s\n','<parameter key ="law" type="','alternative','" optional="','false','">','<section name="damage" />');
            
            % law = 'undamage'
            fprintf(fid,'\t%s%s%s%s%s\n\t\t%s\n','<option value="undamage" type="','string','" default="','true','">','<help> law = undamage </help>');
            % footer for option
            fprintf(fid,'\t%s\n%s\n','</option>');
            
            % law = 'pralong'
            fprintf(fid,'\t%s%s%s%s%s\n\t\t%s\n','<option value="pralong" type="','string','" default="','false','">','<help> law = pralong </help>');
                
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',   '<parameter key ="stress_threshold" type="',class(obj.stress_threshold),'" default="',num2str(obj.stress_threshold),'">','<help> damage stress threshold [Pa] </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',	'<parameter key ="c1" type="', class(obj.c1),'" default="',   num2str(obj.c1),'">',   '<help> damage parameter 1 </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',  	'<parameter key ="c2" type="',          class(obj.c2),'" default="',            num2str(obj.c2),'">','<help> damage parameter 2 </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',  	'<parameter key ="c3" type="',          class(obj.c3),'" default="',            num2str(obj.c3),'">','<help> damage parameter 3 [W/m^2] </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',  	'<parameter key ="c4" type="',            class(obj.c4),'" default="',              num2str(obj.c4),'">','<help> damage parameter 4 </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',  	'<parameter key ="healing" type="', class(obj.healing),'" default="',   num2str(obj.healing),'">','<help> damage healing parameter 1 </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',  	'<parameter key ="equiv_stress" type="',          class(obj.equiv_stress),'" default="',convert2str(obj.equiv_stress),'">','<help> 0: von Mises </help>','</parameter>');
            fprintf(fid,'\t\t%s%s%s%s%s\n\t\t\t%s\n\t\t%s\n',  	'<parameter key ="requested_outputs" type="',          class(obj.requested_outputs),'" default="',convert2str(obj.requested_outputs),'">','<help> additional outputs requested </help>','</parameter>');
             
                
            % footer for option
            fprintf(fid,'\t%s\n%s\n','</option>');
                 
            
            % footer for drop-down
            fprintf(fid,'\t%s\n%s\n%s','<help> damage law (string) from {"undamaged","pralong"} </help>','</parameter>');
            
            
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="D" type="',              class(obj.D),'" default="',                  num2str(obj.D),'">',              '     <section name="damage" />','     <help> damage tensor (scalar) </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="law" type="',            class(obj.law),'" default="',                num2str(obj.law),'">',   '     <section name="damage" />','     <help> damage law (string) from {"undamaged","pralong"} </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="spcdamage" type="',      class(obj.spcdamage),'" default="',          num2str(obj.spcdamage),'">',            '     <section name="damage" />','     <help> damage constraints (NaN means no constraint) </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="max_damage" type="',     class(obj.max_damage),'" default="',         num2str(obj.max_damage),'">',            '     <section name="damage" />','     <help> maximum possible damage (0&amp;lt;=max_damage&amp;lt;1) </help>','</parameter>');
             
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="stabilization" type="',       class(obj.stabilization),'" default="',              num2str(obj.stabilization),'">',              '     <section name="damage" />','     <help> 0: no, 1: artificial_diffusivity, 2: SUPG </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="maxiter" type="',             class(obj.maxiter),'" default="',   num2str(obj.maxiter),'">',   '     <section name="damage" />','     <help> maximum number of non linear iterations </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="penalty_lock" type="',        class(obj.penalty_lock),'" default="',            num2str(obj.penalty_lock),'">',            '     <section name="damage" />','     <help> stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="penalty_threshold" type="',   class(obj.penalty_threshold),'" default="',            num2str(obj.penalty_threshold),'">',            '     <section name="damage" />','     <help> threshold to declare convergence of damage evolution solution (default is 0) </help>','</parameter>');
            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="penalty_factor" type="',      class(obj.penalty_factor),'" default="',            num2str(obj.penalty_factor),'">',            '     <section name="damage" />','     <help> scaling exponent (default is 3) </help>','</parameter>');
 
        end % }}}
        
		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.isdamage=0;
			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;

			%finite element interpolation
			obj.elementinterp='P1';

			%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.isdamage','values',[1,0]);
			if obj.isdamage,
				md = checkfield(md,'fieldname','damage.law','values',{'undamaged','pralong'});
				md = checkfield(md,'fieldname','damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','damage.spcdamage','forcing',1);
				md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0);
				
				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.elementinterp','values',{'P1','P2'});
				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);
			end

			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(domaintype(md.mesh),'2Dhorizontal'),
            list = {'DamageDbar'};
         else
            list = {'DamageD'};
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   Damage:\n'));

			fielddisplay(obj,'isdamage','is damage mechanics being used? {true,false}');
			if obj.isdamage,
				fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}');
				fielddisplay(obj,'D','damage tensor (scalar)');
				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,'elementinterp','interpolation scheme for finite elements {''P1'',''P2''}');
				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)');
			end

			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','isdamage','format','Boolean');
			if obj.isdamage,
				WriteData(fid,'object',obj,'fieldname','law','format','String');
				WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1);
				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','maxiter','format','Integer');
				WriteData(fid,'object',obj,'fieldname','elementinterp','format','String');
				WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
				WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
				WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
			end
	
			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');
			end

			%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
			if obj.isdamage,
				WriteData(fid,'data',outputs,'enum',DamageEvolutionRequestedOutputsEnum,'format','StringArray');
			end

		end % }}}
	end
end