function Rtransient = reconstructSolution(self,u,ut,utt,tlist)
%reconstructSolution - Reconstruct the full-model transient solution
% R = reconstructSolution(rom,u,dudt,d2udt2,tlist) expands the transient
% solution obtained using the reduced-order model, rom, to the full
% solution and returns a pde.TransientStructuralResults object. The inputs
% arguments u, dudt, and d2udt2 are the displacement, velocity, and
% acceleration matrices, respectively. The argument tlist is a vector of
% the time-steps that you used to solve the reduced model. The number of
% rows in the input matrices must match the sum of the numbers of the
% interface degrees of freedom and the number of modes. The number of
% columns in the input matrices must match the number of elements in tlist.
%
%    See also pde.StructuralModel, pde.StructuralModel/structuralProperties
%      pde.StructuralModel/structuralSEInterface,pde.StructuralModel/reduce

% Copyright 2019-2023 The MathWorks, Inc.

narginchk(5,5);

validateattributes(tlist,{'double'},{'real','nonsparse','vector','nonnan','finite','nondecreasing'});
validateattributes(u,{'double'},{'real','nonsparse','2d','nonnan','finite'});
validateattributes(ut,{'double'},{'real','nonsparse','2d','nonnan','finite'});
validateattributes(utt,{'double'},{'real','nonsparse','2d','nonnan','finite'});

    
W = self.TransformationMatrix;
numNodes = size(self.Mesh.Nodes,2);
dim = size(self.Mesh.Nodes,1);
DoFVector = 1:dim*numNodes;

numROMDoF = size(self.K,1);

inputROMResultsSize = [size(u,1), size(ut,1), size(utt,1)];

if any(inputROMResultsSize ~= numROMDoF)
    error(message('pde:structuralROM:incorrectSizeROMResults'));
end

uFull = zeros(numel(DoFVector),numel(tlist));
utFull = zeros(numel(DoFVector),numel(tlist));
uttFull = zeros(numel(DoFVector),numel(tlist));

if isempty(self.ReferenceLocations)
    % No MPC model. Results from user's own workflow.
    uExpanded = W*u;
    utExpanded = W*ut;
    uttExpanded = W*utt;
else
    % MPC model. Results from Simscape simulation.
    T = self.MPCTransformationMatrix;
    W = self.TransformationMatrix;
    uExpanded = T*W*u;
    utExpanded = T*W*ut;
    uttExpanded = T*W*utt;
end

CondensedDoF = setdiff(DoFVector,self.RetainedDoF,'stable')';
numRetainedDoF = numel(self.RetainedDoF);

uFull(self.RetainedDoF,:)   = uExpanded(1:numRetainedDoF,:);
utFull(self.RetainedDoF,:)  = utExpanded(1:numRetainedDoF,:);
uttFull(self.RetainedDoF,:) = uttExpanded(1:numRetainedDoF,:);

uFull(CondensedDoF,:)       = uExpanded(numRetainedDoF+1:end,:);
utFull(CondensedDoF,:)      = utExpanded(numRetainedDoF+1:end,:);
uttFull(CondensedDoF,:)     = uttExpanded(numRetainedDoF+1:end,:);

Rtransient = pde.TransientStructuralResults(self.ModelForReconstruction,uFull,utFull,uttFull,tlist,self.Coefficients);
end
