function hs = internalHeatSource(self,varargin)
%internalHeatSource - Apply internal heat source to a thermal model
%   src = internalHeatSource(thermalmodel,HEATSOURCE) applies heat source,
%   internal heat generation, to a thermal model. HEATSOURCE
%   can be specified as
%        (1) A numerical scalar for a constant heat source.
%        (2) MATLAB function format for a non-constant or nonlinear heat source.
%
%  mtl = internalHeatSource(__,'Face',FACEID) applies heat source
%  only for the specified subdomain of a 2-D geometry with multiple
%  subdomains. FACEID is a positive integer in the range of 1 to
%  thermalmodel.Geometry.NumFaces.
%
%  mtl = internalHeatSource(__,'Cell',CELLID) assigns heat source
%  only for the specified subdomain of a 3-D geometry with multiple
%  subdomains. CELLID is a positive integer in the range of 1 to
%  thermalmodel.Geometry.NumCells.
%
%   Example 1: Specify internal heat source for a steady-state thermal model.
%
%       thermalmodel = createpde('thermal','steadystate');
%       gm = importGeometry(thermalmodel,'SquareBeam.STL');
%       thermalProperties(thermalmodel,'ThermalConductivity',0.2);
%       internalHeatSource(thermalmodel,25);
%
%   Example 2: Specify internal heat source on a specified subdomian.
%
%       thermalmodel = createpde('thermal','steadystate');
%       SQ1 = [3; 4; 0; 3; 3; 0; 0; 0; 3; 3];
%       D1 = [2; 4; 0.5; 1.5; 2.5; 1.5; 1.5; 0.5; 1.5; 2.5];
%       gd = [SQ1,D1];
%       sf = 'SQ1+D1';
%       ns = char('SQ1','D1');
%       ns = ns';
%       dl = decsg(gd,sf,ns);
%       geometryFromEdges(thermalmodel,dl);
%       thermalProperties(thermalmodel,'ThermalConductivity',50);
%       internalHeatSource(thermalmodel,25,'Face',1);
%
%
%    See also pde.ThermalModel, pde.ThermalModel/thermalBC
%             pde.ThermalModel/thermalIC, pde.ThermalModel/solve

%    Copyright 2016-2017 The MathWorks, Inc.

if nargin > 1
    [varargin{:}] = convertStringsToChars(varargin{:});
end

narginchk(2,4);
nargoutchk(0,1);


if isempty(self.Geometry)
    error(message('pde:ThermalModel:thermalModelHasNoGeom'));
end


parser = inputParser;
parser.addParameter('face', [], @pde.EquationModel.isValidEntityID);
parser.addParameter('cell', [], @pde.EquationModel.isValidEntityID);

srcVal = varargin{1};
parser.parse(varargin{2:end});
argsToPass = [varargin(2:end), {'HeatSource', srcVal}];

if ~(isa(srcVal,'double') || isa(srcVal,'function_handle'))
    error(message('pde:ThermalModel:heatSourceNotFuncORNum'))
end

%If a face (for 2-D) or a cell (for 3-D) is not explicitly specified then
%use the same heat source for all faces or cells.
if self.IsTwoD
    if ~isempty(parser.Results.cell)
        error(message('pde:ThermalModel:noHeatSrcOnCellFor2D'));
    end
    if isempty(parser.Results.face)
        argsToPass = [{'face',1:self.Geometry.NumFaces}, argsToPass];
    else
        self.isValidEntityID(parser.Results.face);
        if any(parser.Results.face > self.Geometry.NumFaces)
            error(message('pde:pdeModel:invalidFaceIndex'));
        end
    end
else
    if ~isempty(parser.Results.face)
        error(message('pde:ThermalModel:noHeatSrcOnBoundary'));
    end    
    if isempty(parser.Results.cell)
        argsToPass = [{'cell',1:self.Geometry.NumCells}, argsToPass];
    else
        self.isValidEntityID(parser.Results.cell);
        if any(parser.Results.cell > self.Geometry.NumCells)
            error(message('pde:pdeModel:invalidCellIndex'));
        end
    end
    
end


if isempty(self.HeatSources)
    hsContainer = pde.HeatSourceAssignmentRecords(self);
    hs = pde.HeatSourceAssignment(hsContainer,argsToPass{:});
    hsContainer.HeatSourceAssignments = hs;
    self.HeatSources = hsContainer;
else
    hs = pde.HeatSourceAssignment(self.HeatSources,argsToPass{:});
    self.HeatSources.HeatSourceAssignments(end+1) = hs;
end

% LocalWords:  thermalmodel HEATSOURCE mtl FACEID CELLID createpde steadystate
% LocalWords:  subdomian gd ns dl Func
