function outStruct = reduceSimscape(self,varargin)
% reduceSimscape - Internal function to compute structural ROM for
% Simscape flexible body workflow.

% Copyright 2022 The MathWorks, Inc.

p = varargin{1};
t = varargin{2};
coefstruct = varargin{3};
mpcDefinition = varargin{4};
numModes = varargin{5};
displayVertexIDsGlobal = varargin{6};
numMPCs = size(mpcDefinition.Reference,2);

% Simscape calls reduceSimscape() each time with a new temp model.
self.SuperelementInterfaces = pde.StructuralSEIAssignmentRecords(self);
self.SuperelementInterfaces.StructuralSEIAssignments = pde.StructuralSEIAssignment.empty;
for i = 1:numMPCs
    self.SuperelementInterfaces.StructuralSEIAssignments(end+1) = pde.StructuralSEIAssignment( ...
               self.SuperelementInterfaces,'MPCRequired',true,...
               'Reference',mpcDefinition.Reference(:,i),'MPCNodes',mpcDefinition.NodeIDs{i});
end

femodel = self.assembleSelectedFEMatrices(p, t, coefstruct, [], [], 'MK', false);
Ko = femodel.K;
Mo = femodel.M;
Ko = (Ko+Ko')/2;
Mo = (Mo+Mo')/2;

[AttachmentDoF,Tp, ~] = self.SuperelementInterfaces.packSEI();
[Kr, Mr, W, Tr] = self.craigBamptonImpl(Ko,Mo,AttachmentDoF,Tp,numModes);

% Structural matrix data
StructuralMatrices.NumberOfDofs = size(Kr,1);
StructuralMatrices.StiffnessMatrix = Kr(:);
StructuralMatrices.MassMatrix = Mr(:);

% Recovery data
% Get the mapping of DoFs corresponding displayVertexIDsGlobal to the ROM
% reconstruction matrix rows. ROM ReconMat has attachment DoFs at the
% beginning and then the condensed DoFs.
numNodes = size(self.Mesh.Nodes,2);
DoFVector = 1:self.Dim*numNodes;
CondensedDoF = setdiff(DoFVector,AttachmentDoF,'stable')';
romDoF = [AttachmentDoF.'; CondensedDoF];
[~, globalToROMMap] = ismember(DoFVector,romDoF);
displayVertexIDsROMMatched = [globalToROMMap(displayVertexIDsGlobal); globalToROMMap(displayVertexIDsGlobal+numNodes);  globalToROMMap(displayVertexIDsGlobal+2*numNodes);];

ReconMat = Tr*W;
ReconVec = zeros(self.Dim*size(ReconMat,2),numel(displayVertexIDsGlobal));
for i = 1:numel(displayVertexIDsGlobal)
   rv = ReconMat(displayVertexIDsROMMatched(:,i),:); % Gets three_row x NumROM_DOF matrix, rows correspond to x,y, and z directions.
   ReconVec(:,i) = rv(:);
end
ReconVec = ReconVec(:);

% Pack output
outStruct.StructuralMatrices = StructuralMatrices;
outStruct.Visualization = ReconVec;
end