function R = reduce(self,options)
%reduce - Reduces the thermal model using modes
% R = reduce(modalmodel,'ModalResults',Rmodal) reduces a thermal analysis
% model to the modes in the modal analysis results object Rmodal.

% Copyright 2021-2023 The MathWorks, Inc.

arguments
    self pde.ThermalModel
    options.ModalResults pde.ModalThermalResults;
    options.NumModes (1,1) {mustBeA(options.NumModes,"double"), mustBePositive, mustBeInteger};
end

if ~isfield(options,'ModalResults')
    error(message('pde:ThermalModel:noModalResultsForReduce'));
end

Rm = options.ModalResults;

Np = size(Rm.ModeShapes,2);
if isfield(options,'NumModes')
    validateattributes(options.NumModes,{'double'},{'scalar','<=',Np},'NumModes')
    Np = options.NumModes;
end

coefstruct = packCoefficients(self);
tlist = [0, 1];
u0 = [];
if ~isempty(self.InitialConditions)
    u0 = self.InitialConditions.packInitialConditions();
end

[Kr,Mr,Fr,u0Modal, ~, femodel] = modalProjection(self,Rm,Np,coefstruct,u0,tlist);
Kr = (Kr+Kr')/2;
Mr = (Mr+Mr')/2;

if femodel.vf || femodel.vg
    error(message('pde:ThermalModel:loadsNotLinearForReduce'));
end

modalTransformationMat = [Rm.SnapshotsAverage, Rm.ModeShapes(:,1:Np)];
transientAnalysisType = "transient";
if contains(self.AnalysisType,"axisymmetric")
    transientAnalysisType = "transient-axisymmetric";
end
tempModelForReconstruction = createpde('thermal',transientAnalysisType);
tempModelForReconstruction.Geometry = self.Geometry;
tempModelForReconstruction.Mesh = self.Mesh;
numMtlassignments = numel(self.MaterialProperties.MaterialAssignments);
for i = 1:numMtlassignments
    smai = self.MaterialProperties.MaterialAssignments(i);
    tempModelForReconstruction.thermalProperties(smai.RegionType,smai.RegionID,'ThermalConductivity',smai.ThermalConductivity,'MassDensity',smai.MassDensity,'SpecificHeat',smai.SpecificHeat);
end

R = pde.ReducedThermalModel(self,Kr,Mr,Fr,u0Modal,modalTransformationMat,tempModelForReconstruction,Rm.ModeType,coefstruct);

end

