classdef ElectromagneticModel < pde.EquationModel
    %pde.ElectromagneticModel - Container that represents an electromagnetic analysis model
    %  ElectromagneticModel is a container that stores the geometry,
    %  material properties, sources, boundary conditions, and mesh defining
    %  an electromagnetic problem. You cannot create ElectromagneticModel
    %  directly. Instead, use the CREATEPDE function. ElectromagneticModel
    %  provides functions for defining or importing a geometry, generating
    %  and storing a mesh for the geometry, specifying and storing
    %  properties of the material, applying and storing boundary conditions
    %  that define the behavior at the boundary of the geometry, and
    %  solving the model. The properties of ElectromagneticModel let you
    %  access these parameters of the model.
    %
    % ElectromagneticModel methods:
    %  importGeometry               - Import an STL geometry
    %  geometryFromEdges            - Create a 2-D geometry from DECSG or PDEGEOM
    %  geometryFromMesh             - Create geometry from a mesh
    %  electromagneticProperties    - Specify electromagnetic material properties over a domain or subdomain
    %  electromagneticSource        - Specify an electromagnetic source over a domain or subdomain
    %  electromagneticBC            - Apply electromagnetic boundary conditions to the geometric boundary
    %  generateMesh                 - Generate a mesh from the geometry
    %  solve                        - Solve the model
    %
    % ElectromagneticModel properties:
    %
    %  AnalysisType            - Type of analysis: 'electrostatic','conduction','magnetostatic', 
    %                            'electrostatic-axisymmetric', 'magnetostatic-axisymmetric', or
    %                            'harmonic'
    %  Geometry                - Geometric representation of the domain
    %  MaterialProperties      - Electromagnetic material properties within the domain or subdomain
    %  Sources                 - Electromagnetic source within the domain or subdomain
    %  BoundaryConditions      - Electromagnetic boundary conditions applied to the geometric boundary
    %  VacuumPermeability      - Permeability of vacuum (consistent with the units of the model)
    %  VacuumPermittivity      - Permittivity of vacuum (consistent with the units of the model)
    %  Mesh                    - Discretization of the domain
    %  SolverOptions           - Algorithm options for the solvers
    %  FieldType               - Type of field for harmonic analysis:
    %                            electric or magnetic
    %
    %    See also createpde,
    %    pde.ElectromagneticModel/electromagneticProperties,
    %    pde.ElectromagneticModel/electromagneticBC,
    %    pde.ElectromagneticModel/electromagneticSource,
    %    pde.ElectromagneticModel/electromagneticIC
    %    pde.ElectromagneticModel/solve, pde.PDEModel
    
    % Copyright 2020-2022 The MathWorks, Inc.
    
    properties
        AnalysisType;
    end
    
    properties
        % MaterialProperties - Electromagnetic material properties within the domain or subdomain
        %    An object that contains the electromagnetic material property
        %    assignments within the geometric domain. Use the
        %    electromagneticProperties function to assign material
        %    properties to the geometry. Multiple assignments are all
        %    recorded in this object.
        %
        %    See also pde.ElectromagneticMaterialAssignmentRecords/findElectromagneticProperties
        MaterialProperties;
        
        % Sources - Electromagnetic source within the domain or subdomain
        %    An object that contains the electromagnetic source or
        %    excitation assignments within the geometric domain. Use the
        %    electromagneticSource function to assign electromagnetic
        %    sources to the geometric domains. Multiple assignments are all
        %    recorded in this object.
        %
        %    See also pde.ElectromagneticSourceAssignmentRecords/findElectromagneticSources
        Sources;
        
        % BoundaryConditions - Electromagnetic boundary conditions applied to the geometric boundary
        %    An object that contains the electromagnetic boundary condition
        %    assignments. Use the electromagneticBC function to assign
        %    boundary conditions to the geometric boundary. Multiple
        %    assignments are all recorded in this object.
        %
        %    See also pde.ElectromagneticBCAssignmentRecords/findElectromagneticBC
        BoundaryConditions;
        
        % VacuumPermittivity - Permittivity of vacuum (consistent with the units of the model)
        %  The permittivity of vacuum is a physical constant the applies to
        %  the entire model. Specify a value that is consistent with the
        %  units of the model parameters. For example, if the model
        %  parameters are in the SI system of units, then the value must be
        %  8.8541878128E-12.
        VacuumPermittivity;
        
        % VacuumPermeability - Permeability of vacuum (consistent with the units of the model)
        %  The permeability of vacuum is a physical constant that applies
        %  to the entire model. Specify a value that is consistent with the
        %  units of the model parameters. For example, if the model
        %  parameters are in the SI system of units, then the value must be
        %  1.2566370614E-6.
        VacuumPermeability;

        % FieldType - Choice of electric or magnetic field for harmonic analysis. 
        % The default is electric.
        FieldType = "electric";
    end
    
    methods
        function self = ElectromagneticModel(analysisType)
            arguments
                analysisType char = 'magnetostatic'
            end
            self@pde.EquationModel(1)
            self.AnalysisType = analysisType;
            self.IsTimeDependent = false;
        end
        
        mtl = electromagneticProperties(self,varargin)
        src = electromagneticSource(self,varargin)
        tbc = electromagneticBC(self,varargin)
        eic = electromagneticIC(self,varargin)
        sol = solve(self,varargin);

        function set.AnalysisType(self, inputAnalysis)
            validateattributes(inputAnalysis,{'string','char'},{'nonempty'});
            validAnalysisStrings  = ["electrostatic",...
                "electrostatic-axisymmetric",...
                "magnetostatic",...
                "magnetostatic-axisymmetric",...
                "harmonic",...
                "harmonic-axisymmetric",...
                "conduction"];
            inputAnalysis= pde.EquationModel.validateAnalysisType(inputAnalysis, validAnalysisStrings);
            if ~isempty(self.AnalysisType) %Change analysis
                pde.EquationModel.validateAnalysisChange(self.AnalysisType, inputAnalysis);
            end
            self.AnalysisType = string(inputAnalysis);
        end

        function set.MaterialProperties(self, mpr)
            if(~isempty(mpr) && ~isa(mpr, 'pde.ElectromagneticMaterialAssignmentRecords'))
                error(message('pde:ElectromagneticModel:invalidMaterialProperties'));
            end
            if ~isscalar(mpr) && ~isempty(mpr)
                error(message('pde:ElectromagneticModel:nonScalarMaterialProperties'));
            end
            self.MaterialProperties = mpr;
        end
        
        function set.Sources(self, scr)
            if(~isempty(scr) && ~isa(scr, 'pde.ElectromagneticSourceAssignmentRecords'))
                error(message('pde:ElectromagneticModel:invalidSource'));
            end
            if ~isscalar(scr) && ~isempty(scr)
                error(message('pde:ElectromagneticModel:nonScalarSource'));
            end
            self.Sources = scr;
        end
        
        function set.BoundaryConditions(self, bcr)
            if(~isempty(bcr) && ~isa(bcr, 'pde.ElectromagneticBCAssignmentRecords'))
                error(message('pde:ElectromagneticModel:invalidBoundaryConditions'));
            end
            if ~isscalar(bcr) && ~isempty(bcr)
                error(message('pde:ElectromagneticModel:nonScalarBoundaryConditions'));
            end
            self.BoundaryConditions = bcr;
        end
        
        function set.VacuumPermittivity(self, val)
            self.validateVacuumPerm('VacuumPermittivity',val);
            self.VacuumPermittivity = val;
        end
        
        function set.VacuumPermeability(self, val)
            self.validateVacuumPerm('VacuumPermeability',val);
            self.VacuumPermeability = val;
        end

        function set.FieldType(self,val)
            if ~( isequal(lower(string(val)),"electric")...
                    || isequal(lower(string(val)),"magnetic") )
                error(message("pde:ElectromagneticModel:invalidFieldType"));
            end
            self.FieldType = lower(string(val));
        end

    end
    
    methods (Access = protected,Hidden = true)
        function group = getPropertyGroups(self)
            if contains(self.AnalysisType,'magnetostatic')
            group = matlab.mixin.util.PropertyGroup({'AnalysisType','Geometry', ...
                'MaterialProperties','Sources',...
                'BoundaryConditions','InitialConditions',...
                'VacuumPermeability','Mesh'});
            elseif contains(self.AnalysisType,'electrostatic') 
                 group = matlab.mixin.util.PropertyGroup({'AnalysisType','Geometry', ...
                'MaterialProperties','Sources',...
                'BoundaryConditions','VacuumPermittivity','Mesh'});
            elseif contains(self.AnalysisType,'conduction')
                group = matlab.mixin.util.PropertyGroup({'AnalysisType','Geometry', ...
                'MaterialProperties','BoundaryConditions','Mesh'});
            elseif startsWith(self.AnalysisType,'harmonic','IgnoreCase',true)
                if isequal(self.FieldType,'magnetic')
                    group = matlab.mixin.util.PropertyGroup({'AnalysisType','Geometry', ...
                        'MaterialProperties','Sources',...
                        'BoundaryConditions','VacuumPermittivity',...
                        'VacuumPermeability','Mesh','FieldType'});
                else
                    group = matlab.mixin.util.PropertyGroup({'AnalysisType','Geometry', ...
                        'MaterialProperties','Sources',...
                        'BoundaryConditions','VacuumPermittivity',...
                        'VacuumPermeability','Mesh','FieldType'});
                end
            end
        end
    end
    
    properties (SetAccess='private',Hidden = true)
        IsTimeDependent;
        DCoefficientFlag = true; % EM PDE if first-order in time.
        possibleDomainRegions = ["Face","Cell"]; % Face for 2-D, Cell for 3-D.
        possibleBoundaryRegions = ["Edge", "Face"]; % Edge for 2-D, Face for 3-D.
    end
    
    methods(Hidden)

        function msh = generateMesh(self,varargin)
            % Redirects to base class after setting the mesh to linear if
            % the analysis type is harmonic or 3-D magnetostatics.
            argsToPass = varargin;
            if (contains(self.AnalysisType,'harmonic') && ~isempty(self.Geometry)) ||...
                    ( (contains(self.AnalysisType,'magnetostatic') && ~isempty(self.Geometry) && ~self.IsTwoD) )
                % Get user specified geometric order.
                p = inputParser;
                p.KeepUnmatched = true;
                addParameter(p,'GeometricOrder',[]);
                parse(p,varargin{:});
                inputGeometricOrder = p.Results.GeometricOrder;

                if isempty(inputGeometricOrder)
                    argsToPass = [varargin(:)',{'GeometricOrder'},{'linear'}];
                elseif contains(self.AnalysisType,'magnetostatic')...
                        && startsWith('quadratic',inputGeometricOrder,'IgnoreCase',true)
                    error(message('pde:ElectromagneticModel:meshNotLinearTetFor3DMS'));
                elseif contains(self.AnalysisType,'harmonic')...
                        && startsWith('quadratic',inputGeometricOrder,'IgnoreCase',true)
                    error(message('pde:ElectromagneticModel:meshNotLinearTetHarmonic')); 
                end

            end
            % call base class generateMesh with correct options.
            msh = generateMesh@pde.EquationModel(self,argsToPass{:});
        end

        function geometricConsistencyCheck(self,gm)
            if isempty(gm)
                return;
            end
            if contains(self.AnalysisType,'axisymmetric') && ~gm.IsTwoD
                error(message('pde:ElectromagneticModel:invalid3DGeometryFor2DModel'));
            end
        end

        function ok = interfaceConditionsAllowed(self)
            ok = self.AnalysisType == "conduction";
        end
        
        function [coefstruct,constCoeffs] = packCoefficients(self,modCoef)
            if nargin < 2
                modCoef = [];
            end
            mtlstruct = self.MaterialProperties.packProperties(modCoef);
            emSrcMag = [];
            if ~isempty(self.Sources)
                emSrcStr = self.Sources.packElectromagneticSource();
                if startsWith(string(self.AnalysisType),"magnetostatic",'IgnoreCase',true) && ...
                        any(~cellfun(@isempty,{self.Sources.SourceAssignments.Magnetization}))
                    emSrcMag = self.Sources.packElectromagneticSource("Magnetization");
                end
                constantSource = self.Sources.AllNumericAssignments;
            else % No source was specified
                zeroSource = 0;
                if isequal(self.AnalysisType,"magnetostatic") && ~self.IsTwoD
                    zeroSource = [0;0;0]; % Vector source for 3-D MS.
                end
                if contains(self.AnalysisType,'harmonic')
                    if self.IsTwoD
                        zeroSource = [0 ; 0];
                    else
                        zeroSource = [0 ; 0 ; 0];
                    end
                end
                emSrcStr.Coefficients.f = repmat({zeroSource},1,numel(mtlstruct.Coefficients.c));
                constantSource = true;
            end
            
            coefstruct.ElementsInSubdomain = mtlstruct.ElementsInSubdomain;
            coefstruct.NumElementsInSubdomain = mtlstruct.NumElementsInSubdomain;
            coefstruct.Coefficients.m = repmat({0},1,numel(mtlstruct.Coefficients.c));
            coefstruct.Coefficients.d = repmat({0},1,numel(mtlstruct.Coefficients.c));
            coefstruct.Coefficients.c = mtlstruct.Coefficients.c;
            coefstruct.Coefficients.a = repmat({0},1,numel(mtlstruct.Coefficients.c));
            coefstruct.Coefficients.f = emSrcStr.Coefficients.f;
            if contains(self.AnalysisType,'harmonic')
                coefstruct.Coefficients.m = mtlstruct.Coefficients.m;
            end
            if ~isempty(emSrcMag)
                coefstruct.Coefficients.fc = emSrcMag.Coefficients.f;
            end
            
            coefstruct.IsComplex = [false;false;false;false;false];
            
            constCoeffs = constantSource && self.MaterialProperties.AllNumericAssignments;
        end
        
        function rType = validDomain(self)
            rType = 'Cell';
            if self.IsTwoD
                rType = 'Face';
            end
        end

        function rType = validBoundary(self)
            rType = 'Face';
            if self.IsTwoD
                rType = 'Edge';
            end
        end

        function maxID = maxDomainID(self)
            domainName = self.validDomain;
            maxID = self.Geometry.(['Num',domainName,'s']);
        end

        function maxID = maxBoundaryID(self)
            boundaryName = self.validBoundary;
            maxID = self.Geometry.(['Num',boundaryName,'s']);
        end
        [edges, edge_elem] = constructEdgeElementConnectivity(self);
        [M,K,F,H,R,faceWithBCs] = discretizedModelNedelec(self,coefstruct,edges,edge_elem,numEdges,nodalPotential,computeOnlyNonlinear);
        [M,K,F,H,R,faceWithBCs] = discretizedModelHarmonicNedelec(self,coefstruct,edges,edge_elem,numEdges,frequency);
        [H, R, faceWithBCs] = assembleDirichletBCNedelec(self,edges);
        u = solveNedelecHarmonic(self,coefstruct,frequency);
        modifiedSigma = computeFarFieldBCCoefficients(self);
        F = addCurlComponent(self,p,t,coefstruct,u0,edgeElem,numEdges);
        u = solveNodalESorMS(self,coefstruct,constCoeffsBCs);
        [u, conductorIDs] = solveConduction(self,coefstruct,constCoeffsBCs);
        [u, modelForResults] = solveHarmonic(self,coefstruct,frequency);
        [u, modelForResults] = solveMSNedelec(self,coefstruct);
        Res = solveNedelecMS3D(self,coefstruct);
    end

    methods(Access={?pde.ElectromagneticMaterialAssignment,...
                    ?pde.ElectromagneticSourceAssignment,...
                    ?pde.ElectromagneticBCAssignment,...
                    ?pde.ElectromagneticICAssignment})
        function [complexityInput,sizeInput,res] = validateFcnHdl(self, coef, loc, state)
            if isnumeric(coef)
                complexityInput = ~isreal(coef);
                sizeInput = numel(coef);
                res = coef;
                return
            end

            % %%% Function handle checks %%%

            % 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:pdeCoefficientSpecification:invalidFcnHandleArgs'));
            end

            % Try and rethrow, if user function doesn't run to completion.
            try
                res = coef(loc, state);
            catch ME
                throwAsCaller(ME);
            end

            %User function succeeds, check output.
            complexityInput = ~isreal(res);
            if ~isa(res,'double')
                error(message('pde:pdeCoefficientSpecification:invalidFcnHandleOutputType'));
            end
            sizeInput = size(res,1);

            % Check for nonlinearity
            %state = rmfield(state,'frequency'); % frequency dependency OK.
            for fld = fields(state)'
                if contains(fld{1},'frequency')
                    continue
                end
                state.(fld{1}) = state.(fld{1})*NaN;
            end

            res = coef(loc,state);

            if any(isnan(res))
                error(message('pde:ElectromagneticModel:nonlinearEMModel'));
            end
        end
    end

    methods(Static, Hidden=true)
        femodel = imposeBCOnAxis(femodel);
    end

    methods(Static)
        function ok = validateVacuumPerm(coef,val)
            % validate both vacuum permittivity and permeability
            if isempty(val)
                % allow empty values as a default perm
                ok = true;
                return
            end
            validateattributes(val,{'double'},{'nonnan','nonsparse','real','finite','scalar','positive'},coef);
            ok = true;
        end

        function regionInfo = getRegionText(region,regionID)
            if isscalar(regionID)
                regionInfo = num2str([region ' ' num2str(regionID)]);
            else
                regionInfo = num2str([region ' [' num2str(regionID) ']']);
            end
        end

    end % end of static methods

    methods(Hidden = true, Access = {?pde.ElectromagneticMaterialAssignmentRecords})
        function delistMaterialProperties(self)
            if isvalid(self)
                self.MaterialProperties = [];
            end
        end
    end
    
  methods(Hidden = true, Access = {?pde.ElectromagneticSourceAssignmentRecords})
        function delistElectromagneticSource(self)
            if isvalid(self)
                self.Sources = [];
            end
        end
    end
    
    methods(Hidden = true, Access = {?pde.ElectromagneticBCAssignmentRecords})
        function delistBoundaryConditions(self)
            if isvalid(self)
                self.BoundaryConditions = [];
            end
        end
    end  
end