classdef (Sealed) ElectromagneticMaterialAssignmentRecords < pde.PDECompactDisplay
    % ElectromagneticMaterialAssignmentRecords - Assignments of electromagnetic material properties
    %  ElectromagneticMaterialAssignmentRecords holds a record of all the
    %  assignments of electromagnetic material properties across the geometric
    %  domains. This can 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
    %  ElectromagneticModel class.
    %
    % ElectromagneticMaterialAssignmentRecords methods:
    %    findElectromagneticProperties - Find electromagnetic properties for a geometric region
    %
    % ElectromagneticMaterialAssignmentRecords properties:
    %    MaterialAssignments - Vector of material property assignment objects
    % 
    % See also pde.ElectromagneticModel, pde.ElectromagneticModel/electromagneticProperties
    
    % Copyright 2020-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.ElectromagneticMaterialAssignment object. This object defines the
        %    material properties on a designated domain or subdomain of a
        %    model. Material properties and host domains are defined using the
        %    electromagneticProperties method of the ElectromagneticModel object.
        MaterialAssignments;
    end
    
    methods
        % Method declaration
        mp = findElectromagneticProperties(self, varargin)
    end
    
    methods (Hidden=true, Access={?pde.EquationModel,?pde.ElectrostaticResults,?pde.MagnetostaticResults,?pde.ConductionResults})
        function self = ElectromagneticMaterialAssignmentRecords(empdem)
            self.MaterialAssignments = [];
            self.ParentModel = empdem;
            if self.ParentModel.AnalysisType == "conduction"
                 self.EmptyAllowed = true;
            end
        end
        
        function performSolverPrechecks(self)
            requiredProperty = "RelativePermittivity";
            analysisType = self.ParentModel.AnalysisType;
            if(startsWith(analysisType,"magnetostatic",'IgnoreCase',true))
                requiredProperty = "RelativePermeability";
            end
            
            if contains(string(analysisType),"harmonic")
                requiredProperty = ["RelativePermeability","RelativePermittivity"];
            end
            
            if contains(string(analysisType),"conduction")
                requiredProperty = "Conductivity";
            end

            validDomain = self.ParentModel.validDomain;
            maxDomainID = self.ParentModel.maxDomainID;
            for i = 1:maxDomainID
                for prop = requiredProperty
                    mp = self.findElectromagneticProperties(validDomain,i);
                    
                    if self.EmptyAllowed && isempty(mp)
                        % Skip validating empty mp in case of conduction.
                        continue
                    end

                    if isempty( mp.(prop) )
                        error(message('pde:ElectromagneticModel:missingRequiredMaterialProperties',prop,analysisType));
                    end

                    [loc, state] = self.ParentModel.getSampleLocState(i);

                    complexityAllowed = false;
                    if contains(string(analysisType),"harmonic")
                         %Complex coefficients are allowed only for harmonic analysis.
                        complexityAllowed = true;
                        state.frequency = 1; % just for checks
                    end
                    
                    validateTypeSizeComplexity(mp, loc, state, complexityAllowed);
                    
                    if ~numericMaterialProperties(mp)
                        self.AllNumericAssignments = false;
                    end
                end
            end
        end

        function mtlstruct = packPropertiesHarmonic(self,modCoef)
        % alternative function to pack properties that produces
        % coefficients consistent with the harmonic analysis type
            numelems = size(self.ParentModel.Mesh.Elements,2);
            mtlstruct.ElementsInSubdomain = (1:numelems)';
            eIDtocellID = self.ParentModel.Mesh.MeshAssociation.SubdomainAssociativity;
            mtlstruct.NumElementsInSubdomain = diff(reshape(eIDtocellID,2,[]))+1;
            validDomain = self.ParentModel.validDomain;
            maxDomainID  = self.ParentModel.maxDomainID;
            eps0 = self.ParentModel.VacuumPermittivity;
            mu0 = self.ParentModel.VacuumPermeability;
            if isempty(eps0)
                error(message('pde:ElectromagneticModel:missingVacuumPermittivity'))
            end
            if isempty(mu0)
                error(message('pde:ElectromagneticModel:missingVacuumPermeability'))
            end
            for i = 1:maxDomainID
                ca = self.findElectromagneticProperties(validDomain,i);
                if isequal(string(self.ParentModel.FieldType),"electric")
                    if isempty(ca.Conductivity)
                        ca.Conductivity = 0;
                    end
                    % coefficient for curl-curl term
                    if isa(ca.RelativePermeability,'function_handle')
                        c = @(location,state) 1./(ca.RelativePermeability(location,state)*mu0);
                    else
                        c = 1./(ca.RelativePermeability*mu0);
                    end
                    % coefficient for non-deriv term
                    % permittivity
                    if isa(ca.RelativePermittivity,'function_handle')
                        m = @(location,state) -eps0*ca.RelativePermittivity(location,state)*state.frequency^2;
                    else
                        m = @(location,state) (-eps0*ca.RelativePermittivity*state.frequency^2)*ones(size(location.x));
                    end
                    % conductivity
                    if isa(ca.Conductivity,'function_handle')
                        m = @(location,state) m(location,state) + 1i*ca.Conductivity(location,state)*state.frequency;
                    else
                        m = @(location,state) m(location,state) + 1i*ca.Conductivity*state.frequency;
                    end
                    % PML (use complex to guarantee type, real +
                    % complex(~,0) is real)
                    m = @(location,state) complex(m(location,state) + 1i*modCoef(location,state)*state.frequency);
                else % magnetic
                    % coefficient for curl-curl term
                    if isa(ca.RelativePermittivity,'function_handle')
                        c = @(location,state) 1./(ca.RelativePermittivity(location,state)*eps0);
                    else
                        c = 1./(ca.RelativePermittivity*eps0);
                    end
                    % coefficient for non-deriv term
                    % permeability
                    if isa(ca.RelativePermeability,'function_handle')
                        m = @(location,state) -mu0*ca.RelativePermeability(location,state)*state.frequency^2;
                    else
                        m = @(location,state) (-mu0*ca.RelativePermeability*state.frequency^2)*ones(size(location.x));
                    end
                    % PML (use complex to guarantee type, real +
                    % complex(~,0) is real)
                    m = @(location,state) complex(m(location,state) + 1i*modCoef(location,state)*state.frequency);
                end
                % modify coefs for axisymmetric
                if contains(self.ParentModel.AnalysisType,'axisymmetric')
                    coeffExponent = 1;
                    c = pde.ElectromagneticMaterialAssignmentRecords.axisymmetricC(c,coeffExponent);
                    m = pde.ElectromagneticMaterialAssignmentRecords.axisymmetricC(m,coeffExponent);
                end
                % assign coefficients
                mtlcoefs.c{i} = c;
                mtlcoefs.m{i} = m;
            end
            mtlstruct.Coefficients = mtlcoefs;
        end
        
        function mtlstruct = packProperties(self,modCoef)
            analysisType = self.ParentModel.AnalysisType;
            if contains(analysisType,'harmonic')
                mtlstruct = packPropertiesHarmonic(self,modCoef);
                return;
            end
            numelems = size(self.ParentModel.Mesh.Elements,2);
            mtlstruct.ElementsInSubdomain = (1:numelems)';
            eIDtocellID = self.ParentModel.Mesh.MeshAssociation.SubdomainAssociativity;
            mtlstruct.NumElementsInSubdomain = diff(reshape(eIDtocellID,2,[]))+1;
            validDomain = self.ParentModel.validDomain;
            maxDomainID  = self.ParentModel.maxDomainID;
            for i = 1:maxDomainID
                ca = self.findElectromagneticProperties(validDomain,i);
                mtlcoefs.c{i} = 0;
                if ~isempty(ca)
                    mtlcoefs.c{i} = self.emC(ca,analysisType);
                end
            end
            mtlstruct.Coefficients = mtlcoefs;
        end
        
        function c = emC(self,ca,analysisType)
            coeffExponent  = 1;
            if startsWith(analysisType,"electrostatic",'IgnoreCase',true)
                eps0 = self.ParentModel.VacuumPermittivity;
                if  isa(ca.RelativePermittivity,'function_handle')
                    relCoeff = @(region,state) ca.RelativePermittivity(region,state)*eps0;
                else
                    relCoeff = ca.RelativePermittivity*eps0;
                end
            end

            if startsWith(analysisType,"conduction",'IgnoreCase',true)
                if  isa(ca.Conductivity,'function_handle')
                    relCoeff = @(region,state) ca.Conductivity(region,state);
                else
                    relCoeff = ca.Conductivity;
                end
            end

            if startsWith(analysisType,"magnetostatic",'IgnoreCase',true)
                coeffExponent  = -1;
                curlEquationFactor = 1;
                if ~self.ParentModel.IsTwoD
                    % For 3-D MS PDE is curl(curl(Avector)) = fvector. All other static
                    % analyses are solved using the standard divergence equation: - div(c
                    % grad(u) = f. This flag accounts of the sign difference so that computed
                    % flux direction follows right-hand rule.
                    curlEquationFactor = -1;
                end
                mu0  = self.ParentModel.VacuumPermeability;
                if  isa(ca.RelativePermeability,'function_handle')
                    relCoeff = @(region,state) curlEquationFactor*1./(ca.RelativePermeability(region, state)*mu0);
                else
                    relCoeff = curlEquationFactor*1/(ca.RelativePermeability*mu0);
                end
            end
            
            c = relCoeff;
            if contains(analysisType,'axisymmetric')
                c = pde.ElectromagneticMaterialAssignmentRecords.axisymmetricC(relCoeff,coeffExponent);
            end
        end        
    end
    
    methods(Static, Hidden=true, Access=private)
        function c = axisymmetricC(relCoeff,coeffExponent)
            if isa(relCoeff,'function_handle')
                c = @axisymmetricKFcn;
            else
                c = @(region,state) 2*pi*(region.x.^coeffExponent)*relCoeff;
            end
            
            function c = axisymmetricKFcn(region,state)
                [regionAxisymm, stateAxisymm] = pde.EquationModel.getAxisymmStructs(region,state);               
                c = 2*pi*(region.x.^coeffExponent).*relCoeff(regionAxisymm,stateAxisymm);
            end
        end
    end
    
    methods(Hidden=true, Access={?pde.ElectromagneticMaterialAssignment})
        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
    
    methods(Hidden=true)
        function delete(self)
            nummtl = numel(self.MaterialAssignments);
            for i = 1:nummtl
                if isvalid(self.MaterialAssignments(i))
                    delete(self.MaterialAssignments(i));
                end
            end
            
            if isvalid(self.ParentModel)
                self.ParentModel.delistMaterialProperties();
            end
        end
    end

    properties (Hidden = true, SetAccess='private')
        ParentModel;
        AllNumericAssignments = true;
        EmptyAllowed = false; % Material properties must be specified on all subdomain, only exception being conduction analysis.
    end
end
