function ic = structuralIC(self,varargin)
% structuralIC - Set initial displacement and velocity
% 
% structuralIC(structuralmodel,'Displacement',U0,'Velocity,V0) sets the
% initial displacement U0 and velocity V0 for the entire domain of a
% structural dynamics model. U0 and V0 are vectors of 2 or 3 components of
% displacement and velocity, respectively, or function handles to custom
% functions. These functions must return a 2- or 3-row matrix for a 2-D or
% 3-D model, respectively. Each column of this matrix corresponds to
% initial displacement or velocity at the coordinates provided by the
% solver.
% 
% structuralIC(__,REGIONTYPE,REGIONID) sets the initial conditions for the
% subdomain specified by using REGIONTYPE and REGIONID for a structural
% dynamics model. Here, REGIONTYPE can be 'Cell', 'Face', 'Edge', or
% 'Vertex' for a 3-D model, and 'Face', 'Edge', or 'Vertex' for a 2-D
% model. REGIONID is the ID of the geometric region.
%  
% structuralIC(structuralmodel,SRESULTS) sets the initial conditions using
% the solution SRESULTS obtained by solving a structural model. If SRESULTS
% is a structural dynamics model solution, then structuralIC uses the
% displacement and velocity corresponding to the last time-step.
%  
% structuralIC(structuralmodel,SRESULTS,iT) sets initial conditions using
% the displacement and velocity corresponding to the time-step iT of the
% solution SRESULTS.
% 
% % Example 1: Impulse as an Initial Velocity 
% 
%     % Create a model, and assign a geometry to the model.
%       model = createpde('structural','transient-solid');
%       gm = multicuboid(0.06,0.005,0.01);
%       model.Geometry = gm;
%
%     % Specify zero initial displacement and velocity throughout the domain.
%       structuralIC(model,'Displacement',[0;0;0],'Velocity',[0;0;0]);
%
%     % Update the initial velocity on face 2 to model impulsive excitation.
%       structuralIC(model,'Face',2,'Velocity',[0;60;0]);
% 
% % Example 2: Static Solution as Initial Conditions for Transient Dynamics
%     % Create a static model and assign geometry.
%       staticmodel = createpde('structural','static-solid');
%       gm = multicuboid(0.06,0.005,0.01);
%       staticmodel.Geometry = gm;
%       structuralProperties(staticmodel,'YoungsModulus',210E9, ...
%                                        'PoissonsRatio',0.3,...
%                                        'MassDensity',7800);
%
%     % Apply boundary condition and static load.
%       structuralBC(staticmodel,'Face',5,'Constraint','fixed');
%       structuralBoundaryLoad(staticmodel,'Face',3,'SurfaceTraction',[0;1E6;0]);
%       generateMesh(staticmodel,'Hmax',0.02);
%       Rstatic = solve(staticmodel);
%       
%      % Create a dynamic model, and assign geometry.
%       dynamicmodel = createpde('structural','transient-solid');
%       dynamicmodel.Geometry = gm;
%
%     % Apply boundary condition.
%       structuralBC(dynamicmodel,'Face',5,'Constraint','fixed');
%
%     % Specify initial condition using the static solution
%       generateMesh(dynamicmodel,'Hmax',0.02);
%       structuralIC(dynamicmodel,Rstatic);
%
%    See also pde.StructuralModel, pde.StructuralModel/structuralProperties
%             pde.StructuralModel/structuralBodyLoad, pde.StructuralModel/solve

%    Copyright 2017 The MathWorks, Inc.

narginchk(2,7); 
nargoutchk(0,1);

if isempty(self.Geometry)
    error(message('pde:StructuralModel:structuralModelHasNoGeom'));
end

if ~contains(self.AnalysisType,'transient') 
    error(message('pde:StructuralModel:noICOnNoNTrasnientModel'));
end

isnodal = false;

if (nargin == 2) || (nargin == 3) % Could be nodal
    if isa(varargin{1},'pde.StructuralResults')
        if ~(isa(varargin{1}, 'pde.StaticStructuralResults')|| isa(varargin{1}, 'pde.TransientStructuralResults'))
            error(message('pde:structuralInitialConditions:invalidResultsObject'));
        else
            argsToPass = varargin;
            isnodal = true;
        end
    else 
        if ~(ischar(varargin{1}) || isstring(varargin{1}))
            error(message('pde:StructuralModel:invalidIC'));
        end
        argsToPass = ParseGeometricICs(self, varargin{:});
    end
else
    argsToPass = ParseGeometricICs(self, varargin{:});
end
    
if isempty(self.InitialConditions)
    iccont = pde.StructuralICRecords(self);
    argsToPass = {iccont, argsToPass{1:end}};
    if isnodal
        ic = pde.NodalStructuralICs(argsToPass{:});
    else
        ic = pde.GeometricStructuralICs(argsToPass{:});
    end
    iccont.StructuralICAssignments = ic;
    self.InitialConditions = iccont;
    
else
    argsToPass = {self.InitialConditions, argsToPass{1:end}};
    if isnodal
        ic = pde.NodalStructuralICs(argsToPass{:});
    else
        ic = pde.GeometricStructuralICs(argsToPass{:});
    end
    self.InitialConditions.StructuralICAssignments(end+1) = ic;
end

end


function argsToPass = ParseGeometricICs(self, varargin)
parser = inputParser;
parser.addParameter('Cell', [], @pde.EquationModel.isValidEntityID);
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('Velocity', [], @pde.EquationModel.isValidAttributeVal);
parser.parse(varargin{:});

if (isempty(parser.Results.Cell) && isempty(parser.Results.Face) && ...
           isempty(parser.Results.Edge) && isempty(parser.Results.Vertex)) 
    % Global assignment
    if self.IsTwoD
        argsToPass = [varargin, {'Face', 1:self.Geometry.NumFaces}];
    else
        argsToPass = [varargin, {'Cell', 1:self.Geometry.NumCells}];
    end
else   % Local assignment
    if ~isempty(parser.Results.Cell)
        self.isValidEntityID(parser.Results.Cell);
        if any(parser.Results.Cell > self.Geometry.NumCells)
            error(message('pde:pdeModel:invalidCellIndex'));
        end
    end
    if ~isempty(parser.Results.Face)
        self.isValidEntityID(parser.Results.Face);
        if any(parser.Results.Face > self.Geometry.NumFaces)
            error(message('pde:pdeModel:invalidFaceIndex'));
        end
    end
    if ~isempty(parser.Results.Edge)
        self.isValidEntityID(parser.Results.Edge);
        if any(parser.Results.Edge > self.Geometry.NumEdges)
            error(message('pde:pdeModel:invalidEdgeIndex'));
        end
    end
    if ~isempty(parser.Results.Vertex)
        self.isValidEntityID(parser.Results.Vertex);
        if any(parser.Results.Vertex > self.Geometry.NumVertices)
            error(message('pde:pdeModel:invalidVertexIndex'));
        end
    end
    argsToPass = varargin;
end

end
