classdef (Sealed) GeometricInitialConditions  < pde.PDEInitialConditions &pde.PDECompactDisplay 
% GeometricInitialConditions Initial Conditions over a region or region boundary
%     Initial conditions (ICs) for time dependent problems or the initial guess 
%     for nonlinear stationary problems. The PDE toolbox can solve equations 
%     of the form:
%    
%                m*d^2u/dt^2 + d*du/dt - div(c*grad(u)) + a*u = f
%  
%     The equation to solve is defined in terms of the coefficients m, d, c, 
%     a, f. If the PDE to solve has a non-zero 'm' or 'd' coefficient, then
%     the PDE is time-dependent and Initial Conditions must be specified. 
%     If the 'm' and 'd' coefficients are both zero and the 'c', 'a', or 'f'   
%     coefficient are a function of u or grad(u), then the PDE is a nonlinear 
%     stationary problem. An optional initial guess can be specified, otherwise 
%     a default initial guess of zero is assumed.
%  
%     Instances of this class can only be created by calling the 
%     setInitialConditions method of the PDEModel class. 
%  
%   See also pde.PDEModel, pde.PDEModel/setInitialConditions
 
% Copyright 2015-2022 The MathWorks, Inc.


properties (SetAccess = private)
% RegionType - Type of geometric region to which ICs are specified.
%   The RegionType can be one of the following geometric entities: 'Cell',
%   'Face', 'Edge', or 'Vertex'.
    RegionType;
end

properties
% RegionID - ID of the geometric regions the ICs are assigned to. Where 
%      RegionID satisfies the following:
%        RegionType 'cell'   0 < RegionID(j) < NumCells
%        RegionType 'face'   0 < RegionID(j) < NumFaces
%        RegionType 'edge'   0 < RegionID(j) < NumEdges
%        RegionType 'vertex' 0 < RegionID(j) < NumVertices
%   j is the j'th element of the RegionID column vector. NumFaces and 
%   NumEdges are the properties of the geometric model.
    RegionID;

% InitialValue - The initial value or guess 
%    The initial value for time-dependent problems or the initial-guess for 
%    stationary nonlinear problems. The initial value can be one of the
%    following:
%     (1) A scalar representing a constant initial value/guess throughout 
%         the region.
%     (2) A column vector of length PDESystemSize representing a constant initial 
%         value/guess for each equation across the region.
%     (3) A function handle representing a spatially varying initial value/guess 
%         across the region. 
    InitialValue;

% InitialDerivative	- The initial derivative	
%     The value of the initial derivative for time-dependent problems.
%       The initial derivative can be one of the following:
%      (1) A scalar representing a constant initial derivative value across 
%          the region.
%      (2) A column vector of length PDESystemSize representing a constant initial 
%          derivative value for each equation across the region.
%      (3) A function handle representing a spatially varying initial derivative 
%          value across the region.
    InitialDerivative;
end

methods(Hidden=true, Access={?pde.PDEModel})
    function obj=GeometricInitialConditions(varargin) 
      obj.RecordOwner = varargin{1};
      parser = inputParser;
      obj.InitialValue = [];
      obj.InitialDerivative = [];
      inputargs = varargin;
      if isnumeric(varargin{2}) || isa(varargin{2}, 'function_handle') 
          obj.InitialValue = varargin{2};
          inputargs = {varargin{3:end}};
      end
      if isnumeric(varargin{3}) || isa(varargin{3}, 'function_handle') 
          obj.InitialDerivative = varargin{3};
          inputargs = {varargin{4:end}};
      end            
      parser.addParameter('vertex', []);    
      parser.addParameter('edge', []);     
      parser.addParameter('face', []);     
      parser.addParameter('cell', []);            
      parser.parse(inputargs{:});
     
      if ~isempty(parser.Results.vertex)
              obj.RegionType = 'vertex';
              obj.RegionID = parser.Results.vertex;
      elseif ~isempty(parser.Results.edge)
              obj.RegionType = 'edge';
              obj.RegionID = parser.Results.edge;
       elseif ~isempty(parser.Results.face)
              obj.RegionType = 'face';
              obj.RegionID = parser.Results.face;
      elseif ~isempty(parser.Results.cell)
             obj.RegionType = 'cell';
             obj.RegionID = parser.Results.cell;    
      end            
      systemsize = obj.RecordOwner.ParentPdemodel.PDESystemSize;
      obj.validateICs(systemsize);      
    end
end

methods    
    function set.RegionType(self, rtype)     
       % self.ValidateRegionType(rtype);     
        self.RegionType = rtype;       
    end
    
    function set.RegionID(self, rids)      
        self.ValidateRegionID(rids);      
        self.RegionID = rids;        
    end        
    
    function set.InitialValue(self, val)    
      self.precheckIC(val);     
      self.InitialValue = val;     
    end
    
    function set.InitialDerivative(self, val) 
      self.precheckIC(val);
      self.InitialDerivative = val;   
    end        
end

methods(Hidden=true, Access = {?pde.InitialConditionsRecords})
                      
    function tf=numericIcValue(self)                
         tf = ~isa(self.InitialValue, 'function_handle');                                        
    end
    function tf=numericIcDerivativeValue(self)                
         tf = ~isa(self.InitialDerivative, 'function_handle');                                        
    end         
    function performSolverPrecheck(self, systemsize, numvertices, numedges, numfaces, numcells)  
        if strcmp(self.RegionType, 'vertex')
            if any(self.RegionID > numvertices)
                error(message('pde:pdeInitialConditions:invalidVertexIndexPresolve')); 
            end
        elseif strcmp(self.RegionType, 'edge')
            if any(self.RegionID > numedges)
                error(message('pde:pdeInitialConditions:invalidEdgeIndexPresolve')); 
            end
         elseif strcmp(self.RegionType, 'face')
            if any(self.RegionID > numfaces)
                error(message('pde:pdeInitialConditions:invalidFaceIndexPresolve')); 
            end
         else
            if any(self.RegionID > numcells)
                error(message('pde:pdeInitialConditions:invalidCellIndexPresolve')); 
            end        
         end  
         validateICs(self, systemsize);     
    end         
end

methods(Hidden=true)
    function tf = hasICValue(self) 
        tf = ~isempty(self.InitialValue);
    end
    function tf = hasICDerivativeValue(self)
        tf = ~isempty(self.InitialDerivative);
    end     
end

methods(Hidden=true, Access = private) 
    %
    % validateICs - check the size of the IC matrices
    % 
    function validateICs(self, systemsize)
       if isempty(self.InitialValue)
         error(message('pde:pdeInitialConditions:emptyICValue'));      
       end      
       self.checkIcSize(self.InitialValue, systemsize);
       if hasICDerivativeValue(self)
            self.checkIcSize(self.InitialDerivative, systemsize);           
       end
    end       
end


methods(Static, Access = private)
    function ok=ValidateRegionID(rgnid)
      % Must be real(non-complex), full, natural number.        
      if ~isreal(rgnid) || ~all(rgnid(:) > 0) || issparse(rgnid) || any(mod(rgnid(:),1)~=0)
        error(message('pde:pdeInitialConditions:invalidRegionID'));   
      end      
      ok = true; 
    end
                
    function precheckIC(icval)   
      if isempty(icval)          
          return
      end
      if isfloat(icval)          
        if issparse(icval)
            error(message('pde:pdeInitialConditions:invalidICValueSparse'));   
        elseif any(isnan(icval))
            error(message('pde:pdeInitialConditions:invalidICValueNaN'));
        elseif any(isinf(icval))
            error(message('pde:pdeInitialConditions:invalidICValueInf'));      
        elseif ~isscalar(icval) && ~(iscolumn(icval))
            error(message('pde:pdeInitialConditions:invalidICColVector'));    
        end
      elseif isa(icval, 'function_handle')    
         pde.GeometricInitialConditions.checkIcFcnHdlArgCounts(icval);             
      end          
    end
                    
    function checkIcSize(icval, systemsize)  
         pde.GeometricInitialConditions.precheckIC(icval);
         if (~isscalar(icval) && numel(icval) ~= systemsize)
            error(message('pde:pdeInitialConditions:invalidICColVector'));    
         end                              
    end
   
    
%
%   Coefficient function-handle argument checks
%   
    function checkIcFcnHdlArgCounts(icfcn)        
        location.x = 0;
        location.y = 0;          
        location.z = 0;
        location.subdomain=1;                                 
        try
            icfcn(location);
        catch 
            error(message('pde:pdeInitialConditions:invalidFcnHandleArgs'));   
        end 
        
        try
            res = icfcn(location); %#ok<NASGU>
        catch 
            error(message('pde:pdeInitialConditions:invalidFcnHandleArgs'));
        end 
              
    end  
    
end   

methods(Static, Hidden=true)
    % preDelete - Preserved for backward compatibility. In R2017b and prior
    % releases, the coefficient assignment has an ObjectBeingDestroyed 
    % callback that calls this function.
    function preDelete(self,~)
        if ~isempty(self.RecordOwner) && isvalid(self.RecordOwner)          
            self.RecordOwner.delistPDEInitialConditions(self);
        end
    end
end

methods(Hidden=true)
  function delete(self)
       if ~isempty(self.RecordOwner) && isvalid(self.RecordOwner)
            self.RecordOwner.delistPDEInitialConditions(self);
       end
  end
end


properties (Hidden = true, Access='protected')
    RecordOwner;
end

end
