function u1=parabolic(varargin) % u0,tlist,b,p,e,t,c,a,f,d,varargin
%PARABOLIC Solve parabolic PDE problem
%
%       U1=PARABOLIC(U0,TLIST,B,P,E,T,C,A,F,D) produces the
%       solution to the FEM formulation of the scalar PDE problem
%       d*du/dt-div(c*grad(u))+a*u=f, on a mesh described by P, E,
%       and T, with boundary conditions given by B, and with initial
%       value U0.
%
%       For the scalar case, each row in the solution matrix U1 is the
%       solution at the coordinates given by the corresponding column in
%       P. Each column in U1 is the solution at the time given by the
%       corresponding item in TLIST.  For a system of dimension N
%       with NP node points, the first NP rows of U1 describe the first
%       component of u, the following NP rows of U1 describe the second
%       component of u, and so on.  Thus, the components of u are placed
%       in blocks U as N blocks of node point rows.
%
%       B describes the boundary conditions of the PDE problem.  B
%       can either be a Boundary Condition Matrix or the name of Boundary
%       MATLAB-file. See PDEBOUND for details.
%
%       The coefficients C, A, F, and D of the PDE problem can
%       be given in a wide variety of ways.  See ASSEMPDE for details.
%
%       U1=PARABOLIC(U0,TLIST,B,P,E,T,C,A,F,D,RTOL) and
%       U1=PARABOLIC(U0,TLIST,B,P,E,T,C,A,F,D,RTOL,ATOL) passes
%       absolute and relative tolerances to the ODE solver.
%     
%       U1=PARABOLIC(U0,TLIST,PDEM,...) Allows the boundary conditions B and 
%       the mesh (P,E,T) to be defined by PDEM, a PDEModel object that 
%       contains this data.
%
%       U1=PARABOLIC(U0, ...) can take optional property-value pairs as
%       the last arguments to the function. It is not necessary to include
%       RTOL and ATOL in order to pass optional property-value pairs.
%       Valid property-value pairs include:
% 
%       Property        Value/{Default}         Description
%       -----------------------------------------------------------
%       Stats            {on}|off               Control statistics display 
%                                               from the ode solver.
%
%       U1=PARABOLIC(U0,TLIST,K,F,B,UD,M) produces the solution to
%       the ODE problem B'*M*B*(dui/dt)+K*ui=F, u=B*ui+ud,
%       with initial value U0.
%
%       Note: This function does not support quadratic triangular elements
%
%       See also CREATEPDE, pde.PDEMODEL, ASSEMPDE, HYPERBOLIC

%       Copyright 1994-2020 The MathWorks, Inc.

if nargin > 0
    [varargin{:}] = convertStringsToChars(varargin{:});
end

narginchk(7,14);

global pdeprb
clearN = onCleanup(@() clear('global', 'pdeprb'));

b = varargin{3};
meshedpde=false;  
if isa(b, 'pde.PDEModel') 
   if isa(b.Mesh, 'pde.FEMesh')
      mypde = varargin{3}; 
      meshedpde=true;  
      msh = mypde.Mesh;
      [p,e,t] = msh.meshToPet();    
      b = pde.PDEModel(mypde.PDESystemSize);  
      b.BoundaryConditions = mypde.BoundaryConditions;
      varargin = {varargin{1:2}, varargin{4:end}};
      quadraticTriNotSupported(t);
    elseif ~b.IsTwoD
      error(message('pde:parabolic:pdemodelNoMesh'));  
   end  
end


parser = inputParser;
parser.addRequired('u0', @pdeIsValid.t0);
parser.addRequired('tlist', @pdeIsValid.tlist);
isCoeffForm = ~issparse(b);
if(isCoeffForm)    
  if ~meshedpde
    parser.addRequired('b', @pdeIsValid.b);
    parser.addRequired('p', @pdeIsValid.p2d);
    p = varargin{4};
    e = varargin{5};
    eTest = @(e) pdeIsValid.meshGeomAssoc(e, p);
    parser.addRequired('e', eTest);
    parser.addRequired('t', @pdeIsValid.t2d);  
  end
  parser.addRequired('c', @pdeIsValid.coefficient);
  parser.addRequired('a', @pdeIsValid.coefficient);
  parser.addRequired('f', @pdeIsValid.coefficient);
  parser.addRequired('d', @pdeIsValid.coefficient);
  parser.addOptional('rtol', 1e-3, @isnumeric);
  parser.addOptional('atol', 1e-6, @isnumeric);
  parser.addParamValue('Stats', 'on', @pdeIsValid.onOrOff);
  parser.parse(varargin{:});
  pres = parser.Results;  
  u0 = pres.u0;
  tlist = pres.tlist; 
  c = pres.c;
  a = pres.a;
  f = pres.f;
  d = pres.d;
  if ~meshedpde  
      b = pres.b;
      p = pres.p;
      e = pres.e;
      t = pres.t;
  end
  quadraticTriNotSupported(t);
  pdeprb=pdeParabolicInfo(u0,tlist,b,p,e,t,c,a,f,d);
else
  parser.addRequired('K', @issparse);
  parser.addRequired('F', @isnumeric);
  parser.addRequired('B', @issparse);
  parser.addRequired('ud', @isnumeric);
  parser.addRequired('M', @issparse);
  parser.addOptional('rtol', 1e-3, @isnumeric);
  parser.addOptional('atol', 1e-6, @isnumeric);
  parser.addParamValue('Stats', 'on', @pdeIsValid.onOrOff);  
  parser.parse(varargin{:});
  pres = parser.Results;  
  u0 = pres.u0;
  tlist = pres.tlist; 
  b = pres.K;
  p = pres.F;
  e = pres.B;
  t = pres.ud;
  c = pres.M;
  quadraticTriNotSupported(t);
  pdeprb=pdeParabolicInfo(u0,tlist,b,p,e,t,c);
end


rtol = parser.Results.rtol;
atol = parser.Results.atol;
stats = parser.Results.Stats;

if ischar(u0)
    x=pdeprb.p(1,:)';
    y=pdeprb.p(2,:)';
    if(size(pdeprb.p,1)==2)
        % 2D case
        args = {x,y};
    elseif(size(pdeprb.p,1)==3)
        % 3D case
        z = pdeprb.p(3,:)';
        args = {x,y,z};
    end
    if pdeisfunc(u0)
        u0=feval(u0,args{:}); % supposedly a function name
    else
        if(size(u0,1)==1) % numeric/string from GUI or string from CLI for scalar
            u0Temp = eval(u0);
            if(~isnumeric(u0Temp)) % must be string from GUI,only works for scalar
                u0Temp = eval(u0Temp);
            end
        else                       % CLI systems case only, PDE tool will not work
            u0Temp = zeros(size(pdeprb.p,2),size(u0,1));
            for i=1:size(u0,1)
                u0Temp(:,i) = eval(u0(i,:)); % supposedly an expression of 'x','y' and 'z'
            end
        end
        u0 = u0Temp;
    end
end
u0=u0(:);

if(isCoeffForm)
  nu=size(p,2)*pdeprb.N;
else
  nu=size(pdeprb.ud,1);
end

nuu=nu;
if size(u0,1)==1
  u0=ones(nu,1)*u0;
end

if ~pdeprb.vh
  u0=pdeprb.B'*u0;
else
  [Q,G,H,R]=assemb(b,p,e,u0,tlist(1));
  B=pdenullorth(H);
  u0=B'*u0;
end

% If doJacFiniteDiff is true, ode15s computes the Jacobian by finite
% difference. This is useful for testing purposes.
doJacFinitDiff=false;
if(doJacFinitDiff)
  jac0=pdeprbdf(0, u0);
  [ir,ic]=find(jac0);
  jpat=sparse(ir,ic,ones(size(ir,1),1), size(jac0,1), size(jac0,2));
  options=odeset('JPattern', jpat);
else
  options=odeset('Jacobian',@pdeprbdf);
end
if ~(pdeprb.vc || pdeprb.va || pdeprb.vq || pdeprb.vh || pdeprb.vr)
  options=odeset(options,'JConstant','on');
end
options=odeset(options,'Mass',@pdeprbm);
options=odeset(options,'RelTol',rtol);
options=odeset(options,'AbsTol',atol);
options=odeset(options,'Stats',stats);

[~,uu]=ode15s(@pdeprbf,tlist,u0,options);
numCompTimes = size(uu, 1);
if(numCompTimes < size(tlist, 2))
  warning(message('pde:pdetool:incompSoln'))
end

if(length(tlist)==2)
  % ODE15S behaves differently if length(tlist)==2
  u1=transpose(uu([1 size(uu,1)],:));
  numCompTimes = 2;
else
  u1=transpose(uu);
end

% If Dirichlet constraints are functions of time or u, recovery of the
% full solution requires reevaluating H and R matrices
if pdeprb.vh || pdeprb.vr
    uu1=zeros(nuu,numCompTimes);
    for k=1:numCompTimes
        uu1(:,k)=pdeCalcFullU(b,p,e,u1(:,k),tlist(k),pdeprb.H,pdeprb.R, ...
            rtol, atol(1));
    end
    u1=uu1;
else
  u1=pdeprb.B*u1+pdeprb.ud*ones(1,numCompTimes);
end


end
