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

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

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

		spcepl_head              = NaN;
		mask_eplactive           = NaN;
		epl_compressibility      = 0;
		epl_porosity             = 0;
		epl_thickness            = 0;
		epl_transmitivity        = 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 obj = setdefaultparameters(obj) 
		function obj = setdefaultparameters(obj) 

		%Parameters from de Fleurian 2013
			obj.water_compressibility    = 5.04e-10;
			obj.isefficientlayer         = 1;
			obj.penalty_factor           = 3;
			obj.rel_tol                  = 1.0e-06;
			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_thickness            = 1.0;
			obj.epl_transmitivity        = 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,'hydrology.water_compressibility','>',0,'numel',1);
			md = checkfield(md,'hydrology.isefficientlayer','numel',[1],'values',[0 1]);
			md = checkfield(md,'hydrology.penalty_factor','>',0,'numel',1);
			md = checkfield(md,'hydrology.rel_tol','>',0,'numel',1);
			md = checkfield(md,'hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]);
			md = checkfield(md,'hydrology.transfer_flag','numel',[1],'values',[0 1]);
			if obj.sedimentlimit_flag==1,
				md = checkfield(md,'hydrology.sedimentlimit','>',0,'numel',1);
	    end
			if obj.transfer_flag==1,
				md = checkfield(md,'hydrology.leakage_factor','>',0,'numel',1);
	    end

			md = checkfield(md,'hydrology.spcsediment_head','forcing',1);
			md = checkfield(md,'hydrology.sediment_compressibility','>',0,'numel',1);
			md = checkfield(md,'hydrology.sediment_porosity','>',0,'numel',1);
			md = checkfield(md,'hydrology.sediment_thickness','>',0,'numel',1);
			md = checkfield(md,'hydrology.sediment_transmitivity','>',0,'numel',1);

			if obj.isefficientlayer==1,
				md = checkfield(md,'hydrology.spcepl_head','forcing',1);
				md = checkfield(md,'hydrology.mask_eplactive','size',[md.mesh.numberofvertices 1],'values',[0 1]);
				md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1);
				md = checkfield(md,'hydrology.epl_porosity','>',0,'numel',1);
				md = checkfield(md,'hydrology.epl_thickness','>',0,'numel',1);
				md = checkfield(md,'hydrology.epl_transmitivity','>',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,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]');
			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
			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','active (1) or not (0) EPL');
				fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]');
				fielddisplay(obj,'epl_porosity','epl [dimensionless]');
				fielddisplay(obj,'epl_thickness','epl thickness [m]');
				fielddisplay(obj,'epl_transmitivity','epl transmitivity [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','rel_tol','format','Double');
			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','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','Double');		

			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','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_thickness','format','Double');
				WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double');
			end
		end 
% }}}
  end
end
