function bc = structuralBC(self,varargin)
%structuralBC - Apply boundary conditions to a structural model
%   bc = structuralBC(structuralmodel,'REGIONTYPE',REGIONID,'Name',Value)
%   applies a structural boundary condition (BC) to a structural model. BC
%   is applied on a region of the boundary defined by REGIONTYPE and
%   REGIONID. Here, REGIONTYPE for a 3-D model is 'Face', 'Edge', or
%   'Vertex'. REGIONTYPE for a 2-D model is 'Edge' or 'Vertex'. REGIONID is
%   an integer specifying the ID of a face, edge, or vertex.
%
%   The name-value pairs for various types of structural BCs are shown in
%   the table. Note that only one valid pair of boundary condition can be
%   specified for a geometric boundary.
%
% -------------------------------------------------------------------------
% Name                      Description
% -------------------------------------------------------------------------
%
%'Displacement'             Enforced displacement. It can be specified as
%                           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 enforced displacement vector at
%                           the boundary coordinates provided by the
%                           solver.
%
%'XDisplacement'            x-component of enforced displacement. It can
%                           be a numeric scalar or a function handle to a
%                           custom function. This function must return a
%                           row vector with each column corresponding to
%                           value of x-component of displacement at the
%                           boundary coordinates provided by the solver.
%
%'YDisplacement'            y-component of enforced displacement. It can
%                           be a numeric scalar or a function handle to a
%                           custom function. This function must return a
%                           row vector with each column corresponding to
%                           value of y-component of displacement at the
%                           boundary coordinates provided by the solver.
%
%'ZDisplacement'            z-component of enforced displacement. It can
%                           be a numeric scalar or a function handle to a
%                           custom function. This function must return a
%                           row vector with each column corresponding to
%                           value of z-component of displacement at the
%                           boundary coordinates provided by the solver.
%'RDisplacement'            radial-component of enforced displacement in
%                           an axisymmetric analysis. It can be a numeric
%                           scalar or a function handle to a custom
%                           function. This function must return a row
%                           vector with each column corresponding to value
%                           of radial-component of displacement at the
%                           boundary coordinates provided by the solver.
%
%'Constraint'              Standard structural boundary constraints
%                           'fixed'      - Enforces zero-displacements on the boundary.
%                           'free'       - Stress free boundary.
%                           'roller'     - Prevents displacement perpendicular to the boundary.
%                           'symmetric'  - Enforces planar symmetric boundary condition.
%                           'multipoint' -  Define a multipoint constraint
%                                           using all degrees-of-freedom on
%                                           the specified region and a
%                                           reference point. The default
%                                           reference point is the
%                                           geometric center of all nodes
%                                           on all specified regions. You
%                                           can specify an arbitrary
%                                           reference location using
%                                           'Reference' in conjunction with
%                                           'multipoint'. Further, you can
%                                           restrict the region for
%                                           multipoint constraint to a
%                                           sphere/circle by using
%                                           'Radius'.
%
%                          The 'roller' and 'symmetric' constraints can be
%                          specified only on a face of a 3-D model or an
%                          edge of a 2-D model.
%
%'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 a component of displacement. Use
%                          'StartTime' in conjunction with
%                          'XDisplacement','YDisplacement' or
%                          'ZDisplacement' to define a transient prescribed
%                          displacement in a transient dynamics model.
%
%'EndTime'                 End time for a component of displacement. Use
%                          'EndTime' in conjunction with
%                          'XDisplacement','YDisplacement' or
%                          'ZDisplacement' to define a transient prescribed
%                          displacement in a transient dynamics model.
%
%'RiseTime'                Rise time for a component of displacement. Use
%                          'RiseTime' in conjunction with
%                          'XDisplacement','YDisplacement' or
%                          'ZDisplacement' to define a transient prescribed
%                          displacement in a transient dynamics model.
%
%'FallTime'                Fall time for a component of displacement. Use
%                          'FallTime' in conjunction with
%                          'XDisplacement','YDisplacement' or
%                          'ZDisplacement' to define a transient prescribed
%                          displacement in a transient dynamics model.
%
%
%'Frequency'               Frequency of a sinusoidal displacement, specified
%                          in radians per unit of time. Use 'Frequency' in
%                          conjunction with 'XDisplacement','YDisplacement'
%                          or 'ZDisplacement' to define a sinusoidal
%                          displacement in a transient dynamics model.
%
%'Phase'                   Phase of a sinusoidal displacement, specified in
%                          radians. Use 'Phase' in conjunction with
%                          'XDisplacement','YDisplacement' or
%                          'ZDisplacement'  to define a sinusoidal
%                          displacement in a transient dynamics model.
%
%'Reference'               Coordinates of the reference location for
%                          a multipoint constraint. Set 'Constraint' to
%                          'multipoint' and specify 'Reference' as a
%                          2-element vector for a 2-D model or 3-element
%                          vector for a 3-D model.
%
%'Radius'                  Radius around reference point location for
%                          multipoint constraint. Restricts the region for
%                          multipoint constraint to a sphere (3-D) or
%                          circle (2-D).
%
%'Label'                   Specify a customized label for linearizeInputs
%
%   Example 1: 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]);
%
%    Example 2: Define multipoint constraint interfaces and generate a
%               reduced-order model with six degrees of freedom at each
%               interface.
%
%       structuralmodel = createpde('structural','transient-solid');
%       structuralmodel.Geometry = multicuboid(0.1,0.01,0.01);
%       structuralProperties(structuralmodel,'YoungsModulus',70E9, ...
%                                      'PoissonsRatio',0.3, ...
%                                      'MassDensity',2700);
%       structuralBC(structuralmodel,'Edge',[4,6,9,10],'Constraint','multipoint');
%       structuralBC(structuralmodel,'Edge',[2,8,11,12],'Constraint','multipoint');
%       generateMesh(structuralmodel);
%       R = reduce(structuralmodel,'FrequencyRange',[-Inf,500000])
%
%    See also pde.StructuralModel, pde.StructuralModel/structuralProperties
%             pde.StructuralModel/structuralBodyLoad, 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('Displacement', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('XDisplacement', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('YDisplacement', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('ZDisplacement', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('RDisplacement', [],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Constraints', []);
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('Reference',[],@pde.EquationModel.isValidAttributeVal);
parser.addParameter('Radius',[],@pde.EquationModel.isValidNumericPositiveOnlyVal);
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 && ~isempty(parser.Results.Face)
    error(message('pde:StructuralModel:invalidBCRegionTypeFor2D'))
end

if self.IsTwoD &&  ~contains(self.AnalysisType,'axisymmetric') ...
        && ~isempty(parser.Results.ZDisplacement)
    error(message('pde:StructuralModel:noZDisplacementIn2D'));
end

if ~isempty(parser.Results.Reference) || ~isempty(parser.Results.Radius)
    if isempty(parser.Results.Constraints) || ~startsWith("multipoint",parser.Results.Constraints,'IgnoreCase',true)
        error(message('pde:structuralROM:radiusAndRefOnlyWithConstraint'));
    end
end

if contains(self.AnalysisType,'modal')
    bcOK = checkForNonzeroDisp(parser.Results);
    if ~bcOK
        error(message('pde:StructuralModel:nonZeroDisplacementModal'))
    end
end

if contains(self.AnalysisType,'axisymmetric')

    if ~isempty(parser.Results.XDisplacement) && ~isempty(parser.Results.RDisplacement)
        error(message('pde:StructuralModel:noXandRDisplacement'))
    end

    if ~isempty(parser.Results.YDisplacement) && ~isempty(parser.Results.ZDisplacement)
        error(message('pde:StructuralModel:noYandZDisplacement'))
    end

    % Maps radial displacement specified using RDisplacement to XDisplacement for computation. Passes
    % through if only XDisplacement is specified, treating it as the radial displacement.
    if isempty(parser.Results.XDisplacement) && ~isempty(parser.Results.RDisplacement)
        argsToPass = [{'XDisplacement',parser.Results.RDisplacement}, argsToPass];
    end

    % Maps axial displacement specified using ZDisplacement to YDisplacement for computation. Passes
    % through if only YDisplacement is specified, treating it as the axial displacement.
    if isempty(parser.Results.YDisplacement) && ~isempty(parser.Results.ZDisplacement)
        argsToPass = [{'YDisplacement',parser.Results.ZDisplacement}, argsToPass];
    end
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


function bcOK = checkForNonzeroDisp(userInput)

numericData =   isnumeric(userInput.XDisplacement) && ...
    isnumeric(userInput.YDisplacement) && ...
    isnumeric(userInput.ZDisplacement) && ...
    isnumeric(userInput.Displacement) ;
% Accounts for function handles and other unsupported data formats.
if ~numericData
    bcOK = false;
    return
end

% No non-zero displacement in modal analysis.
nonZeroDisplacement = any(userInput.XDisplacement(:)) ||...
    any(userInput.YDisplacement(:)) || ...
    any(userInput.ZDisplacement(:)) ||...
    any(userInput.Displacement(:));

if nonZeroDisplacement
    bcOK = false;
    return;
end

bcOK = true;
end