classdef (Abstract) EquationModel < handle & matlab.mixin.CustomDisplay & pde.internal.FESolverUtilities
    % pde.EquationModel  An abstract base class for various PDE models.
    %
    %    See also CREATEPDE, pde.PDEModel, pde.ThermalModel, pde.StructuralModel
    
    % Copyright 2016-2022 The MathWorks, Inc.

    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;
    end
  
    methods
        function obj = EquationModel(numEqns)
            validateattributes(numEqns,{'numeric'},...
                {'scalar','integer','positive'});
            obj@pde.internal.FESolverUtilities(numEqns)
        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 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,varargin)
        [g,m] = geometryFromMesh(self, nodes,elems, varargin)
        msh = generateMesh(self, varargin)
    end

    properties (Dependent=true,Hidden=true)
        IsTwoD;
    end

    methods (Hidden = true)
        function ok = geometricConsistencyCheck(~,~)
            % Derived classes can overload this method for additional checks.
            ok = true;
            return
        end

        function ok = interfaceConditionsAllowed(~)
            % Derived classes can overload this method.
            ok = false;
        end

        function varargout = rotate(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = scale(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = translate(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = addFace(~,varargin)
            error(message('pde:pdeModel:functionForNonFunctionHandleGeometry'));
        end

        function varargout = addVertex(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end

        function varargout = nearestFace(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = nearestEdge(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = facesAttachedToEdges(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = cellFaces(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = cellEdges(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end

        function varargout = faceEdges(~,varargin)
            error(message('pde:pdeModel:functionForGeomOnly'));
        end
    end

    methods (Hidden = true)
        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)
        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

        function analysisType = validateAnalysisType(analysisTypeInput,validStrings)
            analysisTypeInput = lower(analysisTypeInput);
            exactMatches = validStrings == analysisTypeInput;
            if ~any(exactMatches)
                partialMatches = startsWith(validStrings,analysisTypeInput);
                if nnz(partialMatches) > 1
                    error(message('pde:createpde:invalidAnalysisType',analysisTypeInput))
                end
                analysisTypeInput = validatestring(analysisTypeInput,validStrings); % complete to shortest
            end
            analysisType = analysisTypeInput;
        end

        function  validateAnalysisChange(currentAnalysis,inputAnalysis)
            currentAnalysisParts = [];
            if ~isempty(currentAnalysis)
                currentAnalysisParts = split(currentAnalysis,"-");
            end
            inputAnalysisParts = split(inputAnalysis,"-");

            if numel(currentAnalysisParts) ~= numel(inputAnalysisParts)
                error(message('pde:createpde:invalidChangeOfAnalysisType',currentAnalysis,inputAnalysis));
            end

            if numel(currentAnalysisParts) == 2
                currentSpatialDomain = currentAnalysisParts(2);
                if inputAnalysisParts(2) ~= currentSpatialDomain
                    error(message('pde:createpde:invalidChangeOfAnalysisType',currentAnalysis,inputAnalysis));
                end
            end
        end

        function [regionAxisymm, stateAxisymm] = getAxisymmStructs(region,state)
            regionAxisymm = struct;
            stateAxisymm = struct;
            if ~isempty(region)
                regionAxisymm = region;
                regionAxisymm.r = region.x;
                regionAxisymm.z = region.y;
            end

            if ~isempty(state)
                stateAxisymm = state;
                stateAxisymm.ur = state.ux;
                stateAxisymm.uz = state.uy;
            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
        function ok = isValidHface(hf)
            if ~iscell(hf) || bitget(numel(hf),1) || isempty(hf)
                error(message('pde:pdeModel:invalidHface'));
            end
            ok = true;
        end
        function ok = isValidHedge(he)
            if ~iscell(he) || bitget(numel(he),1) || isempty(he)
                error(message('pde:pdeModel:invalidHedge'));
            end
            ok = true;
        end
        function ok = isValidHvertex(hv)
            if ~iscell(hv) || bitget(numel(hv),1) || isempty(hv)
                error(message('pde:pdeModel:invalidHvertex'));
            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 ok=ValidateLabelVal(labelVal)
            % Label Val must be string or char
            if ~isstring(labelVal) && ~ischar(labelVal)
                error(message('pde:linearize:InvalidType'));
            end

	    %Label val can't string vector
            if isstring(labelVal) && ~isscalar(labelVal)
                error(message('pde:linearize:invalidStringArray'));
            end

            % Label value must not be empty
            if isempty(strtrim(convertStringsToChars(labelVal)))
                error(message('pde:StructuralModel:NoLabelVal'));
            end            

            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

