%Hydrologydc class definition
%
%   Usage:
%      hydrologydc=hydrologydc();

classdef hydrologydc
	properties (SetAccess=public) 
		water_compressibility    = 0;
		isefficientlayer         = 0;
		penalty_factor           = 0;
		penalty_lock             = 0;
		rel_tol                  = 0;
		max_iter                 = 0;
		sedimentlimit_flag       = 0;
		sedimentlimit            = 0;
		transfer_flag            = 0;
		leakage_factor           = 0;
		basal_moulin_input       = NaN;

		spcsediment_head         = NaN;
		sediment_transmitivity   = NaN;
		sediment_compressibility = 0;
		sediment_porosity        = 0;
		sediment_thickness       = 0;


		spcepl_head              = NaN;
		mask_eplactive_node      = NaN;
		epl_compressibility      = 0;
		epl_porosity             = 0;
		epl_initial_thickness    = 0;
		epl_max_thickness    = 0;
		epl_conductivity         = 0;
  end
	methods
		% {{{ function obj = hydrologydc(varargin) 
		function obj = hydrologydc(varargin) 
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				otherwise
					error('constructor not supported');
	    end
		end 
		% }}}
		function self = initialize(self,md) % {{{
			if isnan(self.basal_moulin_input),
				self.basal_moulin_input=zeros(md.mesh.numberofvertices,1);
				disp('      no hydrology.basal_moulin_input specified: values set as zero');
			end

		end % }}}
		% {{{ function obj = setdefaultparameters(obj) 
		function obj = setdefaultparameters(obj) 

		%Parameters from de Fleurian 2014
			obj.water_compressibility    = 5.04e-10;
			obj.isefficientlayer         = 1;
			obj.penalty_factor           = 3;
			obj.rel_tol                  = 1.0e-06;
			obj.max_iter                 = 100;
			obj.sedimentlimit_flag       = 0;
			obj.sedimentlimit            = 0;
			obj.transfer_flag            = 0;
			obj.leakage_factor           = 10.0;

			obj.sediment_compressibility = 1.0e-08;
			obj.sediment_porosity        = 0.4;
			obj.sediment_thickness       = 20.0;
			obj.sediment_transmitivity   = 8.0e-04;

			obj.epl_compressibility      = 1.0e-08;
			obj.epl_porosity             = 0.4;
			obj.epl_initial_thickness    = 1.0;
			obj.epl_max_thickness    = 5.0;
			obj.epl_conductivity         = 8.0e-02;

		end 
		% }}}
		% {{{ function md = checkconsistency(obj,md,solution,analyses) 
		function md = checkconsistency(obj,md,solution,analyses) 
		%Early return
			if ~ismember(HydrologyDCInefficientAnalysisEnum(),analyses) & ~ismember(HydrologyDCEfficientAnalysisEnum(),analyses),
				return;
			end

			md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0 1]);
			md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.rel_tol','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]);
			md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0 1]);
			if obj.sedimentlimit_flag==1,
				md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1);
	    end
			if obj.transfer_flag==1,
				md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1);
	    end
			md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1);

			md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1);
			md = checkfield(md,'fieldname','hydrology.sediment_compressibility','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0,'numel',1);
			md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices 1]);

			if obj.isefficientlayer==1,
				md = checkfield(md,'fieldname','hydrology.spcepl_head','forcing',1);
				md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices 1],'values',[0 1]);
				md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0,'numel',1);
				md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0,'numel',1);
				md = checkfield(md,'fieldname','hydrology.epl_initial_thickness','>',0,'numel',1);
				md = checkfield(md,'fieldname','hydrology.epl_max_thickness','>',0,'numel',1);
				md = checkfield(md,'fieldname','hydrology.epl_conductivity','>',0,'numel',1);
	    end
		end 
		% }}}
		% {{{ function disp(obj)
		function disp(obj) 
			disp(sprintf('   hydrology Dual Porous Continuum Equivalent parameters:'));
			disp(sprintf('   - general parameters'));
			fielddisplay(obj,'water_compressibility','compressibility of water [Pa^-1]');
			fielddisplay(obj,'isefficientlayer','do we use an efficient drainage system [1: true; 0: false]');
			fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]');
			fielddisplay(obj,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
			fielddisplay(obj,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]');
			fielddisplay(obj,'max_iter','maximum number of nonlinear iteration');
			fielddisplay(obj,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer');
			disp(sprintf('%55s  0: no limit',' '));
			disp(sprintf('%55s  1: user defined: %s',' ','sedimentlimit'));
			disp(sprintf('%55s  2: hydrostatic pressure',' '));
			disp(sprintf('%55s  3: normal stress',' '));
			if obj.sedimentlimit_flag==1,
				fielddisplay(obj,'sedimentlimit','user defined upper limit for the inefficient layer [m]');
	    end
			fielddisplay(obj,'transfer_flag',['what kind of transfer method is applied between the layers']);
			disp(sprintf('%55s  0: no transfer',' '));
			disp(sprintf('%55s  1: constant leakage factor: %s',' ','leakage_factor'));
			if obj.transfer_flag==1,
				fielddisplay(obj,'leakage_factor','user defined leakage factor [m]');
	    end
			fielddisplay(obj,'basal_moulin_input','Figure out what it is');
			disp(sprintf('   - for the sediment layer'));
			fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]');
			fielddisplay(obj,'sediment_compressibility','sediment compressibility [Pa^-1]');
			fielddisplay(obj,'sediment_porosity','sediment [dimensionless]');
			fielddisplay(obj,'sediment_thickness','sediment thickness [m]');
			fielddisplay(obj,'sediment_transmitivity','sediment transmitivity [m^2/s]');

			if obj.isefficientlayer==1,
				disp(sprintf('   - for the epl layer'));
				fielddisplay(obj,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]');
				fielddisplay(obj,'mask_eplactive_node','active (1) or not (0) EPL');
				fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]');
				fielddisplay(obj,'epl_porosity','epl [dimensionless]');
				fielddisplay(obj,'epl_initial_thickness','epl initial thickness [m]');
				fielddisplay(obj,'epl_max_thickness','epl maximal thickness [m]');
				fielddisplay(obj,'epl_conductivity','epl conductivity [m^2/s]');
	    end

		end 
		% }}}
		% {{{ function marshall(obj,md,fid) 
		function marshall(obj,md,fid) 
			WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer');
			WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double');
			WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
			WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
			WriteData(fid,'object',obj,'fieldname','rel_tol','format','Double');
			WriteData(fid,'object',obj,'fieldname','max_iter','format','Integer');
			WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer');
			WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer');
			if obj.sedimentlimit_flag==1,
				WriteData(fid,'object',obj,'fieldname','sedimentlimit','format','Double');
	    end
			if obj.transfer_flag==1,
				WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double');
	    end
			WriteData(fid,'object',obj,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)

			WriteData(fid,'object',obj,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'fieldname','sediment_compressibility','format','Double');			
			WriteData(fid,'object',obj,'fieldname','sediment_porosity','format','Double');			
			WriteData(fid,'object',obj,'fieldname','sediment_thickness','format','Double');
			WriteData(fid,'object',obj,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1');		

			if obj.isefficientlayer==1,	
				WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);	
				WriteData(fid,'object',obj,'fieldname','mask_eplactive_node','format','DoubleMat','mattype',1);
				WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double');			
				WriteData(fid,'object',obj,'fieldname','epl_porosity','format','Double');			
				WriteData(fid,'object',obj,'fieldname','epl_initial_thickness','format','Double');
				WriteData(fid,'object',obj,'fieldname','epl_max_thickness','format','Double');
				WriteData(fid,'object',obj,'fieldname','epl_conductivity','format','Double');
			end
		end 
% }}}
  end
end
