classdef (Sealed) MaterialAssignmentRecords < pde.PDECompactDisplay
% MaterialAssignmentRecords - Assignments of the material properties
%    The MaterialAssignmentRecords holds a record of all the assignments of
%    material properties across the geometric domain. This may be a single
%    assignment throughout the domain or multiple assignments that vary
%    across various subdomains. This object cannot be created directly, it
%    can only be created by the ThermalModel class.
%
% MaterialAssignmentRecords methods:
%    findThermalProperties - Find thermal properties for a geometric region
%
% MaterialAssignmentRecords properties:
%    MaterialAssignments - Vector of material property assignments
%
% See also pde.ThermalModel, pde.ThermalModel/thermalProperties

% Copyright 2016-2022 The MathWorks, Inc.

properties
    % MaterialAssignments - Vector of material property assignments
    %    A vector containing the material property assignments to the
    %    geometric domain or subdomains. Each entry in the vector is a
    %    pde.ThermalMaterialAssignment object. This object defines the
    %    material properties on a designated domain or subdomain of a
    %    thermal model. Material properties and host domains are defined
    %    using the thermalProperties method of the ThermalModel.

    MaterialAssignments;
end

methods
    % Method declaration
    mp = findThermalProperties(self, varargin)
    tf = materialPropertiesSpecified(self, varargin)
    tf = globalMaterialProperties(self, varargin)
end

methods (Hidden=true)
    function self = MaterialAssignmentRecords(tpdem)
        self.MaterialAssignments = [];
        self.ParentModel = tpdem;
    end

    function tf = allMaterialPropertiesNumeric(self)
        tf = true;
        if ~self.materialPropertiesSpecified()
            error(message('pde:materialPropertySpecification:subdomainsWithoutMtlProps'));
        end

        for i = 1:self.nRegion
            mp = self.findThermalProperties(self.qregiontype,i);
            if ~numericMaterialProperties(mp)
                tf = false;
                return
            end
        end

    end


    function performSolverPrechecks(self)
        if ~self.materialPropertiesSpecified()
            error(message('pde:materialPropertySpecification:subdomainsWithoutMtlProps'));
        end
        if ~consistentMaterialPropertyDefinition(self)
            error(message('pde:materialPropertySpecification:inconsistentMtlProps'));
        end

        syssize = 1;

        nf = self.ParentModel.Geometry.NumFaces;
        nc = self.ParentModel.Geometry.NumCells;


        for i = 1:self.nRegion
            ca = self.findThermalProperties(self.qregiontype,i);
            ca.performSolverPrecheck(syssize, nf, nc);
        end
    end

    function coefstruct = packProperties(self)

        if ~isempty(self.ParentModel.HeatSources)
            isHeatSourceGlobal = self.ParentModel.HeatSources.globalHeatSource();
        else
            isHeatSourceGlobal = true;
        end

        msh = self.ParentModel.Mesh;
        ma = msh.MeshAssociation;
        numelems = size(msh.Elements,2);
        if self.globalMaterialProperties() && isHeatSourceGlobal
            [loc, state] = self.ParentModel.getSampleLocState(1);
            coefstruct.ElementsInSubdomain = (1:numelems)';
            coefstruct.NumElementsInSubdomain = numelems;
            ca = self.MaterialAssignments(end);
            thisstruct.m{1} = 0;
            thisstruct.d{1} = pde.MaterialAssignmentRecords.computeDcoefficient(ca.MassDensity,ca.SpecificHeat,self.ParentModel.AnalysisType);
            thisstruct.c{1} = pde.MaterialAssignmentRecords.conductivityCCoefficient(ca.ThermalConductivity,self.ParentModel.AnalysisType);
            thisstruct.a{1} = 0;
            thisstruct.f{1} = 0;
            coefstruct.Coefficients = thisstruct;
            coefstruct.IsComplex = hasComplexCoefficient(ca, loc, state);
        else

            eIDtocellID = ma.SubdomainAssociativity;
            coefstruct.ElementsInSubdomain = zeros(numelems,1)';
            coefstruct.NumElementsInSubdomain = zeros(self.nRegion,1);
            complexcoef = false(3,1);
            endidx = 0;

            for i = 1:self.nRegion
                [loc, state] = self.ParentModel.getSampleLocState(i);
                celems = eIDtocellID(2,i)-eIDtocellID(1,i)+1;
                numcelems = celems;
                startidx = endidx;
                endidx = startidx+numcelems;
                startidx = startidx+1;
                coefstruct.NumElementsInSubdomain(i) = numcelems;
                coefstruct.ElementsInSubdomain(startidx:endidx)=eIDtocellID(1,i):eIDtocellID(2,i);
                ca = self.findThermalProperties(self.qregiontype,i);
                thisstruct.m{i} = 0;
                thisstruct.d{i} = pde.MaterialAssignmentRecords.computeDcoefficient(ca.MassDensity,ca.SpecificHeat,self.ParentModel.AnalysisType);
                thisstruct.c{i} = pde.MaterialAssignmentRecords.conductivityCCoefficient(ca.ThermalConductivity,self.ParentModel.AnalysisType);
                thisstruct.a{i} = 0;
                thisstruct.f{i} = 0;
                coefstruct.Coefficients = thisstruct;
                complexcoef = complexcoef | hasComplexCoefficient(ca, loc, state);
            end
            coefstruct.IsComplex = complexcoef;
        end
    end
end

methods(Hidden=true, Access=private)
    % Returns true if each material assignment has consistent definition across all
    % subdomains. For example, mass density defined in one region is
    % also defined in all other regions.
    function tf = consistentMaterialPropertyDefinition(self)
        tf = globalMaterialProperties(self);
        if tf
            return
        end
        tf = true;

        thiscoef = self.findThermalProperties(self.qregiontype, 1);
        for i = 2:self.nRegion
            othercoef = self.findThermalProperties(self.qregiontype, i);
            if ~thiscoef.coefficientsMatch(othercoef)
                tf = false;
                return;
            end
        end
    end

    function rType = qregiontype(self)
        if self.ParentModel.IsTwoD
            rType = 'face';
        else
            rType = 'cell';
        end
    end

    function n = nRegion(self)
        if self.ParentModel.IsTwoD
            n = self.ParentModel.Geometry.NumFaces;
        else
            n = self.ParentModel.Geometry.NumCells;
        end
    end
end

methods(Static, Hidden=true, Access=private)
    function c = conductivityCCoefficient(k,analysisType)
        if contains(analysisType,'axisymmetric')
            c = pde.MaterialAssignmentRecords.axisymmetricK(k);
        else
            c = k;
        end
    end

    function c = axisymmetricK(k)
        if isa(k,'function_handle')
            c = @axisymmetricKFcn;
        else
            c = @(region,state) 2*pi*region.x.*k;
        end
        
        function c = axisymmetricKFcn(region,state)
            regionAxsym = region;
            regionAxsym.r = region.x;
            regionAxsym.z = region.y;
            
            stateAxsym = state;
            stateAxsym.ur = state.ux;
            stateAxsym.uz = state.uy;

            c = 2*pi*region.x.*k(regionAxsym,stateAxsym);
        end
    end
    
    function dcoef = computeDcoefficient(massDensity,specificHeat,analysisType)
        if ~isa(massDensity, 'function_handle') && ~isa(specificHeat, 'function_handle')
            d = massDensity.*specificHeat;
        elseif ~isa(massDensity, 'function_handle') && isa(specificHeat, 'function_handle')
            d = @(region,state) massDensity.*specificHeat(region,state);
        elseif isa(massDensity, 'function_handle') && ~isa(specificHeat, 'function_handle')
            d = @(region,state) massDensity(region,state).*specificHeat;
        else
            d = @(region,state) massDensity(region,state).*specificHeat(region,state);
        end
        
        if ~contains(analysisType,'axisymmetric')
            dcoef = d;
            return;
        else
            dcoef = pde.MaterialAssignmentRecords.axisymmetricMass(d);
        end
    end
    
    function dcoef = axisymmetricMass(d)
        if isa(d,'function_handle')
            dcoef = @axisymmetricDFcn;
        else
            dcoef = @(region,state) 2*pi*region.x.*d;
        end
        function dcoef = axisymmetricDFcn(region,state)
            regionAxsym = region;
            regionAxsym.r = region.x;
            regionAxsym.z = region.y;
            
            stateAxsym = state;
            stateAxsym.ur = state.ux;
            stateAxsym.uz = state.uy;
            
            dcoef = 2*pi*region.x.*d(regionAxsym,stateAxsym);
        end
    end
        
    
end

methods(Hidden=true, Access={?pde.ThermalMaterialAssignment})
    function delistMaterialAssignments(self, mtltodelete)
        nummtl = numel(self.MaterialAssignments);
        for i = 1:nummtl
            thismtl = self.MaterialAssignments(i);
            if thismtl == mtltodelete
                self.MaterialAssignments(i) = [];
                break
            end
        end
        nummtl = numel(self.MaterialAssignments);
        if nummtl == 0 && isvalid(self.ParentModel)
            self.ParentModel.delistMaterialProperties();
        end
    end
end


properties (Hidden = true, SetAccess='private')
    ParentModel;
end



end
