function R = reduce(self,varargin)
%reduce - Reduces the structural model using component mode synthesis
% R = reduce(modalmodel,'FrequencyRange',[omega1,omega2]) reduces a
% structural analysis model to the fixed interface modes in the frequency
% range [omega1,omega2] and the boundary interface degrees-of-freedom.
%
%   Example 1: Reduce a beam model to superelement, retaining only end faces
%              degrees of freedom.
%
%       structuralmodel = createpde('structural','transient-solid');    
%       structuralmodel.Geometry = multicuboid(0.1,0.01,0.01);
%       structuralProperties(structuralmodel,'YoungsModulus',70E9, ...
%                                      'PoissonsRatio',0.3, ...
%                                      'MassDensity',2700);
%       structuralSEInterface(structuralmodel,'Edge',[4,6,9,10]);
%       structuralSEInterface(structuralmodel,'Edge',[2,8,11,12]);
%       generateMesh(structuralmodel);
%       R = reduce(structuralmodel,'FrequencyRange',[-Inf,500000])
%
%    Example 2: Define multipoint constraint interfaces and generate a
%               reduced-order model with six degrees of freedom at each 
%               interface.
%
%       structuralmodel = createpde('structural','transient-solid');
%       structuralmodel.Geometry = multicuboid(0.1,0.01,0.01);
%       structuralProperties(structuralmodel,'YoungsModulus',70E9, ...
%                                      'PoissonsRatio',0.3, ...
%                                      'MassDensity',2700);
%       structuralBC(structuralmodel,'Edge',[4,6,9,10],'Constraint','multipoint');
%       structuralBC(structuralmodel,'Edge',[2,8,11,12],'Constraint','multipoint');
%       generateMesh(structuralmodel);
%       R = reduce(structuralmodel,'FrequencyRange',[-Inf,500000])
%
%    See also pde.StructuralModel, pde.StructuralModel/structuralProperties
%             pde.StructuralModel/structuralSEInterface

% Copyright 2019-2023 The MathWorks, Inc.

narginchk(3,3);
nargoutchk(0,1);

parser = inputParser;
parser.addParameter('FrequencyRange',@pde.EquationModel.isValidAttributeVal);
parser.parse(varargin{:});

validateattributes(parser.Results.FrequencyRange,{'double'},...
    {'real', 'nonsparse','vector','numel',2,'nondecreasing', 'nonnan'},mfilename,'FrequencyRange');

%Convert the input frequency range to eigenvalue range.
frequencyRangeSign = sign(parser.Results.FrequencyRange);
evalinterval = frequencyRangeSign.*parser.Results.FrequencyRange.^2;

performSolverPrecheck(self, false);

if ~isprop(self,'PackedBCs') 
    p = self.addprop('PackedBCs');
    p.Transient = true;
    p.Hidden = true;    
end
clearN = onCleanup(@() solverCleanup(p));

mtlstruct = self.MaterialProperties.packProperties();

if any(cellfun(@isempty,mtlstruct.Coefficients.m))
    error(message('pde:structuralROM:massRequiredForROM'));
end

coefstruct.ElementsInSubdomain = mtlstruct.ElementsInSubdomain;
coefstruct.NumElementsInSubdomain = mtlstruct.NumElementsInSubdomain;
coefstruct.Coefficients.m = mtlstruct.Coefficients.m;
coefstruct.Coefficients.d = mtlstruct.Coefficients.d;
coefstruct.Coefficients.c = mtlstruct.Coefficients.c;
coefstruct.Coefficients.a = mtlstruct.Coefficients.a;
coefstruct.IsComplex = [false;false;false;false;false];

%Assign zeros f-coefs for all subdomains to ignore loading.
for sd = 1:length(coefstruct.Coefficients.c)
    coefstruct.Coefficients.f{sd} = zeros(self.PDESystemSize,1);
end

% Compute BC matrices, this would account for zero displacements.
[p,e,t] = self.Mesh.meshToPet();
u0 = pdeuxpd(p,0,self.PDESystemSize);
tdummy = 0:0.1:1;

femodel = pde.DynamicDiscretizedPDEModel(self,p,e,t,coefstruct,u0,tdummy,true);
Ko = femodel.K;
Mo = femodel.Mass;

Ko = (Ko+Ko')/2;
Mo = (Mo+Mo')/2;

if isempty(self.SuperelementInterfaces)
    error(message('pde:structuralROM:noSEInterfaces'));
end
self.SuperelementInterfaces.performSolverPrechecks();

seiCleanObj = onCleanup(@() cleanupAutoSEIs(self.SuperelementInterfaces)); % Deletes auto SEs due to MPC BCs/
[AttachmentDoF,Tp, MPCReference] = self.SuperelementInterfaces.packSEI();

[Kr, Mr, W, Tr, Np] = craigBamptonImpl(self,Ko,Mo,AttachmentDoF,Tp,evalinterval);

elasticityType = split(self.AnalysisType,'-');
transientAnalysisType = join(["transient",elasticityType(2)],"-");
tempModelForReconstruction = createpde('structural',transientAnalysisType);
tempModelForReconstruction.Geometry = self.Geometry;
tempModelForReconstruction.Mesh = self.Mesh;
numSMassignments = numel(self.MaterialProperties.MaterialAssignments);
for i = 1:numSMassignments
    smai = self.MaterialProperties.MaterialAssignments(i);
    tempModelForReconstruction.structuralProperties(smai.RegionType,smai.RegionID,'YoungsModulus',smai.YoungsModulus,'PoissonsRatio',smai.PoissonsRatio,'MassDensity',smai.MassDensity,'CTE', smai.CTE);
end

R = pde.ReducedStructuralModel(self,Kr,Mr,AttachmentDoF,W,Np,MPCReference',Tr,tempModelForReconstruction,coefstruct);
end

function solverCleanup(prop)
% Clear the packed BCs for the current invocation of solver.
delete(prop);
end

function cleanupAutoSEIs(seiRecords)
% Clear the auto SEIs that were created at solve time.
for i = numel(seiRecords.StructuralSEIAssignments):-1:1
    if (seiRecords.StructuralSEIAssignments(i).MPCRequired)
        delete(seiRecords.StructuralSEIAssignments(i))
    end
end

numsei = numel(seiRecords.StructuralSEIAssignments);
if numsei == 0 && isvalid(seiRecords.ParentModel)
    seiRecords.ParentModel.delistSEAssignment();
end

end

