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 structural models.
%
%   LINEARIZEOUTPUT(structuralmodel,'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 a positive
%   integer specifying the ID of a face, edge, or vertex.
%
%   LINEARIZEOUTPUT(structuralmodel,...,'Component',XYZ)
%   further specifies which of the X,Y,Z coordinates to include. For
%   example,
%      linearizeOutput(pdemodel,'Face',10,'Component','xz')
%   only selects the X and Z DoFs for Face 10.
%
%   OUT = LINEARIZEOUTPUT(...) returns a structure OUT with the
%   linearization output description.
%
%   See also pde.StructuralModel/linearizeInput, pde.StructuralModel/linearize.

%    Copyright 2021-2022 The MathWorks, Inc.

arguments
    self pde.StructuralModel
    op.Vertex (1,1) double
    op.Edge (1,1) double
    op.Face (1,1) double
    op.Cell (1,1) double
    op.Component (1,:) char {mustBeMember(op.Component,{'x','y','z','xy','xz','yz','xyz'})}
end

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

if ~isfield(op,'Component') %When 'Component' is not specified
    regionOfInterest= struct('RegionType',Region,...
        'RegionID',ID,'Component','xyz');

    if isempty(self.LinearizeOutputs)
        self.LinearizeOutputs=regionOfInterest;
    else
        if outputIsInvalid(self,regionOfInterest)
            error(message('pde:linearize:RepeatedOutput',regionOfInterest.RegionType,regionOfInterest.RegionID));
        end
        self.LinearizeOutputs(end+1)=regionOfInterest;
    end
else
    %When 'Component' is specified
    validateComponent(self,op.Component);

    regionOfInterest= struct('RegionType',Region,...
        'RegionID',ID,'Component',op.Component);

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

end
end

function Region=validateRegionType(self,op)
fn=fieldnames(op);
p={'Vertex','Edge','Face','Cell','Component'};
r=startsWith(p,fn,'ignoreCase',true);
idx=find(r);
if numel(idx)>2 && r(end) || numel(idx)>1 && ~r(end)
    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 validateComponent(self,component)
if self.IsTwoD && any(strfind(component,'z'))
    error(message('pde:linearize:InvalidCompVal'));
end
end

function invalid=outputIsInvalid(self,req)
invalid=false;
op=self.LinearizeOutputs;

if any(ismember({op.RegionType},req.RegionType))
    indx=ismember({op.RegionType},req.RegionType);
    IDs=[op.RegionID];
    if any(ismember(IDs(indx),req.RegionID))
        invalid=true;
    end
end
end
