function regionOfInterest = linearizeOutput(self,op)
%LINEARIZEOUTPUT  Specify outputs of linearized model.
%
%   LINEARIZEINPUT and LINEARIZEOUTPUT are used in conjunction with
%   LINEARIZE to extract sparse linear models from PDE models.
%
%   LINEARIZEOUTPUT(thermalmodel,'REGIONTYPE',REGIONID) selects all degrees
%   of freedom (DoF) associated with the region defined by REGIONTYPE and
%   REGIONID. For a 3-D model, REGIONTYPE is 'Face', 'Edge', or 'Vertex'.
%   For a 2-D model, REGIONTYPE is 'Edge' or 'Vertex'. REGIONID is
%   an integer specifying the ID of a face, edge, or vertex.
%
%   OUT = LINEARIZEOUTPUT(...) returns a structure OUT with the
%   linearization output description.
%
%   Example:
%      tmodel = createpde('thermal','transient');
%      geometryFromEdges(tmodel,@crackg);
%      thermalProperties(tmodel,'ThermalConductivity',1,...
%                            'MassDensity',1,'SpecificHeat',1);
%
%      thermalBC(tmodel,'Edge',1,'HeatFlux',-10,'Label','Flux');
%      internalHeatSource(tmodel,25,'Label','HeatSource');
%      % Create outputs for temperature of Edge 2
%      linearizeOutput(tmodel,'Edge',2);
%
%   See also pde.ThermalModel/linearizeInput, pde.ThermalModel/linearize.

%    Copyright 2021 The MathWorks, Inc.

arguments
    self pde.ThermalModel
    op.Vertex (1,1) double
    op.Edge (1,1) double
    op.Face (1,1) double
    op.Cell (1,1) double
end

Region=validateRegionType(self,op);
ID=validateRegionID(self,Region,op);

regionOfInterest= struct('RegionType',Region,'RegionID',ID);
if isempty(self.LinearizeOutputs)
    self.LinearizeOutputs=regionOfInterest;
else
    if requestIsInvalid(self,regionOfInterest)
        error(message('pde:linearize:RepeatedOutput',regionOfInterest.RegionType,regionOfInterest.RegionID));
    end
    self.LinearizeOutputs(end+1)=regionOfInterest;
end
end

function Region=validateRegionType(self,op)
fn=fieldnames(op);
p={'Vertex','Edge','Face','Cell'};
r=startsWith(p,fn,'ignoreCase',true);
idx=find(r);
if numel(idx)>1
    error(message('pde:linearize:TooManyRegions'));
end

if ~self.IsTwoD || self.IsTwoD && ~r(4)
    Region=p(idx(1));
elseif self.IsTwoD && r(4)
    error(message('pde:pdeSanityChecker:NoCellGeometry'));
end
end

function ID=validateRegionID(self,Region,op)
Region=Region{:};

if strcmp(Region,'Vertex')
    numG='NumVertices';
else
    numG=strcat('Num',Region,'s');
end

if ~any(op.(Region)>self.Geometry.(numG))
    ID=op.(Region);
else
    if strcmp(Region,'Vertex')
        error(message('pde:pdeModel:invalidVertexIndex'));
    elseif strcmp(Region,'Edge')
        error(message('pde:pdeModel:invalidEdgeIndex'));
    elseif strcmp(Region,'Face')
        error(message('pde:pdeModel:invalidFaceIndex'));
    else
        error(message('pde:pdeModel:invalidCellIndex'));
    end
end
end

function invalid=requestIsInvalid(self,req)
invalid=false;
op=self.LinearizeOutputs;
for i=1:numel(op)
    if isequaln(op(i),req)
        invalid=true;
    end
end
end
