classdef pdeODEInfo
  %pdeODEInfo Internal class used to communicate data to ode15s helper functions.
  % This undocumented function may be removed in a future release.
  
  %       Copyright 2012-2015 The MathWorks, Inc.
 
  properties
    % input parameters
    p, e, t, b, c, a, f, d;
    % Current values of global matrices
    K, M, F, Q, G, H, R, MM, B, ud, Nu, Or;
    % Flags indicating whether coefficient is dependent on time or u
    vq, vg, vh, vr, vc, va, vf, vd;
    N, ts;
    numNodes;
    totalNumEqns, numConstrainedEqns;
    isMatrixMode;
    thePde;
  end
  
  properties(Access=private)
    % Flags indicating whether coefficients are functions of time
    qt, gt, ht, rt, ct, at, ft, dt;
    % Flags indicating whether coefficients are functions of solution
    qu, gu, hu, ru, cu, au, fu, du;
  end
  
  methods(Access=protected)
    function obj=pdeODEInfo(tlist, isMatMode)
      obj.qt=false; obj.gt=false; obj.ht=false; obj.rt=false; obj.ct=false; 
      obj.at=false; obj.ft=false; obj.dt=false;
      obj.qu=false; obj.gu=false; obj.hu=false; obj.ru=false; obj.cu=false; 
      obj.au=false; obj.fu=false; obj.du=false;
      obj.vq=false; obj.vg=false; obj.vh=false; obj.vr=false; obj.vc=false; 
      obj.va=false; obj.vf=false; obj.vd=false;
      obj.isMatrixMode = isMatMode;
      obj.checkTlist(tlist);
    end
  end
  
  methods(Access=protected)
    
    function self=checkFuncDepen(self, u0, tlist)
      %checkFuncDepen determine whether coefficients are functions of t or u
      
      self.numNodes = size(self.p,2);
      
      nt=size(tlist,1);
      t0 = min(tlist);
      if(t0 == 0)
        t0 = eps;
      end
      t1 = NaN;
      
      u0 = evalU0Expr(self, u0);
      self.thePde = getPde(self.p, self.e, self.t, self.b, self.f, u0, t0);
      self.N = self.thePde.PDESystemSize;
      if(size(u0,1) == 1)
        u0=pdeuxpd(self.p, u0, self.N);
      end
      [MM0,K0,M0,F0,Q0,G0,H0,R0] = self.getMats(u0, t0);
      [MM1,K1,M1,F1,Q1,G1,H1,R1] = self.getMats(u0, t1);
        
      u1 = NaN*ones(size(u0,1), 1);   
      [MM2,K2,M2,F2,Q2,G2,H2,R2] = self.getMats(u1, t0);
      
      % check if coefficients are functions of time    
      if(any(isnan(MM1(:))))
        self.dt=true;
      end
      if(any(isnan(K1(:))))
        self.ct=true;
      end
      if(any(isnan(M1(:))))
        self.at=true;
      end
      if(any(isnan(F1(:))))
        self.ft=true;
      end
      if(any(isnan(Q1(:))))
        self.qt=true;
      end
      if(any(isnan(G1(:))))
        self.gt=true;
      end
      if(any(isnan(H1(:))))
        self.ht=true;
      end
      if(any(isnan(R1(:))))
        self.rt=true;
      end

      % check if coefficients are functions of solution
      % In the current code, u-dependence and time-dependence are
      % handled the same. So check again only if no time-dependence. 
      if(~self.dt && any(isnan(MM2(:))))
        self.du=true;
      end
      if(~self.ct && any(isnan(K2(:))))
        self.cu=true;
      end
      if(~self.at && any(isnan(M2(:))))
        self.au=true;
      end
      if(~self.ft && any(isnan(F2(:))))
        self.fu=true;
      end
      if(~self.qt && any(isnan(Q2(:))))
        self.qu=true;
      end
      if(~self.gt && any(isnan(G2(:))))
        self.gu=true;
      end
      if(~self.ht && any(isnan(H2(:))))
        self.hu=true;
      end
      if(~self.rt && any(isnan(R2(:))))
        self.ru=true;
      end

      self.vq = self.qt || self.qu; 
      self.vg = self.gt || self.gu; 
      self.vh = self.ht || self.hu; 
      self.vr = self.rt || self.ru; 
      self.vc = self.ct || self.cu; 
      self.va = self.at || self.au; 
      self.vf = self.ft || self.fu; 
      self.vd = self.dt || self.du;
      
      self.MM = MM0;
      self.K = K0;
      self.M = M0;
      self.F = F0;
      self.Q = Q0;
      self.G = G0;
      self.H = H0;
      self.R = R0;
      
      [KK,FF,B,ud]=assempde(K0,M0,F0,Q0,G0,H0,R0);
      [self.totalNumEqns, self.numConstrainedEqns]=size(B);
      [self.Nu,self.Or]=pdenullorth(H0);
      self.B=B;
      self.ud=ud;
      self.ts=max(tlist)-min(tlist);

    end
    
  end
  
  methods(Access=private)
    
    function [MM,K,M,F,Q,G,H,R] = getMats(self, u, time)
      if(self.isMatrixMode)
        MM=self.MM;
        K=self.K;
        M=self.M;
        F=self.F;
        Q=self.Q;
        G=self.G;
        H=self.H;
        R=self.R;
      else
        [~,MM]=self.thePde.assema(self.p,self.t,0,self.d,self.f,u,time);
        [K,M,F]=self.thePde.assema(self.p,self.t,self.c,self.a,self.f,u,time);
        [Q,G,H,R]=assemb(self.b,self.p,self.e,u,time);
      end
    end
    
    function u0=evalU0Expr(self, u0In)
        if(ischar(u0In))
            % String case
            x=self.p(1,:)';
            y=self.p(2,:)';
            if(size(self.p,1)==2)
                % 2D case
                args = {x,y};
            elseif(size(self.p,1)==3)
                % 3D case
                z = self.p(3,:)';
                args = {x,y,z};
            end
            if pdeisfunc(u0In)
                u0=feval(u0In,args{:}); % supposedly a function name
            else
                if(size(u0In,1)==1) % numeric/string from GUI or string from CLI for scalar
                    u0 = eval(u0In);
                    if(~isnumeric(u0)) % must be string from GUI,only works for scalar
                        u0 = eval(u0);
                    end
                else                   % CLI systems case only, PDE tool will not work
                    u0 = zeros(size(self.p,2),size(u0In,1));
                    for i=1:size(u0In,1)
                        u0(:,i) = eval(u0In(i,:)); % supposedly an expression of 'x','y' and 'z'
                    end
                    u0 = u0(:);
                end
            end
        else
            % Numeric case
            u0 = u0In;
        end      
    end
    
  end
  
  methods(Access=protected, Static)
    
    function eq=eqMat(m1, m2)
      mm1 = max(abs(m1(:)));
      mm2 = max(abs(m2(:)));
      maxVal = max(mm1, mm2);
      if(any(abs(m1(:)-m2(:)) > eps(maxVal)))
        eq=false;
      else
        eq=true;
      end
    end
    
    function checkTlist(tlist)
      if length(tlist) < 2
        error(message('pde:pdeODEInfo:tlistTooShort'));
      end
      t0 = tlist(1);
      tfinal = tlist(end);
      if(t0 == tfinal)
        error(message('pde:pdeODEInfo:tlistEndpointsNotDistinct'));
      end
      tdir = sign(tfinal - t0);
      if any( tdir*diff(tlist) <= 0 )
        error(message('pde:pdeODEInfo:tlistUnordered'));
      end
    end
    
  end
  
end

