classdef (Sealed) StructuralMaterialAssignmentRecords < pde.PDECompactDisplay
    % StructuralMaterialAssignmentRecords - Assignments of structural material properties
    %  StructuralMaterialAssignmentRecords holds a record of all the
    %  assignments of structural material properties across the geometric
    %  domain. This may be a single assignment throughout the domain or
    %  multiple assignments that vary across various subdomains. This
    %  object cannot be created directly; it can only be created by the
    %  StructuralModel class.
    %
    % StructuralMaterialAssignmentRecords methods:
    %    findStructuralProperties - Find structural properties for a geometric region
    %
    % StructuralMaterialAssignmentRecords properties:
    %    MaterialAssignments - Vector of material property assignments
    %
    % See also pde.StructuralModel, pde.StructuralModel/structuralProperties
    
    % Copyright 2017-2022 The MathWorks, Inc.
    
    properties
        % MaterialAssignments - Vector of material property assignments
        %    A vector containing the material property assignments to the
        %    geometric domain or subdomains. Each entry in the vector is a
        %    pde.StructuralMaterialAssignment object. This object defines the
        %    material properties on a designated domain or subdomain of a
        %    model. Material properties and host domains are defined using the
        %    structuralProperties method of the StructuralModel object.
        MaterialAssignments;
    end
    
    methods
        % Method declaration
        mp = findStructuralProperties(self, varargin)
        tf = materialPropertiesSpecified(self, varargin)
        tf = globalMaterialProperties(self, varargin)
    end
    
    methods (Hidden=true, Access={?pde.EquationModel, ?pde.StructuralResults})
        function self = StructuralMaterialAssignmentRecords(spdem)
            self.MaterialAssignments = [];
            self.ParentModel = spdem;
        end
        
        function tf = allMaterialPropertiesNumeric(self)
            tf = true;
            if ~self.materialPropertiesSpecified()
                error(message('pde:materialPropertySpecification:subdomainsWithoutMtlProps'));
            end
            
            for i = 1:self.nRegion
                mp = self.findStructuralProperties(self.qregiontype,i);
                if ~numericMaterialProperties(mp)
                    tf = false;
                    return
                end
            end
        end
        
        function performSolverPrechecks(self)
            if ~self.materialPropertiesSpecified()
                error(message('pde:materialPropertySpecification:subdomainsWithoutStrMtlProps'));
            end
            if ~consistentMaterialPropertyDefinition(self)
                error(message('pde:materialPropertySpecification:inconsistentMtlProps'));
            end
            
            syssize = self.ParentModel.PDESystemSize;
            nf = self.ParentModel.Geometry.NumFaces;
            nc = self.ParentModel.Geometry.NumCells;
            
            for i = 1:self.nRegion
                ca = self.findStructuralProperties(self.qregiontype,i);
                ca.performSolverPrecheck(syssize, nf, nc);
            end
        end
        
        function mtlstruct = packProperties(self)
            
            numdim = self.ParentModel.Dim;
            [loc, state] = getSampleLocState(self);
            ElasticityType = self.ParentModel.ElasticityType;
            msh = self.ParentModel.Mesh;
            ma = msh.MeshAssociation;
            numelems = size(msh.Elements,2);
            mtlstruct.ElementsInSubdomain = (1:numelems)';
            complexcoef = false(5,1);
            
            isBodyLoadGlobal = true;
            if ~isempty(self.ParentModel.BodyLoads)
                isBodyLoadGlobal = self.ParentModel.BodyLoads.globalBodyLoad();
            end
            
            numTrueSubdomain = self.nRegion;
            eIDtocellID = ma.SubdomainAssociativity;
            if self.globalMaterialProperties() && isBodyLoadGlobal
                numTrueSubdomain = 1;
                eIDtocellID = [1;numelems];
            end
            
            mtlstruct.NumElementsInSubdomain = zeros(numTrueSubdomain,1);
            
            for i = 1:numTrueSubdomain
                numcelems = eIDtocellID(2,i)-eIDtocellID(1,i)+1;
                mtlstruct.NumElementsInSubdomain(i) = numcelems;
                ca = self.findStructuralProperties(self.qregiontype,i);
                mtlcoefs.m{i} = pde.internal.FESolverUtilities.elasticityM(ca.MassDensity,ElasticityType);
                mtlcoefs.c{i} = pde.StructuralMaterialAssignmentRecords.elasticityC(ca.YoungsModulus,ca.PoissonsRatio,numdim,ElasticityType);
                mtlcoefs.a{i} = pde.internal.FESolverUtilities.elasticityA(ca.YoungsModulus,ca.PoissonsRatio,ElasticityType);
                mtlcoefs.t{i} = pde.internal.FESolverUtilities.elasticityT(ca.YoungsModulus,ca.PoissonsRatio,ElasticityType);
                mtlcoefs.d{i} = 0;
                mtlcoefs.f{i} = 0;
                mtlcoefs.CTE{i} = ca.CTE;
                mtlcoefs.HystereticDamping{i} = ca.HystereticDamping;
                mtlcoefs.YoungsModulus{i} = ca.YoungsModulus;
                mtlcoefs.PoissonsRatio{i} = ca.PoissonsRatio;
                mtlcoefs.MassDensity{i} = ca.MassDensity;
                
                mtlstruct.Coefficients = mtlcoefs;
                complexcoef = complexcoef | hasComplexCoefficient(ca, loc, state);
            end
            
            mtlstruct.IsComplex = complexcoef;
            if any(mtlstruct.IsComplex)
                error(message('pde:materialPropertySpecification:invalidMtlPropValueComplex'));
            end
        end
    end
    
    methods(Hidden=true, Access=private)
        % Returns true if each material assignment has consistent definition across all
        % subdomains. For example, mass density defined in one region is
        % also defined in all other regions.
        function tf = consistentMaterialPropertyDefinition(self)
            tf = globalMaterialProperties(self);
            if tf
                return
            end
            tf = true;
            
            thiscoef = self.findStructuralProperties(self.qregiontype, 1);
            for i = 2:self.nRegion
                othercoef = self.findStructuralProperties(self.qregiontype, i);
                if ~thiscoef.coefficientsMatch(othercoef)
                    tf = false;
                    return;
                end
            end
        end
        
        function [location, state] = getSampleLocState(self)
            % This function is not necessary for isotropic material,
            % retaining it for future expansion purpose.
            msh = self.ParentModel.Mesh;
            nodecoords =msh.Nodes(:,1);
            location.x = nodecoords(1);
            location.y = nodecoords(2);
            if numel(nodecoords) == 3
                location.z = nodecoords(3);
            else
                location.z = 0;
            end
            location.subdomain=1;
            N = self.ParentModel.PDESystemSize;
            state.u(1:N,1)=0;   % state.u(1:NSystemSize, NumLocations)
            state.ux(1:N,1)=0;
            state.uy(1:N,1)=0;
            state.uz(1:N,1)=0;
            state.time=0;
        end
    end
    
    methods(Static, Hidden=true, Access=private)
        function c = elasticityC(E,nu,numdim,TwoDType)
            if (nu >= 0.5)
                error(message('pde:materialPropertySpecification:invalidIsotropicNu'))
            end
            
            if numdim == 2
                c = pde.internal.FESolverUtilities.elasticityC2D(E,nu,TwoDType);
            else
                c = pde.StructuralMaterialAssignmentRecords.elasticityC3D(E,nu);
            end
        end
        
        function cmat = elasticityC3D(E,nu)
            %   3-D solid elasticity c-matrix, symmetric form.
            
            G = E/(2*(1+nu));
            c1 = E*(1-nu)/((1+nu)*(1-2*nu));
            c12 = c1*nu/(1-nu);
            C11 = [c1 0 G 0 0 G];
            C12 = [0 G 0 c12 0 0 0 0 0];
            C13 = [0 0 G 0 0 0 c12 0 0];
            C22 = [G 0 c1 0 0 G];
            C23 = [0 0 0 0 0 G 0 c12 0];
            C33 = [G 0 G 0 0 c1];
            cmat  = [C11 C12 C22 C13 C23 C33]';
        end
 
    end
    
    methods(Hidden=true, Access={?pde.StructuralMaterialAssignment})
        function rType = qregiontype(self)
            if self.ParentModel.IsTwoD
                rType = 'face';
            else
                rType = 'cell';
            end
        end
        
        function n = nRegion(self)
            if ~isempty(self.ParentModel.Geometry)
                if self.ParentModel.IsTwoD
                    n = self.ParentModel.Geometry.NumFaces;
                else
                    n = self.ParentModel.Geometry.NumCells;
                end
            else
                n = [];
            end
        end
        
        function delistMaterialAssignments(self, mtltodelete)
            nummtl = numel(self.MaterialAssignments);
            for i = 1:nummtl
                thismtl = self.MaterialAssignments(i);
                if thismtl == mtltodelete
                    self.MaterialAssignments(i) = [];
                    break
                end
            end
            nummtl = numel(self.MaterialAssignments);
            if nummtl == 0 && isvalid(self.ParentModel)
                self.ParentModel.delistMaterialProperties();
            end
        end
    end
    
    methods(Static, Hidden=true)
        function delete(self,~)
            nummtl = numel(self.MaterialAssignments);
            for i = 1:nummtl
                if isvalid(self.MaterialAssignments(i))
                    delete(self.MaterialAssignments(i));
                end
            end
            
            if isvalid(self.ParentModel)
                self.ParentModel.delistMaterialProperties();
            end
        end
    end


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