classdef (Sealed) ElectromagneticSourceAssignmentRecords < pde.PDECompactDisplay
    % ElectromagneticSourceAssignmentRecords - Assignments of electromagnetic sources
    % ElectromagneticSourceAssignmentRecords holds a record of all the
    %  assignments of electromagnetic sources 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.
    %
    % ElectromagneticSourceAssignmentRecords methods:
    %    findElectromagneticSource - Find electromagnetic source for a geometric region
    %
    % ElectromagneticSourceAssignmentRecords properties:
    %    SourceAssignments - Vector of electromagnetic source assignments
    %
    % See also pde.ElectromagneticModel, pde.ElectromagneticModel/electromagneticSource
    
    % Copyright 2020-2022 The MathWorks, Inc.
    
    properties
        % SourceAssignments - Vector of source assignments 
        % A vector containing the source assignments to the geometric
        % domain or subdomains. Each entry in the vector is a
        % pde.ElectromagneticSourceAssignment object. This object defines
        % the source per unit volume for the entire domain or subdomain of
        % a model. Source values and host domains are defined using the
        % electromagneticSource method of the ElectromagneticModel object.
        SourceAssignments;
    end
    
    methods
        % Method declaration
        mp = findElectromagneticSource(self, varargin)
    end
    
    methods(Hidden=true, Access={?pde.EquationModel, ?pde.ElectromagneticResults})
        function self = ElectromagneticSourceAssignmentRecords(spdem)
            self.SourceAssignments = [];
            self.ParentModel = spdem;
        end
        
        function performSolverPrechecks(self)
            validDomain = self.ParentModel.validDomain;
            maxDomainID = self.ParentModel.maxDomainID;
            for i = 1:maxDomainID
                ca = self.findElectromagneticSource(validDomain,i);
                if isempty(ca)
                    continue
                end

                [loc, state] = self.ParentModel.getSampleLocState(i);
                
                complexAllowed = false;
                if contains(self.ParentModel.AnalysisType,'harmonic','IgnoreCase',true)
                    % can be complex for harmonic
                    complexAllowed = true;
                end
                validateTypeSizeComplexity(ca, loc, state, complexAllowed);
                if ~numericSource(ca)
                    self.AllNumericAssignments = false;
                end
            end
        end
        
        function bodyloadstruct = packElectromagneticSource(self,varargin)

            validSourceInput = 'ChargeDensity';
            analysisType = self.ParentModel.AnalysisType;
            numVecComponents = 1;
            if(startsWith(analysisType,"magnetostatic",'IgnoreCase',true))
                validSourceInput = 'CurrentDensity';
                if ~self.ParentModel.IsTwoD
                    numVecComponents = 3; % One vector equation with three components.
                end
            end
            if contains(string(analysisType),"harmonic")
                validSourceInput = 'CurrentDensity';
                numVecComponents = 2;
                if ~self.ParentModel.IsTwoD
                    numVecComponents = 3;
                end
            end
            % overwrite valid input if requesting magnetization
            if (nargin > 1 && isequal(string(varargin{1}),"Magnetization"))
                validSourceInput = "Magnetization";
                numVecComponents = 2;
                if ~self.ParentModel.IsTwoD
                    numVecComponents = 3;
                end
            end
            validDomain = self.ParentModel.validDomain;
            maxDomainID = self.ParentModel.maxDomainID;
            for i = 1:maxDomainID
                ca = self.findElectromagneticSource(validDomain,i);
                if isempty(ca)
                    blcoefs.f{i}  = zeros(numVecComponents,1);
                    continue
                end
                if isempty(ca.(validSourceInput))
                    blcoefs.f{i}  = zeros(numVecComponents,1);
                    continue
                end

                analysisType = self.ParentModel.AnalysisType;

                if contains(analysisType,'harmonic')
                    % check size of density. 
                    % this is done here and not with the coefficient
                    % prechecks because the field type must be defined.
                    % This check must happen after solve.
                    CurrentDensity = ca.CurrentDensity;
                    if isa(ca.CurrentDensity,'function_handle')
                        [location, state] = self.ParentModel.getSampleLocState(1);
                        res = ca.CurrentDensity(location, state);
                    else
                        res = ca.CurrentDensity;
                    end

                    % check if input is vector with correct number of
                    % components
                    if ~(isequal(size(res),[numVecComponents,1]) || isequal(size(res),[1,numVecComponents])) && ...
                            isequal(string(self.ParentModel.FieldType),"electric")
                        error(message('pde:ElectromagneticModel:invalidSourceDimensionsElectric'))
                    elseif ( (numVecComponents == 2 && length(res) ~= 1) ||... % 2-D, f is scalar
                            (numVecComponents == 3 &&  ~(isequal(size(res),[numVecComponents,1]) || isequal(size(res),[1,numVecComponents]))) )... % 3-D, f has dim 3
                            && isequal(string(self.ParentModel.FieldType),"magnetic")
                        error(message('pde:ElectromagneticModel:invalidSourceDimensionsMagnetic'))
                    end

                    % transpose input if need be
                    if size(res,2) > 1
                        if ~isnumeric(ca.CurrentDensity)
                            CurrentDensity = @(location,state) ca.CurrentDensity(structfun(@transpose,location,'UniformOutput',false),state).';
                        else
                            CurrentDensity = CurrentDensity.';
                        end
                    end

                    % set up the f coefficients
                    if isequal(string(self.ParentModel.FieldType),"electric")
                        if ~isnumeric(CurrentDensity)
                            f = @(region,state) -1i*state.frequency*CurrentDensity(region,state);
                        else
                            f = @(region,state) -1i*state.frequency*CurrentDensity*ones(size(region.x));
                        end
                    else % magnetic
                        % add on fictitious zero in 2D (don't change
                        % computation, but have the correct num equations
                        % in 2D for passing to C++ code)
                        matProps = self.ParentModel.MaterialProperties.findElectromagneticProperties(validDomain,i); 
                        % no explicit error check here, findElectromagneticProperties does its own
                        
                        % if either density or permittivity is function
                        % handle, treat both as handles
                        if isa(CurrentDensity,'function_handle') || isa(matProps.RelativePermittivity,'function_handle')
                            if isa(matProps.RelativePermittivity,'function_handle')
                                permittivity = @(region,state) self.ParentModel.VacuumPermittivity.*matProps.RelativePermittivity(region,state);
                            else
                                permittivity = @(region,state) self.ParentModel.VacuumPermittivity.*matProps.RelativePermittivity;
                            end
                            if isa(CurrentDensity,'function_handle')
                                density = @(region,state) CurrentDensity(region,state);
                            else
                                density = @(region,state) CurrentDensity;
                            end
                        else
                            permittivity = self.ParentModel.VacuumPermittivity.*matProps.RelativePermittivity;
                            density = CurrentDensity;
                            if isequal(permittivity,0)
                                error(message('pde:ElectromagneticModel:invalidPermittivity'))
                            end
                        end
                        if self.ParentModel.IsTwoD
                            if ~isnumeric(CurrentDensity) || ~isnumeric(permittivity)
                                f = @(region,state) [density(region,state)./permittivity(region,state);0*region.x];
                            else
                                f = [density./permittivity;0];
                            end
                        else % 3D
                            if ~isnumeric(CurrentDensity) || ~isnumeric(permittivity)
                                f = @(region,state) density(region,state)./permittivity(region,state);
                            else
                                f = density./permittivity;
                            end
                        end
                    end
                    if contains(analysisType,'axisymmetric')
                        f = pde.ElectromagneticSourceAssignmentRecords.axisymmSource(f);
                    end
                end
                
                if (startsWith(analysisType,"electrostatic",'IgnoreCase',true))
                    if ~isnumeric(ca.ChargeDensity)
                        f = @(region,state) ca.ChargeDensity(region,state);
                    else
                        f = ca.ChargeDensity;
                    end
                    
                    if contains(analysisType,'axisymmetric')
                        f = pde.ElectromagneticSourceAssignmentRecords.axisymmSource(f);
                    end
                end

                if(startsWith(analysisType,"magnetostatic",'IgnoreCase',true) && ...
                        isequal(string(validSourceInput),"CurrentDensity"))
                    axisymmFactor = 1;
                    if contains(analysisType,'axisymmetric')
                        axisymmFactor = 2*pi;
                    end
                    if isnumeric(ca.CurrentDensity)
                        f = ca.CurrentDensity(:)*axisymmFactor;
                    elseif isa(ca.CurrentDensity,'function_handle')
                        f = @(r,s)ca.CurrentDensity(r,s)*axisymmFactor;
                    elseif isa(ca.CurrentDensity,'pde.ConductionResults')
                        f = @(r,s) conductorFlux(ca,r,s);
                    end
                end

                if(startsWith(analysisType,"magnetostatic",'IgnoreCase',true) && ...
                        isequal(string(validSourceInput),"Magnetization"))
                    % axisymmetric set up
                    axisymmFactor = 1;
                    if contains(analysisType,'axisymmetric')
                        axisymmFactor = 2*pi;
                    end
                    % incorporate permeability
                    matProps = self.ParentModel.MaterialProperties.findElectromagneticProperties(validDomain,i);
                    permeability = matProps.RelativePermeability;
                    if ~isnumeric(ca.Magnetization) || ~isnumeric(permeability)
                        if ~isnumeric(ca.Magnetization)
                            magnetization = @(r,s) ca.Magnetization(r,s);
                        else
                            magnetization = @(r,s) ca.Magnetization;
                        end
                        if ~isnumeric(permeability)
                            perm = @(r,s) permeability(r,s);
                        else
                            perm = @(r,s) permeability;
                        end
                        source = @(r,s) magnetization(r,s)./perm(r,s);
                    else % both permeability and magnetization are constant for this subdomain
                        source = ca.Magnetization(:)./permeability(:);
                    end
                    % define source
                    if ~isnumeric(source)
                        f = @(r,s)source(r,s)*axisymmFactor;
                    else
                        f = source(:)*axisymmFactor;
                    end
                end
                
                blcoefs.f{i} = f;
            end
            bodyloadstruct.Coefficients = blcoefs;
        end
    end
    
    methods(Static, Hidden=true, Access=private)
        function q = axisymmSource(fin)
            if isa(fin,'function_handle')
                q = @axisymmetricQFcn;
            else
                q = @(region,state) 2*pi*region.x.*fin;
            end
            function q = axisymmetricQFcn(region,state)
                [regionAxisymm, stateAxisymm] = pde.EquationModel.getAxisymmStructs(region,state);                
                q = 2*pi*region.x.*fin(regionAxisymm,stateAxisymm);
            end
        end
    end

    methods(Hidden=true, Access={?pde.ElectromagneticSourceAssignment})
        function delistSourceAssignments(self, bltodelete)
            numhs = numel(self.SourceAssignments);
            for i = 1:numhs
                thishs = self.SourceAssignments(i);
                if thishs == bltodelete
                    self.SourceAssignments(i) = [];
                    break
                end
            end
            numhs = numel(self.SourceAssignments);
            if numhs == 0 && isvalid(self.ParentModel)
                self.ParentModel.delistElectromagneticSource();
            end
        end
    end

    methods(Hidden=true)
        function delete(self)
            nummtl = numel(self.SourceAssignments);
            for i = 1:nummtl
                if isvalid(self.SourceAssignments(i))
                    delete(self.SourceAssignments(i));
                end
            end

            if isvalid(self.ParentModel)
                self.ParentModel.delistElectromagneticSource();
            end
        end
    end


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

function fcond = conductorFlux(ca, r,~)
intrpI = ca.CurrentDensity.interpolateCurrentDensity(r.x,r.y,r.z);
fcond = [intrpI.Jx';intrpI.Jy'; intrpI.Jz'];
end
