classdef AnalyticGeometry < pde.GeometricModel & matlab.mixin.internal.CompactDisplay
    % pde.AnalyticGeometry 2D geometry object from PDEGEOM or DECSG geometry matrix
    %   G = AnalyticGeometry(PDEGEOM) Creates an AnalyticGeometry object
    %   from PDEGEOM, a geometry representation defined in a MATLAB-file, see
    %   PDEGEOM for details.
    %
    %   G = AnalyticGeometry(DECSG) Creates an AnalyticGeometry object
    %   from DECSG, a decomposed geometry matrix that represents the boundary
    %   edge segments of a 2D domain, see DECSG for details.
    %
    %   An AnalyticGeometry is typically created using the geometryFromEdges
    %   method of the pde.PDEModel class with similar arguments. The PDE toolbox
    %   supports different geometry formats; namely, DECSG and PDEGEOM for
    %   2D, and STL geometry for 3D. The object-based representation unifies
    %   these different formats to promote consistency across the toolbox.
    %
    % AnalyticGeometry properties:
    %   NumCells - The number of volumetric cells in the geometry
    %   NumFaces - The number of faces in the geometry
    %   NumEdges - The number of edges in the geometry
    %   NumVertices - The number of vertices in the geometry
    %
    % See also CREATEPDE, DECSG, PDEGEOM, pde.PDEModel, pde.PDEModel/geometryFromEdges,
    %          pde.GeometricModel
    
    % Copyright 2014-2020 The MathWorks, Inc.
    
    properties(SetAccess=private)
        NumCells;
        NumFaces;
        NumEdges;
        NumVertices;
    end
    
    properties(Hidden = true,SetAccess=protected)
        % Edges - array of pdeEdge objects defining the boundary of this geometry
        Edges = pde.internal.pdeEdge.empty;
    end
    
    properties(Hidden = true)
        geom;
    end
    
    methods
        function self=AnalyticGeometry(geom)
            % geom must be a handle to a user function
            if(isa(geom, 'function_handle'))
                self.geom = geom;
                self.NumCells = 0;
                self = self.setPdegeomEdges(geom());
                self = self.setPdegeomFaceCount();
                self = self.setPdegeomEdgeCount();
                self = self.setPdegeomVertexCount();
            elseif(isnumeric(geom))
                if isscalar(geom)
                    error(message('pde:pdeGeometry:scalarGeomEdges'))
                end
                if size(geom,1) < 7 || size(geom,1) > 12
                    error(message('pde:pdeGeometry:invalidDecsgMatrixSize'))
                end
                if(~(all(geom(1,:)==1 | geom(1,:)==2 | geom(1,:)==4)))
                    error(message('pde:pdeGeometry:invalidDecsgRowOne'))
                end
                validateattributes(geom,{'double'},{'real', 'finite','nonsparse', 'nonnan'});
                faceids = geom(6:7,:);
                faceids = unique(faceids(:));
                fcheck = diff(faceids)-1;
                if(any(fcheck))
                    error(message('pde:pdeGeometry:invalidFaceNumbering'))
                end
                % first row of dgm must be one of these integers
                self.geom = geom;
                self.NumCells = 0;
                self = self.setDecsgEdges(size(geom, 2));
                self = self.setDecsgFaceCount();
                self = self.setDecsgEdgeCount();
                self = self.setDecsgVertexCount();
            else
                error(message('pde:pdeGeometry:invalidArg'))
            end
        end
        
        function newgm = extrude(self,height)
            mdl = createpde;
            facets = facetAnalyticGeometry(mdl,self,0);
            newgm = analyticToDiscrete(facets);
            newgm = extrude(newgm, height);
        end
    end
    
    methods(Hidden=true)
        function vxl = vertexLabelLocations(self)
            [vxl, ~] = DeriveVertexData(self);
            % Make consistent
            vxl = {vxl(:,1), vxl(:,2)};
        end
        function videid = vertexIdToEdgeId(self)
            [~, videid] = DeriveVertexData(self);
        end
        function tf = IsTwoD(self) %#ok<MANU>
            tf = true;
        end
        
        function eaf=edgeAdjacentFaces(self, edis)
            if(isa(self.geom, 'function_handle'))
                dsgm=PdegeomToDecsg(self);
            else
                dsgm = self.geom;
            end
            lf = dsgm(6,:)';
            idx = (lf == 0);
            lf(idx) = NaN;
            rf = dsgm(7,:)';
            idx = (rf == 0);
            rf(idx) = NaN;
            eaf = [lf, rf];
            eaf = eaf(edis,:);
        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(Access=private)
        function self=setPdegeomEdges(self, numEdges)
            % Create an array of pdeEdge objects of length numEdges
            for i=1:numEdges
                self.Edges(i) = pde.internal.pdeEdge(self, i);
            end
        end
        function self=setPdegeomFaceCount(self)
            ne = self.geom();
            selr = self.geom(1:ne);
            faceids = selr(3:4,:);
            faceids = unique(faceids(:));
            zid = (faceids == 0);
            faceids(zid) = [];
            self.NumFaces=numel(faceids);
        end
        
        function self=setPdegeomEdgeCount(self)
            self.NumEdges=self.geom();
        end
        
        function self=setPdegeomVertexCount(self)
            ne = self.geom();
            selr = self.geom(1:ne);
            sp = selr(1,:);
            eidx = 1:ne;
            [sx, sy] = self.geom(eidx(:),sp(:));
            ep = selr(2,:);
            [ex, ey] = self.geom(eidx(:),ep(:));
            v = [sx(:) sy(:); ex(:) ey(:)];
            comptol = max(max(abs(v(:))));
            comptol = max(comptol,1)*eps*1000;
            v = uniquetol(v, comptol, 'ByRows',true);
            self.NumVertices = size(v,1);
        end
        
        function self=setDecsgEdges(self, numEdges)
            % Create an array of pdeEdge objects of length numEdges
            for i=1:numEdges
                self.Edges(i) = pde.internal.pdeEdge(self, i);
            end
        end
        
        function self=setDecsgFaceCount(self)
            faceids = self.geom(6:7,:);
            faceids = unique(faceids(:));
            zid = (faceids == 0);
            faceids(zid) = [];
            self.NumFaces=numel(faceids);
        end
        
        function self=setDecsgEdgeCount(self)
            self.NumEdges=size(self.geom, 2);
        end
        
        function self=setDecsgVertexCount(self)
            sx = self.geom(2,:)';
            ex = self.geom(3,:)';
            sy = self.geom(4,:)';
            ey = self.geom(5,:)';
            v = [sx sy; ex ey];
            comptol = max(max(abs(v(:))));
            comptol = max(comptol,1)*eps*1000;
            v = uniquetol(v, comptol, 'ByRows',true);
            self.NumVertices = size(v,1);
        end
        
        % Derive vertex coordinates and
        % the vertex to edge mapping.
        function [vxc, vxtoedge] = DeriveVertexData(self)
            if(isa(self.geom, 'function_handle'))
                dsgm=PdegeomToDecsg(self);
            else
                dsgm = self.geom;
            end
            sx = dsgm(2,:)';
            ex = dsgm(3,:)';
            sy = dsgm(4,:)';
            ey = dsgm(5,:)';
            vxc = [sx sy; ex ey];
            comptol = max(max(abs(vxc(:))));
            comptol = max(comptol,1)*eps*1000;
            [~, I, ~] = uniquetol(vxc, comptol, 'ByRows',true);
            vxc = vxc(sort(I),:);
            eid = [(1:numel(sx))';(1:numel(sx))';];
            vxtoedge = eid(sort(I));
        end
        
        function dgm=PdegeomToDecsg(self)
            geomFunc = self.geom;
            ne = geomFunc();
            edgeParams=geomFunc(1:ne);
            s = edgeParams(1,:);
            [xs,ys] = geomFunc(1:ne, s);
            s = edgeParams(2,:);
            [xe,ye] = geomFunc(1:ne, s);
            
            dgm = zeros(10,ne);
            dgm(1,:) = 2;
            dgm(2,:) = xs;
            dgm(3,:) = xe;
            dgm(4,:) = ys;
            dgm(5,:) = ye;
            dgm(6,:) = edgeParams(3,:);
            dgm(7,:) = edgeParams(4,:);
        end
        
    end
end


