classdef StructuralModel < pde.EquationModel
    %pde.StructuralModel  A container that represents a structural analysis model
    %  A structural model is a container that stores the geometry, material
    %  properties, body and boundary loads, boundary conditions, and mesh that
    %  define a structural problem. To create a StructuralModel, use the
    %  CREATEPDE function. The StructuralModel provides functions that allow
    %  you to: define 2D geometry or import 2D and 3D geometry, generate and
    %  store a mesh for these geometries, specify and store structural
    %  properties of the material, apply and store loads and boundary
    %  conditions that define the behavior at the boundary of the geometry, and
    %  solve the structural analysis model. The properties of the
    %  StructuralModel allow you to access this information easily.
    %
    % StructuralModel methods:
    %  geometryFromEdges       - Create 2D geometry from DECSG or PDEGEOM
    %  importGeometry          - Create 2D and 3D geometry from file (STL only)
    %  geometryFromMesh        - Create 3D geometry from a triangulated mesh
    %  generateMesh            - Generate a mesh from the geometry
    %  structuralProperties    - Specify structural material properties over a domain or subdomain
    %  structuralBodyLoad      - Specify body loads that act over a domain or subdomain
    %  structuralDamping       - Specify proportional damping parameters for transient dynamic analysis
    %  structuralBoundaryLoad  - Specify boundary loads that act over a geometric boundary
    %  structuralBC            - Apply a structural boundary conditions to the geometric boundary
    %  structuralIC            - Set initial displacement and velocity
    %  solve                   - Solve the model
    %
    % StructuralModel properties:
    %  AnalysisType            - Static analysis: 'static-solid','static-planestress', or 'static-planestrain'
    %                          - Modal analysis: 'modal-solid','modal-planestress', or 'modal-planestrain'
    %                          - Transient dynamics: 'transient-solid','transient-planestress', or 'transient-planestrain'
    %  Geometry                - Geometric representation of the domain
    %  MaterialProperties      - Structural material properties within the domain or subdomain
    %  BodyLoads               - Loads that act the domain or subdomain
    %  BoundaryConditions      - Structural loads and boundary conditions applied to the geometric boundary
    %  ReferenceTemperature    - Reference temperature for thermal loading
    %  InitialConditions       - Initial displacement and velocity
    %  SuperelementInterfaces  - Superelement interfaces
    %  Mesh                    - Discretization of the domain
    %  SolverOptions           - Algorithm options for the solvers
    %
    %    See also createpde, pde.StructuralModel/structuralProperties,
    %    pde.StructuralModel/structuralBodyLoad, pde.StructuralModel/structuralBC,
    %    pde.StructuralModel/structuralBoundaryLoad, pde.StructuralModel/structuralIC,
    %    pde.StructuralModel/solve, pde.PDEModel

    % Copyright 2017-2023 The MathWorks, Inc.

    properties
        %AnalysisType - Analysis attribute
        AnalysisType;
    end

    properties
        % MaterialProperties - Structural material properties within the domain or subdomain
        %    An object that contains the structural material property
        %    assignments within the geometric domain. Material properties are
        %    specified and assigned to the geometry using the
        %    structuralProperties function. If multiple assignments are made,
        %    then they are all recorded within this object. Material properties
        %    must be specified in units consistent with the geometry and other
        %    model attributes.
        %
        %    See also pde.StructuralMaterialAssignmentRecords/findStructuralProperties,
        %             pde.StructuralMaterialAssignmentRecords/globalMaterialProperties
        MaterialProperties;

        % BodyLoads - Loads that act on complete domain or a subdomain
        %    An object that contains the load assignments
        %    for the geometric domain. Body loads are specified and
        %    assigned to the geometry using the structuralBodyLoad function. If
        %    multiple assignments are made, then they are all recorded within
        %    this object. Body loads must be specified in units consistent
        %    with geometry and other model attributes.
        %
        %    See also pde.BodyLoadAssignmentRecords/findBodyLoad,
        %             pde.BodyLoadAssignmentRecords/globalBodyLoad
        BodyLoads;

        % BoundaryConditions - Structural boundary conditions applied to the geometric boundary
        %    An object that contains the structural boundary condition assignments
        %    on the boundary of the geometry. Boundary conditions are specified
        %    and assigned to the geometry using the structuralBC function. If
        %    multiple assignments are made, then they are all recorded within
        %    this object.
        %
        %    See also pde.StructuralBCRecords/findStructuralBC
        BoundaryConditions;

        %DampingModels - Damping model for transient dynamic analysis
        %    An object that contains the damping assignments for the whole
        %    domain. Damping model and its values specified and assigned to the
        %    model using the structuralDamping function. If multiple
        %    assignments are made, then they are all recorded within this
        %    object. Damping values must be specified in units consistent with
        %    geometry and other model attributes.
        %
        %    See also pde.StructuralDampingAssignmentRecords/findStructuralDamping
        DampingModels;

        %InitialConditions - Initial displacement and velocity
        %    An object that contains the initial condition assignments within
        %    the geometric domain. Initial condition is specified and assigned
        %    to the geometry using the structuralIC function. If multiple
        %    assignments are made, then they are all recorded within this
        %    object.
        %
        %    See also pde.StructuralICRecords/findStructuralIC,
        %             pde.GeometricStructuralIC, pde.NodalStructuralIC
        InitialConditions;

        %ReferenceTemperature - Reference temperature for thermal loading
        %    A reference temperature that corresponds to zero thermal-stress
        %    state of the model. The default zero temperature implies that the
        %    thermal load is specified in terms of the temperature change and
        %    its derivatives. Reference temperature must be a numeric scalar.
        %    To specify the thermal load itself, use the structuralBodyLoad
        %    function.
        %
        %    See also pde.StructuralModel/structuralBodyLoad.
        ReferenceTemperature;

        %SuperelementInterfaces - Superelement interfaces
        %    An object that contains the superelement interface specifications.
        %    Superelement interfaces can be specified using
        %    structuralSEInterface or structuralBC. Use structuralSEInterface
        %    to retain all the degrees-of-freedom associated with the specified
        %    geometric region. Alternatively, specify a 'multipoint' constraint
        %    using the structuralBC function to define a superelement interface
        %    with six degrees of freedom at the reference point of the
        %    multipoint constraint.
        SuperelementInterfaces;

        LinearizeInputs;
        LinearizeOutputs;
    end


    methods
        function self = StructuralModel(analysisType)
            arguments
                analysisType (1,:) char = "static-solid"
            end
            self@pde.EquationModel(3);
            self.AnalysisType = analysisType;
            self.PDESystemSize = self.Dim;
            self.SolverOptions.RigidMotionCheck = true;
        end

        %methods declaration.
        mtl = structuralProperties(self,varargin)
        sbl = structuralBodyLoad(self,varargin)
        ssl = structuralBoundaryLoad(self,varargin)
        sbc = structuralBC(self,varargin)
        sdm = structuralDamping(self,varargin)
        sol = solve(self,varargin);
        sol = structuralSEInterface(self,varargin);
        rom = reduce(self,varargin);
        li = linearizeInput(self,varargin);
        lo = linearizeOutput(self,varargin);
        output=linearize(self,varargin);

        function set.AnalysisType(self, inputAnalysis)
            validateattributes(inputAnalysis,{'string','char'},{'nonempty'});
            validAnalysisStrings  = ["static-solid","static-planestress","static-planestrain","static-axisymmetric",...
                "modal-solid","modal-planestress","modal-planestrain","modal-axisymmetric",...
                "transient-solid","transient-planestress","transient-planestrain","transient-axisymmetric",...
                "frequency-solid","frequency-planestress","frequency-planestrain","frequency-axisymmetric"];

            try
                inputAnalysis= pde.EquationModel.validateAnalysisType(inputAnalysis, validAnalysisStrings);
            catch ME
                if (strcmp(ME.identifier,'MATLAB:unrecognizedStringChoice'))
                    error(message('pde:StructuralModel:invalidAnalysisType'));
                else
                    rethrow(ME)
                end
            end

            if ~isempty(self.AnalysisType) %Change analysis
                pde.EquationModel.validateAnalysisChange(self.AnalysisType, inputAnalysis);
            end
            
            %Tighten tolerance for transient analysis, if they are not user
            % defined.
            if contains(inputAnalysis,"transient", 'IgnoreCase',true) && ~isempty(self.SolverOptions)
                % Get default tolerance.
                defaultSolverOption = pde.PDESolverOptions;

                % Tighten tolerance if they are equal to default.
                if (self.SolverOptions.AbsoluteTolerance == defaultSolverOption.AbsoluteTolerance)
                    self.SolverOptions.AbsoluteTolerance = 1E-7;
                end
                if (self.SolverOptions.RelativeTolerance == defaultSolverOption.RelativeTolerance)
                    self.SolverOptions.RelativeTolerance = 1E-5;
                end
            end

            self.AnalysisType = string(inputAnalysis);
        end

        function set.ReferenceTemperature(self, t)
            validateattributes(t,{'double'},{'scalar','real', 'nonsparse','finite'},'ReferenceTemperature');
            if (contains(self.AnalysisType,{'modal','transient'})) && ~isempty(t) %#ok<MCSUP>
                warning(message('pde:StructuralModel:thermalLoadingOnlyInStatic'));
            end
            self.ReferenceTemperature =t;
        end

        function set.BoundaryConditions(self, bc)

            if isempty(bc) && ~isempty(self.BoundaryConditions)
                delete(self.BoundaryConditions);
                self.BoundaryConditions = [];
                return
            end

            if ~isempty(bc) && ~isa(bc, 'pde.StructuralBCRecords')
                error(message('pde:pdeModel:invalidBoundaryConditions'));
            end

            if ~isempty(bc) &&  ~isscalar(bc)
                error(message('pde:pdeModel:nonScalarBoundaryConditions'));
            end

            if ~isempty(bc) && (self ~= bc.ParentModel)
                if self.Geometry ~= bc.ParentModel.Geometry
                    error(message('pde:StructuralModel:inconsistentBCandGeometry'));
                end
                BCcont = pde.StructuralBCRecords(self);
                BCcont.StructuralBCAssignments =  bc.StructuralBCAssignments;
                self.BoundaryConditions = BCcont;
            else
                self.BoundaryConditions = bc;
            end

        end


        function set.MaterialProperties(self, mp)

            if(~isempty(mp) && ~isa(mp, 'pde.StructuralMaterialAssignmentRecords'))
                error(message('pde:StructuralModel:invalidMaterialProperties'));
            end

            if ~isscalar(mp) && ~isempty(mp)
                error(message('pde:StructuralModel:nonScalarMaterialProperties'));
            end

            if ~isempty(mp) && (self ~= mp.ParentModel)
                if self.Geometry ~= mp.ParentModel.Geometry
                    error(message('pde:StructuralModel:inconsistentMtlandGeometry'));
                end
                mtlContainer = pde.StructuralMaterialAssignmentRecords(self);
                mtlContainer.MaterialAssignments = mp.MaterialAssignments;
                self.MaterialProperties =  mtlContainer;
            else
                self.MaterialProperties = mp;
            end
        end

        function set.BodyLoads(self, bl)

            if(~isempty(bl) && ~isa(bl, 'pde.BodyLoadAssignmentRecords'))
                error(message('pde:StructuralModel:invalidBodyLoad'));
            end

            if ~isscalar(bl) && ~isempty(bl)
                error(message('pde:StructuralModel:nonScalarBodyLoad'));
            end

            if ~isempty(bl) && (self ~= bl.ParentModel)
                if self.Geometry ~= bl.ParentModel.Geometry
                    error(message('pde:StructuralModel:inconsistentBodyLoadandGeometry'));
                end
                blContainer = pde.BodyLoadAssignmentRecords(self);
                blContainer.BodyLoadAssignments = bl.BodyLoadAssignments;
                self.BodyLoads =  blContainer;
            else
                self.BodyLoads = bl;
            end
        end


        function set.DampingModels(self, dmpr)

            if(~isempty(dmpr) && ~isa(dmpr, 'pde.StructuralDampingAssignmentRecords'))
                error(message('pde:StructuralModel:invalidDampingAssignment'));
            end

            if ~isscalar(dmpr) && ~isempty(dmpr)
                error(message('pde:StructuralModel:nonScalarDampingAssignment'));
            end

            if ~isempty(dmpr) && (self ~= dmpr.ParentModel)
                if self.Geometry ~= dmpr.ParentModel.Geometry
                    error(message('pde:StructuralModel:inconsistentDampingModelandGeometry'));
                end
                dmContainer = pde.StructuralDampingAssignmentRecords(self);
                dmContainer.StructuralDampingAssignments = dmpr.StructuralDampingAssignments;
                self.DampingModels = dmContainer;
            else
                self.DampingModels = dmpr;
            end
        end


        function set.InitialConditions(self, icr)

            if(~isempty(icr) && ~isa(icr, 'pde.StructuralICRecords'))
                error(message('pde:pdeModel:invalidInitialConditions'));
            end

            if ~isscalar(icr) && ~isempty(icr)
                error(message('pde:pdeModel:nonScalarInitialConditions'));
            end

            if ~isempty(icr) && (self ~= icr.ParentPdemodel)
                if self.Geometry ~= icr.ParentPdemodel.Geometry
                    error(message('pde:StructuralModel:inconsistentICandGeometry'));
                end
                dmContainer = pde.StructuralICRecords(self);
                dmContainer.StructuralICAssignments = icr.StructuralICAssignments;
                self.InitialConditions = dmContainer;
            else
                self.InitialConditions = icr;
            end
        end

        function set.LinearizeInputs(self,sli)
            self.LinearizeInputs=sli;
        end

        function  out = get.IsTimeDependent(self)
            out = false;
            if  contains(self.AnalysisType,"transient")
                out = true;
            end
        end   

        function  out = get.Dim(self)
            out = 2;
            if  contains(self.AnalysisType,"solid")
                out = 3;
            end
        end   

        function  out = get.ElasticityType(self)
            if  contains(self.AnalysisType,"solid")
                out = '3D';
            elseif contains(self.AnalysisType,"planestress")
                out = 'ps';
            elseif contains(self.AnalysisType,"planestrain")
                out = 'pn';
            elseif contains(self.AnalysisType,"axisymmetric")
                out = 'as';
            end
        end


    end

    methods(Hidden)

        function geometricConsistencyCheck(self,gm)

            if isempty(gm)
                return;
            end

            if strcmpi(self.ElasticityType,'3D') && gm.IsTwoD
                error(message('pde:StructuralModel:invalid2DGeometryFor3DModel'));
            end

            if ~strcmpi(self.ElasticityType,'3D') && ~gm.IsTwoD
                error(message('pde:StructuralModel:invalid3DGeometryFor2DModel'));
            end
        end

        function ok = interfaceConditionsAllowed(self)
            ok =  isprop(self,'IsThermalStress');
        end

        function validateRegionID(self,inputVertex, inputEdge, inputFace)
            regionsLogical = [~isempty(inputVertex), ~isempty(inputEdge),  ~isempty(inputFace)];
            regionsString = ["Vertex","Edge","Face"];
            regionsStringPlural = ["vertices","edges","faces"];

            if nnz(regionsLogical) ~= 1 % Specify only one region type in a call.
                error(message('pde:StructuralModel:multipleRegionTypes'));
            end

            regionIDs = [inputVertex, inputEdge, inputFace]; % Only one type of region would be non-empty.
            maxRegionID = [self.Geometry.NumVertices, self.Geometry.NumEdges, self.Geometry.NumFaces];

            if any(regionIDs > maxRegionID(regionsLogical))
                error(message('pde:StructuralModel:invalidRegionIndex',regionsString(regionsLogical),regionsStringPlural(regionsLogical)));
            end
        end

        function [coefstruct,constCoeffs,coefstructForThermalStress, hystereticDamping] = packCoefficients(self)

            dim = self.Dim;
            mtlstruct = self.MaterialProperties.packProperties();
            hystereticDamping = mtlstruct.Coefficients.HystereticDamping;
            if ~isempty(self.BodyLoads)
                massDensity = mtlstruct.Coefficients.MassDensity;
                bodyLoadStr = self.BodyLoads.packBodyLoad(massDensity);
                isBodyLoadGlobal = self.BodyLoads.globalBodyLoad();
                bodyLoadStr.IsGlobal = isBodyLoadGlobal;
                self.IsThermalStress = contains(self.AnalysisType,'static') && any([mtlstruct.Coefficients.CTE{:}]) && ~isempty([bodyLoadStr.Coefficients.Temperature{:}]);
            else
                bodyLoadStr.Coefficients.f = ...
                    repmat({zeros(dim,1)},1,numel(mtlstruct.Coefficients.m));
                bodyLoadStr.IsComplex =false;
                bodyLoadStr.Coefficients.Temperature = ...
                    repmat([],1,numel(mtlstruct.Coefficients.m));
                bodyLoadStr.Coefficients.TimeStep = [];
                bodyLoadStr.Coefficients.needsInterpolation = {false};
                bodyLoadStr.IsGlobal = true;
                self.IsThermalStress = false;
            end

            coefstruct.ElementsInSubdomain = mtlstruct.ElementsInSubdomain;
            coefstruct.NumElementsInSubdomain = mtlstruct.NumElementsInSubdomain;
            coefstruct.Coefficients.m = mtlstruct.Coefficients.m;
            if isempty(mtlstruct.Coefficients.m{1})
                coefstruct.Coefficients.m = mtlstruct.Coefficients.d; % Static-cases
            end
            coefstruct.Coefficients.d = mtlstruct.Coefficients.d;
            coefstruct.Coefficients.c = mtlstruct.Coefficients.c;
            coefstruct.Coefficients.a = mtlstruct.Coefficients.a;
            coefstruct.Coefficients.f = bodyLoadStr.Coefficients.f;
            coefstruct.IsComplex = [false;false;false;false;false];

            isMaterialAssignmentGlobal = self.MaterialProperties.globalMaterialProperties();
            coefstructForThermalStress = struct;
            if self.IsThermalStress
                if bodyLoadStr.Coefficients.needsInterpolation{1}
                    bodyLoadStr.Coefficients.Tintrp{1} = self.interpolateThermalBodyLoad(bodyLoadStr);
                else
                    bodyLoadStr.Coefficients.Tintrp{1} = bodyLoadStr.Coefficients.Temperature{1};
                end

                if self.IsTwoD
                   numSubDomains = self.Geometry.NumFaces;
                else
                   numSubDomains = self.Geometry.NumCells;
                end
                coefstruct.Coefficients.f = self.constructFCoefficient(mtlstruct,bodyLoadStr,isMaterialAssignmentGlobal,numSubDomains,self.ElasticityType);
                coefstructForThermalStress.Coefficients.PoissonsRatio = mtlstruct.Coefficients.PoissonsRatio;
                coefstructForThermalStress.Coefficients.YoungsModulus = mtlstruct.Coefficients.YoungsModulus;
                coefstructForThermalStress.Coefficients.CTE = mtlstruct.Coefficients.CTE;
                coefstructForThermalStress.Coefficients.T = bodyLoadStr.Coefficients.Temperature;
                coefstructForThermalStress.Coefficients.Tintrp = bodyLoadStr.Coefficients.Tintrp;
                coefstructForThermalStress.Coefficients.needsInterpolation = bodyLoadStr.Coefficients.needsInterpolation;
                coefstructForThermalStress.Coefficients.TimeStep = bodyLoadStr.Coefficients.TimeStep;
                coefstructForThermalStress.IsGlobal =   bodyLoadStr.IsGlobal;
            end

            % Internal coefficient for axisymmetric analysis.
            if contains(self.AnalysisType,"axisymmetric")
                coefstruct.Coefficients.t = mtlstruct.Coefficients.t;
                coefstruct.IsComplex(6) = false;
            end
            allbodyLoadsnumeric = true; % Functional form of body load is not supported.
            constCoeffs = self.MaterialProperties.allMaterialPropertiesNumeric() && allbodyLoadsnumeric;
        end

        %methods declaration.
        romStruct = reduceSimscape(self,varargin);
    end

    methods (Access = protected,Hidden = true)
        function group = getPropertyGroups(self)
            if contains(self.AnalysisType,"frequency")
                groupNames = {'AnalysisType','Geometry', ...
                    'MaterialProperties','BodyLoads','BoundaryConditions',...
                    'DampingModels','SuperelementInterfaces','Mesh','SolverOptions'};
            elseif contains(self.AnalysisType,"modal")
                groupNames = {'AnalysisType','Geometry', ...
                    'MaterialProperties',...
                    'BoundaryConditions','SuperelementInterfaces','Mesh','SolverOptions'};
            elseif contains(self.AnalysisType,"transient")
                groupNames = {'AnalysisType','Geometry', ...
                    'MaterialProperties','BodyLoads','BoundaryConditions',...
                    'DampingModels','InitialConditions','SuperelementInterfaces','Mesh','SolverOptions'};
            else
                groupNames = {'AnalysisType','Geometry', ...
                    'MaterialProperties','BodyLoads',...
                    'BoundaryConditions','ReferenceTemperature','SuperelementInterfaces','Mesh','SolverOptions'};
            end

            group = matlab.mixin.util.PropertyGroup(groupNames);

        end
    end

    properties (SetAccess='private',Hidden = true)
        IsThermalStress;
        DCoefficientFlag = false; % Second-order PDE, m-coeff is used.
    end

    properties (SetAccess='private',Hidden = true, Dependent = true)
        IsTimeDependent;
        ElasticityType;
        Dim;
    end

    methods(Hidden = true, Static  = true)
        function ok = isValidDampingModel(dm)
            if isempty(dm)
                return;
            end

            if ~any(contains(["proportional","modal"],dm,'IgnoreCase',true))
                error(message('pde:StructuralModel:invalidDampingModel'))
            end
            ok=true;
        end
    end

    methods(Hidden = true, Access = {?pde.StructuralMaterialAssignmentRecords})
        function delistMaterialProperties(self)
            if isvalid(self)
                self.MaterialProperties = [];
            end
        end
    end

    methods(Hidden = true, Access = {?pde.BodyLoadAssignmentRecords})
        function delistBodyLoads(self)
            if isvalid(self)
                self.BodyLoads = [];
            end
        end
    end


    methods(Hidden = true, Access = {?pde.StructuralDampingAssignmentRecords})
        function delistDampingAssignment(self)
            if isvalid(self)
                self.DampingModels = [];
            end
        end
    end

    methods(Hidden = true, Access = {?pde.StructuralBCRecords})
        function delistBoundaryConditions(self)
            if isvalid(self)
                self.BoundaryConditions = [];
            end
        end
    end

    methods(Hidden = true, Access = {?pde.StructuralICRecords})
        function delistInitialConditions(self)
            if isvalid(self)
                self.InitialConditions = [];
            end
        end
    end

    methods(Hidden = true, Access = {?pde.StructuralSEIAssignmentRecords})
        function delistSEAssignment(self)
            if isvalid(self)
                self.SuperelementInterfaces = [];
            end
        end
    end


end


