classdef AnalyticGeometry < pde.GeometricModel &pde.PDECompactDisplay
    % 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 methods:
    %    facesAttachedToEdges - Find faces attached to the specified edges.
    %    nearestFace - Find faces nearest to the specified point.
    %    faceEdges - Find edges attached to the specified faces. 
    %    nearestEdge - Find edges nearest to the specified point.
    %
    % AnalyticGeometry properties:
    %   NumCells - The number of volumetric cells in the geometry. For an 
    %     AnalyticGeometry object, this number is always 0.
    %   NumFaces - The number of faces in the geometry
    %   NumEdges - The number of edges in the geometry
    %   NumVertices - The number of vertices in the geometry
    %   Vertices - The coordinates of the vertices in the geometry
    %
    % See also CREATEPDE, DECSG, PDEGEOM, pde.PDEModel, pde.PDEModel/geometryFromEdges,
    %          pde.GeometricModel
    
    % Copyright 2014-2023 The MathWorks, Inc.
    
    properties(SetAccess=private)
        NumCells;
        NumFaces;
        NumEdges;
        NumVertices;
        Vertices;
    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
    
    % Vertices are packed when the object is loaded. Needed when loading
    % objects pre-R2021a.
    methods (Hidden = true, Static = true, Access = protected)
        function obj = loadobj(obj)
            % function called during loading an object of this type from
            % a MAT-file, vertices property need to be created
            obj.Vertices = getVertices(obj);
         end
    end
    
    % Object method to get the vertices
    methods (Hidden = true)
        function v = getVertices(obj)
            if(isa(obj.geom, 'function_handle'))
                v = getPdegeomVertices(obj);
            else
                v = getDecsgVertices(obj);
            end
            [~,I,~] = uniquetol(v,'ByRows',true);
            v = v(sort(I),:);
        end
    end
    
    methods
        function self=AnalyticGeometry(geom)
            % geom must be a handle to a user function
            if(isa(geom, 'function_handle'))
                checkFunctionHandle(self,geom);
                self.geom        = geom;
                self.NumCells    = 0;
                self             = self.setPdegeomEdges(geom());
                self             = self.setPdegeomFaceCount();
                self.NumEdges    = self.geom();
                vertices         = getPdegeomVertices(self);
            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'});
                % check that unique sorted face ids have consecutive values
                % without gaps (a gap in ids indicates missing face)
                if(any(diff(unique([geom(6,:) geom(7,:)]))-1))
                    error(message('pde:pdeGeometry:invalidFaceNumbering'))
                end
                self.geom        = geom;
                self.NumCells    = 0;
                self             = self.setDecsgEdges(size(geom, 2));
                self             = self.setDecsgFaceCount();
                self.NumEdges    = size(self.geom, 2);
                vertices         = getDecsgVertices(self);
            else
                error(message('pde:pdeGeometry:invalidArg'))
            end
            % set unique vertices sorted by index
            [~,I,~] = uniquetol(vertices,'ByRows',true);
            self.Vertices    = vertices(sort(I),:);
            self.NumVertices = size(self.Vertices,1);
        end

        function [varargout] = addFace(self,varargin)
            % addFace - Adds an internal face to geometry
            %
            %   NewGM = addFace(GM,EDGES) adds new geometric faces bounded by the
            %   specified edges to geometry GM.
            %   EDGES is a vector of IDs of edges or a cell array of such vectors.
            %   Edges specified in each vector must form a unique closed contour.
            %   Each contour in the cell array must be unique.
            %   NewGM is a handle to the geometry object GM.
            %
            %   [NewGM,FaceID] = addFace(GM, EDGES) also returns a row vector FaceID
            %   containing IDs of added faces.
            %
            %   A new face in a 2-D geometry must fill an internal void region in the
            %   geometry. A new face in a 3-D geometry splits a geometric cell and
            %   becomes a boundary between new cells. For both 2-D and 3-D geometries,
            %   you cannot add a new face on an exterior geometric boundary.
            %   
            %   New faces modify a geometry, but they do not modify the mesh. After
            %   adding faces, always call generateMesh. Creating a new mesh ensures
            %   a proper mesh association with the new geometry.
            %
            %   Example: Add a face to a 2-D geometry to fill an internal void region.
            %
            %       model = createpde;
            %       gm = importGeometry(model,'PlateHolePlanar.stl')
            %       newgm = addFace(gm,5)
            %       pdegplot(gm,'FaceLabels','on')
            %
            %   See also pde.DiscreteGeometry, pde.EquationModel/importGeometry,
            %            pde.DiscreteGeometry/addVertex,
            %            pde.EquationModel/generateMesh

            nargoutchk(0,2);
            narginchk(2,2);
            if ~isscalar(self)
                error(message('pde:pdeGeometryClass:GeometryMustBeScalar'))
            end
            if (isa(self.geom, 'function_handle'))
                error(message('pde:pdeModel:functionForNonFunctionHandleGeometry'))
            end
            contours = varargin{1};
            
            % emulate addFace on discrete geometry to check for errors
            facets = facetAnalyticGeometry(createpde(),self,0);
            [~,fids] = addFace(analyticToDiscrete(facets),contours);
            clearvars facets;

            % replace 0-regions by new face for all input edges
            dl = self.geom;
            if iscell(contours)
                for i = 1:numel(fids)
                    edges = contours{i};
                    dl(6,edges(dl(6,edges) == 0)) = fids(i);
                    dl(7,edges(dl(7,edges) == 0)) = fids(i);
                end
            else
                dl(6,contours(dl(6,contours) == 0)) = fids;
                dl(7,contours(dl(7,contours) == 0)) = fids;
            end
            self.geom = dl;
            self = setDecsgFaceCount(self);
            varargout{1} = self;
            if nargout == 2
                varargout{2} = fids;
            end
        end

        function newgm = extrude(self,varargin)
            % extrude - Create 3-D geometry by vertically extruding 2-D geometry
            %
            %   NewGM = extrude(GM,HEIGHT) creates a 3-D discrete geometry by
            %   extruding a 2-D geometry by HEIGHT along the z-direction.
            %   You can create a stacked multilayered 3-D discrete geometry by
            %   specifying a vector HEIGHT of thicknesses of the layers.
            %   The geometry GM is a 2-D DiscreteGeometry or AnalyticGeometry
            %   object.
            %   HEIGHT is a vector of thicknesses.
            %   NewGM is a handle to the geometry object GM if GM is a 2-D
            %   DiscreteGeometry object.
            %   NewGM is a new geometry object if GM is an AnalyticGeometry
            %   object.
            %
            %   Example: Extrude a 2-D discrete geometry.
            %
            %       model = createpde;
            %       gm = importGeometry(model,'PlateHolePlanar.stl')
            %       newgm = extrude(gm,1)
            %       pdegplot(newgm,'FaceLabels','on')
            %
            %   Example: Extrude an analytic geometry.
            %
            %       gd = [1,0.5,-0.25,0.25]';
            %       dl = decsg(gd);
            %       model = createpde;
            %       gm = geometryFromEdges(model,dl);
            %       newgm = extrude(gm,1)
            %       pdegplot(newgm,'FaceLabels','on')
            %
            %   See also pde.DiscreteGeometry, pde.AnalyticGeometry
            %            pde.EquationModel/importGeometry, geometryFromEdges,
            %            pde.EquationModel/generateMesh
            %
            newgm = transformToDiscrete(self);
            newgm = extrude(newgm, varargin{:});
        end
        
        function newgm = rotate(self,varargin)
            % rotate - Rotates a geometry
            % 
            % ROTATE(geom,theta) rotates the geometry geom around the Z-axis by the 
            % angle theta. Theta must be specified in degrees. The direction of 
            % rotation is based on the right-hand rule.
            % 
            % ROTATE(geom,theta,refpoint) uses the rotation axis specified by the 
            % reference point, refpoint. the axis of rotation is the line in the 
            % Z-direction passing through the reference point.
            % 
            % ROTATE modifies a geometry, but it does not modify a mesh. After 
            % modifying a geometry, re-generate the mesh to ensure a proper 
            % mesh association with the modified geometry.
            %
            % Example 1: Rotate a L-shape geometry around the Z-axis by 45 degrees.
            %     model = createpde();
            %     geometryFromEdges(model,@lshapeg);
            %     model.Geometry = rotate(model.Geometry,45);
            %     pdegplot(model)
            %
            % Example 2: Specify the rotation axis for a decomposed geometry 
            % matrix of a square geometry.
            %     model = createpde();
            %     dl = [2 2 2 2; 0 1 1 0; 1 1 0 0; 0 0 1 1;
            %           0 1 1 0; 1 1 1 1; 0 0 0 0];
            %     geometryFromEdges(model,dl);
            %     model.Geometry = rotate(model.Geometry,45,[1 1]);
            %     pdegplot(model)
            %
            if ~isscalar(self)
                error(message('pde:pdeGeometryClass:GeometryMustBeScalar'));            
            end
            newgm = transformToDiscrete(self);
            newgm = rotate(newgm, varargin{:});
        end
        
        function newgm = scale(self,varargin)
            % scale - Scales a geometry
            % 
            % SCALE(geom,s) scales the geometry geom by the factor s with respect to 
            % the origin. s can be a scalar or a vector. Use a scalar for 
            % uniform scaling in all directions. Use a vector of 2 elements 
            % to specify different scaling factors along the x-, and y-axis.
            % 
            % SCALE(geom,s,refpoint) scales the geometry with respect to the reference 
            % point, refpoint. 
            % 
            % SCALE modifies a geometry, but it does not modify a mesh. After modifying 
            % a geometry, re-generate the mesh to ensure proper mesh association with 
            % the modified geometry.
            % 
            % Example 1: Scale a 2-D geometry by a factor of 2.
            %     model = createpde();
            %     geometryFromEdges(model,@lshapeg);
            %     model.Geometry = scale(model.Geometry,2);
            %     pdegplot(model)
            %
            % Example 2:  Scale a square geometry by factors of 2 and 3 along 
            % the x- and y- respectively
            %     model = createpde();
            %     geometryFromEdges(model,@squareg);
            %     model.Geometry = scale(model.Geometry,45,[2 3]);
            %     pdegplot(model)
            %
            % Example 3: Scale a decomposed geometry matrix of a circle by a 
            % factor of 2 around the point [1 1].
            %     model = createpde();
            %     dl = [1  1 1 1; -1 0 1 0; 0  1 0 -1; 0 -1 0 1;-1 0 1 0;
            %           1  1 1 1; 0  0 0 0; 0  0 0 0; 0  0 0 0; 1  1 1 1];
            %     geometryFromEdges(model,dl);
            %     model.Geometry = scale(model.Geometry,2, [1 1]);
            %     pdegplot(model)
            %
            if ~isscalar(self)
                error(message('pde:pdeGeometryClass:GeometryMustBeScalar'));            
            end
            newgm = transformToDiscrete(self);
            newgm = scale(newgm, varargin{:});
        end
        
        function newgm = translate(self,varargin)
            % translate - Translates a geometry
            %
            % TRANSLATE(geom,s) translates the geometry geom by the distance s.  
            % The distance s is a vector of 2 elements for a 2-D geometry. The 
            % elements specify the distance along the x-, and y-axis.
            %
            % TRANSLATE modifies a geometry, but it does not modify the mesh. After 
            % modifying a geometry, re-generate the mesh to ensure proper mesh 
            % association with the modified geometry.
            %
            % Example 1: Translate a 2-D geometry by 1 along the X-axis and by 2 along 
            % the Y-axis 
            %     model = createpde();
            %     geometryFromEdges(model,@lshapeg);
            %     model.Geometry = translate(model.Geometry,[1 2]);
            %     pdegplot(model) 
            %
            % Example 2: Translate a decomposed geometry matrix of a circle by 
            % 2 along the X-axis and by 3 along the Y-axis.
            %     model = createpde();
            %     dl = [1 1 1 1; -1 0 1 0; 0 1 0 -1; 0 -1 0 1;-1 0 1 0;
            %           1 1 1 1;  0 0 0 0;  0 0 0 0; 0  0 0 0; 1 1 1 1];
            %     geometryFromEdges(model,dl);
            %     model.Geometry = translate(model.Geometry,[2 3]);
            %     pdegplot(model);
            %
            if ~isscalar(self)
                error(message('pde:pdeGeometryClass:GeometryMustBeScalar'));            
            end
            newgm = transformToDiscrete(self);
            newgm = translate(newgm, varargin{:});
        end
                
        function fae=facesAttachedToEdges(self, eids,filterType)
            % facesAttachedToEdges - Find faces attached to specified edges
            % 
            % facesAttachedToEdges(GEOM,REGIONID,FILTERTYPE) finds faces 
            % attached to the edges with ID numbers listed in REGIONID. Here, 
            % GEOM is a DiscreteGeometry or AnalyticGeometry object. 
            % REGIONID is a single edge ID or a vector of Edge IDs.
            % 
            % For a 3-D geometry, FILTERTYPE is the optional argument 
            % specifying the type of the returned faces:
            % 'internal' - return faces shared between multiple cells.
            % 'external' - return faces not shared between multiple cells.
            % 'all' - return all faces attached to the specified edges. 
            % This is the default option.
            % 
            % 
            % Example 1: Find faces attached to edges 1 and 2 in the L-shape 
            % membrane.
            % model = createpde;
            % geom = geometryFromEdges(model,@lshapeg);
            % facesAttachedToEdges(geom,[1 2])
            %
            arguments
                self (1,1) pde.AnalyticGeometry
                eids (1,:) {mustBeInteger,mustBeNumeric,mustBeNonsparse,mustBePositive}
                filterType (1,1) string = "all"
            end
            filterType = partialMatch(filterType);
            if (filterType == "internal")
                fae = [];
            else  % all edges in a 2D geom are external
                if (isa(self.geom,'function_handle'))
                    nedges = self.geom();
                    if (any(eids > nedges))
                        error(message('pde:pdeGeometryClass:InvalidEdgeIndex'));
                    else
                        selr1 = self.geom(1:nedges);
                        fae = nonzeros(unique((selr1(3:4,eids))))';
                    end
                elseif (isnumeric(self.geom))
                    nedges = size(self.geom(),2);
                    if (any(eids > nedges))
                        error(message('pde:pdeGeometryClass:InvalidEdgeIndex'));
                    else
                        fae = nonzeros(unique((self.geom(6:7,eids))))'; 
                    end                                     
                end
            end            
        end

        function nf = nearestFace(self,coords)
            % nearestFace - Find faces nearest to specified point
            %
            % nearestFace(GEOM,COORDS) finds faces of the geometry GEOM 
            % nearest to the point with the coordinates COORDS. A coordinate 
            % vector COORDS has 2 or 3 elements for a 2-D or 3-D geometry, 
            % respectively.
            %
            % Example 1: Find external faces attached to an analytic geometry
            %
            %   pdem=createpde;
            %   geom = pdem.geometryFromEdges(@lshapeg); 
            %   nearestFace(geom,[-0.1 0.6]) 
            %   pdegplot(geom,'EdgeLabels','on','FaceLabels','on','FaceAlpha',0.2);
            %
            arguments
                self (1,1) pde.AnalyticGeometry
                coords
            end
            model = createpde;
            facets = facetAnalyticGeometry(model,self,0);
            newgm = analyticToDiscrete(facets);
            nf = nearestFace(newgm,coords);
        end
                
        function eaf=faceEdges(self,fids,filterType)
            % faceEdges - Find edges attached to specified faces
            %
            % faceEdges(GEOM,REGIONID,FILTERTYPE) finds edges of the 
            % geometry GEOM attached to the faces with ID numbers listed in 
            % REGIONID. REGIONID is a single face ID or a vector of face IDs.
            %
            % For a 3-D geometry, FILTERTYPE is the optional argument specifying
            % the type of the returned edges:
            % 'internal' - return edges attached to internal faces. Internal
            % faces are faces shared between multiple cells.
            % 'external' - return edges attached to external faces. External
            % faces are faces not shared between multiple cells.
            % 'all' - return all edges attached to the specified faces. This is
            % the default option.
            %
            % Example 1: Find edges attached to a face in an L-shape membrane.
            % model = createpde;
            % geom = geometryFromEdges(model,@lshapeg);
            % faceEdges(geom,1)
            %       
            arguments
                self (1,1) pde.AnalyticGeometry
                fids (1,:) {mustBeInteger,mustBeNumeric,mustBeNonsparse,mustBePositive}
                filterType (1,1) string = "all"
            end
            
            filterType = partialMatch(filterType);
            if (filterType == "internal")
                eaf = [];
            else  % all edges in a 2D geom are external
                if (isa(self.geom,'function_handle'))
                    nedges = self.geom();
                    selr1 = self.geom(1:nedges);
                    ef = selr1(3:4,:);
                    nfaces = max(max(ef));
                    if (any(fids > nfaces))
                        error(message('pde:pdeGeometryClass:InvalidFaceIndex'));
                    else
                        eaf = [];
                        for i=1:numel(fids)                            
                            ef2 = find(any(ef == fids(i)));
                            eaf = [eaf ef2];
                        end                        
                    end
                    eaf = unique(eaf);  
                elseif (isnumeric(self.geom))
                    ef = self.geom(6:7,:);
                    nfaces = max(max(ef));
                    if (any(fids > nfaces))
                        error(message('pde:pdeGeometryClass:InvalidFaceIndex'));
                    else
                        eaf = [];
                        for i=1:numel(fids)                            
                            ef2 = find(any(ef == fids(i)));
                            eaf = [eaf ef2];
                        end                        
                    end
                    eaf = unique(eaf);                  
                end
            end
        end       
        
        function ne = nearestEdge(self,coords)
            % nearestEdge - Find edges nearest to specified point
            %
            % nearestEdge(GEOM,COORDS) finds edges of the geometry GEOM 
            % nearest to the point with the coordinates COORDS. A coordinate 
            % vector COORDS has 2 or 3 elements for a 2-D or 3-D geometry,
            % respectively.
            %
            % Example 1: Find edges nearest to the specified point.
            % model = createpde;
            % geom = geometryFromEdges(model,@lshapeg);
            % nearestEdge(geom,[-0.9 -0.5])
            %
            arguments
                self (1,1) pde.AnalyticGeometry
                coords
            end
            model = createpde;
            facets = facetAnalyticGeometry(model,self,0);
            newgm = analyticToDiscrete(facets);
            ne = nearestEdge(newgm,coords);
        end
    end
    
    methods(Hidden=true)
        function vxl = vertexLabelLocations(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,:)';
            vxl = [sx sy; ex ey];
            [~,I,~] = uniquetol(vxl,'ByRows',true);
            vxl = vxl(sort(I),:);
            vxl = {vxl(:,1), vxl(:,2)};
        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 = addVertex(~,varargin)
            error(message('pde:pdeModel:functionForDiscreteGeomOnly'));
        end

        function eac=cellEdges(~,~,~)
            error(message('pde:pdeGeometryClass:GeometryMustBe3D'));
        end

        function fac = cellFaces(~,~,~)
            error(message('pde:pdeGeometryClass:GeometryMustBe3D'));
        end
    end
    methods(Access=private)
        function self=setPdegeomEdges(self, numEdges)
            for i=1:numEdges
                self.Edges(i) = pde.internal.pdeEdge(self, i);
            end
        end
        function self=setPdegeomFaceCount(self)
            nedges = self.geom();
            selr = self.geom(1:nedges);
            faceids = selr(3:4,:);
            faceids = unique(faceids(:));
            zid = (faceids == 0);
            faceids(zid) = [];
            self.NumFaces=numel(faceids);
        end

        function v=getPdegeomVertices(self)
            nedges = self.geom();
            selr = self.geom(1:nedges);
            sp = selr(1,:);
            eidx = 1:nedges;
            [sx, sy] = self.geom(eidx(:),sp(:));
            ep = selr(2,:);
            [ex, ey] = self.geom(eidx(:),ep(:));
            v = [sx(:) sy(:); ex(:) ey(:)];
        end

        function self=setDecsgEdges(self, 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 v = getDecsgVertices(self)
            sx = self.geom(2,:)';
            ex = self.geom(3,:)';
            sy = self.geom(4,:)';
            ey = self.geom(5,:)';
            v = [sx sy; ex ey]; 
        end
 
        function dgm=PdegeomToDecsg(self)
            geomFunc = self.geom;
            nedges = geomFunc();
            edgeParams=geomFunc(1:nedges);
            s = edgeParams(1,:);
            [xs,ys] = geomFunc(1:nedges, s);
            s = edgeParams(2,:);
            [xe,ye] = geomFunc(1:nedges, s);
            
            dgm = zeros(10,nedges);
            dgm(1,:) = 2;
            dgm(2,:) = xs;
            dgm(3,:) = xe;
            dgm(4,:) = ys;
            dgm(5,:) = ye;
            dgm(6,:) = edgeParams(3,:);
            dgm(7,:) = edgeParams(4,:);
        end
        
        function checkFunctionHandle(~,geom)
            nedges = geom();
            if (nedges ~= geom())
                error(message('pde:pdeGeometry:invalidGeomFunction'))
            end
            atol = 1e-8;
            selr1 = geom(1:nedges);
            selr2 = geom(1:nedges);
            if (max(abs(selr1(:) - selr2(:))) > atol * max(abs(selr1(:))))
                error(message('pde:pdeGeometry:invalidGeomFunction'))
            end
        end
    end
end

% Helper functions

function str = partialMatch(str)
    validStrings = ["all","external","internal"];
    str = validatestring(str,validStrings);
end

function newgm = transformToDiscrete(geom)
    mdl = createpde;
    facets = facetAnalyticGeometry(mdl,geom,0);
    newgm = analyticToDiscrete(facets);
end
