function FEMatricesOut = assembleFEMatrices(pdem, varargin)
% assembleFEMatrices - Assemble intermediate FE matrices 
%    FEM = assembleFEMatrices(model) assembles all FE matrices. 
%    FEM = assembleFEMatrices(model, matrixoption) assembles a subset
%          of FE matrices based on the specified matrixoption.
%    FEM = assembleFEMatrices(model, bcmethod) assembles all FE matrices
%          and returns the resulting system matrices after imposing BCs
%          using the approach specified in bcmethod.
%    FEM = assembleFEMatrices(__,state) assembles finite element matrices
%          using the time or solution at evaluation points specified in the
%          state structure array. The function uses the 'time' field of the
%          structure for time-dependent models and the solution field 'u'
%          for nonlinear models. The 'time' field must contain a scalar
%          real number and the 'u' field must contain a valid solution
%          matrix. The number of rows in the matrix must equal the number
%          of PDEs. The number of columns must equal the number of nodes.
%          Use any arguments from the previous syntaxes.
%
%    model is one of the following model objects:
%       PDEModel        - Generic PDE model created by using createpde(N),
%                         where N is the PDE system size.
%       ThermalModel    - Thermal analysis model, created using
%                         createpde('thermal').
%       StructuralModel - Structural analysis model, created using
%                         createpde('structural').
%
%    bcmethod option is one of the following:
%       'none'         - (default) 
%       'nullspace'    - elimination boundary condition method 
%       'stiff-spring' - penalty boundary condition method
%
%    matrixoption is one of the following:
%       'domain'       - Assemble all domain-related matrices. 
%       'boundary'     - Assemble all matrices related to geometry boundaries.
%       '<identifier>' - FE matrix identifiers: 'K','A', 'F', 'Q', 'G',
%                       'H','R', 'M', or any combination of these, such as,
%                       'KF','FG','MKAR'.
%
%    assembleFEMatrices assembles the intermediate FE matrices for generic
%    PDE equations of the following form:
%  
%              m*d^2u/dt^2 + d*du/dt - div(c*grad(u)) + a*u = f
%
%    ThermalModel maps to this equation with m = 0 and a = 0. Here, the
%    thermal conductivity maps to the c-coefficient. The product of the
%    mass density and the specific heat maps to the d-coefficient. The
%    internal heat source maps to the f-coefficient.
%
%    StructuralModel maps to this equation with a = 0. The Young's modulus
%    and Poisson's ratio map to the c-coefficient. The mass density maps to
%    the m-coefficient. The body loads map to the f-coefficient. When you
%    specify damping model by using Rayleigh damping parameters Alpha and
%    Beta, the discretized damping matrix C is computed by using the mass
%    matrix M and the stiffness matrix K as C = AlphaM+BetaK.
%
%    The intermediate FE matrices can subsequently be used to compute a 
%    solution to the PDE using linear algebra methods.
%
%    Note: If both m- and d-coefficients are nonzero, then the
%    d-coefficient must be a matrix of the same size the discretized mass
%    matrix M, or the stiffness matrix Kc mentioned below.
%
%    When bcmethod is set to 'none', the output argument FEM is a
%    structural array with the following fields:
%           K: Stiffness matrix 
%           A: Absorption/Reaction matrix
%           F: Volumetric load vector  
%           Q: Neumann boundary condition matrix  
%           G: Neumann boundary condition vector 
%           H: Dirichlet boundary condition matrix  
%           R: Dirichlet boundary condition vector 
%           M: Mass matrix  
%
%    When bcmethod is set to 'nullspace', the output argument FEM is a
%    structural array with the following fields:
%           Kc: Stiffness matrix  
%           Fc: Load vector
%            B: Dirichlet nullspace matrix
%           ud: Dirichlet vector
%            M: Mass matrix  
%
%    When bcmethod is set to 'stiff-spring', the output argument FEM is a
%    structural array with the following fields:
%           Ks: Stiffness matrix  
%           Fs: Load vector
%            M: Mass matrix  
%
% See also pde.PDEModel, pde.StructuralModel, pde.ThermalModel
%

% Copyright 2015-2020 The MathWorks, Inc.

if nargin > 1
    [varargin{:}] = convertStringsToChars(varargin{:});
end

narginchk(1, 3);
BCEnforcementOption = 'none';
matrix = [];
state =[];

if ~isa(pdem, 'pde.EquationModel')    
    error(message('pde:assembleFEMatrices:invalidModel'));
end

if nargin == 1
    FEMatricesOut = assembleFEMatricesInternal(pdem,BCEnforcementOption,matrix,state);
    return
end

if nargin == 2
   if  ~(isstruct(varargin{1}) || isstring(varargin{1}) || ischar(varargin{1}))
       error(message('pde:assembleFEMatrices:invalidSecondArg'));
   end
end

if nargin == 3
   if  ~(isstring(varargin{1}) || ischar(varargin{1})) 
       error(message('pde:assembleFEMatrices:invalidSecondArg'));
   end
   
   if ~isstruct(varargin{2})
       error(message('pde:assembleFEMatrices:invalidThirdArg'));
   end
   state = varargin{2};
end

if isstring(varargin{1}) || ischar(varargin{1})
    nc = numel(varargin{1});
    if (startsWith("none",varargin{1},'IgnoreCase',true) ||...
            startsWith("nullspace",varargin{1},'IgnoreCase',true) ||...
            startsWith("stiff-spring",varargin{1},'IgnoreCase',true))
        if nc == 1 && strncmpi(varargin{1},'n',nc)
            error(message('pde:assembleFEMatrices:ambiguousBcmethodType'));
        end
        BCEnforcementOption= varargin{1};
    elseif startsWith("domain",varargin{1},'IgnoreCase',true)
        matrix = 'MKAF';
    elseif startsWith("boundary",varargin{1},'IgnoreCase',true)
        matrix = 'HRGQ';
    else
        inputStr = sort(unique(upper(char(varargin{1}))));
        validMatIDs = intersect('AFGHKMQRT',inputStr);
        if isequal(inputStr,validMatIDs)
            matrix = upper(varargin{1});
        else
            error(message('pde:assembleFEMatrices:invalidSecondArg'));
        end
    end
else
    state = varargin{1};
end

FEMatricesOut = assembleFEMatricesInternal(pdem,BCEnforcementOption,matrix,state);

end

