function [u,dudt,dudt2] = solveTimeDependent(self, coefstruct, u0, ut0, tlist, tsecondOrder)
% solveTimeDependent - Internal function for solving a time-dependent PDE
%

% Copyright 2015-2022 The MathWorks, Inc.

global femodel


p = self.addprop('PackedBCs');
p.Transient = true;
p.Hidden = true;

dbc = self.addprop('DiscretizedBCs');
dbc.Transient = true;
dbc.Hidden = true;


clearN = onCleanup(@() temporalSolverCleanup(p,dbc));
import pde.internal.odefunchandles.*

rtol =  self.SolverOptions.RelativeTolerance;
atol =  self.SolverOptions.AbsoluteTolerance;
stats = self.SolverOptions.ReportStatistics;
[p,e,t] = self.Mesh.meshToPet();

femodel=pde.DynamicDiscretizedPDEModel(self,p,e,t,coefstruct,u0,tlist,tsecondOrder);

if (isa(self,'pde.StructuralModel') || isa(self,'pde.internal.FESolverUtilities')) && (isprop(self,'DampingModels') && ~isempty(self.DampingModels))
    if isa(self.DampingModels,'pde.StructuralDampingAssignmentRecords') && ...
            contains(self.DampingModels.StructuralDampingAssignments(end).DampingModel,"modal")
        warning(message('pde:dampingModelSpecification:noModalDampingInDirectSolver'));
    else
        femodel.C = evaluateDampingMatrix(self,femodel,coefstruct);
    end
end

if(~femodel.IsSpatialCoefficientsNonlinear)
    % Spatial coefficients are independent of solution and time.
    % Impose BC on the system matrices only once in the beginning.
    nun=femodel.numConstrainedEqns;
    [null,orth]=pdenullorth(femodel.H);
    if size(orth,2)==0
        ud=zeros(size(femodel.K,2),1);
    else
        ud=full(orth*((femodel.H*orth)\femodel.R));
    end
    KK=femodel.K+femodel.A+femodel.Q;
    FF=null'*((femodel.F+femodel.G)-KK*ud);
    KK=null'*KK*null;
    
    KK = self.checkForSymmetry(KK);
    if femodel.IsSecondOrderODE
        femodel.K=[sparse(nun,nun) -speye(nun,nun); KK sparse(nun,nun)];
        femodel.F=[zeros(nun,1);FF];
    else
        femodel.K=KK;
        femodel.F=FF;
    end
    
end

femodel.K = self.checkForSymmetry(femodel.K);
femodel.Mass = self.checkForSymmetry(femodel.Mass);

massNonSingular = 'maybe';
if (isa(self,'pde.internal.FESolverUtilities') && isprop(self,'AnalysisType')) ||...
        (~isempty(self.ParentModel) && isprop(self.ParentModel,'AnalysisType'))
    massNonSingular = 'no';
end

unConstrainedDoF = size(p,2)*self.PDESystemSize;
totalDoF = unConstrainedDoF;

if tsecondOrder
    [u0,ut0] = imposeBConICsecondOrderODE(self,u0,ut0,tlist);
    unConstrainedDoF=size(u0,1);
    uu0=[u0;ut0];
    if isa(self,'pde.PDEModel') && (self.EquationCoefficients.dDefined())
        femodel.C = coefstruct.Coefficients.d{1};
    end
    %Following pragma, %#function, are needed for deployment.
    %#function pde.internal.odefunchandles.secondOrderODEf
    fcn  = str2func('secondOrderODEf');
    %#function pde.internal.odefunchandles.secondOrderODEdf
    dfcn = str2func('secondOrderODEdf');
    %#function pde.internal.odefunchandles.secondOrderODEm
    mfcn = str2func('secondOrderODEm');
else
    uu0 = imposeBConICfirstOrderODE(self,u0,tlist);
    unConstrainedDoF=size(uu0,1);
    %Following pragma, %#function, are needed for deployment.
    %#function pde.internal.odefunchandles.firstOrderODEf
    fcn  = str2func('firstOrderODEf');
    %#function pde.internal.odefunchandles.firstOrderODEdf
    dfcn = str2func('firstOrderODEdf');
    %#function pde.internal.odefunchandles.firstOrderODEm
    mfcn = str2func('firstOrderODEm');
end
odeoptions = constructODEoptions(uu0,rtol,atol,stats,unConstrainedDoF,dfcn,mfcn,tsecondOrder,massNonSingular);

sol=ode15s(fcn,tlist,uu0,odeoptions);
if(sol.x(end) < tlist(end))
    error(message('pde:pdetool:incompSoln'))
end
[y,yp] = deval(sol,tlist);
numCompTimes = size(y, 2);

DoFIndex=1:unConstrainedDoF; % vector of 1 to number of degree of freedom
u  = y(DoFIndex,:);
dudt = zeros(size(u));
dudt2 = zeros(size(u));
    
if tsecondOrder
    velIndex=unConstrainedDoF+1:size(y,1);
    dudt(DoFIndex,:) = y(velIndex,:);
    dudt2(DoFIndex,:) = yp(velIndex,:);
end

% If Dirichlet constraints are functions of time or u, recovery of the
% full solution requires reevaluating H and R matrices
if (femodel.vh || femodel.vr)
    % new constrained values calc
    uFull=zeros(totalDoF,numCompTimes);
    dudtFull=zeros(totalDoF,numCompTimes);
    dudt2Full=zeros(totalDoF,numCompTimes);
    
    if (~femodel.hu || femodel.ru)
        for k=1:numCompTimes
            [uFull(:,k), dudtFull(:,k), dudt2Full(:,k)] = recoverFullSolutionLinear(self,tsecondOrder,u(:,k),dudt(:,k),dudt2(:,k),tlist(k),femodel.tspan,femodel.H);
        end
    else
        for k=1:numCompTimes
            [uFull(:,k), dudtFull(:,k), dudt2Full(:,k)] = recoverFullSolution(self,tsecondOrder,u(:,k),dudt(:,k),dudt2(:,k),tlist(k),femodel.tspan,femodel.H,femodel.R, ...
                rtol, atol(1));
        end
    end
else
    uFull=femodel.B*u+femodel.ud*ones(1,numCompTimes);
    if tsecondOrder
        dudtFull=femodel.B*dudt+femodel.dudt*ones(1,numCompTimes);
        dudt2Full=femodel.B*dudt2+femodel.dudt*zeros(1,numCompTimes); % No support for acceleration BC.
    end
end

u = uFull;
if tsecondOrder
    dudt = dudtFull;
    dudt2 = dudt2Full;
end

end


function temporalSolverCleanup(p,dbc)
% Clear the packed BCs for the current invocation of solver.
delete(p);
delete(dbc);
%Clear global data.
clear('global', 'femodel')
end


function u0 = imposeBConICfirstOrderODE(thePde, u0,tlist)
global femodel
if ~(femodel.vh || femodel.vr)
    u0=femodel.B'*u0;
else
    [~,~,H,~]=thePde.assembleBoundary(u0,tlist(1));
    B=pdenullorth(H);
    u0=B'*u0;
end

end

function [u0, ut0] = imposeBConICsecondOrderODE(thePde, u0,ut0,tlist)
global femodel
if ~(femodel.vh || femodel.vr)
    u0=femodel.B'*u0;
    ut0=femodel.B'*ut0;
else
    t0=tlist(1);
    t1=t0+femodel.tspan*sqrt(eps);
    dt=t1-t0;
    [~,~,H,R]=thePde.assembleBoundary(u0,t0);
    [N,O]=pdenullorth(H);
    ud=O*((H*O)\R);
    [~,~,H,R]=thePde.assembleBoundary(u0,t1);
    [N1,O]=pdenullorth(H);
    ud1=O*((H*O)\R);
    u0=N'*u0;
    ut0=N'*(ut0-((N1-N)*u0-(ud1-ud))/dt);
end
end


function odeoptions = constructODEoptions(uu0,rtol,atol,stats,unConstrainedDoF,dfcn,mfcn,tsecondOrder,massNonSingular)
global femodel

odeoptions = odeset;
odeoptions=odeset(odeoptions, 'Jacobian',dfcn);
if(~(femodel.vc || femodel.va || femodel.vq || femodel.vh || femodel.vr))
    odeoptions=odeset(odeoptions,'JConstant','on');
end

if(~(femodel.vd || femodel.vh))
    odeoptions=odeset(odeoptions,'Mass', mfcn(0, uu0));
else
    odeoptions=odeset(odeoptions,'Mass',mfcn);
end
odeoptions=odeset(odeoptions,'RelTol',rtol);

if (tsecondOrder)
    % Mask out dudt
    atol=atol*ones(2*unConstrainedDoF,1);
    atol(unConstrainedDoF+1:unConstrainedDoF+unConstrainedDoF)=Inf*ones(unConstrainedDoF,1);
    odeoptions=odeset(odeoptions,'MaxOrder',2);
end
odeoptions=odeset(odeoptions,'MassSingular',massNonSingular);
odeoptions=odeset(odeoptions,'AbsTol',atol);
odeoptions=odeset(odeoptions,'Stats',stats);

end


function [uf, dudtf, dudt2f] = recoverFullSolutionLinear(thePde,tsecondOrder,ur,durdt,dur2dt2,t,tspan,H)
%This function is used to transform the reduced solution with constrained DOFs
% removed to the full solution with all DOFs at all FE nodes for linear
% system.

dt=tspan*eps^(1/3);
t1 = t+dt;
tm1 = t-dt;

ufi = zeros(size(H,2),1);
[~,~,H,R]=thePde.assembleBoundary(ufi,t);

[null,orth]=pdenullorth(H);
ud=full(orth*((H*orth)\R));
uf=null*ur + ud;

dudtf = zeros(size(uf));
dudt2f = zeros(size(uf));

if tsecondOrder
% Compute time-derivative using a central difference approximation.
[~,~,H,R]=thePde.assembleBoundary(uf,t1);
ud1=full(orth*((H*orth)\R));
[~,~,H,R]=thePde.assembleBoundary(uf,tm1);
udm1=full(orth*((H*orth)\R));
udDot = (ud1 - udm1)./(2*dt);
udDDot = (ud1 -2*ud+ udm1)./dt.^2;
dudtf = null*durdt + udDot;
dudt2f = null*dur2dt2 + udDDot;
end
end



function [uf, dudtf, dudt2f, isConverged] = recoverFullSolution(thePde,tsecondOrder,ur,durdt,dur2dt2,t,tspan,H,R, rtol, atol)
%This function is used to transform the reduced solution with constrained DOFs
% removed to the full solution with all DOFs at all FE nodes.
%
% When the H or R matrices are functions of the dependent variables
% (i.e. a nonlinear system), multiple iterations may be required to
% obtain a converged solution.

maxIter = 50;
ufi = zeros(size(H,2),1);
dt=tspan*eps^(1/3);
t1 = t+dt;
tm1 = t-dt;

for i=1:maxIter
    [null,orth]=pdenullorth(H);
    ud=full(orth*((H*orth)\R));
    uf=null*ur + ud;
    
    dudtf = zeros(size(uf));
    dudt2f = zeros(size(uf));
    
    if tsecondOrder
        % Compute time-derivative using a central difference approximation.
        [~,~,H,R]=thePde.assembleBoundary(uf,t1);
        ud1=full(orth*((H*orth)\R));
        [~,~,H,R]=thePde.assembleBoundary(uf,tm1);
        udm1=full(orth*((H*orth)\R));
        udDot = (ud1 - udm1)./(2*dt);
        udDDot = (ud1 -2*ud+ udm1)./dt.^2;
        dudtf = null*durdt + udDot;
        dudt2f = null*dur2dt2 + udDDot;
    end
    
    err=max(abs(ufi-uf)); % Check convergence only on displacements
    isConverged = all(abs(ufi-uf) < max(rtol*abs(uf), atol));
    if(isConverged)
        break;
    end
    ufi=uf;
    [~,~,H,R]=thePde.assembleBoundary(uf,t);
end
end