classdef (Sealed) StructuralICRecords  < pde.PDECompactDisplay 
%StructuralICRecords - Records of assignment of structural initial conditions
% The StructuralICRecords holds a record of all the assignments of
% structural initial condition (IC), displacement and velocity, across the
% geometric domain. The records can contain a single assignment that
% represents a single IC throughout the domain or multiple assignments,
% where the ICs differ across geometric regions. This object cannot be
% created directly, it can only be created by the StructuralModel class.
% 
% StructuralICRecords methods:
%  findStructuralIC - Find IC assignment for a geometric region
% 
% StructuralICRecords properties:
%  StructuralICAssignments - Vector of IC assignments
% 
% See also pde.StructuralModel, pde.StructuralModel/structuralIC

% Copyright 2017-2022 The MathWorks, Inc.

properties
        
% StructuralICAssignments - Vector of structural initial condition assignments
%    A vector containing the IC assignments to the geometry or geometric
%    boundaries. Each entry in the vector is either a
%    GeometricStructuralICs or NodalStructuralICs object. This object
%    defines the ICs that should be set on a designated region or region
%    boundary. ICs and host regions are set using the structuralIC method
%    of the StructuralModel.
    StructuralICAssignments;
  
end

methods  
   ca = findStructuralIC(self, varargin)   
end

methods(Hidden=true, Access={?pde.EquationModel})
    function obj=StructuralICRecords(spdem) 
       obj.StructuralICAssignments = [];
       obj.ParentPdemodel = spdem;
    end         
    
    % Check that the domain or each subdomain has initial conditions set if
    % they are required. Also perform a basic sanity check of the initial 
    % conditions in the event that they were edited since assigned.
  
    function performSolverPrechecks(self)
   
        if (self.ParentPdemodel.IsTimeDependent)
            if ~self.icValueDefined() && ~self.icDerivativeValueDefined
                error(message('pde:structuralInitialConditions:transientNeedsIcValue')); 
            end
        end
        
        icassigntype = classifyICAssignment(self);
        if(strcmp(icassigntype, 'GEOMETRIC'))
            performGeometricPrechecks(self);
        elseif(strcmp(icassigntype, 'NODAL'))
            performNodalPrechecks(self);
        else
            performHybridPrechecks(self);
        end
       
               
    end
    
    % Check Geometric IC assignments - no nodal assignments
    % This check ensures all Faces(2-D) and Cells(3-D) have an IC.
    % It then checks all ICs to ensure entities on Edges/vertices
    % are valid.
    function performGeometricPrechecks(self)
        systemsize = self.ParentPdemodel.PDESystemSize;
        nv = self.ParentPdemodel.Geometry.NumVertices;
        ne = self.ParentPdemodel.Geometry.NumEdges;
        nf = self.ParentPdemodel.Geometry.NumFaces;
        nc = self.ParentPdemodel.Geometry.NumCells;
        numic = numel(self.StructuralICAssignments);
        
        if self.ParentPdemodel.IsTwoD
            %Check for valid global assignment
            for i = 1:nf
                thisic = self.findStructuralIC('Face',i);
                thisic.performSolverPrecheck(systemsize, nv, ne, nf, nc);
            end
            %Check for valid local assignments
            for idx = 1:numic
                thisic = self.StructuralICAssignments(idx);
                if ( strcmp(thisic.RegionType, 'Edge') || ...
                        strcmp(thisic.RegionType, 'Vertex') )
                    thisic.performSolverPrecheck(systemsize, nv, ne, nf, nc);
                end
            end
        else
            %Check for valid global assignment
            for i = 1:nc
                thisic = self.findStructuralIC('Cell',i);
                thisic.performSolverPrecheck(systemsize, nv, ne, nf, nc);
            end
            %Check for valid local assignments
            for idx = 1:numic
                thisic = self.StructuralICAssignments(idx);
                if (strcmp(thisic.RegionType, 'Face') || ...
                        strcmp(thisic.RegionType, 'Edge')|| ...
                        strcmp(thisic.RegionType, 'Vertex'))
                    thisic.performSolverPrecheck(systemsize, nv, ne, nf, nc);
                end
            end
            
            
        end
    end
   
    % Check Nodal IC assignments - no geometric assignments
    function performNodalPrechecks(self)
        % The last assignment must be Nodal        
        systemsize = self.ParentPdemodel.PDESystemSize;
        numnodes = size(self.ParentPdemodel.Mesh.Nodes,2);
        thisic = self.StructuralICAssignments(end);
        thisic.performSolverPrecheck(systemsize, numnodes);    
    end
    
    % performHybridPrechecks, find the most recent nodal assignment
    % and sanity check that assignment, then check all subsequent
    % geometric assignments
    function performHybridPrechecks(self)        
        numic = numel(self.StructuralICAssignments);                
        idx = findIndexOfNodalIC(self);
        systemsize = self.ParentPdemodel.PDESystemSize;
        numnodes = size(self.ParentPdemodel.Mesh.Nodes,2);
        thisic = self.StructuralICAssignments(idx);
        thisic.performSolverPrecheck(systemsize, numnodes);    
        nv = self.ParentPdemodel.Geometry.NumVertices;
        ne = self.ParentPdemodel.Geometry.NumEdges;    
        nf = self.ParentPdemodel.Geometry.NumFaces;
        nc = self.ParentPdemodel.Geometry.NumCells; 
        idx = idx+1;
        for i = idx:numic
            thisic = self.StructuralICAssignments(i);
            thisic.performSolverPrecheck(systemsize, nv, ne, nf, nc);
        end 
    end
    
    % Returns a string 'GEOMETRIC', 'NODAL' or 'MIXED'
    % to classify the nature of the IC assignment
    function icassigntype = classifyICAssignment(self)
        thisic = self.StructuralICAssignments(end);
        if isa(thisic, 'pde.NodalStructuralICs')
            icassigntype = 'NODAL';
            return;
        end
        numic = numel(self.StructuralICAssignments);
        for i = 1:numic
            thisic = self.StructuralICAssignments(i);
            if isa(thisic, 'pde.NodalStructuralICs')
                icassigntype = 'MIXED';
                return;
            end 
        end           
        icassigntype = 'GEOMETRIC';
    end
    
    function [u0, ut0] = packInitialConditions(self)
       icassigntype = classifyICAssignment(self);
        if(strcmp(icassigntype, 'GEOMETRIC'))
            [u0, ut0] = packGeometricInitialConditions(self);
        elseif(strcmp(icassigntype, 'NODAL'))
            [u0, ut0] = packNodalInitialConditions(self);
        else
            [u0, ut0] = packMixedInitialConditions(self);
        end
    end
end
  
methods(Hidden=true, Access = private) 
% Declaration
tf = activeGeometricIC(self, varargin);

function entcnt = CountEntities(self, entitytype)
    switch entitytype
        case 'Cell'
            entcnt = self.ParentPdemodel.Geometry.NumCells; 
        case 'Face'
            entcnt = self.ParentPdemodel.Geometry.NumFaces; 
        case 'Edge'
            entcnt = self.ParentPdemodel.Geometry.NumEdges; 
        case 'Vertex'
            entcnt = self.ParentPdemodel.Geometry.NumVertices;                            
    end        
end

% Finds the index of the last nodal-IC.
function idx = findIndexOfNodalIC(self)
    numic = numel(self.StructuralICAssignments);                
    for idx = numic:-1:1
        thisic = self.StructuralICAssignments(idx);
        if isa(thisic, 'pde.NodalStructuralICs')
          return;
        end 
    end     
end


function [u0, ut0] = packGeometricInitialConditions(self)
    nodes = self.ParentPdemodel.Mesh.Nodes;
    syssize = self.ParentPdemodel.PDESystemSize;
    numnodes = size(nodes,2);
    u0 = ones(numnodes,syssize);
    ut0 = ones(numnodes,syssize);

    % Fill the array of IC values         
    [u0, ut0] = packEntityInitialConditions(self, 'Cell', u0, ut0);
    [u0, ut0] = packEntityInitialConditions(self, 'Face', u0, ut0);
    [u0, ut0] = packEntityInitialConditions(self, 'Edge', u0, ut0);
    [u0, ut0] = packEntityInitialConditions(self, 'Vertex', u0, ut0);                           
    u0 = u0(:);
    ut0 = ut0(:);
end

function [u0, ut0] = packMixedInitialConditions(self)
    [u0, ut0] = packNodalInitialConditions(self);
    u0 = reshape(u0,[],self.ParentPdemodel.PDESystemSize);
    ut0 = reshape(ut0,[],self.ParentPdemodel.PDESystemSize);
    [u0, ut0] = packEntityInitialConditions(self, 'Cell', u0, ut0);
    [u0, ut0] = packEntityInitialConditions(self, 'Face', u0, ut0);
    [u0, ut0] = packEntityInitialConditions(self, 'Edge', u0, ut0);
    [u0, ut0] = packEntityInitialConditions(self, 'Vertex', u0, ut0);                           
    u0 = u0(:);
    ut0 = ut0(:);    
end

function [u0, ut0] = packNodalInitialConditions(self)
   idx = findIndexOfNodalIC(self);
   thisic = self.StructuralICAssignments(idx);
   u0 = thisic.InitialDisplacement;
   ut0 = thisic.InitialVelocity;
   u0 = u0(:);
   ut0 = ut0(:);    
end


function [u0, ut0] = packEntityInitialConditions(self, entitytype, u0, ut0)
    entitycount = CountEntities(self, entitytype);
    if entitycount == 0
        return
    end
    is3D = ~self.ParentPdemodel.IsTwoD;
    msh = self.ParentPdemodel.Mesh;
    ma = msh.MeshAssociation;
    hasic = self.activeGeometricIC(entitytype,1:entitycount);
    icentities = find(hasic);
    numicentities = numel(icentities);
    for i = 1:numicentities
        thisentid = icentities(i);
        thisic = self.findStructuralIC(entitytype,thisentid);
        fnid = ma.getNodes(entitytype,thisentid);
        if thisic.numericIcValue()
            if thisic.defaultDisplacement()
                u0Input = zeros(1, self.ParentPdemodel.PDESystemSize);
            else
                u0Input = reshape(thisic.InitialDisplacement,1,self.ParentPdemodel.PDESystemSize);
            end
            u0(fnid,:) = bsxfun(@times,ones(size(u0(fnid,:))),(u0Input));
        else
            fhdl = thisic.InitialDisplacement;
            locations.x = msh.Nodes(1,fnid);
            locations.y = msh.Nodes(2,fnid);
            if is3D
                locations.z = msh.Nodes(3,fnid);
            end
            u0(fnid,:) = (fhdl(locations))';
        end
        
        if thisic.numericIcDerivativeValue()
            if thisic.defaultVelocity()
                ut0Input = zeros(1, self.ParentPdemodel.PDESystemSize);
            else
                ut0Input = reshape(thisic.InitialVelocity,1,self.ParentPdemodel.PDESystemSize);
            end
            ut0(fnid,:) = bsxfun(@times,ones(size(ut0(fnid,:))),(ut0Input));
        else
            fhdl = thisic.InitialVelocity;
            locations.x = msh.Nodes(1,fnid);
            locations.y = msh.Nodes(2,fnid);
            if is3D
                locations.z = msh.Nodes(3,fnid);
            end
            ut0(fnid,:) = (fhdl(locations))';
        end
        
    end
end

end


methods(Hidden=true)     
    function tf = icValueDefined(self)
       tf = false;
       numassigns = numel(self.StructuralICAssignments);
       if numassigns
         ic = self.StructuralICAssignments(end);
         tf = ic.hasICValue();
       end   
    end
    
    function [tf,ic] = icDerivativeValueDefined(self)
       tf = false;
       numassigns = numel(self.StructuralICAssignments);
       if numassigns
         ic = self.StructuralICAssignments(end);
         tf = ic.hasICDerivativeValue();
       end   
    end
end

    
methods(Hidden=true, Access={?pde.GeometricStructuralICs, ?pde.NodalStructuralICs, ?pde.PDEInitialConditions})
    function delistStructuralICAssignments(self, ictodelete)
        numic = numel(self.StructuralICAssignments);
        for i = 1:numic
            thisic = self.StructuralICAssignments(i);
            if thisic == ictodelete
                self.StructuralICAssignments(i) = [];
                break
            end
        end   
        numic = numel(self.StructuralICAssignments);
        if numic == 0 && isvalid(self.ParentPdemodel)
           self.ParentPdemodel.delistInitialConditions(); 
        end
    end 
end


methods(Hidden=true)
    function delete(self,~)
        
        numic = numel(self.StructuralICAssignments);
        for i = numic:-1:1
            if isvalid(self.StructuralICAssignments(i))
                delete(self.StructuralICAssignments(i));
            end
        end
        
        if isvalid(self.ParentPdemodel)
            self.ParentPdemodel.delistInitialConditions();
        end
    end
end


properties (Hidden = true, SetAccess='private')
    ParentPdemodel;
end  
  
end
