function bl = findBodyLoad(self,varargin)
% findBodyLoad - Find body load assigned to a geometric region 
% BL = findBodyLoad(BLR,REGIONTYPE, REGIONID) returns the
% BodyLoadAssignment object, BL, that defines the body load assigned to the
% specified REGIONTYPE and REGIONID. Here BLR is the BodyLoads property of
% the structural model object, REGIONTYPE is 'Cell' for a 3-D model or
% 'Face' for a 2-D model, and REGIONID is the cell of face ID.
%
% % Example: Find body load for cells of a 3-D geometry.
%   structuralmodel = createpde('structural','static-solid');
%   importGeometry(structuralmodel,'Block.stl');
%   structuralBodyLoad(structuralmodel,'GravitationalAcceleration',[0,-9.8,0]);
%   findBodyLoad(structuralmodel.BodyLoads,'Cell',1)
%
%   See also pde.StructuralModel/structuralBodyLoad

% Copyright 2017-2021 The MathWorks, Inc.

narginchk(3,3);  
if isempty(self.ParentModel.Geometry)
    return;
end

parser = inputParser;
parser.addParameter('cell', [], @pde.EquationModel.isValidEntityID); 
parser.addParameter('face', [], @pde.EquationModel.isValidEntityID);  
parser.addParameter('label', [], @pde.EquationModel.ValidateLabelVal);
parser.parse(varargin{:});

if ~isempty(parser.Results.label)
    queryLabel = string(parser.Results.label);
    allExistingLabels={self.BodyLoadAssignments.Label};
    allExistingLabels(cellfun(@isempty,allExistingLabels))=[];
    allExistingLabels = string(allExistingLabels);
    
    if ~ismember(queryLabel,allExistingLabels)
    error(message('pde:linearize:NonExistingLabel',queryLabel));
    end
    numassigns = numel(self.BodyLoadAssignments);
    for j = numassigns:-1:1
       thisbc =  self.BodyLoadAssignments(j);       
       if queryLabel == string(thisbc.Label)       
          bl = thisbc;
          break
       end
    end
    return;
end

if ~ismember('face', parser.UsingDefaults)
    if ~self.ParentModel.IsTwoD
      error(message('pde:bodyLoadSpecification:noBoundaryBodyLoadQueryRegion')); 
   end    
    qregiontype = 'Face';
    qregionid = parser.Results.face;
    maxid = self.ParentModel.Geometry.NumFaces;       
else
    qregiontype = 'Cell';
    qregionid = parser.Results.cell;
    maxid = self.ParentModel.Geometry.NumCells;   
end


if any(qregionid > maxid)
    error(message('pde:bodyLoadSpecification:invalidRegionID'));
end

numassigns = numel(self.BodyLoadAssignments);
numqueries = numel(qregionid);

if ~bodyLoadSpecified(self,qregiontype,qregionid)
bl = [];
return;
end
       
bl = repmat(self.BodyLoadAssignments(1),size(qregionid));
for i = 1:numqueries
    rid = qregionid(i);
    %mpfound = false;
    for j = numassigns:-1:1
       thisca =  self.BodyLoadAssignments(j);   
       if ~bodyLoadSpecified(self,qregiontype,rid)
           continue
       end
       if strcmp(qregiontype,thisca.RegionType) && ismember(rid, thisca.RegionID)        
          bl(i) = thisca;
     %     mpfound = true;
          break
       end
    end
% It is not necessary to have body load assignment to all region.
%     if ~mpfound
%       error(message('pde:bodyLoadSpecification:entitiesWithoutCoefficients'));
%     end   
end

end

