function u1=hyperbolic(varargin) % u0,ut0,tlist,b,p,e,t,c,a,f,d,varargin
%HYPERBOLIC Solve hyperbolic PDE problem.
%
%       U1=HYPERBOLIC(U0,UT0,TLIST,B,P,E,T,C,A,F,D) produces the
%       solution to the FEM formulation of the scalar PDE problem
%       d*d^2u/dt^2-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 and initial derivative UT0.
%
%       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=HYPERBOLIC(U0,UT0,TLIST,B,P,E,T,C,A,F,D,RTOL) and
%       U1=HYPERBOLIC(U0,UT0,TLIST,B,P,E,T,C,A,F,D,RTOL,ATOL) passes
%       absolute and relative tolerances to the ODE solver.
%
%       U1=HYPERBOLIC(U0,UT0,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=HYPERBOLIC(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=HYPERBOLIC(U0,UT0,TLIST,K,F,B,UD,M) or
%       U1=HYPERBOLIC(U0,UT0,TLIST,K,F,B,UD,M,RTOL) or
%       U1=HYPERBOLIC(U0,UT0,TLIST,K,F,B,UD,M,RTOL,ATOL) 
%       produces the solution to the following ODE
%       B'*M*B*(d^2ui/dt^2) + K*ui = F, 
%       u = B*ui + ud, with initial values U0 and UT0.
%
%       In addition to the optional property-value pairs listed above for 
%       the coefficient-form, this matrix-form of HYPERBOLIC can
%       take the following optional, property-value pair.
%
%       Property        Value/{Default}         Description
%       -----------------------------------------------------------
%       DampingMatrix     D               Sparse damping matrix with
%                                         the same size as the M-matrix.
%        
%       If the DampingMatrix option is included, the following ODE
%       problem is solved B'*M*B*(d^2ui/dt^2) + B'*D*B*(dui/dt) + K*ui = F, 
%       u = B*ui + ud, with initial values U0 and UT0.
%
%       Note: This function does not support quadratic triangular elements
%
%       See also CREATEPDE, pde.PDEMODEL, ASSEMPDE, PARABOLIC

%       Copyright 1994-2020 The MathWorks, Inc.

    if nargin > 0
        [varargin{:}] = convertStringsToChars(varargin{:});
    end
    
    narginchk(8,17);
    global pdehyp
    clearN = onCleanup(@() clear('global', 'pdehyp'));

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

    parser = inputParser;
    parser.addRequired('u0', @pdeIsValid.t0);
    parser.addRequired('ut0', @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{5};
        e = varargin{6};
        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;
      ut0 = pres.ut0;
      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);
      pdehyp = pdeHyperbolicInfo(u0,ut0,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.addParamValue('DampingMatrix', [], @issparse);
      parser.parse(varargin{:});
      pres = parser.Results;  
      u0 = pres.u0;
      ut0 = pres.ut0;
      tlist = pres.tlist; 
      b = pres.K;
      p = pres.F;
      e = pres.B;
      t = pres.ud;
      c = pres.M;
      quadraticTriNotSupported(t);
      pdehyp = pdeHyperbolicInfo(u0,ut0,tlist,b,p,e,t,c);
      pdehyp.C = pres.DampingMatrix;
      % Make sure damping matrix and mass matrix (c argument) have same size
      if(~isempty(pdehyp.C) && ~all(size(c) == size(pdehyp.C)))
          error(message('pde:hyperbolic:dampsize'));
      end
    end
  
 
  rtol = parser.Results.rtol;
  atol = parser.Results.atol;
  stats = parser.Results.Stats;
    
  % Expand initial condition
  if ischar(u0) || ischar(ut0)
      x=pdehyp.p(1,:)';
      y=pdehyp.p(2,:)';
      
      if(size(pdehyp.p,1)==2)
          % 2D case
          args = {x,y};
      elseif(size(pdehyp.p,1)==3)
          % 3D case
          z = pdehyp.p(3,:)';
          args = {x,y,z};
      end
  end
  
  if ischar(u0)
      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(pdehyp.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
  if ischar(ut0)
      if pdeisfunc(ut0)
          ut0=feval(ut0,args{:}); % supposedly a function name
      else
          if(size(u0,1)==1) % numeric/string from GUI or string from CLI for scalar
              ut0Temp = eval(ut0);
              if(~isnumeric(ut0Temp)) % must be string from GUI,only works for scalar
                  ut0Temp = eval(ut0Temp);
              end
          else                        % CLI systems case only, PDE tool will not work
              ut0Temp = zeros(size(pdehyp.p,2),size(ut0,1));
              for i=1:size(ut0,1)
                  ut0Temp(:,i) = eval(ut0(i,:)); % supposedly an expression of 'x','y' and 'z'
              end
          end
          ut0 = ut0Temp;
      end
  end
  
  u0=u0(:);
  ut0=ut0(:);
  
  if(isCoeffForm)
    nu=size(p,2)*pdehyp.N;
  else
    nu=size(pdehyp.ud,1);
  end
  
  nuu=nu;

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

if ~(pdehyp.vh || pdehyp.vr)
  u0=pdehyp.B'*u0;
  ut0=pdehyp.B'*ut0;
else
  t0=tlist(1);
  t1=t0+pdehyp.ts*sqrt(eps);
  dt=t1-t0;
  [Q,G,H,R]=assemb(b,p,e,u0,t0);
  [N,O]=pdenullorth(H);
  ud=O*((H*O)\R);
  [Q,G,H,R]=assemb(b,p,e,u0,t1);
  [N1,O]=pdenullorth(H);
  ud1=O*((H*O)\R);
  u0=N'*u0;
  ut0=N'*(ut0-((N1-N)*u0-(ud1-ud))/dt);
end

nu=size(u0,1);

uu0=[u0;ut0];

options=odeset('MaxOrder',2);
% If doJacFiniteDiff is true, ode15s computes the Jacobian by finite
% difference. This is useful for testing purposes.
doJacFiniteDiff=false;
if(doJacFiniteDiff)
  jac0=pdehypdf(0, uu0);
  [ir,ic]=find(jac0);
  jpat=sparse(ir,ic,ones(size(ir,1),1), size(jac0,1), size(jac0,2));
  options=odeset(options, 'JPattern', jpat);
else
  options=odeset(options, 'Jacobian',@pdehypdf);
end
if(~(pdehyp.vc || pdehyp.va || pdehyp.vq || pdehyp.vh || pdehyp.vr))
  options=odeset(options,'JConstant','on');
end

if(~(pdehyp.vd || pdehyp.vh))
  options=odeset(options,'Mass', pdehypm(0, uu0));
else
  options=odeset(options,'Mass',@pdehypm);
end
options=odeset(options,'RelTol',rtol);

% Mask out ut
atol=atol*ones(2*nu,1);
atol(nu+1:nu+nu)=Inf*ones(nu,1);
options=odeset(options,'AbsTol',atol);
options=odeset(options,'Stats',stats);

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

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

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

end

