function ic = thermalIC(self,varargin)
% thermalIC - Set initial temperature or initial guess
%
% thermalIC(thermalmodel,T0) sets the initial temperature T0 for the entire
% domain of a transient model. T0 is a scalar or a function handle.
%
% thermalIC(thermalmodel,T0,REGIONTYPE,REGIONID) sets the initial
% temperature T0 for the subdomain specified by using REGIONTYPE and
% REGIONID for a transient thermal model. Here, REGIONTYPE can be 'Cell',
% 'Face', 'Edge', or 'Vertex' for a 3-D model and 'Face', 'Edge', or 'Vertex'
% for a 2-D model, and REGIONID is the ID of the geometric region.
%
% thermalIC(thermalmodel,TRESULTS) sets the initial temperature using the
% solution TRESULTS obtained by solving a thermal model. If TRESULTS is a
% transient thermal solution, then the temperature corresponding to the
% last time-step is chosen.
%
% thermalIC(thermalmodel,TRESULTS,iT) applies IC from of solution obtained
% using a transient thermal model using the time-step specified by iT.
%
%  Example 1: Set initial temperature and solve.
%
%         thermalmodel = createpde('thermal','transient');
%         gm = @squareg;
%         geometryFromEdges(thermalmodel,gm);
%         generateMesh(thermalmodel);
%         thermalProperties(thermalmodel,'ThermalConductivity',50,...
%                                     'MassDensity',2500,...
%                                     'SpecificHeat',600);
%         bc = thermalBC(thermalmodel,'Edge',[1,3],'Temperature',100);
%         qs = internalHeatSource(thermalmodel,2);
%         ic = thermalIC(thermalmodel,0);
%         tlist = 0:0.1:1;
%         result1 = solve(thermalmodel,tlist)
%
%
%  Example 2: Restart thermal analysis using the solution from Example 1.
%
%         thermalmodelR1 = createpde('thermal','transient');
%         gm = @squareg;
%         geometryFromEdges(thermalmodelR1,gm);
%         generateMesh(thermalmodelR1);
%         thermalProperties(thermalmodelR1,'ThermalConductivity',50,...
%                                          'MassDensity',2500,...
%                                          'SpecificHeat',600);
%         bc = thermalBC(thermalmodelR1,'Edge',[1,3],'Temperature',100);
%         qs = internalHeatSource(thermalmodelR1,2);
%         ic = thermalIC(thermalmodelR1,result1);
%         tlist = 1:0.1:2;
%         result2 = solve(thermalmodelR1,tlist)
%
%
%  Example 3: Restart thermal analysis using the solution from Example 2
%             initializing from a specified time.
%
%         thermalmodelR1Ti = createpde('thermal','transient');
%         gm = @squareg;
%         geometryFromEdges(thermalmodelR1Ti,gm);
%         generateMesh(thermalmodelR1Ti);
%         thermalProperties(thermalmodelR1Ti,'ThermalConductivity',50,...
%                                            'MassDensity',2500,...
%                                            'SpecificHeat',600);
%         bc = thermalBC(thermalmodelR1Ti,'Edge',[1,3],'Temperature',100);
%         qs = internalHeatSource(thermalmodelR1Ti,2);
%         ic = thermalIC(thermalmodelR1Ti,result2,5);
%         tlist = 0.4:0.1:1;
%         result3 = solve(thermalmodelR1Ti,tlist)
%
%    See also pde.ThermalModel, pde.ThermalModel/thermalBC
%             pde.ThermalModel/internalHeatSource,
%             pde.ThermalModel/thermalProperties, pde.ThermalModel/solve

%    Copyright 2016-2017 The MathWorks, Inc.

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


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

isnodal = false;
if isa(varargin{1}, 'pde.SteadyStateThermalResults')||isa(varargin{1}, 'pde.TransientThermalResults')
    if numel(varargin) == 1
        argsToPass = {varargin{1}};
    elseif numel(varargin) == 2
        argsToPass = {varargin{1},  varargin{2}};
    else
        error(message('pde:pdeModel:setICMultipleVarArgs'));
    end
    isnodal = true;
else
    argsToPass = ParseGeometricICs(self, varargin{:});
end

if isempty(self.InitialConditions)
    iccont = pde.ThermalICRecords(self);
    argsToPass = {iccont, argsToPass{1:end}};
    if isnodal
        ic = pde.NodalThermalICs(argsToPass{:});
    else
        ic = pde.GeometricThermalICs(argsToPass{:});
    end
    iccont.ThermalICAssignments = ic;
    self.InitialConditions = iccont;
    
else
    argsToPass = {self.InitialConditions, argsToPass{1:end}};
    if isnodal
        ic = pde.NodalThermalICs(argsToPass{:});
    else
        ic = pde.GeometricThermalICs(argsToPass{:});
    end
    self.InitialConditions.ThermalICAssignments(end+1) = ic;
end

end




function argsToPass = ParseGeometricICs(self, varargin)
if ~(isnumeric(varargin{1}) || isa(varargin{1}, 'function_handle'))
    error(message('pde:pdeModel:initialCondBadValue')); % Invalid U0
end

inputargs = {};
if numel(varargin) == 3
    inputargs = varargin(2:end);
elseif numel(varargin) == 4
    inputargs = varargin(3:end);
end

parser = inputParser;
parser.KeepUnmatched=true;
parser.addParameter('cell', [], @pde.EquationModel.isValidEntityID);
parser.addParameter('face', [], @pde.EquationModel.isValidEntityID);
parser.addParameter('edge', [], @pde.EquationModel.isValidEntityID);
parser.addParameter('vertex', [], @pde.EquationModel.isValidEntityID);
parser.parse(inputargs{:});

if numel(parser.UsingDefaults) == 4 % Global assignment
    if self.IsTwoD
        argsToPass = [varargin, {'face', 1:self.Geometry.NumFaces}];
    else
        argsToPass = [varargin, {'cell', 1:self.Geometry.NumCells}];
    end
else   % Local assignment
    if ~isempty(parser.Results.cell)
        self.isValidEntityID(parser.Results.cell);
        if any(parser.Results.cell > self.Geometry.NumCells)
            error(message('pde:pdeModel:invalidCellIndex'));
        end
    end
    if ~isempty(parser.Results.face)
        self.isValidEntityID(parser.Results.face);
        if any(parser.Results.face > self.Geometry.NumFaces)
            error(message('pde:pdeModel:invalidFaceIndex'));
        end
    end
    if ~isempty(parser.Results.edge)
        self.isValidEntityID(parser.Results.edge);
        if any(parser.Results.edge > self.Geometry.NumEdges)
            error(message('pde:pdeModel:invalidEdgeIndex'));
        end
    end
    if ~isempty(parser.Results.vertex)
        self.isValidEntityID(parser.Results.vertex);
        if any(parser.Results.vertex > self.Geometry.NumVertices)
            error(message('pde:pdeModel:invalidVertexIndex'));
        end
    end
    argsToPass = varargin;
end

end
