function bc = structuralBoundaryLoad(self,varargin)
%structuralBoundaryLoad - Apply structural loads on the geometric boundary.
%   sbl = structuralBoundaryLoad(structuralmodel,'REGIONTYPE',REGIONID,'Name',Value)
%   applies a structural boundary load to a structural model. Load is
%   applied on a region of the boundary defined by REGIONTYPE and REGIONID.
%  When specifying a surface traction, pressure, or translational
%  stiffness, REGIONTYPE is 'Face' for a 3-D model and 'Edge' for a 2-D
%  model. When specifying force, REGIONTYPE must always be 'Vertex'.
%  REGIONID is an integer specifying the ID of a face, edge, or vertex.
%
%  The name-value pairs for various types of structural boundary loads are:
%
% -------------------------------------------------------------------------
% Name			      	    Description
% -------------------------------------------------------------------------
%'SurfaceTraction'          Normal and tangential distributed forces on a
%                           boundary, resolved along global Cartesian
%                           coordinates system. Surface traction must be a
%                           numeric vector or a function handle to a custom
%                           function. A numeric vector must contain 2
%                           elements for a 2-D model and 3 elements for a
%                           3-D model. The function must return a 2- or
%                           3-row matrix, for a 2-D or 3-D model,
%                           respectively. Each column of the matrix
%                           corresponds to the surface traction vector at
%                           the boundary coordinates provided by the
%                           solver.
%
%'Pressure'                 Pressure acting normal to the boundary. It can 
%                           be specified as a numeric scalar or a function
%                           handle to a custom function. The function must
%                           return a row vector with each column
%                           corresponding to value of pressure at the
%                           boundary coordinates provided by the solver.
%                           Pressure acts in the direction normal to the
%                           boundary. A positive value of pressure acts
%                           into the boundary (for example, compression),
%                           while a negative value of pressure acts away
%                           from the boundary (for example, suction).
%
%'Force'                    Concentrated force acting at a point. Force 
%                           must be a numeric vector or a function handle
%                           to a custom function. A numeric vector must
%                           contain 2 elements for a 2-D model and 3
%                           elements for a 3-D model. The function must
%                           return a 2- or 3-element vector, for a 2-D or
%                           3-D model, respectively. 
%
%'TranslationalStiffness'   Distributed spring stiffness for each 
%                           translational direction used to model elastic
%                           foundation. It can be a numeric vector or a
%                           function handle to a custom function. A numeric
%                           vector must contain 2 elements for a 2-D model
%                           and 3 elements for a 3-D model. The custom
%                           function must return a 2- or 3-row matrix, for
%                           a 2-D or 3-D model, respectively. Each column
%                           of this matrix corresponds to the stiffness
%                           vector at the boundary coordinates provided by
%                           the solver.
%
%'Vectorized'              If set to 'on', indicates that the custom
%                          function can accept input arguments at multiple
%                          locations and returns boundary condition
%                          matrices defined at these multiple locations.
%                          The default is 'off'. This setting is useful
%                          when one of the other parameters are defined as
%                          a function handle.
%
%'StartTime'               Start time for pressure load. Use 'StartTime' in
%                          conjunction with 'Pressure' or 'Force to define
%                          a transient load in a transient dynamics model.
%
%'EndTime'                 End time for pressure load. Use 'EndTime' in
%                          conjunction with 'Pressure' or 'Force to define
%                          a transient load in a transient dynamics model.
%
%'RiseTime'                Rise time for pressure load. Use 'RiseTime' in
%                          conjunction with 'Pressure' or 'Force to define
%                          a transient load in a transient dynamics model.
%
%'FallTime'                Fall time for pressure load. Use 'FallTime' in
%                          conjunction with 'Pressure' or 'Force to define
%                          a transient load in a transient dynamics model.
%
%'Frequency'               Frequency of a sinusoidal displacement, specified 
%                          in radians per unit of time. Use 'Frequency' in
%                          conjunction with 'Pressure' or 'Force' to define
%                          a sinusoidal amplitude in a transient dynamics
%                          model.
%
%'Phase'                   Phase of a sinusoidal displacement, specified in
%                          radians. Use 'Phase' in conjunction with
%                          'Pressure' or 'Force' to define a sinusoidal
%                           amplitude in a transient dynamics model.
%
%'Label'                   Specify a customized label for linearizeInputs
%
%   The default boundary load is stress-free BC.
%
%   Example: Apply fixed boundary and traction on two ends of a bimetallic cable.
%       structuralmodel = createpde('structural','static-solid');
%       structuralmodel.Geometry = multicylinder([0.01,0.015],0.05);
%       pdegplot(structuralmodel,'FaceLabels','on')
%       structuralProperties(structuralmodel,'Cell',1,'YoungsModulus',110E9, ...
%                                                  'PoissonsRatio',0.28);
%       structuralProperties(structuralmodel,'Cell',2,'YoungsModulus',210E9, ...
%                                                  'PoissonsRatio',0.3);
%
%       structuralBC(structuralmodel,'Face',[1,4],'Constraint','fixed');
%       structuralBoundaryLoad(structuralmodel,'Face',[2,5],'SurfaceTraction',[0;0;100]);
%
%    See also pde.StructuralModel, pde.StructuralModel/structuralProperties
%             pde.StructuralModel/structuralBodyLoad,pde.StructuralModel/structuralBC,
%             pde.StructuralModel/solve, pde.DiscreteGeometry/addVertex

%    Copyright 2017-2021 The MathWorks, Inc.

narginchk(5,15); 
nargoutchk(0,1);
    
if isempty(self.Geometry)
    error(message('pde:StructuralModel:structuralModelHasNoGeom'));
end

argsToPass = varargin;

parser = inputParser;
parser.addParameter('Face', [],@pde.EquationModel.isValidEntityID);
parser.addParameter('Edge', [],@pde.EquationModel.isValidEntityID);
parser.addParameter('Vertex', [],@pde.EquationModel.isValidEntityID);
parser.addParameter('SurfaceTraction', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Pressure', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('TranslationalStiffness',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Force',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('StartTime',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('EndTime',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('RiseTime',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('FallTime',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Frequency',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Phase',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Vectorized', []);
parser.addParameter('Label', [],@pde.EquationModel.ValidateLabelVal);
parser.parse(varargin{:});

validateRegionID(self,parser.Results.Vertex, parser.Results.Edge, parser.Results.Face);

if self.IsTwoD
    if ~isempty(parser.Results.Face)
        error(message('pde:StructuralModel:invalidBCRegionTypeFor2D'))
    end
else
    if ~isempty(parser.Results.Edge) % No boundary load supported on an edge of 3-D.
        error(message('pde:StructuralModel:invalidCombinationRegionTypeAndLoad'))
    end
end

if ~isempty(parser.Results.Force) && isempty(parser.Results.Vertex)
    error(message('pde:StructuralModel:forceMustbeOnVertex'));
end

if ~isempty(parser.Results.Vertex) &&( ~isempty(parser.Results.SurfaceTraction) || ...
                                              ~isempty(parser.Results.Pressure) || ...
                                    ~isempty(parser.Results.TranslationalStiffness) )
    if self.IsTwoD
        error(message('pde:StructuralModel:invalidBCRegionTypeFor2D'));
    else
        error(message('pde:StructuralModel:invalidCombinationRegionTypeAndLoad'))
    end
end

if contains(self.AnalysisType,'modal') && ...
        (~isempty(parser.Results.SurfaceTraction) || ~isempty(parser.Results.Pressure))
    warning(message('pde:StructuralModel:boundaryLoadsIgnoredInModal'))
end

if isempty(self.BoundaryConditions)
    BCcont = pde.StructuralBCRecords(self);
    bc = pde.StructuralBC(BCcont,argsToPass{:});
    BCcont.StructuralBCAssignments = bc;
    self.BoundaryConditions = BCcont;
else
    bc = pde.StructuralBC(self.BoundaryConditions,argsToPass{:});
    self.BoundaryConditions.StructuralBCAssignments(end+1) = bc;
end
end
