classdef (Abstract) EquationModel < handle & matlab.mixin.CustomDisplay & dynamicprops
    % pde.EquationModel  An abstract base class for various PDE models.
    %
    %    See also CREATEPDE, pde.PDEModel, pde.ThermalModel, pde.StructuralModel
    
    % Copyright 2016-2020 The MathWorks, Inc.
    
    properties
        % PDESystemSize - Number of PDEs in the system
        PDESystemSize;
    end
    
    properties (SetObservable)
        % Geometry - Geometric representation of the domain
        %    A scalar object that defines the geometric domain. The geometry object
        %    is created by the geometryFromEdges function for 2D or the importGeometry
        %    function for 3D. Assemblies of geometric models are currently not supported.
        Geometry;
        
        % Mesh - Discretization of the domain
        %    An object that defines the meshed domain. The mesh object stores the
        %    node and element data. The mesh is created using the generateMesh function.
        Mesh;
    end
    
    properties (SetAccess='protected')
        % SolverOptions - Algorithm options for the PDE solvers
        %
        %    See also pde.PDESolverOptions
        SolverOptions;
    end
    
    methods
        function obj = EquationModel(numEqns)
            validateattributes(numEqns,{'numeric'},...
                {'scalar','integer','positive'});
            obj.PDESystemSize = numEqns;
            obj.SolverOptions = pde.PDESolverOptions();
        end
        
        function set.Geometry(self, gm)
            if( ~isempty(gm) && (~isa(gm, 'pde.DiscreteGeometry') && ~isa(gm, 'pde.AnalyticGeometry')) )
                error(message('pde:pdeModel:invalidGeometry'));
            end
            if ~isscalar(gm) && ~isempty(gm)
                error(message('pde:pdeModel:noAssemSupport'));
            end
            self.geometricConsistencyCheck(gm);
            
            self.Geometry = gm;
        end
        
        function set.Mesh(self, msh)
            if(~isempty(msh) && ~isa(msh, 'pde.FEMesh'))
                error(message('pde:pdeModel:invalidMesh'));
            end
            if ~isscalar(msh) && ~isempty(msh)
                error(message('pde:pdeModel:noAssemSupport'));
            end
            self.Mesh = msh;
        end
        
        function set.PDESystemSize(self, sz)
            validateattributes(sz,{'numeric'},{'scalar','integer','positive'});
            self.PDESystemSize = sz;
        end
        
        function tf = get.IsTwoD(self)
            if isempty(self.Geometry)
                tf = [];
            else
                tf = self.Geometry.IsTwoD();
            end
        end
        
        % Method declaration
        gm = geometryFromEdges(self, gedges)
        gm = importGeometry(self, geofilename)
        [g,m] = geometryFromMesh(self, nodes,elems, varargin)
        msh = generateMesh(self, varargin)
    end
    
    properties (Dependent=true,Hidden=true)
        IsTwoD;
    end
    
    methods (Hidden = true)
        %Internal methods declarations
        [a,b,c] = assema(self, varargin)
        sol = solveStationary(self, varargin)
        [sol,dsol,ddsol] = solveTimeDependent(self, varargin)
        [u, res] = solveStationaryNonlinear(self,varargin);
        [evec, evalue] = solveEigenvalue(self,varargin);
        femstruct = assembleFEMatricesInternal(self, varargin)
        varargout = assembleBoundary(self,varargin)
        femat=assembleSelectedFEMatrices(self,varargin)
        msh = generateMeshFacettedAnalytic(self,Hmax,Hmin,Hgrad,geomOrder)
        gm = analyticToDiscrete(self)
        F = facetAnalyticGeometry(self, g, Hmax)
        
        function res = addprop(varargin)
            %Dynamic properties are used for internal purposes by
            %inheriting from dynamicprops class. This method makes the
            %public method of dynamicprops, addprop, hidden.
            res = addprop@dynamicprops(varargin{:});
        end
        
        function ok = geometricConsistencyCheck(~,~)
            % Derived classes can overload this method for additional checks.
            ok = true;
            return
        end
        
        function varargout = rotate(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end
        
        function varargout = scale(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end
        
        function varargout = translate(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end
        
        function varargout = addFace(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end
        
        function varargout = addVertex(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end
        
    end
    
    methods (Hidden = true, Access = protected)
        performSolverPrecheck(self,checkics)
    end
    
    methods (Hidden = true)
        function [location, state] = getSampleLocState(self,sid)
            nodecoords = zeros(3,1);
            msh = self.Mesh;
            if ~isempty(msh) && isvalid(msh)
                %get coordinate of first node in the subdomain
                subdomainNodes = msh.findNodes('region','Face',sid);
                nodecoords = msh.Nodes(:,subdomainNodes(1));
            end
            systemsize = self.PDESystemSize;
            location.x =  nodecoords(1);
            location.y =  nodecoords(2);
            location.subdomain=sid;
            state.u = ones(systemsize, 1);
            state.ux = ones(systemsize, 1);
            state.uy = ones(systemsize, 1);
            state.time = 1;
            
            if ~self.IsTwoD
                location.z =  nodecoords(3);
                state.uz = ones(systemsize, 1);
            end
            
            if isprop(self,'AnalysisType') && contains(self.AnalysisType,'axisymmetric')
                location.r = location.x;
                location.z = location.y;
                state.ur = state.ux;
                state.uz = state.uy;
            end
        end
    end
    
    methods (Hidden = true, Static = true)
        K  =  checkForSymmetry(K);
    
        function checkTlist(tlist)
            if length(tlist) < 2
                error(message('pde:pdeODEInfo:tlistTooShort'));
            end
            t0 = tlist(1);
            tfinal = tlist(end);
            if(t0 == tfinal)
                error(message('pde:pdeODEInfo:tlistEndpointsNotDistinct'));
            end
            tdir = sign(tfinal - t0);
            if any( tdir*diff(tlist) <= 0 )
                error(message('pde:pdeODEInfo:tlistUnordered'));
            end
        end
    end
    
    methods(Hidden = true, Static  = true, Access = protected)
        function ok = isValidHmax(hval)
            if ~isreal(hval) || ~isscalar(hval) || ischar(hval) || hval < 0 || issparse(hval) || ~isfinite(hval)
                error(message('pde:pdeModel:invalidHmax'));
            end
            ok = true;
        end
        function ok = isValidHmin(hval)
            if ~isreal(hval) || ~isscalar(hval) || ischar(hval) || hval < 0 || issparse(hval) || ~isfinite(hval)
                error(message('pde:pdeModel:invalidHmin'));
            end
            ok = true;
        end
        function ok=isValidHgrad(hgrad)
            if ~isreal(hgrad) || ~isscalar(hgrad) || ischar(hgrad) || issparse(hgrad) || ~isfinite(hgrad)
                error(message('pde:pdeModel:invalidHgrad'));
            elseif hgrad<1 || hgrad>2
                error(message('pde:pdeModel:invalidHgrad'));
            end
            ok=true;
        end
        function ok = isValidGeomOrder(go)
            go = convertStringsToChars(go);
            if(~ischar(go) || isempty(go))
                error(message('pde:pdeModel:invalidGeomOrder'))
            end
            nc = numel(go);
            if ~(strncmpi(go,'linear',nc) || strncmpi(go,'quadratic',nc))
                error(message('pde:pdeModel:invalidGeomOrder'))
            end
            ok=true;
        end
    end
    
    
    methods(Static, Hidden = true)
        function ok = isValidEntityID(entityids)
            if strcmp(entityids,'default')
                ok=true;
                return
            end
            validateattributes(entityids,{'numeric'},{'integer','positive', 'nonzero','real', 'nonsparse','vector'});
            ok=true;
        end
        
        function ok = isValidAttributeVal(attributeVal)
            validateattributes(attributeVal,{'double','function_handle'},{})
            ok=true;
        end
        
        function ok = isValidNumericOnlyVal(attributeVal)
            validateattributes(attributeVal,{'double'},{})
            ok=true;
        end
        
        function ok = isValidNumericPositiveOnlyVal(attributeVal)
            validateattributes(attributeVal,{'double'},{'positive'})
            ok=true;
        end
        
        function delistMesh(~,~)
            %         if isvalid(self.ParentPdemodel)
            %             self.ParentPdemodel.Mesh = [];
            %         end
        end
        function delistGeometry(~,~)
            %         if isvalid(self.ParentPdemodel)
            %             self.ParentPdemodel.Geometry = [];
            %         end
        end
    end
    
    
end

