function [Kr,Mr,Fr,u0Modal,Phi,femodel] = modalProjection(self,Rm,Np,coefstruct,u0,tlist)
% modalProjection - Internal function to compute modal projections.

% Copyright 2021-2022 The MathWorks, Inc.

[p,e,t] = self.Mesh.meshToPet();

% Discretize the model.
femodel=pde.DynamicDiscretizedPDEModel(self,p,e,t,coefstruct,u0,tlist,false);

%Verify that constraints remain unchanged between modal analysis and
%transient simulation.
if ~isequal(femodel.B,Rm.ModalConstraintsMap)
    error(message('pde:StructuralModel:inconsistentConstraintsModalAndTransient'));
end

isModelNonlinear = femodel.vc || femodel.va || femodel.vd ||femodel.vq;
if isModelNonlinear
    error(message('pde:ThermalModel:modelNotLinearForModalProjection'))
end

% Construct modal matrix.
Phi = Rm.ModeShapes(:,1:Np);
if Rm.ModeType == "PODModes"
    Phi = [Rm.SnapshotsAverage,Phi];
end
% Impose constraints on the modes.
Phi = femodel.B'*Phi;

% Construct projected system of matrices
Kr = Phi.'*femodel.B'*(femodel.K+femodel.Q)*femodel.B*Phi;
Mr = Phi.'*femodel.B'*femodel.Mass*femodel.B*Phi;
Fr = Phi.'*femodel.B'*(femodel.F +femodel.G- (femodel.K+femodel.Q)*femodel.ud);

Kr = (Kr+Kr.')/2;
Mr = (Mr+Mr.')/2;

% Project IC on to modal basis.
u0Modal = [];
if ~isempty(u0)
    if Rm.ModeType == "PODModes"
        u0Modal = [1; Phi(:,2:end).'*femodel.B'*femodel.Mass*(u0-Rm.SnapshotsAverage)];
    else
        u0Modal = Phi.'*femodel.B'*femodel.Mass*u0;
    end
end

end

