classdef pdeEquations < handle
% pdeEquations Calculate the K and M matrices and the F vector for tetrahedral models.
% Models may be composed of 4-node and 10-node tetrahedron elements.
% This undocumented function may be changed or removed in a future release.

%       Copyright 2014-2016 The MathWorks, Inc.

  properties(SetAccess = private)
    BoundaryConditions, Points, MeshAssoc, ElemConn;
    N, Dimensionality;
  end
  
  properties(SetAccess = private, GetAccess=public)
    K, F, M, Q, G, H, R;
  end
  
  methods(Access=public)
    
    function obj = pdeEquations(b,p,gm,t,c,a,f, varargin)
      import pde.internal.*;
      parser = inputParser;
      parser.addRequired('b', @pdeEquations.isValidB);
      parser.addRequired('p', @pdeEquations.isValidP);
      parser.addRequired('gm', @pdeEquations.isValidGm);
      parser.addRequired('t', @pdeEquations.isValidT);
      parser.addRequired('c', @pdeEquations.isValidC);
      parser.addRequired('a', @pdeEquations.isValidA);
      parser.addRequired('f', @pdeEquations.isValidF);
      parser.addOptional('u', [], @pdeEquations.isValidU);
      parser.addOptional('time', [], @pdeEquations.isValidTime);
      parser.addOptional('sdl', [], @isnumeric);
      parser.parse(b,p,gm,t,c,a,f, varargin{:});
      
      obj.BoundaryConditions = b;
      obj.thePde = b;
      obj.Points = p;
      obj.MeshAssoc = gm;
      obj.ElemConn = t;
      obj.Dimensionality = size(p,1);
      [obj.isFuncC, obj.c] = pdeEquations.setCoeff(c);
      [obj.isFuncA, obj.a] = pdeEquations.setCoeff(a);
      [obj.isFuncF, obj.f] = pdeEquations.setCoeff(f);
      obj.u = parser.Results.u;
      if(length(obj.u) == 1)
        obj.u = obj.u*ones(size(p,2)*obj.N, 1);
      end
      obj.time = parser.Results.time;
	  obj.N = b.PDESystemSize;
    end
    
  end
  
  methods
    [ K,M,F ] = getKMF(self);
  end
  
  methods(Access=private, Static)
    
    function [isFuncHandle, c]=setCoeff(coef)
      if(isa(coef, 'function_handle'))
        isFuncHandle = true;
        c = coef;
      elseif(ischar(coef))
        isFuncHandle = true;
        % string. may be a function name or string expression
        if(pdeisfunc(coef))
          c = str2func(coef);
        else
          c = @(thePde, loc, state) pde.internal.pdeEquations.stringExprCoefEval(coef, thePde, loc, state);
        end
      else
        isFuncHandle = false;
        c = coef;
      end
    end
    
    function c=stringExprCoefEval(strExpr, thePde, loc, state)
	  % Handle the case when the coefficient is defined as a matrix of string expressions
      x = loc.x; y = loc.y; z = loc.z;
      sd = loc.subdomain;
      stateArgs = {};
      if(~ isempty(state))
        if(isfield(state, 'time'))
          time = state.time;
          stateArgs = [stateArgs state.time];
        else
          time = [];
        end
        if(isfield(state, 'u'))
          u = state.u;
          ux = state.ux;
          uy = state.uy;
          uz = state.uz;
        else
          u = []; ux = []; uy = []; uz = [];
        end
      else
        time = []; u = []; ux = []; uy = []; uz = [];
      end
      ncx = size(x,2);
      nrc = size(strExpr, 1);
      c = zeros(nrc, ncx);
      for i=1:nrc
        %ci = eval(strExpr(i,:));
        ci = pde.internal.pdeEquations.stringExprCoefEntryEval(x,y,z,sd,strExpr(i,:), u,ux,uy,uz,time);
        %ci = pde.internal.pdeEquations.stringExprCoefEntryEval(x,y,z,sd,strExpr(i,:), stateArgs{:});
        if(~(isscalar(ci) || all(size(ci)==[1,ncx])))
          error(message('pde:pdeEquations:invalidNumEntriesStringFunc', strExpr(i,:), ...
            ncx, size(ci,1), size(ci,2)));
        end
        c(i,:) = ci;
      end
    end
    
    function ci = stringExprCoefEntryEval(x,y,z,sd,f,u,ux,uy,uz,t)
	% Evaluate a single coefficient entry defined by a string expression.
      if(isempty(u))
        clear('u', 'ux', 'uy', 'uz');
      end
      if(isempty(t))
        clear('t');
      end
      ncx = size(x,2);
      ic=strfind(f,'!');
      nf=length(ic)+1;
      if nf>1
        X=x;
        Y=y;
        Z=z;
        SD=sd;
        fff=zeros(1,ncx);
        lf=length(f);
        ic=[0 ic lf+1];
        for ii=1:nf
          ff=f(ic(ii)+1:ic(ii+1)-1);
          iii=find(SD==ii);         
          ffold=ff;
          try
            ff=eval(ffold);
          catch ex
            error(message('pde:pdeEquations:evalSubdomainStringExpr', ii, ...
              ffold, ex.message));
          end
          fff(iii)=ff.*ones(size(iii));
        end
        ci = fff;
      else
        try
          ci = eval(f);
        catch ex
          error(message('pde:pdeEquations:evalStringExpr', f, ex.message));
        end
      end
      if(~(isscalar(ci) || all(size(ci)==[1,ncx])))
          error(message('pde:pdeEquations:invalidNumEntriesStringFunc', f, ...
            ncx, size(ci,1), size(ci,2)));
      end
    end
    
    function ok = isValidB(b)
      ok = isa(b, 'pde.PDEModel');
    end
    
    function ok = isValidP(v)
      ok = isnumeric(v);
    end
    
    function ok = isValidGm(v)
      ok = true;
    end
    
    function ok = isValidT(v)
      ok = isnumeric(v);
    end
    
    function ok = isValidC(v)
      ok = pdeIsValid.coefficient(v);
    end
    
    function ok = isValidA(v)
      ok = pdeIsValid.coefficient(v);
    end
    
    function ok = isValidF(v)
      ok = pdeIsValid.coefficient(v);
    end
    
    function isValidU(v)
      validateattributes(v,{'numeric'}, {});
    end
    
    function isValidTime(v)
      if ~isempty(v)
        validateattributes(v,{'numeric'},...
          {'scalar','real','nonnegative'});
      end
    end
    
  end %methods
  
  properties(Access = private)
    c,a,f,u,time;
    isFuncC, isFuncA, isFuncF, isFuncD;
    thePde;
  end
  
end

