classdef (Sealed) BoundaryCondition < pde.PDECompactDisplay 
%BoundaryCondition Defines a boundary condition (BC) for the PDE
%
%     The PDE toolbox supports Dirichlet BC of the form:
%    
%                h*u = r
%      
%     and Generalized Neumann BC of the form:
% 							
%                n*(c*grad(u)) + q*u = g
%
%     BoundaryCondition objects is defined for each region of the boundary 
%     where boundary conditions are required. If multiple regions have identical 
%     boundary conditions, they can be defined in a single BoundaryCondition. 
%     The default boundary condition for a region of boundary is zero Neumann, 
%     with g = 0. A BoundaryCondition is generally not created directly, it 
%     is created using the applyBoundaryCondition method of the PDEModel class. 
%     The BC object is appended to the BoundaryConditionsRecords. 
% 
%     BoundaryCondition properties:
%          BCType - Type of BC is assigned to 
%          RegionType - Type of region the BC is assigned to 
%          RegionID - ID of the region the BC is assigned to
%          r - r-coefficient in the Dirichlet boundary condition equation  
%          h - h-coefficient in the Dirichlet boundary condition equation
%          g - g-coefficient in the Neumann boundary condition equation
%          q - q-coefficient in the Neumann boundary condition equation
%          u - prescribed values of dependent variables on a boundary
%          EquationIndex - equation number (1-N) for the prescribed values in u
%          Vectorized - Indicates whether a user-defined function accepts vector
%                       arguments.
% 
%     See also pde.PDEModel, pde.PDEModel/applyBoundaryCondition

%     Copyright 2014-2022 The MathWorks, Inc.
  
  properties
    % BCType - Type of the BC is assigned to 
    BCType;
    % RegionType - Type of region the BC is assigned to 
    RegionType;
    % RegionID - ID of the region the BC is assigned to
    RegionID;       
    % r - r-coefficient in the Dirichlet boundary condition equation
    r;
    % h - h-coefficient in the Dirichlet boundary condition equation
    h;
    % g - g-coefficient in the Neumann boundary condition equation
    g;
    % q - q-coefficient in the Neumann boundary condition equation
    q;
    % u - prescribed values of dependent variables on a boundary
    u;
    % EquationIndex - equation number (1-N) for the prescribed values in u
    EquationIndex;
    % Vectorized - Indicates whether a user-written function will return
    %              boundary condition matrices at multiple locations.
    Vectorized;
  end
  
  properties(Hidden = true, Transient = true)  
    ApplicationRegion;  
  end
  
  properties (Hidden = true, Access='private')
      RecordOwner;
  end  
  
  methods
    function obj=BoundaryCondition(bcr,varargin)
      obj.RecordOwner = bcr;
      import pde.*  
      import pde.internal.*
      narginchk(1,16);
      obj.CoeffsSet = BoundaryCondition.noneSet;
      parser = inputParser;
      if nargin > 0 && (isa(varargin{1}, 'pdeEdge') || isa(varargin{1}, 'pdeFace'))
        parser.addRequired('ApplicationRegion', @BoundaryCondition.ValidateAppRegion);
      else
        parser.addRequired('BCType', @BoundaryCondition.ValidateBCType);
        parser.addParameter('Face', [], @BoundaryCondition.ValidateRegionID);
        parser.addParameter('Edge', [], @BoundaryCondition.ValidateRegionID);        
      end              
      parser.addParameter('Vectorized', 'off', @pdeIsValid.onOrOff);
      parser.addParameter('r', [], @BoundaryCondition.ValidateBCVal);
      parser.addParameter('h', [], @BoundaryCondition.ValidateBCVal);
      parser.addParameter('g', [], @BoundaryCondition.ValidateBCVal);
      parser.addParameter('q', [], @BoundaryCondition.ValidateBCVal);
      parser.addParameter('EquationIndex', [], @BoundaryCondition.ValidateBCVal);
      parser.addParameter('u', [], @BoundaryCondition.ValidateBCVal);
      parser.addParameter('SystemSize', 1);
      
      parser.parse(varargin{:});
      if isa(varargin{1}, 'pdeEdge') ||  isa(varargin{1}, 'pdeFace')
        obj.ApplicationRegion = varargin{1};
      else
          obj.BCType = varargin{1};
          if ~isempty(parser.Results.Face) &&  ~isempty(parser.Results.Edge)
             error(message('pde:pdeBoundaryConditions:singleRegionType'));
          elseif isempty(parser.Results.Face) && isempty(parser.Results.Edge)
             error(message('pde:pdeBoundaryConditions:noRegionType'));
          end
          if ~isempty(parser.Results.Face)
             obj.RegionType = 'Face';
             obj.RegionID = parser.Results.Face;
          elseif ~isempty(parser.Results.Edge)
             obj.RegionType = 'Edge';
             obj.RegionID = parser.Results.Edge;                
          end          
      end
      % convert old BC to having a BCType
      systemsize = parser.Results.SystemSize;     
      nc = numel(obj.BCType);
      if strncmpi(obj.BCType,'unknown',nc)
          if (~isempty(parser.Results.EquationIndex) && (length(parser.Results.EquationIndex)~=systemsize))
              obj.BCType = 'mixed';
          elseif ~isempty(parser.Results.g) || ~isempty(parser.Results.q)
              obj.BCType = 'neumann';
          else
              obj.BCType = 'dirichlet';
          end
      end
       
      if strncmpi(obj.BCType,'mixed',nc) 
          if ~isempty(parser.Results.u)
            obj.Type = 'mixedvalue';
          else
             obj.Type = 'mixeddirichlet';
          end
      elseif strncmpi(obj.BCType,'dirichlet',nc) 
          if ~isempty(parser.Results.u)
            obj.Type = 'value';
          else
             obj.Type = 'dirichlet';
          end
      else
          obj.Type = 'neumann';
      end
      
      obj.Vectorized = parser.Results.Vectorized;
      obj.r = parser.Results.r;
      obj.h = parser.Results.h;
      obj.g = parser.Results.g;   
      obj.q = parser.Results.q;
      obj.EquationIndex = parser.Results.EquationIndex;
      obj.u = parser.Results.u;
      
      obj.ValidateBCTypeCoeff(systemsize);
          
    end
    
    function set.BCType(self, bctype)     
        self.ValidateBCType(bctype);
        nc = numel(bctype);
        if strncmpi(bctype,'mixed',nc)
            self.BCType = 'mixed';       
        elseif strncmpi(bctype,'neumann',nc)
            self.BCType = 'neumann';  
        elseif strncmpi(bctype,'dirichlet',nc)
            self.BCType = 'dirichlet';      
        else
            self.BCType = bctype; 
        end
    end
    
    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.ApplicationRegion(self, appReg)
      import pde.internal.*   
      self.ValidateAppRegion(appReg);     
      self.ApplicationRegion = appReg;
      if isa(appReg, 'pdeEdge')
        self.RegionType = 'Edge';
      else
        self.RegionType = 'Face';
      end
      eids = ones(size(appReg));
      nume = numel(appReg);
      for i = 1:nume
        eids(i) = appReg(i).ID;  
      end
       self.RegionID = eids;   
    end
    
    function set.r(self, coef)    
      self.ValidateCoeff(coef, self.rSet, self.gSet|self.qSet|self.rSet|self.hSet);     
      self.r = coef;
    end
    
    function set.h(self, coef) 
      self.ValidateCoeff(coef, self.hSet, self.gSet|self.qSet|self.rSet|self.hSet);     
      self.h = coef;
    end
    
    function set.g(self, coef)  
      self.ValidateCoeff(coef, self.gSet, ~self.noneSet);      
      self.g = coef;
    end
    
    function set.q(self, coef)    
      self.ValidateCoeff(coef, self.qSet, ~self.noneSet);     
      self.q = coef;
    end
    
    function set.u(self, coef)  
      self.ValidateCoeff(coef, self.uSet, self.gSet|self.qSet|self.uSet|self.iSet);      
      self.u = coef;
    end
    
    function set.EquationIndex(self, coef)    
      self.ValidateCoeff(coef, self.iSet, self.gSet|self.qSet|self.uSet|self.iSet);     
      self.EquationIndex = coef;
    end
    
    function set.Vectorized(self, vec)
      try
        pdeIsValid.onOrOff(vec);
      catch exc
        throwAsCaller(exc);
      end
      self.Vectorized = vec;
    end
    
  end % methods
  
  methods(Static, Hidden=true)
      function obj = loadobj(A)
          if ~(isprop(A(1),'RecordOwner') || isfield(A(1),'RecordOwner'))
              mypde = createpde();
              obj = pde.BoundaryConditionRecords(mypde);
              for i = 1:length(A)
                  if ~isfield(A(i), 'BCType')
                      % Old way never allows 'mixed' case BC
                      nc = numel(A(i).Type);
                      if strncmpi(A(i).Type,'neumann',nc)
                          A(i).BCType = A(i).Type;
                      else
                          A(i).BCType = 'dirichlet';
                      end
                  end
                  bc = pde.BoundaryCondition(obj,A(i).BCType,A(i).RegionType,A(i).RegionID);
                  bc.Vectorized = A(i).Vectorized;
                  bc.r = A(i).r;
                  bc.h = A(i).h;
                  bc.g = A(i).g;
                  bc.q = A(i).q;
                  bc.EquationIndex = A(i).EquationIndex;
                  bc.u = A(i).u;
                  bc.Type = A(i).Type;
                  if isempty(obj.BoundaryConditionAssignments)
                    obj.BoundaryConditionAssignments = bc;
                  else
                    obj.BoundaryConditionAssignments(end+1) = bc;
                  end
              end
          else
              if ~isfield(A, 'BCType')
                  % Old way never allows 'mixed' case BC
                  nc = numel(A.Type);
                  if strncmpi(A.Type,'neumann',nc)
                      A.BCType = A.Type;
                  else
                      A.BCType = 'dirichlet';
                  end
              end
              obj = pde.BoundaryCondition(A.RecordOwner,A.BCType,A.RegionType,A.RegionID);
              obj.Vectorized = A.Vectorized;
              obj.r = A.r;
              obj.h = A.h;
              obj.g = A.g;
              obj.q = A.q;
              obj.EquationIndex = A.EquationIndex;
              obj.u = A.u;
              obj.Type = A.Type;
          end
      end
  end
  
  properties(Access = {?pde.internal.pde3DBCImpl})
    % Type - string indicating type of boundary condition
    %        Allowable values are 'Value', 'Neumann', 'Dirichlet',
    %        'mixedvalue' and 'mixeddirichlet'
    Type;
    % CoeffsSet - Boolean array indicating which coefficients have been set
    CoeffsSet;
  end
  
  methods(Access = private)
    
      function ValidateCoeff(self, coeff, cInd, cSet)
          if(~isempty(coeff))
              self.ValidateBCVal(coeff);
              % Is this coefficient compatible with the ones already set?
              if(any((~cSet & self.CoeffsSet)))
                  msg = '';
                  startMsg = true;
                  for i=1:length(self.CoeffsSet)
                      ci = self.CoeffsSet(i);
                      if(ci)
                          if ~startMsg
                              % separate parameter names with commas
                              msg = [msg ',']; %#ok
                          end
                          msg = [msg ' ' self.coefNames{i}]; %#ok
                          startMsg = false;
                      end
                  end
                  error(message('pde:pdeBoundaryConditions:invalidParameter', ...
                      self.coefNames{cInd}, msg));
              end
              % Flag this coefficient as having been set.
              self.CoeffsSet(cInd) = true;
          else
              % Coefficient is empty so flag it as not set.
              self.CoeffsSet(cInd) = false;
          end
      end
      
      function ValidateBCTypeCoeff(self,systemsize)
          % identify conflict BC coefficient assignment with BC type
          nc = numel(self.BCType);
          if strncmpi(self.BCType,'dirichlet',nc)
              if ~(isempty(self.g) && isempty(self.q))
                  error(message('pde:pdeBoundaryConditions:invalidDirichlet'));
              elseif (~isempty(self.EquationIndex)) && (length(self.EquationIndex)<systemsize)
                  error(message('pde:pdeBoundaryConditions:nonDirichlet'));
              end
          elseif strncmpi(self.BCType,'neumann',nc)
              if ~(isempty(self.u) && isempty(self.EquationIndex) ...
                      && isempty(self.h) && isempty(self.r))
                  error(message('pde:pdeBoundaryConditions:invalidNeumann'));
              end
          else % mixed case
              if ~isempty(self.u)
                  if isempty(self.EquationIndex)
                      error(message('pde:pdeBoundaryConditions:nonMixed'));
                  elseif length(self.EquationIndex)==systemsize
                      error(message('pde:pdeBoundaryConditions:nonMixed'));
                  end
              end
              if ~(isempty(self.q) && isempty(self.g))
                  if isempty(self.EquationIndex) && isempty(self.u) && isempty(self.r) && isempty(self.h)
                      error(message('pde:pdeBoundaryConditions:nonMixed'));
                  end                  
              end              
          end
      end
       
  end
  
  methods(Hidden=true, Access = {?pde.PDEModel})      
      function tf=numericCoefficients(self,varargin) 
          tf = ~(isa(self.r, 'function_handle')  || ...
                isa(self.h, 'function_handle')  || ...
                isa(self.g, 'function_handle')  || ...
                isa(self.q, 'function_handle')  || ...
                isa(self.u, 'function_handle'));
      end
  end
  
  methods(Static, Access = private)
    function ok=ValidateBCVal(v)
      ok = isa(v, 'function_handle') || isnumeric(v);
      if(~ok)
        error(message('pde:pdeBoundaryConditions:invalidTermType'));
      end
      if isa(v, 'function_handle') && ~(nargin(v) == 2 || nargin(v) == 3)
         error(message('pde:pdeBoundaryConditions:invalidFunctionHandle'));
      end
    end
    
    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:pdeBoundaryConditions:invalidRegionID'));   
      end      
      ok = true; 
    end
    
     function ok=ValidateRegionType(rgntype)   
      nc = numel(rgntype);
      if ~(strncmpi(rgntype,'Face',nc) || strncmpi(rgntype,'Edge',nc))
        error(message('pde:pdeBoundaryConditions:invalidRegionType'));   
      end
      ok = true;      
     end
     
     function ok=ValidateBCType(bctype)   
      nc = numel(bctype);
      if ~(strncmpi(bctype,'dirichlet',nc) || strncmpi(bctype,'neumann',nc)...
              || strncmpi(bctype,'mixed',nc) || strncmpi(bctype,'unknown',nc))
        error(message('pde:pdeBoundaryConditions:invalidBCType'));   
      end
      ok = true;      
     end
    
    function ok=ValidateAppRegion(appReg)
      % In general, application region must be an instance of the abstract
      % pdeBoundaryGeom class. For now, only pdeEdge qualifies.
      import pde.internal.*
      if(~ isa(appReg, 'pdeBoundaryGeom'))
        error(message('pde:pdeBoundaryConditions:invalidAppRegion'));
      end
      ok = true;
    end
    
  end % static methods
    
  methods(Hidden=true)
      function delete(self)
          if ~isempty(self.RecordOwner) && isvalid(self.RecordOwner) 
              self.RecordOwner.delistBoundaryConditions(self);
          end
      end
  end
  
  properties(Constant, Access=private)
    coefNames = {'r', 'h', 'g', 'q', 'u', 'EquationIndex'};
    % Define Boolean arrays for the different option names
    % These will be used for checking valid combinations of options
    % using Boolean ops in the "set" functions defined above.
    noneSet = logical([0 0 0 0 0 0]);
    rSet = logical([1 0 0 0 0 0]);
    hSet = logical([0 1 0 0 0 0]);
    gSet = logical([0 0 1 0 0 0]);
    qSet = logical([0 0 0 1 0 0]);
    uSet = logical([0 0 0 0 1 0]);
    iSet = logical([0 0 0 0 0 1]);
  end

end

