function R = solve(self,options)
% solve - Solve electromagnetic model
% R = solve(emagmodel) solves an electromagnetic model and returns the
% solution as a results object that corresponds to the analysis type.
%
%  Example: Find the electrostatic potential in an air-filled annular quadrilateral frame.
%
%      % Create an electromagnetic model for electrostatic analysis.
%      emagmodel = createpde('electromagnetic','electrostatic');
%
%      % Import and plot the geometry.
%      importGeometry(emagmodel,'Frame.STL');
%      pdegplot(emagmodel,'EdgeLabels','on')
%
%      % Specify the relative permittivity of the material.
%      electromagneticProperties(emagmodel,'RelativePermittivity',1.00059);
%
%      % Specify the electrostatic potential at the inner boundary.
%      electromagneticBC(emagmodel,'Voltage',1000,'Edge',[1 2 4 6]);
%
%      % Specify the electrostatic potential at the outer boundary.
%      electromagneticBC(emagmodel,'Voltage',0,'Edge',[3 5 7 8]);
%
%      % Generate the mesh.
%      generateMesh(emagmodel);
%
%      % Specify the vacuum permittivity value in the SI system of units and solve.
%      emagmodel.VacuumPermittivity = 8.8541878128E-12;
%      R = solve(emagmodel);
%
%      % Plot the electric potential using the Contour parameter to display equipotential lines.
%      pdeplot(emagmodel,'XYData', R.ElectricPotential,'Contour','on');
%
% See also pde.ElectromagneticModel, pde.ElectromagneticModel/electromagneticProperties,
%          pde.ElectromagneticModel/electromagneticBC,pde.ElectromagneticModel/electromagneticSource

% Copyright 2020-2022 The MathWorks, Inc.

arguments
    self pde.ElectromagneticModel
    options.Frequency (1,:) {mustBeNumeric,mustBeReal}
end

narginchk(1,3);
checkics = false;
performSolverPrecheck(self, checkics);

frequency = [];
if isfield(options,'Frequency')
    frequency = options.Frequency;
end

if contains(self.AnalysisType,'harmonic') && isempty(frequency)
    error(message('pde:ElectromagneticModel:HarmonicNoFrequency'));
end

allbcsnumeric = true;
if ~isempty(self.BoundaryConditions)
    numbcs = numel(self.BoundaryConditions.BCAssignments);
    for i = 1:numbcs
        thisbc = self.BoundaryConditions.BCAssignments(i);
        if ~thisbc.numericCoefficients()
            allbcsnumeric = false;
            break;
        end
    end
end

% compute info for farfield BCs if present
modCoef = @(loc,state) 0;
if ~isempty(self.BoundaryConditions) && ~any(isempty([self.BoundaryConditions.BCAssignments.FarField]))
    modCoef = computeFarFieldBCCoefficients(self);
end

[coefstruct,constCoeffs] = packCoefficients(self,modCoef);
constCoeffsBCs = constCoeffs && allbcsnumeric;

analysisType = self.AnalysisType;
if startsWith(analysisType,"electrostatic",'IgnoreCase',true)
    u = solveNodalESorMS(self,coefstruct,constCoeffsBCs);
    R = pde.ElectrostaticResults(self,u);
elseif startsWith(analysisType,"conduction",'IgnoreCase',true)
    [u,conductorIDs] = solveConduction(self,coefstruct,constCoeffsBCs);
    R = pde.ConductionResults(self,u,conductorIDs);
elseif startsWith(analysisType,"magnetostatic",'IgnoreCase',true)
    if self.IsTwoD
        u = solveNodalESorMS(self,coefstruct,constCoeffsBCs);
        modelForResults = self;
    else
        [u, modelForResults] = solveMSNedelec(self,coefstruct);
    end
    R = pde.MagnetostaticResults(modelForResults,u);
elseif startsWith(analysisType,"harmonic",'IgnoreCase',true)
    [u, modelForResults] = solveHarmonic(self,coefstruct,frequency);
    R = pde.HarmonicResults(modelForResults,u,frequency);
end
end % End of solve method.
