classdef (Sealed) ThermalMaterialAssignment < pde.PDECompactDisplay
% ThermalMaterialAssignment - Specify material properties over domain or subdomain
%  For a thermal model, the toolbox solves equations of the following form:
%
%  (MassDensity*SpecificHeat)*dT/dt - div(ThermalConductivity*grad(T)) = HeatSource
%
%     This method creates an object representing the thermal properties:
%     MassDensity, SpecificHeat, and ThermalConductivity in a domain or
%     subdomain and appends the object to the MaterialProperties property.
%
%     Instances of this class can only be created by calling the
%     thermalProperties method of the ThermalModel class.
%
% See also pde.ThermalModel, pde.ThermalModel/thermalProperties

% Copyright 2016-2022 The MathWorks, Inc.

properties (SetAccess = private)   
    
    % RegionType - Type of geometric region the material is assigned to.
    %    The RegionType can be one of the following geometric entities: 'cell'
    %    or 'face'.
    RegionType;
end

properties
    % RegionID - ID of the geometric regions the coefficients are assigned to.
    %    For 2-D each RegionID satisfies 0 < RegionID(j) < NumFaces in the
    %    geometry. For 3-D each RegionID satisfies 0 < RegionID(j) < NumCells
    %    in the geometry.
    RegionID;

    %ThermalConductivity of the material
    %   A numerical scalar, vector or a matrix, or a function handle that
    %   compute thermal conductivity for a give spatial location at a
    %   given temperature.
    ThermalConductivity;

    %MassDensity of the material
    %   A numerical scalar, vector or a matrix, or a function handle that
    %   compute mass density for a give spatial location at a
    %   given temperature.
    MassDensity;

    %SpecificHeat of the material
    %   A numerical scalar, vector or a matrix, or a function handle that
    %   compute specific heat capacity for a give spatial location at a
    %   given temperature.
    SpecificHeat;
end

methods %ctor
    function self = ThermalMaterialAssignment(tmar,varargin)
        self.RecordOwner = tmar;
        parser = inputParser;
        parser.addParameter('face', []);
        parser.addParameter('cell', []);
        parser.addParameter('ThermalConductivity', []);
        parser.addParameter('MassDensity', []);
        parser.addParameter('SpecificHeat', []);
        parser.parse(varargin{:});

        numdims = 2;
        if ~isempty(parser.Results.face)
            self.RegionType = 'face';
            self.RegionID = parser.Results.face;
        elseif ~isempty(parser.Results.cell)
            self.RegionType = 'cell';
            self.RegionID = parser.Results.cell;
            numdims = 3;
        end
        systemsize = 1;

        self.ThermalConductivity = parser.Results.ThermalConductivity;
        if ~isempty(parser.Results.MassDensity)
            self.MassDensity = parser.Results.MassDensity;
        end

        if ~isempty(parser.Results.SpecificHeat)
            self.SpecificHeat = parser.Results.SpecificHeat;
        end

        self.checkAllMatrixMtlSizes(systemsize, numdims);
        self.checkFcnHdlArgCounts();
    end
end


methods % Setter methods

    function set.RegionType(self, rtype)
        self.RegionType = rtype;
    end

    function set.RegionID(self, rids)
        self.ValidateRegionID(rids);
        self.RegionID = rids;
    end

    function set.MassDensity(self, coef)
        self.CoefPrecheck(coef,'MassDensity');
        self.MassDensity = coef;
    end

    function set.SpecificHeat(self, coef)
        self.CoefPrecheck(coef,'SpecificHeat');
        self.SpecificHeat = coef;
    end

    function set.ThermalConductivity(self, coef)
        self.CoefPrecheck(coef,'ThermalConductivity');
        self.ThermalConductivity = coef;
    end


end


methods(Hidden=true, Access = {?pde.MaterialAssignmentRecords})
    function tf=sameCoefficients(self,other,varargin)
        if ~coefficientsMatch(self, other)
            tf = false;
            return
        end
        if isempty(varargin)
            tf = (isequal(self.MassDensity, other.MassDensity) && ...
                isequal(self.SpecificHeat, other.SpecificHeat) && ...
                isequal(self.ThermalConductivity, other.ThermalConductivity));
        else
            cname = varargin{1};
            switch cname
                case 'ThermalConductivity'
                    tf = isequal(self.ThermalConductivity, other.ThermalConductivity);
                case 'MassDensity'
                    tf = isequal(self.MassDensity, other.MassDensity);
                case 'SpecificHeat'
                    tf = isequal(self.SpecificHeat, other.SpecificHeat);
            end
        end
    end

    function tf=numericMaterialProperties(self,varargin)
        if isempty(varargin)

            tf = ~any(isa(self.ThermalConductivity, 'function_handle')  || ...
                isa(self.MassDensity, 'function_handle')  || ...
                isa(self.SpecificHeat, 'function_handle') );
        else
            cname = varargin{1};
            switch cname
                case 'ThermalConductivity'
                    tf = ~isa(self.ThermalConductivity, 'function_handle');
                case 'MassDensity'
                    tf = ~isa(self.MassDensity, 'function_handle');
                case 'SpecificHeat'
                    tf = ~isa(self.SpecificHeat, 'function_handle');
            end
        end
    end


    function tf = coefficientsMatch(self, other)
        tf = true;
        tf = tf & (self.ThermalConductivityDefined() == other.ThermalConductivityDefined());
        tf = tf & (self.MassDensityDefined() == other.MassDensityDefined());
        tf = tf & (self.SpecificHeatDefined() == other.SpecificHeatDefined());
    end

    function performSolverPrecheck(self, systemsize, numfaces, numcells)
        ndims = 2;
        if strcmp(self.RegionType, 'face')
            if any(self.RegionID > numfaces)
                error(message('pde:materialPropertySpecification:invalidFaceIndexPresolve'));
            end
        else
            if any(self.RegionID > numcells)
                error(message('pde:materialPropertySpecification:invalidCellIndexPresolve'));
            end
            ndims = 3;
        end
        checkAllMatrixMtlSizes(self, systemsize, ndims);
        checkFcnHdlArgCounts(self)
    end

end


methods(Hidden=true, Access = private)
    function checkAllMatrixMtlSizes(self, systemsize, ndims)
        if ~isempty(self.MassDensity)
            self.checkMassDensitySize(self.MassDensity, systemsize);
        end
        if ~isempty(self.SpecificHeat)
            self.checkSpecificHeatSize(self.SpecificHeat, systemsize);
        end
        %ThermalConductivity is a must define property for thermal model.
        self.checkThermalConductivitySize(self.ThermalConductivity, systemsize,ndims);
    end
    %

    function checkFcnHdlArgCounts(self)
        if self.ThermalConductivityDefined()
            self.checkCoefFcnHdlArgCounts(self.ThermalConductivity);
        end
        if self.MassDensityDefined()
            self.checkCoefFcnHdlArgCounts(self.MassDensity);
        end
        if self.SpecificHeatDefined()
            self.checkCoefFcnHdlArgCounts(self.SpecificHeat);
        end
    end
    
    function checkCoefFcnHdlArgCounts(self,coef)
        if ~isa(coef, 'function_handle')
            return
        end
        % check number of input and output
        numFcnInputs = nargin(coef);
        numFcnOutputs = abs(nargout(coef)); % abs for anonymous function handle, which returns -1.
        if (numFcnInputs ~= 2) || (numFcnOutputs ~= 1)
            error(message('pde:materialPropertySpecification:invalidFcnHandleArgs'));
        end
        % Try function with a sample input to check if it runs through and
        % returns a numeric output.
        [location, state] = self.RecordOwner.ParentModel.getSampleLocState(1);
        try
        res = coef(location, state);
        catch ME
            throwAsCaller(ME);
        end
        if ~isa(res,'double')
            error(message('pde:materialPropertySpecification:invalidFcnHandleOutputType'));
        end
    end

end % Private methods



methods(Hidden=true, Access = {?pde.MaterialAssignmentRecords})
    function tf = ThermalConductivityDefined(self)
        tf = self.coefDefined(self.ThermalConductivity);
    end

    function tf = MassDensityDefined(self)
        tf = self.coefDefined(self.MassDensity);
    end

    function tf = SpecificHeatDefined(self)
        tf = self.coefDefined(self.SpecificHeat);
    end


    function tf = hasComplexCoefficient(self, loc, state)
        tf = false(3,1);
        tf(1) = pde.ThermalMaterialAssignment.coefIsComplexNumericOrFcnHdl(self.ThermalConductivity, loc, state);
        tf(2) = pde.ThermalMaterialAssignment.coefIsComplexNumericOrFcnHdl(self.MassDensity, loc, state);
        tf(3) = pde.ThermalMaterialAssignment.coefIsComplexNumericOrFcnHdl(self.SpecificHeat, loc, state);
        
        if any(tf)
            error(message('pde:materialPropertySpecification:invalidMtlPropValueComplex'));
        end
        
    end

end % end of methods


methods(Static, Access = private)
    function tf = coefDefined(coef)
        tf = (isnumeric(coef) && ~(isscalar(coef) && coef == 0) || isa(coef,'function_handle'));
    end

    function tf = coefIsComplexNumericOrFcnHdl(coef, loc, state)
        tf = false;
        if ~pde.ThermalMaterialAssignment.coefDefined(coef)
            return
        end
        if isnumeric(coef)
            if ~isreal(coef)
                tf = true;
            end
        else % isa(coef, 'function_handle')
            res = coef(loc, state);
            if ~isreal(res)
                tf = true;
            end
        end
    end

    
    function ok=ValidateRegionID(rgnid)
        % Must be real(non-complex), full, natural number.
        if ~isreal(rgnid) || ~all(rgnid(:) > 0) || issparse(rgnid) || any(mod(rgnid(:),1)~=0)
            error(message('pde:materialPropertySpecification:invalidRegionID'));
        end
        ok = true;
    end

    function ok=CoefPrecheck(coef,coefName)
        
      validateattributes(coef,{'double','function_handle'},{},'',coefName)
        if isa(coef,'function_handle')
            validateattributes(coef,{'function_handle'},{'scalar'},'',coefName);
        else
            validateattributes(coef,{'double'},{'nonnan','nonempty','nonsparse','real','finite','nonnegative'},coefName);
            if all(coef==0) % Matrix form of coef can have some zeros, but not all.
                error(message('pde:materialPropertySpecification:invalidMtlPropNonPositive'));
            end
        end

        ok = true;
    end

    function checkMassDensitySize(dcoef, systemsize)
        if isscalar(dcoef) && ~isa(dcoef, 'function_handle') && isempty(dcoef)
            return
        elseif isa(dcoef, 'function_handle')
            return
        end
        dveclen = numel(dcoef);
        lengthok = (dveclen == 0 || dveclen == 1 || dveclen == systemsize || ...
            dveclen == (systemsize*(systemsize+1)/2) || dveclen == systemsize*systemsize);
        if ~(isvector(dcoef)  && iscolumn(dcoef) && lengthok)
            error(message('pde:materialPropertySpecification:invalidMassDensitySize'));
        end
    end

    function checkSpecificHeatSize(cpcoef,systemsize)
        %specific heat should have same size as mass density,
        %Use mass density check
        try
        pde.ThermalMaterialAssignment.checkMassDensitySize(cpcoef,systemsize)
        catch
            error(message('pde:materialPropertySpecification:invalidSpecificHeatSize'));
        end
    end

    function checkThermalConductivitySize(ccoef, systemsize, ndims)
        if isscalar(ccoef) && ~isa(ccoef, 'function_handle') && isempty(ccoef)
            return
        elseif isa(ccoef, 'function_handle')
            return
        end
        cveclen = numel(ccoef);
        if ndims == 2
            lengthok = (cveclen == 1 || cveclen == 2 || cveclen == 3 || ...
                cveclen == 4 || cveclen == systemsize || ...
                cveclen == 2*systemsize || cveclen == 3*systemsize || ...
                cveclen == 4*systemsize || cveclen == 2*systemsize*((2*systemsize)+1)/2 || ...
                cveclen == 4*(systemsize^2));
            if ~(isvector(ccoef)  && iscolumn(ccoef) && lengthok)
                error(message('pde:materialPropertySpecification:invalidThermConductivitySize'));
            end
        else
            lengthok = (cveclen == 1 || cveclen == 3 || cveclen == 6 || ...
                cveclen == 9 || cveclen == systemsize || ...
                cveclen == 3*systemsize || cveclen == 6*systemsize || ...
                cveclen == 9*systemsize || cveclen == 3*systemsize*((3*systemsize)+1)/2 || ...
                cveclen == 9*(systemsize^2));
            if ~(isvector(ccoef)  && iscolumn(ccoef) && lengthok)
                error(message('pde:materialPropertySpecification:invalidThermConductivitySize'));
            end
        end
    end
end

methods(Hidden=true)
    function delete(self)
        if ~isempty(self.RecordOwner) && isvalid(self.RecordOwner)
            self.RecordOwner.delistMaterialAssignments(self);
        end
    end
end


properties (Hidden = true, Access='private')
    RecordOwner;
end



end
